match Namespace Reference

Functions

bool matchPrism (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, const Eigen::MatrixXd &refPoints, std::vector< int > *basal1, std::vector< int > *basal2, std::vector< double > *rmsdPerAtom, bool isPerfect=true)
 
bool matchUntetheredPrism (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, const Eigen::MatrixXd &refPoints, std::vector< int > *basal1, std::vector< int > *basal2, std::vector< double > *rmsdPerAtom)
 
bool matchPrismBlock (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, const Eigen::MatrixXd &refPoints, std::vector< int > *basal1, std::vector< int > *basal2, int *beginIndex)
 
int updatePerAtomRMSDRing (std::vector< int > basalRing, int startingIndex, std::vector< double > rmsdFromMatch, std::vector< double > *rmsdPerAtom)
 Update the per-particle RMSD for a prism block basal ring. More...
 
int updateRMSDRing (std::vector< int > basalRing, int startingIndex, double rmsdVal, std::vector< double > *rmsdPerAtom)
 

Function Documentation

◆ matchPrism()

bool match::matchPrism ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< std::vector< int >>  nList,
const Eigen::MatrixXd &  refPoints,
std::vector< int > *  basal1,
std::vector< int > *  basal2,
std::vector< double > *  rmsdPerAtom,
bool  isPerfect = true 
)

Shape-matching for a pair of polygon basal rings. Returns true if the pair of basal rings form a prism block.

Definition at line 17 of file shapeMatch.cpp.

21  {
22  //
23  int ringSize = (*basal1).size(); // Number of nodes in each basal ring
24  std::vector<int> matchedBasal1,
25  matchedBasal2; // Re-ordered basal rings 1 and 2
26  int dim = 3; // Number of dimensions
27  // Matrices for the point sets of the target basal rings
28  Eigen::MatrixXd basal1Set(ringSize,
29  dim); // Point set for the first basal ring
30  Eigen::MatrixXd basal2Set(ringSize,
31  dim); // Point set for the second basal ring
32  int startingIndex; // Index in the basal rings from which the reference point
33  // set is matched
34  double cutoffAngle = 15; // Tolerance for the angular distance between the
35  // quaternions of the basal rings
36  double angDist; // Calculated angular distance between the two basal rings
37  // Variables for the absolute orientation
38  std::vector<double> quat1, quat2; // quaternion rotation
39  double rmsd1, rmsd2; // least total RMSD
40  std::vector<double> rmsdList1,
41  rmsdList2; // List of RMSD per atom in the order fed in
42  double scale1, scale2; // Scale factor
43  // Deformed block classification
44  bool doAngleCriterion = false; // For deformed blocks
45  // -----------------------
46  // Getting the target Eigen vectors
47  // Get the re-ordered matched basal rings, ordered with respect to each other
48  pntToPnt::relOrderPrismBlock(yCloud, *basal1, *basal2, nList, &matchedBasal1,
49  &matchedBasal2);
50  // -----------------------
51  // Match the basal rings with a complete prism block, given the relatively
52  // ordered basal rings This actually only needs to be done for deformed prism
53  // blocks
54  bool blockMatch = match::matchPrismBlock(
55  yCloud, nList, refPoints, &matchedBasal1, &matchedBasal2, &startingIndex);
56  // -----------------------
57  // Check for deformed prisms
58  if (!isPerfect) {
59  if (!blockMatch) {
60  return blockMatch;
61  }
62  }
63  // If the deformed prism does not fulfil the shape matching criterion, then do
64  // not calculate the RMSD per atom
65  // -----------------------
66  // Section for RMSD per atom, based on basal ring matching
67 
68  basal1Set =
69  pntToPnt::fillPointSetPrismRing(yCloud, matchedBasal1, startingIndex);
70  // Fill up the point set for basal2
71  basal2Set =
72  pntToPnt::fillPointSetPrismRing(yCloud, matchedBasal2, startingIndex);
73  // Use Horn's algorithm to calculate the absolute orientation and RMSD etc.
74  absor::hornAbsOrientation(refPoints, basal1Set, &quat1, &rmsd1, &rmsdList1,
75  &scale1); // basal1
76  absor::hornAbsOrientation(refPoints, basal2Set, &quat2, &rmsd2, &rmsdList2,
77  &scale2); // basal2
78  // // ------------
79  // // Update the per-atom RMSD for each particle
80  // // Basal1
81  // match::updatePerAtomRMSDRing(matchedBasal1, startingIndex, rmsdList1,
82  // rmsdPerAtom);
83  // // Basal2
84  // match::updatePerAtomRMSDRing(matchedBasal2, startingIndex, rmsdList2,
85  // rmsdPerAtom);
86  // // ------------
87  // ------------
88  // Update the RMSD (obtained for each ring) for each particle
89  // Basal1
90  match::updateRMSDRing(matchedBasal1, startingIndex, rmsd1, rmsdPerAtom);
91  // Basal2
92  match::updateRMSDRing(matchedBasal2, startingIndex, rmsd2, rmsdPerAtom);
93  // ------------
94 
95  return true;
96 } // end of function
int hornAbsOrientation(const Eigen::MatrixXd &refPoints, const Eigen::MatrixXd &targetPoints, std::vector< double > *quat, double *rmsd, std::vector< double > *rmsdList, double *scale)
Get the absolute orientation using Horn's algorithm (with quaternions)
bool matchPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, const Eigen::MatrixXd &refPoints, std::vector< int > *basal1, std::vector< int > *basal2, int *beginIndex)
Definition: shapeMatch.cpp:247
int updateRMSDRing(std::vector< int > basalRing, int startingIndex, double rmsdVal, std::vector< double > *rmsdPerAtom)
Definition: shapeMatch.cpp:216
Eigen::MatrixXd fillPointSetPrismRing(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basalRing, int startingIndex)
int relOrderPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, std::vector< std::vector< int >> nList, std::vector< int > *outBasal1, std::vector< int > *outBasal2)

