Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

1201 lines (966 sloc) 35.262 kb
/****************************************************************************
* Copyright (C) 2009-2011 GGA Software Services LLC
*
* This file is part of Indigo toolkit.
*
* This file may be distributed and/or modified under the terms of the
* GNU General Public License version 3 as published by the Free Software
* Foundation and appearing in the file LICENSE.GPL included in the
* packaging of this file.
*
* This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
* WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
***************************************************************************/
#include "layout/molecule_layout_graph.h"
#include "layout/attachment_layout.h"
#include "graph/biconnected_decomposer.h"
#include "graph/cycle_enumerator.h"
#include "graph/embedding_enumerator.h"
#include "graph/morgan_code.h"
using namespace indigo;
enum
{
QUERY_BOND_SINGLE_OR_DOUBLE = 5,
QUERY_BOND_SINGLE_OR_AROMATIC = 6,
QUERY_BOND_DOUBLE_OR_AROMATIC = 7,
QUERY_BOND_ANY = 8
};
// Make relative coordinates of a component absolute
void MoleculeLayoutGraph::_copyLayout (MoleculeLayoutGraph &component)
{
int i;
for (i = component.vertexBegin(); i < component.vertexEnd(); i = component.vertexNext(i))
{
LayoutVertex &vert = component._layout_vertices[i];
_layout_vertices[vert.ext_idx].pos.copy(vert.pos);
_layout_vertices[vert.ext_idx].type = vert.type;
}
for (i = component.edgeBegin(); i < component.edgeEnd(); i = component.edgeNext(i))
{
LayoutEdge &edge = component._layout_edges[i];
_layout_edges[edge.ext_idx].type = edge.type;
}
}
static int _vertex_cmp (int &n1, int &n2, void *context)
{
const MoleculeLayoutGraph &graph = *(MoleculeLayoutGraph *)context;
const LayoutVertex &v1 = graph.getLayoutVertex(n1);
const LayoutVertex &v2 = graph.getLayoutVertex(n2);
if (v1.is_cyclic != v2.is_cyclic)
{
if (v1.is_cyclic == true)
return 1;
return -1;
}
return v1.morgan_code - v2.morgan_code;
}
void MoleculeLayoutGraph::_assignAbsoluteCoordinates (float bond_length)
{
BiconnectedDecomposer bc_decom(*this);
QS_DEF(Array<int>, bc_tree);
QS_DEF(ObjArray<MoleculeLayoutGraph>, bc_components);
QS_DEF(Array<int>, fixed_components);
bool all_trivial = true;
int i, j = -1;
int n_comp = bc_decom.decompose();
fixed_components.clear_resize(n_comp);
fixed_components.zerofill();
bc_components.clear();
for (i = 0; i < n_comp; i++)
{
Filter comp;
bc_decom.getComponent(i, comp);
bc_components.push().makeLayoutSubgraph(*this, comp);
}
bc_tree.clear_resize(vertexEnd());
_makeComponentsTree(bc_decom, bc_components, bc_tree);
// 1. Find biconnected components forming connected subgraph from fixed vertices
if (_n_fixed > 0)
{
QS_DEF(Array<int>, fixed_count);
fixed_count.clear_resize(n_comp);
fixed_count.zerofill();
// calculate number of fixed vertices in each component
for (i = 0; i < n_comp; i++)
{
Filter filter;
bc_decom.getComponent(i, filter);
for (j = vertexBegin(); j < vertexEnd(); j = vertexNext(j))
if (filter.valid(j) && _fixed_vertices[j])
fixed_count[i]++;
}
// keep only with fixed number greater than a half
for (i = 0; i < n_comp; i++)
{
Filter filter;
bc_decom.getComponent(i, filter);
if (fixed_count[i] > filter.count(*this) / 2)
fixed_components[i] = 1;
}
_fixed_vertices.zerofill();
// update fixed vertices
for (i = 0; i < n_comp; i++)
{
if (!fixed_components[i])
continue;
MoleculeLayoutGraph &component = bc_components[i];
for (j = component.vertexBegin(); j < component.vertexEnd(); j = component.vertexNext(j))
_fixed_vertices[component.getVertexExtIdx(j)] = 1;
}
Filter fixed_filter(_fixed_vertices.ptr(), Filter::EQ, 1);
Graph fixed_graph;
QS_DEF(Array<int>, fixed_mapping);
fixed_graph.makeSubgraph(*this, fixed_filter, &fixed_mapping, 0);
if (Graph::isConnected(fixed_graph))
_n_fixed = fixed_filter.count(*this);
else
{ // fixed subgraph is not connected - choose its greatest component
int n = fixed_graph.countComponents();
const Array<int> &decomposition = fixed_graph.getDecomposition();
fixed_count.clear_resize(n);
fixed_count.zerofill();
for (i = fixed_graph.vertexBegin(); i < fixed_graph.vertexEnd(); i = fixed_graph.vertexNext(i))
fixed_count[decomposition[i]]++;
j = 0;
for (i = 1; i < n; i++)
if (fixed_count[i] > fixed_count[j])
j = i;
Filter max_filter(decomposition.ptr(), Filter::EQ, j);
// update fixed vertices
_fixed_vertices.zerofill();
_n_fixed = 0;
for (i = fixed_graph.vertexBegin(); i < fixed_graph.vertexEnd(); i = fixed_graph.vertexNext(i))
{
if (max_filter.valid(i))
{
_fixed_vertices[fixed_mapping[i]] = 1;
_n_fixed++;
}
}
for (i = 0; i < n_comp; i++)
{
if (!fixed_components[i])
continue;
MoleculeLayoutGraph &component = bc_components[i];
if (!max_filter.valid(component.getVertexExtIdx(component.vertexBegin())))
fixed_components[i] = 0;
}
}
}
// Possible solutions:
// 1. a) vertex code is calculated inside component (doesn't depend on neighbors) or
// b) vertex code is calculated respecting whole graph
// 2. a) component code is the sum of 1a codes
// b) component code is the sum of 1b codes
// Initially was 1a and 2b then changed to 1b and 2b
for (i = 0; i < n_comp; i++)
{
MoleculeLayoutGraph &component = bc_components[i];
component.max_iterations = max_iterations;
//component._calcMorganCodes();
component._total_morgan_code = 0;
for (j = component.vertexBegin(); j < component.vertexEnd(); j = component.vertexNext(j))
component._total_morgan_code += _layout_vertices[component.getLayoutVertex(j).ext_idx].morgan_code;
// Mark cyclic atoms
if (!component.isSingleEdge())
{
all_trivial = false;
for (j = component.vertexBegin(); j < component.vertexEnd(); j = component.vertexNext(j))
{
component._layout_vertices[j].is_cyclic = true;
_layout_vertices[component._layout_vertices[j].ext_idx].is_cyclic = true;
}
for (j = component.edgeBegin(); j < component.edgeEnd(); j = component.edgeNext(j))
{
component._layout_edges[j].is_cyclic = true;
_layout_edges[component._layout_edges[j].ext_idx].is_cyclic = true;
}
}
int fixed = fixed_components[i];
component._assignRelativeCoordinates(fixed, *this);
if (fixed != fixed_components[i])
{
fixed_components[i] = fixed;
// update fixed vertices
_fixed_vertices.resize(vertexEnd());
_fixed_vertices.zerofill();
_n_fixed = 0;
for (j = 0; j < n_comp; j++)
{
if (!fixed_components[j])
continue;
Filter fix_filter;
bc_decom.getComponent(j, fix_filter);
for (int k = vertexBegin(); k < vertexEnd(); k = vertexNext(k))
if (!_fixed_vertices[k])
{
_fixed_vertices[k] = 1;
_n_fixed++;
}
}
}
}
if (_n_fixed > 0)
{
for (i = 0; i < n_comp; i++)
if (fixed_components[i])
{
_copyLayout(bc_components[i]);
j = i;
}
MoleculeLayoutGraph &component = bc_components[j];
_first_vertex_idx = component._layout_vertices[component.vertexBegin()].ext_idx;
} else
{
// ( 0]. Nucleus.;
// Begin from nontrivial component with maximum code
// if there's no then begin from vertex with maximum code and its neighbor with maximum code too
int nucleus_idx = 0;
if (!all_trivial)
{
nucleus_idx = -1;
for (i = 0; i < n_comp; i++)
{
MoleculeLayoutGraph &component = bc_components[i];
if (!component.isSingleEdge())
{
if (nucleus_idx == -1 || component._total_morgan_code > bc_components[nucleus_idx]._total_morgan_code)
nucleus_idx = i;
}
}
if (nucleus_idx < 0)
throw Error("Internal error: cannot find nontrivial component");
MoleculeLayoutGraph &nucleus = bc_components[nucleus_idx];
_copyLayout(nucleus);
_first_vertex_idx = nucleus._layout_vertices[nucleus._first_vertex_idx].ext_idx;
} else
{
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
if (_layout_vertices[i].morgan_code > _layout_vertices[nucleus_idx].morgan_code)
nucleus_idx = i;
const Vertex &nucleus = getVertex(nucleus_idx);
int nucleus_idx2 = nucleus.neiBegin();
for (j = nucleus.neiNext(nucleus_idx2); j < nucleus.neiEnd(); j = nucleus.neiNext(j))
if (_layout_vertices[nucleus.neiVertex(j)].morgan_code > _layout_vertices[nucleus.neiVertex(nucleus_idx2)].morgan_code)
nucleus_idx2 = j;
int nucleus_edge = nucleus.neiEdge(nucleus_idx2);
nucleus_idx2 = nucleus.neiVertex(nucleus_idx2);
_first_vertex_idx = nucleus_idx;
_layout_vertices[nucleus_idx].type = ELEMENT_BOUNDARY;
_layout_vertices[nucleus_idx].pos.set(0.f, 0.f);
_layout_vertices[nucleus_idx2].type = ELEMENT_BOUNDARY;
_layout_vertices[nucleus_idx2].pos.set(1.f, 0.f);
_layout_edges[nucleus_edge].type = ELEMENT_BOUNDARY;
}
}
// ( 1] atoms assigned absolute coordinates and adjacent to atoms not;
// assigned coordinates are put on a list;
QS_DEF(Array<int>, assigned_list);
QS_DEF(Array<int>, adjacent_list);
while (true)
{
assigned_list.clear();
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
{
if (_layout_vertices[i].type == ELEMENT_NOT_DRAWN)
continue;
const Vertex &vert = getVertex(i);
for (j = vert.neiBegin(); j < vert.neiEnd(); j = vert.neiNext(j))
{
if (_layout_vertices[vert.neiVertex(j)].type == ELEMENT_NOT_DRAWN)
{
assigned_list.push(i);
break;
}
}
}
if (assigned_list.size() == 0)
{
_refineCoordinates(bc_decom, bc_components, bc_tree);
return;
}
// ( 2] the list is ordered with cyclic atoms at the top of the list;
// with descending ATCD numbers and acyclic atoms at the bottom;
// of the list with descending ATCD numbers;;
assigned_list.qsort(_vertex_cmp, this);
// ( 3.i] let k = 0 ( top of the list];;
while (assigned_list.size() != 0)
{
int k = assigned_list.pop();
const Vertex &vert_k = getVertex(k);
// ( 3.ii] a list of atoms adjacent to atom Uzel and not previously;
// assigned coordinates is created and ordered with cyclic atoms;
// at the top of the list with descending ATCD numbers and acyclic atoms;
// at the bottom of the list with descending ATCD numbers;;
adjacent_list.clear();
for (i = vert_k.neiBegin(); i < vert_k.neiEnd(); i = vert_k.neiNext(i))
if (_layout_vertices[vert_k.neiVertex(i)].type == ELEMENT_NOT_DRAWN)
adjacent_list.push(vert_k.neiVertex(i));
if (adjacent_list.size() == 0)
break;
// When all components outgoing from vertex are trivial (edges) then use tree algorithm
all_trivial = true;
for (i = 0; i < bc_decom.getIncomingCount(k); i++)
if (!bc_components[bc_decom.getIncomingComponents(k)[i]].isSingleEdge())
{
all_trivial = false;
break;
}
if (all_trivial && bc_tree[k] != -1 && !bc_components[bc_tree[k]].isSingleEdge())
all_trivial = false;
if (all_trivial)
{
adjacent_list.qsort(_vertex_cmp, this);
_attachDandlingVertices(k, adjacent_list);
} else
{
// Component layout in current vertex should have the same angles between components.
// So it depends on component order and their flipping (for nontrivial components)
AttachmentLayout att_layout(bc_decom, bc_components, bc_tree, *this, k);
// ( 3.iii] Look over all possible orders of component layouts
// (vertex itself is already drawn means one component is already drawn)
// ( 3.iv] Choose layout with minimal energy
LayoutChooser layout_chooser(att_layout);
layout_chooser.perform();
att_layout.markDrawnVertices();
}
// ( 3.v] let k = k + 1;;
// ( 3.vi] repeat steps 3.ii-3.v until all atoms in the list have been processed;;
}
// ( 4] repeat steps 1-3 until all atoms have been assigned absolute coordinates.;
}
}
bool MoleculeLayoutGraph::_match_pattern_bond (Graph &subgraph, Graph &supergraph, int self_idx, int other_idx, void *userdata)
{
if (userdata == 0 || ((MoleculeLayoutGraph *)userdata)->_molecule == 0)
return true;
BaseMolecule &mol = *((MoleculeLayoutGraph *)userdata)->_molecule;
const int *mapping = ((MoleculeLayoutGraph *)userdata)->_molecule_edge_mapping;
int layout_idx = ((const MoleculeLayoutGraph &)supergraph).getLayoutEdge(other_idx).ext_idx;
const PatternBond &pattern_bond = ((const PatternLayout &)subgraph).getBond(self_idx);
switch (pattern_bond.type)
{
case BOND_SINGLE:
case BOND_DOUBLE:
case BOND_TRIPLE:
case BOND_AROMATIC:
if (!mol.possibleBondOrder(mapping[layout_idx], pattern_bond.type))
return false;
break;
case QUERY_BOND_SINGLE_OR_DOUBLE:
if (!mol.possibleBondOrder(mapping[layout_idx], BOND_SINGLE) &&
!mol.possibleBondOrder(mapping[layout_idx], BOND_DOUBLE))
return false;
break;
case QUERY_BOND_SINGLE_OR_AROMATIC:
if (!mol.possibleBondOrder(mapping[layout_idx], BOND_SINGLE) &&
!mol.possibleBondOrder(mapping[layout_idx], BOND_AROMATIC))
return false;
break;
case QUERY_BOND_DOUBLE_OR_AROMATIC:
if (!mol.possibleBondOrder(mapping[layout_idx], BOND_DOUBLE) &&
!mol.possibleBondOrder(mapping[layout_idx], BOND_AROMATIC))
return false;
break;
}
int parity = mol.cis_trans.getParity(mapping[layout_idx]);
if (parity != 0 && parity != pattern_bond.parity)
return false;
return true;
}
int MoleculeLayoutGraph::_pattern_embedding (Graph &subgraph, Graph &supergraph, int *core_sub, int *core_super, void *userdata)
{
if (userdata == 0)
return 1;
MoleculeLayoutGraph &layout_graph = *(MoleculeLayoutGraph *)userdata;
const PatternLayout &pattern_graph = (const PatternLayout &)subgraph;
int i;
// TODO: correct element marking (internal and non-planar)?
for (i = layout_graph.vertexBegin(); i < layout_graph.vertexEnd(); i = layout_graph.vertexNext(i))
{
layout_graph._layout_vertices[i].pos = pattern_graph.getAtom(core_super[i]).pos;
layout_graph._layout_vertices[i].type = ELEMENT_BOUNDARY;
}
for (i = layout_graph.edgeBegin(); i < layout_graph.edgeEnd(); i = layout_graph.edgeNext(i))
layout_graph._layout_edges[i].type = ELEMENT_BOUNDARY;
layout_graph._first_vertex_idx = layout_graph.vertexBegin();
if (layout_graph._outline.get() == 0)
layout_graph._outline.create();
layout_graph._outline->copy(pattern_graph.getOutline());
return 0;
}
void MoleculeLayoutGraph::_assignRelativeCoordinates (int &fixed_component, const MoleculeLayoutGraph &supergraph)
{
int i;
if (isSingleEdge())
{
// Trivial component layout
int idx1 = vertexBegin();
int idx2 = vertexNext(idx1);
_layout_vertices[idx1].type = ELEMENT_BOUNDARY;
_layout_vertices[idx2].type = ELEMENT_BOUNDARY;
if (fixed_component)
{
_layout_vertices[idx1].pos = supergraph.getPos(getVertexExtIdx(idx1));
_layout_vertices[idx2].pos = supergraph.getPos(getVertexExtIdx(idx2));
} else
{
_layout_vertices[idx1].pos.set(0.f, 0.f);
_layout_vertices[idx2].pos.set(0.f, 1.f);
}
_layout_edges[edgeBegin()].type = ELEMENT_BOUNDARY;
return;
}
// 2.1. Use layout of fixed components and find border edges and vertices
if (fixed_component)
{
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos = supergraph.getPos(getVertexExtIdx(i));
CycleEnumerator ce(*this);
ce.context = this;
ce.cb_handle_cycle = _border_cb;
if (ce.process())
return;
fixed_component = 0;
} else
{ // try to find pattern
TL_GET(ObjArray<PatternLayout>, _patterns);
MorganCode morgan(*this);
QS_DEF(Array<long>, morgan_codes);
morgan.calculate(morgan_codes, 3, 7);
long morgan_code = 0;
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
morgan_code += morgan_codes[i];
int left = 0;
int right = _patterns.size() - 1;
int pat_idx = 0;
int cmp;
while (left < right)
{
if (right - left == 1)
{
if (_pattern_cmp2(_patterns[left], vertexCount(), edgeCount(), morgan_code) == 0)
pat_idx = left;
else if (_pattern_cmp2(_patterns[right], vertexCount(), edgeCount(), morgan_code) == 0)
pat_idx = right;
break;
}
pat_idx = (right + left) / 2;
cmp = _pattern_cmp2(_patterns[pat_idx], vertexCount(), edgeCount(), morgan_code);
if (cmp < 0)
left = pat_idx;
else
right = pat_idx;
}
while (pat_idx > 0 && _pattern_cmp2(_patterns[pat_idx - 1], vertexCount(), edgeCount(), morgan_code) == 0)
pat_idx--;
while (pat_idx < _patterns.size() && _pattern_cmp2(_patterns[pat_idx], vertexCount(), edgeCount(), morgan_code) == 0)
{
// Match pattern
// TODO: check different attachment points
PatternLayout &pattern = _patterns[pat_idx];
EmbeddingEnumerator ee(*this);
ee.setSubgraph(pattern);
ee.cb_match_edge = _match_pattern_bond;
ee.cb_embedding = _pattern_embedding;
ee.userdata = this;
if (!ee.process())
{
if (pattern.isFixed())
fixed_component = 1;
return;
}
pat_idx++;
}
}
//TODO: repair exception with vec2f
QS_DEF(ObjPool<Cycle>, cycles);
QS_DEF(Array<int>, sorted_cycles);
cycles.clear();
int n_cycles = sssrCount();
for (i = 0; i < n_cycles; i++)
{
int cycle_idx = cycles.add(sssrEdges(i), *this);
cycles[cycle_idx].canonize();
}
sorted_cycles.clear();
for (i = cycles.begin(); i < cycles.end(); i = cycles.next(i))
{
cycles[i].calcMorganCode(*this);
sorted_cycles.push(i);
}
sorted_cycles.qsort(Cycle::compare_cb, &cycles);
_assignFirstCycle(cycles[sorted_cycles[0]]);
cycles.remove(sorted_cycles[0]);
sorted_cycles.remove(0);
bool chain_attached;
// Try to attach chains with one, two or more common edges outside drawn part
do
{
chain_attached = false;
for (i = 0; !chain_attached && i < sorted_cycles.size(); )
{
if (_attachCycleOutside(cycles[sorted_cycles[i]], 1.f, 1))
{
cycles.remove(sorted_cycles[i]);
sorted_cycles.remove(i);
chain_attached = true;
} else
i++;
}
for (i = 0; !chain_attached && i < sorted_cycles.size(); )
{
if (_attachCycleOutside(cycles[sorted_cycles[i]], 1.f, 2))
{
cycles.remove(sorted_cycles[i]);
sorted_cycles.remove(i);
chain_attached = true;
} else
i++;
}
for (i = 0; !chain_attached && i < sorted_cycles.size(); )
{
if (_attachCycleOutside(cycles[sorted_cycles[i]], 1.f, 0))
{
cycles.remove(sorted_cycles[i]);
sorted_cycles.remove(i);
chain_attached = true;
} else
i++;
}
} while (chain_attached);
// Try to attach chains inside
for (i = 0; i < sorted_cycles.size(); )
{
if (_attachCycleInside(cycles[sorted_cycles[i]], 1.f))
{
cycles.remove(sorted_cycles[i]);
sorted_cycles.remove(i);
} else
i++;
}
// Try to attach chains inside with lower edge length
for (i = 0; i < sorted_cycles.size(); )
{
if (_attachCycleInside(cycles[sorted_cycles[i]], 0.75f))
{
cycles.remove(sorted_cycles[i]);
sorted_cycles.remove(i);
} else
i++;
}
do
{
chain_attached = false;
for (i = 0; !chain_attached && i < sorted_cycles.size(); )
{
// 1.5f (> 1) means to calculate new length;
if (_attachCycleOutside(cycles[sorted_cycles[i]], 1.5f, 0))
{
cycles.remove(sorted_cycles[i]);
sorted_cycles.remove(i);
chain_attached = true;
} else
i++;
}
} while (chain_attached);
do
{
chain_attached = false;
for (i = 0; !chain_attached && i < sorted_cycles.size(); )
{
if (_attachCycleWithIntersections(cycles[sorted_cycles[i]], 1.f))
{
cycles.remove(sorted_cycles[i]);
sorted_cycles.remove(i);
chain_attached = true;
} else
i++;
}
} while (chain_attached);
_attachCrossingEdges();
for (i = edgeBegin(); i < edgeEnd(); i = edgeNext(i))
{
if (_layout_edges[i].type == ELEMENT_NOT_PLANAR)
{
_buildOutline();
break;
}
}
}
void MoleculeLayoutGraph::_assignFirstCycle (const Cycle &cycle)
{
// TODO: Start drawing from vertex with maximum code and continue to the right with one of two which has maximum code
int i, n;
float phi;
n = cycle.vertexCount();
for (i = 0; i < n; i++)
{
_layout_vertices[cycle.getVertex(i)].type = ELEMENT_BOUNDARY;
_layout_edges[cycle.getEdge(i)].type = ELEMENT_BOUNDARY;
}
_first_vertex_idx = cycle.getVertex(0);
_layout_vertices[cycle.getVertex(0)].pos.set(0.f, 0.f);
_layout_vertices[cycle.getVertex(1)].pos.set(1.f, 0.f);
phi = (float)M_PI * (n - 2) / n;
for (i = 1; i < n - 1; i++)
{
const Vec2f &v1 = _layout_vertices[cycle.getVertex(i - 1)].pos;
const Vec2f &v2 = _layout_vertices[cycle.getVertex(i)].pos;
_layout_vertices[cycle.getVertex(i + 1)].pos.rotateAroundSegmentEnd(v1, v2, phi);
}
}
// If vertices are already drawn
// draw edges with intersections
void MoleculeLayoutGraph::_attachCrossingEdges ()
{
int i, j, pr;
bool intersection;
for (i = edgeBegin(); i < edgeEnd(); i = edgeNext(i))
{
const Edge &edge_i = getEdge(i);
if (_layout_vertices[edge_i.beg].type != ELEMENT_NOT_DRAWN &&
_layout_vertices[edge_i.end].type != ELEMENT_NOT_DRAWN &&
_layout_edges[i].type == ELEMENT_NOT_DRAWN)
{
intersection = true;
while (intersection)
{
intersection = false;
for (j = edgeBegin(); j < edgeEnd(); j = edgeNext(j))
{
if (_layout_edges[j].type != ELEMENT_NOT_DRAWN)
{
pr = _calcIntersection(i, j);
// 1. If the edge 1 ends on the edge 2
// then shift edge 2 from edge 1 by epsilon orthogonally to the edge 2
if (pr == 222 || pr == 223)
{
_shiftEdge(j, 0.2f);
intersection = true;
break;
}
if (pr == 224 || pr == 225)
{
_shiftEdge(i, 0.2f);
intersection = true;
break;
}
// 2. If the edge 1 overlaps some other edge shift it by epsilon orthogonally
if (pr == 3 || pr == 4)
{
_shiftEdge(i, 0.2f);
intersection = true;
break;
}
}
}
}
_layout_edges[i].type = ELEMENT_NOT_PLANAR;
}
}
}
void MoleculeLayoutGraph::_buildOutline (void)
{
Vec2f v, inter;
Vec2f pos_i;
int i, j;
int first_idx = vertexBegin();
float min_y = getPos(first_idx).y;
const float EPS = 0.0001f;
for (i = vertexNext(first_idx); i < vertexEnd(); i = vertexNext(i))
{
if (getPos(i).y < min_y)
{
min_y = getPos(i).y;
first_idx = i;
}
}
i = first_idx;
float max_angle, cur_angle;
float i_angle = 0;
int next_nei = 0;
pos_i = getPos(i);
if (_outline.get() == 0)
_outline.create();
else
_outline->clear();
while (true)
{
const Vertex &vert = getVertex(i);
if (i != first_idx)
{
v = pos_i;
pos_i = getPos(i);
v.sub(pos_i);
i_angle = v.tiltAngle2();
} else if (_outline->size() > 0)
break;
_outline->push(pos_i);
max_angle = 0.f;
for (j = vert.neiBegin(); j < vert.neiEnd(); j = vert.neiNext(j))
{
const Vec2f &pos_nei = getPos(vert.neiVertex(j));
v.diff(pos_nei, pos_i);
cur_angle = v.tiltAngle2() - i_angle;
if (cur_angle < 0.f)
cur_angle += 2 * PI;
if (max_angle < cur_angle)
{
max_angle = cur_angle;
next_nei = j;
}
}
i = vert.neiVertex(next_nei);
float dist, min_dist = 0.f;
int int_edge = -1;
Vec2f cur_v1 = pos_i;
Vec2f cur_v2 = getPos(i);
int prev_edge = -1;
int cur_edge = vert.neiEdge(next_nei);
while (min_dist < 10000.f)
{
min_dist = 10001.f;
for (j = edgeBegin(); j < edgeEnd(); j = edgeNext(j))
{
const Edge &edge = getEdge(j);
const Vec2f &cur_v3 = getPos(edge.beg);
const Vec2f &cur_v4 = getPos(edge.end);
if (Vec2f::intersection(cur_v1, cur_v2, cur_v3, cur_v4, v))
if ((dist = Vec2f::dist(cur_v1, v)) < min_dist)
{
if (dist > EPS || j != prev_edge)
{
inter = v;
min_dist = dist;
int_edge = j;
}
}
}
if (min_dist < 10000.f)
{
if (min_dist > EPSILON)
_outline->push(v);
const Edge &edge = getEdge(int_edge);
const Vec2f &cur_v3 = getPos(edge.beg);
const Vec2f &cur_v4 = getPos(edge.end);
Vec2f cur_v1v;
Vec2f cur_v3v;
Vec2f cur_v4v;
cur_v1v.diff(cur_v1, inter);
cur_v3v.diff(cur_v3, inter);
cur_v4v.diff(cur_v4, inter);
float angle1 = cur_v1v.tiltAngle2();
float angle3 = cur_v3v.tiltAngle2() - angle1;
float angle4 = cur_v4v.tiltAngle2() - angle1;
if (angle3 < 0)
angle3 += 2 * PI;
if (angle4 < 0)
angle4 += 2 * PI;
cur_v1 = inter;
if (angle3 > angle4)
{
cur_v2 = cur_v3;
i = edge.beg;
} else
{
cur_v2 = cur_v4;
i = edge.end;
}
prev_edge = cur_edge;
cur_edge = int_edge;
}
}
}
}
// Return 1 - with maximum code, 2 - neighbor of the 1 with maximum code
// 3 - neighbor with maximum code from the rest or -1 if it doesn't exist.
void MoleculeLayoutGraph::_getAnchor (int &v1, int &v2, int &v3) const
{
int i;
if (vertexCount() ==1)
{
v1 = v2 = vertexBegin();
v3 = -1;
return;
}
if (vertexCount() == 2)
{
v1 = vertexBegin();
v2 = vertexNext(v1);
v3 = -1;
if (_layout_vertices[v1].morgan_code < _layout_vertices[v2].morgan_code)
{
v2 = vertexBegin();
v1 = vertexNext(v2);
}
return;
}
v1 = vertexBegin();
for (i = vertexNext(v1); i < vertexEnd(); i = vertexNext(i))
if (_layout_vertices[i].morgan_code > _layout_vertices[v1].morgan_code)
v1 = i;
const Vertex &vert = getVertex(v1);
v2 = vert.neiBegin();
for (i = vert.neiNext(v2); i < vert.neiEnd(); i = vert.neiNext(i))
if(_layout_vertices[vert.neiVertex(i)].morgan_code > _layout_vertices[vert.neiVertex(v2)].morgan_code)
v2 = i;
if (vert.degree() < 2)
{
v2 = vert.neiVertex(v2);
v3 = -1;
return;
}
v3 = vert.neiBegin();
if (v3 == v2)
v3 = vert.neiNext(v2);
for (i = vert.neiBegin(); i < vert.neiEnd(); i = vert.neiNext(i))
if (i != v2 && _layout_vertices[vert.neiVertex(i)].morgan_code > _layout_vertices[vert.neiVertex(v3)].morgan_code)
v3 = i;
v2 = vert.neiVertex(v2);
v3 = vert.neiVertex(v3);
}
// Scale and transform
void MoleculeLayoutGraph::_assignFinalCoordinates (float bond_length, const Array<Vec2f> &src_layout)
{
int i;
if (_n_fixed > 0)
{
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos.scale(bond_length);
return;
}
if (vertexCount() == 1)
{
getPos(vertexBegin()).set(0.f, 0.f);
return;
}
// flip molecule vertically if R1 is not above other R-groups
// flip molecule horizontally if R1 is not on the left
if (_molecule != 0 && _n_fixed == 0 && _molecule->countRSites() > 1)
{
QS_DEF(Array<int>, rgroup_list);
Vec2f r1_pos, highest_pos(0.f, -1000.f);
bool r1_exist = false;
float center_x = 0.f;
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
{
if (_molecule->isRSite(_layout_vertices[i].ext_idx))
{
_molecule->getAllowedRGroups(_layout_vertices[i].ext_idx, rgroup_list);
if (rgroup_list.size() == 1 && rgroup_list[0] == 1)
{
r1_pos = _layout_vertices[i].pos;
r1_exist = true;
} else if (_layout_vertices[i].pos.y > highest_pos.y)
{
highest_pos = _layout_vertices[i].pos;
}
}
center_x += _layout_vertices[i].pos.x;
}
center_x /= vertexCount();
if (r1_exist)
{
if (r1_pos.y < highest_pos.y)
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos.y *= -1;
if (r1_pos.x > center_x)
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos.x *= -1;
}
}
// 1. Choose scale ratio and first edge to match
float scale = bond_length, src_norm, norm;
int v1, v2, v3;
Vec2f p1, p2, p;
_getAnchor(v1, v2, v3);
p1.diff(src_layout[v2], src_layout[v1]);
p2.diff(getPos(v2), getPos(v1));
src_norm = p1.length();
norm = p2.length();
if (norm < 0.0001)
throw("too small edge");
// 2.1. If matching edge has zero length - just move to this point and scale
if (src_norm < 0.001)
{
p1 = src_layout[v1];
p2 = getPos(v1);
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
{
p.diff(getPos(i), p2);
p.scale(scale);
_layout_vertices[i].pos.sum(p1, p);
}
return;
}
// 2.2. If it has length from L/2 to 2L - scale by it, otherwise by L
if (src_norm >= bond_length / 2 && src_norm <= 2 * bond_length)
scale = src_norm / norm;
// 3. Move first vertex to (0,0)
p = getPos(v1);
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos.sub(p);
// 4. Rotate CCW on Alpha angle between (first, second) edge and (first, second) edge in source graph
float phi1, phi2, alpha, sina, cosa;
phi1 = p1.tiltAngle();
phi2 = p2.tiltAngle();
alpha = phi1 - phi2;
sina = sin(alpha);
cosa = cos(alpha);
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos.rotate(sina, cosa);
// 5. Scale
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos.scale(scale);
// 6. Match first vertices - shift by vector Pos(first)
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos.add(src_layout[v1]);
// 7. If needed turn around (first, second)
if (vertexCount() > 2)
{
float crit = 0.f;
// If v3 lays on the other side of line (first, second) - turn
p1 = getPos(v1);
p.diff(getPos(v2), p1);
if (v3 == -1)
{
for (v3 = vertexBegin(); v3 < vertexEnd(); v3 = vertexNext(v3))
{
if (fabs(p.x) < 0.001f)
crit = (src_layout[v3].x - p1.x) * (getPos(v3).x - p1.x);
else if (fabs(p.y) < 0.001f)
crit = (src_layout[v3].y - p1.y) * (getPos(v3).y - p1.y);
else
{
crit = (p.y * (src_layout[v3].x - p1.x) - p.x * (src_layout[v3].y - p1.y)) *
(p.y * (getPos(v3).x - p1.x) - p.x * (getPos(v3).y - p1.y));
}
if (fabs(crit) > 0.001)
break;
}
} else
crit = -1.0;
if (crit < 0 && v3 < vertexEnd())
{
// Move first vertex to (0,0)
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
_layout_vertices[i].pos.sub(p1);
// Turn by -phi1 and flip vertically
sina = -sin(phi1);
cosa = cos(phi1);
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
{
_layout_vertices[i].pos.rotate(sina, cosa);
_layout_vertices[i].pos.y *= -1;
}
// Turn by phi1 and translate back
sina = -sina;
for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
{
_layout_vertices[i].pos.rotate(sina, cosa);
_layout_vertices[i].pos.add(p1);
}
}
}
}
Jump to Line
Something went wrong with that request. Please try again.