◆ matchPrismBlock()

bool match::matchPrismBlock ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< std::vector< int >>  nList,
const Eigen::MatrixXd &  refPoints,
std::vector< int > *  basal1,
std::vector< int > *  basal2,
int *  beginIndex 
)

Shape-matching for a pair of polygon basal rings, comparing with a complete prism block. Returns true if the pair of basal rings form a prism block.

Definition at line 247 of file shapeMatch.cpp.

250  {
251  //
252  int ringSize = (*basal1).size(); // Number of nodes in the basal rings
253  bool isMatch; // Qualifier for whether the prism block matches or not
254  int dim = 3; // Number of dimensions
255  int nop = 2 * ringSize; // Number of particles in the prism block
256  Eigen::MatrixXd refPrismBlock(
257  nop, dim); // eigen matrix for the reference prism block
258  Eigen::MatrixXd targetPrismBlock(
259  nop, dim); // eigen matrix for the target prism block
260  int startingIndex; // Index from which the basal rings will be ordered
261  // (permutations)
262  // variables for shape-matching
263  std::vector<double> quat; // Quaternion for the rotation
264  double rmsd; // RMSD value for the entire shape
265  std::vector<double> rmsdList; // RMSD value for each point
266  double scale; // Scale factor
267  // ----------------------------------------------------
268  // Get the reference prism block
269  refPrismBlock =
270  pntToPnt::createPrismBlock(yCloud, refPoints, ringSize, *basal1, *basal2);
271  // ----------------------------------------------------
272 
273  // Loop through possible point-to-point correspondences
274  if (ringSize % 2 == 0 || ringSize == 3) {
275  startingIndex = 0;
276  // Fill up the point set for the target prism block
277  targetPrismBlock =
278  pntToPnt::fillPointSetPrismBlock(yCloud, *basal1, *basal2, 0);
279  // Shape-matching
280  absor::hornAbsOrientation(refPrismBlock, targetPrismBlock, &quat, &rmsd,
281  &rmsdList,
282  &scale); // basal2
283  } // even or if there are 3 nodes
284  else {
285  // Define the vector, RMSD etc:
286  std::vector<double> currentQuat; // quaternion rotation
287  double currentRmsd; // least total RMSD
289  currentRmsdList; // List of RMSD per atom in the order fed in
290  double currentScale;
291  // Loop through all possible startingIndex
292  for (int i = 0; i < ringSize; i++) {
293  //
294  // Fill up the point set for basal1
295  targetPrismBlock =
296  pntToPnt::fillPointSetPrismBlock(yCloud, *basal1, *basal2, i);
297  // Use Horn's algorithm to calculate the absolute orientation and RMSD
298  // etc. for basal1
299  absor::hornAbsOrientation(refPrismBlock, targetPrismBlock, &currentQuat,
300  &currentRmsd, &currentRmsdList, &currentScale);
301  // Comparison to get the least RMSD for the correct mapping
302  if (i == 0) {
303  // Init
304  quat = currentQuat;
305  rmsd = currentRmsd;
306  rmsdList = currentRmsdList;
307  scale = currentScale;
308  startingIndex = i;
309  } // init
310  else {
311  // Check to see if the calculated RMSD is less than the RMSD already
312  // saved
313  if (currentRmsd < rmsd) {
314  quat = currentQuat;
315  rmsd = currentRmsd;
316  rmsdList = currentRmsdList;
317  scale = currentScale;
318  startingIndex = i;
319  } // end of check to see if the current RMSD is smaller
320  } // end of comparison and filling
321  } // end of loop through startingIndex
322  } // ringSize is odd, so every point must be tried
323 
324  // Update the index from which matching is started
325  *beginIndex = startingIndex;
326 
327  // // TEST DELETE BLOCK LATER
328  // if (ringSize == 6) {
329  // std::fstream rmsdFile;
330  // rmsdFile.open("../runOne/topoINT/rmsd6.dat",
331  // std::fstream::in | std::fstream::out | std::fstream::app);
332  // rmsdFile << rmsd << "\n";
333  // rmsdFile.close();
334  // } else if (ringSize == 4) {
335  // std::fstream rmsdFile;
336  // rmsdFile.open("../runOne/topoINT/rmsd4.dat",
337  // std::fstream::in | std::fstream::out | std::fstream::app);
338  // rmsdFile << rmsd << "\n";
339  // rmsdFile.close();
340  // }
341  // // END OF BLOCK TO BE DELETED
342 
343  // Condition for shape-matching
344  if (rmsd <= 6) {
345  isMatch = true; // Is a prism block
346  } // within RMSD range
347  else {
348  isMatch = false;
349  }
350 
351  // Return
352  return isMatch;
353 }
Eigen::MatrixXd createPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, int ringSize, std::vector< int > basal1, std::vector< int > basal2)
Eigen::MatrixXd fillPointSetPrismBlock(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< int > basal1, std::vector< int > basal2, int startingIndex)

◆ matchUntetheredPrism()

bool match::matchUntetheredPrism ( molSys::PointCloud< molSys::Point< double >, double > *  yCloud,
std::vector< std::vector< int >>  nList,
const Eigen::MatrixXd &  refPoints,
std::vector< int > *  basal1,
std::vector< int > *  basal2,
std::vector< double > *  rmsdPerAtom 
)

Shape-matching for a pair of polygon basal rings. Returns true if the pair of basal rings form a prism block.

For the pentagonal nanochannels in amorphous ice.

These must be aligned.

Shape-matching for a pair of polygon basal rings. Returns true if the pair of basal rings form a prism block.

Definition at line 106 of file shapeMatch.cpp.

110  {
111  //
112  int ringSize = (*basal1).size(); // Number of nodes in each basal ring
113  std::vector<int> matchedBasal1,
114  matchedBasal2; // Re-ordered basal rings 1 and 2
115  int dim = 3; // Number of dimensions
116  // Matrices for the point sets of the target basal rings
117  Eigen::MatrixXd basal1Set(ringSize,
118  dim); // Point set for the first basal ring
119  Eigen::MatrixXd basal2Set(ringSize,
120  dim); // Point set for the second basal ring
121  int startingIndex; // Index in the basal rings from which the reference point
122  // set is matched
123  double angDist; // Calculated angular distance between the two basal rings
124  // Variables for the absolute orientation
125  std::vector<double> quat1, quat2; // quaternion rotation
126  double rmsd1, rmsd2; // least total RMSD
127  std::vector<double> rmsdList1,
128  rmsdList2; // List of RMSD per atom in the order fed in
129  double scale1, scale2; // Scale factor
130  // Deformed block classification
131  bool doAngleCriterion = false; // For deformed blocks
132 
133  // Getting the target Eigen vectors
134  // Get the re-ordered matched basal rings, ordered with respect to each other
135  pntToPnt::relOrderPrismBlock(yCloud, *basal1, *basal2, &matchedBasal1,
136  &matchedBasal2);
137  // -----------------------
138  // Match the basal rings with a complete prism block, given the relatively
139  // ordered basal rings This actually only needs to be done for deformed prism
140  // blocks
141  bool blockMatch = match::matchPrismBlock(
142  yCloud, nList, refPoints, &matchedBasal1, &matchedBasal2, &startingIndex);
143  // -----------------------
144  // Check to see if the prism matches the reference prism block
145  if (!blockMatch) {
146  return blockMatch;
147  }
148 
149  // If the deformed prism does not fulfil the shape matching criterion, then do
150  // not calculate the RMSD per atom
151  // -----------------------
152  // Section for RMSD per atom, based on basal ring matching
153 
154  basal1Set =
155  pntToPnt::fillPointSetPrismRing(yCloud, matchedBasal1, startingIndex);
156  // Fill up the point set for basal2
157  basal2Set =
158  pntToPnt::fillPointSetPrismRing(yCloud, matchedBasal2, startingIndex);
159  // Use Horn's algorithm to calculate the absolute orientation and RMSD etc.
160  absor::hornAbsOrientation(refPoints, basal1Set, &quat1, &rmsd1, &rmsdList1,
161  &scale1); // basal1
162  absor::hornAbsOrientation(refPoints, basal2Set, &quat2, &rmsd2, &rmsdList2,
163  &scale2); // basal2
164  // // Calculate the angular distance between basal1 and basal2
165  // angDist = gen::angDistDegQuaternions(quat1, quat2);
166  // // Check if the shapes are aligned
167  // if (angDist > cutoffAngle) {
168  // return false;
169  // } // If not aligned, it is not a prism block
170  // ------------
171  // Update the RMSD (obtained for each ring) for each particle
172  // Basal1
173  match::updateRMSDRing(matchedBasal1, startingIndex, rmsd1, rmsdPerAtom);
174  // Basal2
175  match::updateRMSDRing(matchedBasal2, startingIndex, rmsd2, rmsdPerAtom);
176  // ------------
177 
178  return true;
179 } // end of function

◆ updatePerAtomRMSDRing()

int match::updatePerAtomRMSDRing ( std::vector< int >  basalRing,
int  startingIndex,
std::vector< double >  rmsdFromMatch,
std::vector< double > *  rmsdPerAtom 
)

Update the per-particle RMSD for a prism block basal ring.

Definition at line 182 of file shapeMatch.cpp.

184  {
185  //
186  int atomIndex; // Atom particle index, in rmsdPerAtom
187  int ringSize = basalRing.size(); // Size of the basal ring
188  int index; // Corresponding index in the basal ring
189  double iRMSD; // Per-particle RMSD
190 
191  // -----------------
192  // Update the RMSD per particle
193  for (int i = 0; i < ringSize; i++) {
194  index = i + startingIndex; // Corresponding index in the basal ring
195  // Wrap-around
196  if (index >= ringSize) {
197  index -= ringSize;
198  } // end of wrap-around
199  //
200  // Get the atom index
201  atomIndex = basalRing[index];
202  // Get and update the RMSD per particle
203  iRMSD = rmsdFromMatch[i];
204  if ((*rmsdPerAtom)[atomIndex] == -1) {
205  (*rmsdPerAtom)[atomIndex] = iRMSD;
206  } // end of update of RMSD
207 
208  } // end of updating the RMSD per particle
209  // -----------------
210  // finito
211  return 0;
212 } // end of function
T size(T... args)

◆ updateRMSDRing()

int match::updateRMSDRing ( std::vector< int >  basalRing,
int  startingIndex,
double  rmsdVal,
std::vector< double > *  rmsdPerAtom 
)

Update the RMSD of each particle in a prism block basal ring with the RMSD of the ring.

Update the RMSD for each particle with the RMSD of each ring for a prism block basal ring.

Definition at line 216 of file shapeMatch.cpp.

217  {
218  //
219  int atomIndex; // Atom particle index, in rmsdPerAtom
220  int ringSize = basalRing.size(); // Size of the basal ring
221  int index; // Corresponding index in the basal ring
222 
223  // -----------------
224  // Update the RMSD per particle
225  for (int i = 0; i < ringSize; i++) {
226  index = i + startingIndex; // Corresponding index in the basal ring
227  // Wrap-around
228  if (index >= ringSize) {
229  index -= ringSize;
230  } // end of wrap-around
231  //
232  // Get the atom index
233  atomIndex = basalRing[index];
234  // The RMSD value to update with is rmsdVal
235  if ((*rmsdPerAtom)[atomIndex] == -1) {
236  (*rmsdPerAtom)[atomIndex] = rmsdVal;
237  } // end of update of RMSD
238 
239  } // end of updating the RMSD for each particle
240  // -----------------
241  // finito
242  return 0;
243 } // end of function