source: src/Fragmentation/Exporters/SphericalPointDistribution.cpp@ 1ae9aa

Last change on this file since 1ae9aa was 1ae9aa, checked in by Frederik Heber <heber@…>, 11 years ago

Extracted joinPoints() function to make it accessible to tests.

  • added unit test.
  • moved some function definitions around and added documentation.
  • Property mode set to 100644
File size: 26.9 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2014 Frederik Heber. All rights reserved.
5 *
6 *
7 * This file is part of MoleCuilder.
8 *
9 * MoleCuilder is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * MoleCuilder is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
21 */
22
23/*
24 * SphericalPointDistribution.cpp
25 *
26 * Created on: May 30, 2014
27 * Author: heber
28 */
29
30// include config.h
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35#include "CodePatterns/MemDebug.hpp"
36
37#include "SphericalPointDistribution.hpp"
38
39#include "CodePatterns/Assert.hpp"
40#include "CodePatterns/IteratorAdaptors.hpp"
41#include "CodePatterns/Log.hpp"
42#include "CodePatterns/toString.hpp"
43
44#include <algorithm>
45#include <boost/math/quaternion.hpp>
46#include <cmath>
47#include <functional>
48#include <iterator>
49#include <limits>
50#include <list>
51#include <vector>
52#include <map>
53
54#include "LinearAlgebra/Line.hpp"
55#include "LinearAlgebra/Plane.hpp"
56#include "LinearAlgebra/RealSpaceMatrix.hpp"
57#include "LinearAlgebra/Vector.hpp"
58
59// static entities
60const double SphericalPointDistribution::warn_amplitude = 1e-2;
61
62typedef std::vector<double> DistanceArray_t;
63
64// class generator: taken from www.cplusplus.com example std::generate
65struct c_unique {
66 int current;
67 c_unique() {current=0;}
68 int operator()() {return current++;}
69} UniqueNumber;
70
71inline
72DistanceArray_t calculatePairwiseDistances(
73 const std::vector<Vector> &_points,
74 const SphericalPointDistribution::IndexList_t &_indices
75 )
76{
77 DistanceArray_t result;
78 for (SphericalPointDistribution::IndexList_t::const_iterator firstiter = _indices.begin();
79 firstiter != _indices.end(); ++firstiter) {
80 for (SphericalPointDistribution::IndexList_t::const_iterator seconditer = firstiter;
81 seconditer != _indices.end(); ++seconditer) {
82 if (firstiter == seconditer)
83 continue;
84 const double distance = (_points[*firstiter] - _points[*seconditer]).NormSquared();
85 result.push_back(distance);
86 }
87 }
88 return result;
89}
90
91/** Calculate the center of a given set of points in \a _positions but only
92 * for those indicated by \a _indices.
93 *
94 */
95inline
96Vector calculateCenter(
97 const SphericalPointDistribution::VectorArray_t &_positions,
98 const SphericalPointDistribution::IndexList_t &_indices)
99{
100 Vector Center;
101 Center.Zero();
102 for (SphericalPointDistribution::IndexList_t::const_iterator iter = _indices.begin();
103 iter != _indices.end(); ++iter)
104 Center += _positions[*iter];
105 if (!_indices.empty())
106 Center *= 1./(double)_indices.size();
107
108 return Center;
109}
110
111/** Decides by an orthonormal third vector whether the sign of the rotation
112 * angle should be negative or positive.
113 *
114 * \return -1 or 1
115 */
116inline
117double determineSignOfRotation(
118 const Vector &_oldPosition,
119 const Vector &_newPosition,
120 const Vector &_RotationAxis
121 )
122{
123 Vector dreiBein(_oldPosition);
124 dreiBein.VectorProduct(_RotationAxis);
125 dreiBein.Normalize();
126 const double sign =
127 (dreiBein.ScalarProduct(_newPosition) < 0.) ? -1. : +1.;
128 LOG(6, "DEBUG: oldCenter on plane is " << _oldPosition
129 << ", newCenter in plane is " << _newPosition
130 << ", and dreiBein is " << dreiBein);
131 return sign;
132}
133
134/** Convenience function to recalculate old and new center for determining plane
135 * rotation.
136 */
137inline
138void calculateOldAndNewCenters(
139 Vector &_oldCenter,
140 Vector &_newCenter,
141 const SphericalPointDistribution::VectorArray_t &_referencepositions,
142 const SphericalPointDistribution::VectorArray_t &_currentpositions,
143 const SphericalPointDistribution::IndexList_t &_bestmatching)
144{
145 const size_t NumberIds = std::min(_bestmatching.size(), (size_t)3);
146 SphericalPointDistribution::IndexList_t continuousIds(NumberIds, -1);
147 std::generate(continuousIds.begin(), continuousIds.end(), UniqueNumber);
148 _oldCenter = calculateCenter(_referencepositions, continuousIds);
149 // C++11 defines a copy_n function ...
150 SphericalPointDistribution::IndexList_t::const_iterator enditer = _bestmatching.begin();
151 std::advance(enditer, NumberIds);
152 SphericalPointDistribution::IndexList_t firstbestmatchingIds(NumberIds, -1);
153 std::copy(_bestmatching.begin(), enditer, firstbestmatchingIds.begin());
154 _newCenter = calculateCenter( _currentpositions, firstbestmatchingIds);
155}
156/** Returns squared L2 error of the given \a _Matching.
157 *
158 * We compare the pair-wise distances of each associated matching
159 * and check whether these distances each match between \a _old and
160 * \a _new.
161 *
162 * \param _old first set of points (fewer or equal to \a _new)
163 * \param _new second set of points
164 * \param _Matching matching between the two sets
165 * \return pair with L1 and squared L2 error
166 */
167std::pair<double, double> SphericalPointDistribution::calculateErrorOfMatching(
168 const std::vector<Vector> &_old,
169 const std::vector<Vector> &_new,
170 const IndexList_t &_Matching)
171{
172 std::pair<double, double> errors( std::make_pair( 0., 0. ) );
173
174 if (_Matching.size() > 1) {
175 LOG(3, "INFO: Matching is " << _Matching);
176
177 // calculate all pair-wise distances
178 IndexList_t keys(_Matching.size());
179 std::generate (keys.begin(), keys.end(), UniqueNumber);
180 const DistanceArray_t firstdistances = calculatePairwiseDistances(_old, keys);
181 const DistanceArray_t seconddistances = calculatePairwiseDistances(_new, _Matching);
182
183 ASSERT( firstdistances.size() == seconddistances.size(),
184 "calculateL2ErrorOfMatching() - mismatch in pair-wise distance array sizes.");
185 DistanceArray_t::const_iterator firstiter = firstdistances.begin();
186 DistanceArray_t::const_iterator seconditer = seconddistances.begin();
187 for (;(firstiter != firstdistances.end()) && (seconditer != seconddistances.end());
188 ++firstiter, ++seconditer) {
189 const double gap = *firstiter - *seconditer;
190 // L1 error
191 if (errors.first < gap)
192 errors.first = gap;
193 // L2 error
194 errors.second += gap*gap;
195 }
196 } else
197 ELOG(3, "calculateErrorOfMatching() - Given matching's size is less than 2.");
198 LOG(3, "INFO: Resulting errors for matching (L1, L2): "
199 << errors.first << "," << errors.second << ".");
200
201 return errors;
202}
203
204SphericalPointDistribution::Polygon_t SphericalPointDistribution::removeMatchingPoints(
205 const VectorArray_t &_points,
206 const IndexList_t &_matchingindices
207 )
208{
209 SphericalPointDistribution::Polygon_t remainingpoints;
210 IndexArray_t indices(_matchingindices.begin(), _matchingindices.end());
211 std::sort(indices.begin(), indices.end());
212 LOG(4, "DEBUG: sorted matching is " << indices);
213 IndexArray_t remainingindices(_points.size(), -1);
214 std::generate(remainingindices.begin(), remainingindices.end(), UniqueNumber);
215 IndexArray_t::iterator remainiter = std::set_difference(
216 remainingindices.begin(), remainingindices.end(),
217 indices.begin(), indices.end(),
218 remainingindices.begin());
219 remainingindices.erase(remainiter, remainingindices.end());
220 LOG(4, "DEBUG: remaining indices are " << remainingindices);
221 for (IndexArray_t::const_iterator iter = remainingindices.begin();
222 iter != remainingindices.end(); ++iter) {
223 remainingpoints.push_back(_points[*iter]);
224 }
225
226 return remainingpoints;
227}
228
229/** Recursive function to go through all possible matchings.
230 *
231 * \param _MCS structure holding global information to the recursion
232 * \param _matching current matching being build up
233 * \param _indices contains still available indices
234 * \param _matchingsize
235 */
236void SphericalPointDistribution::recurseMatchings(
237 MatchingControlStructure &_MCS,
238 IndexList_t &_matching,
239 IndexList_t _indices,
240 unsigned int _matchingsize)
241{
242 LOG(4, "DEBUG: Recursing with current matching " << _matching
243 << ", remaining indices " << _indices
244 << ", and sought size " << _matching.size()+_matchingsize);
245 //!> threshold for L1 error below which matching is immediately acceptable
246 const double L1THRESHOLD = 1e-2;
247 if (!_MCS.foundflag) {
248 LOG(4, "DEBUG: Current matching has size " << _matching.size() << ", places left " << _matchingsize);
249 if (_matchingsize > 0) {
250 // go through all indices
251 for (IndexList_t::iterator iter = _indices.begin();
252 (iter != _indices.end()) && (!_MCS.foundflag);) {
253 // add index to matching
254 _matching.push_back(*iter);
255 LOG(5, "DEBUG: Adding " << *iter << " to matching.");
256 // remove index but keep iterator to position (is the next to erase element)
257 IndexList_t::iterator backupiter = _indices.erase(iter);
258 // recurse with decreased _matchingsize
259 recurseMatchings(_MCS, _matching, _indices, _matchingsize-1);
260 // re-add chosen index and reset index to new position
261 _indices.insert(backupiter, _matching.back());
262 iter = backupiter;
263 // remove index from _matching to make space for the next one
264 _matching.pop_back();
265 }
266 // gone through all indices then exit recursion
267 if (_matching.empty())
268 _MCS.foundflag = true;
269 } else {
270 LOG(3, "INFO: Found matching " << _matching);
271 // calculate errors
272 std::pair<double, double> errors = calculateErrorOfMatching(
273 _MCS.oldpoints, _MCS.newpoints, _matching);
274 if (errors.first < L1THRESHOLD) {
275 _MCS.bestmatching = _matching;
276 _MCS.foundflag = true;
277 } else if (_MCS.bestL2 > errors.second) {
278 _MCS.bestmatching = _matching;
279 _MCS.bestL2 = errors.second;
280 }
281 }
282 }
283}
284
285/** Finds combinatorially the best matching between points in \a _polygon
286 * and \a _newpolygon.
287 *
288 * We find the matching with the smallest L2 error, where we break when we stumble
289 * upon a matching with zero error.
290 *
291 * As points in \a _polygon may be have a weight greater 1 we have to match it to
292 * multiple points in \a _newpolygon. Eventually, these multiple points are combined
293 * for their center of weight, which is the only thing follow-up function see of
294 * these multiple points. Hence, we actually modify \a _newpolygon in the process
295 * such that the returned IndexList_t indicates a bijective mapping in the end.
296 *
297 * \sa recurseMatchings() for going through all matchings
298 *
299 * \param _polygon here, we have indices 0,1,2,...
300 * \param _newpolygon and here we need to find the correct indices
301 * \return list of indices: first in \a _polygon goes to first index for \a _newpolygon
302 */
303SphericalPointDistribution::IndexList_t SphericalPointDistribution::findBestMatching(
304 const WeightedPolygon_t &_polygon,
305 Polygon_t &_newpolygon
306 )
307{
308 MatchingControlStructure MCS;
309 MCS.foundflag = false;
310 MCS.bestL2 = std::numeric_limits<double>::max();
311 for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
312 iter != _polygon.end(); ++iter)
313 MCS.oldpoints.push_back(iter->first);
314 MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
315
316 // search for bestmatching combinatorially
317 {
318 // translate polygon into vector to enable index addressing
319 IndexList_t indices(_newpolygon.size());
320 std::generate(indices.begin(), indices.end(), UniqueNumber);
321 IndexList_t matching;
322
323 // walk through all matchings
324 const unsigned int matchingsize = _polygon.size();
325 ASSERT( matchingsize <= indices.size(),
326 "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new points to choose for matching to old ones.");
327 recurseMatchings(MCS, matching, indices, matchingsize);
328 }
329
330 // combine multiple points and create simple IndexList from IndexTupleList
331 IndexTupleList_t bestmatching;
332 for (IndexList_t::const_iterator iter = MCS.bestmatching.begin();
333 iter != MCS.bestmatching.end(); ++iter)
334 bestmatching.push_back(IndexList_t(1, *iter));
335 const SphericalPointDistribution::IndexList_t IndexList =
336 joinPoints(_newpolygon, MCS.newpoints, bestmatching);
337
338 return IndexList;
339}
340
341SphericalPointDistribution::IndexList_t SphericalPointDistribution::joinPoints(
342 Polygon_t &_newpolygon,
343 const VectorArray_t &_newpoints,
344 const IndexTupleList_t &_bestmatching
345 )
346{
347 // combine all multiple points
348 IndexList_t IndexList;
349 IndexArray_t removalpoints;
350 unsigned int UniqueIndex = _newpolygon.size(); // all indices up to size are used right now
351 VectorArray_t newCenters;
352 newCenters.reserve(_bestmatching.size());
353 for (IndexTupleList_t::const_iterator tupleiter = _bestmatching.begin();
354 tupleiter != _bestmatching.end(); ++tupleiter) {
355 ASSERT (tupleiter->size() > 0,
356 "findBestMatching() - encountered tuple in bestmatching with size 0.");
357 if (tupleiter->size() == 1) {
358 // add point and index
359 IndexList.push_back(*tupleiter->begin());
360 } else {
361 // combine into weighted and normalized center
362 Vector Center = calculateCenter(_newpoints, *tupleiter);
363 Center.Normalize();
364 _newpolygon.push_back(Center);
365 LOG(5, "DEBUG: Combining " << tupleiter->size() << "points to weighted center "
366 << Center << " with new index " << UniqueIndex);
367 // mark for removal
368 removalpoints.insert(removalpoints.end(), tupleiter->begin(), tupleiter->end());
369 // add new index
370 IndexList.push_back(UniqueIndex++);
371 }
372 }
373 // IndexList is now our new bestmatching (that is bijective)
374 LOG(4, "DEBUG: Our new bijective IndexList reads as " << IndexList);
375
376 // modifying _newpolygon: remove all points in removalpoints, add those in newCenters
377 Polygon_t allnewpoints = _newpolygon;
378 {
379 _newpolygon.clear();
380 std::sort(removalpoints.begin(), removalpoints.end());
381 size_t i = 0;
382 IndexArray_t::const_iterator removeiter = removalpoints.begin();
383 for (Polygon_t::iterator iter = allnewpoints.begin();
384 iter != allnewpoints.end(); ++iter, ++i) {
385 if ((removeiter != removalpoints.end()) && (i == *removeiter)) {
386 // don't add, go to next remove index
387 ++removeiter;
388 } else {
389 // otherwise add points
390 _newpolygon.push_back(*iter);
391 }
392 }
393 }
394 LOG(4, "DEBUG: The polygon with recentered points removed is " << _newpolygon);
395
396 // map IndexList to new shrinked _newpolygon
397 typedef std::set<unsigned int> IndexSet_t;
398 IndexSet_t SortedIndexList(IndexList.begin(), IndexList.end());
399 IndexList.clear();
400 {
401 size_t offset = 0;
402 IndexSet_t::const_iterator listiter = SortedIndexList.begin();
403 IndexArray_t::const_iterator removeiter = removalpoints.begin();
404 for (size_t i = 0; i < allnewpoints.size(); ++i) {
405 if ((removeiter != removalpoints.end()) && (i == *removeiter)) {
406 ++offset;
407 ++removeiter;
408 } else if ((listiter != SortedIndexList.end()) && (i == *listiter)) {
409 IndexList.push_back(*listiter - offset);
410 ++listiter;
411 }
412 }
413 }
414 LOG(4, "DEBUG: Our new bijective IndexList corrected for removed points reads as "
415 << IndexList);
416
417 return IndexList;
418}
419
420SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPlaneAligningRotation(
421 const VectorArray_t &_referencepositions,
422 const VectorArray_t &_currentpositions,
423 const IndexList_t &_bestmatching
424 )
425{
426 bool dontcheck = false;
427 // initialize to no rotation
428 Rotation_t Rotation;
429 Rotation.first.Zero();
430 Rotation.first[0] = 1.;
431 Rotation.second = 0.;
432
433 // calculate center of triangle/line/point consisting of first points of matching
434 Vector oldCenter;
435 Vector newCenter;
436 calculateOldAndNewCenters(
437 oldCenter, newCenter,
438 _referencepositions, _currentpositions, _bestmatching);
439
440 if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
441 LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
442 oldCenter.Normalize();
443 newCenter.Normalize();
444 if (!oldCenter.IsEqualTo(newCenter)) {
445 // calculate rotation axis and angle
446 Rotation.first = oldCenter;
447 Rotation.first.VectorProduct(newCenter);
448 Rotation.second = oldCenter.Angle(newCenter); // /(M_PI/2.);
449 } else {
450 // no rotation required anymore
451 }
452 } else {
453 LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
454 if ((oldCenter.IsZero()) && (newCenter.IsZero())) {
455 // either oldCenter or newCenter (or both) is directly at origin
456 if (_bestmatching.size() == 2) {
457 // line case
458 Vector oldPosition = _currentpositions[*_bestmatching.begin()];
459 Vector newPosition = _referencepositions[0];
460 // check whether we need to rotate at all
461 if (!oldPosition.IsEqualTo(newPosition)) {
462 Rotation.first = oldPosition;
463 Rotation.first.VectorProduct(newPosition);
464 // orientation will fix the sign here eventually
465 Rotation.second = oldPosition.Angle(newPosition);
466 } else {
467 // no rotation required anymore
468 }
469 } else {
470 // triangle case
471 // both triangles/planes have same center, hence get axis by
472 // VectorProduct of Normals
473 Plane newplane(_referencepositions[0], _referencepositions[1], _referencepositions[2]);
474 VectorArray_t vectors;
475 for (IndexList_t::const_iterator iter = _bestmatching.begin();
476 iter != _bestmatching.end(); ++iter)
477 vectors.push_back(_currentpositions[*iter]);
478 Plane oldplane(vectors[0], vectors[1], vectors[2]);
479 Vector oldPosition = oldplane.getNormal();
480 Vector newPosition = newplane.getNormal();
481 // check whether we need to rotate at all
482 if (!oldPosition.IsEqualTo(newPosition)) {
483 Rotation.first = oldPosition;
484 Rotation.first.VectorProduct(newPosition);
485 Rotation.first.Normalize();
486
487 // construct reference vector to determine direction of rotation
488 const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
489 Rotation.second = sign * oldPosition.Angle(newPosition);
490 LOG(5, "DEBUG: Rotating plane normals by " << Rotation.second
491 << " around axis " << Rotation.first);
492 } else {
493 // else do nothing
494 }
495 }
496 } else {
497 // TODO: we can't do anything here, but this case needs to be dealt with when
498 // we have no ideal geometries anymore
499 if ((oldCenter-newCenter).Norm() > warn_amplitude)
500 ELOG(2, "oldCenter is " << oldCenter << ", yet newCenter is " << newCenter);
501 // else they are considered close enough
502 dontcheck = true;
503 }
504 }
505
506#ifndef NDEBUG
507 // check: rotation brings newCenter onto oldCenter position
508 if (!dontcheck) {
509 Line Axis(zeroVec, Rotation.first);
510 Vector test = Axis.rotateVector(newCenter, Rotation.second);
511 LOG(4, "CHECK: rotated newCenter is " << test
512 << ", oldCenter is " << oldCenter);
513 ASSERT( (test - oldCenter).NormSquared() < std::numeric_limits<double>::epsilon()*1e4,
514 "matchSphericalPointDistributions() - rotation does not work as expected by "
515 +toString((test - oldCenter).NormSquared())+".");
516 }
517#endif
518
519 return Rotation;
520}
521
522SphericalPointDistribution::Rotation_t SphericalPointDistribution::findPointAligningRotation(
523 const VectorArray_t &remainingold,
524 const VectorArray_t &remainingnew,
525 const IndexList_t &_bestmatching)
526{
527 // initialize rotation to zero
528 Rotation_t Rotation;
529 Rotation.first.Zero();
530 Rotation.first[0] = 1.;
531 Rotation.second = 0.;
532
533 // recalculate center
534 Vector oldCenter;
535 Vector newCenter;
536 calculateOldAndNewCenters(
537 oldCenter, newCenter,
538 remainingold, remainingnew, _bestmatching);
539
540 Vector oldPosition = remainingnew[*_bestmatching.begin()];
541 Vector newPosition = remainingold[0];
542 LOG(6, "DEBUG: oldPosition is " << oldPosition << " and newPosition is " << newPosition);
543 if (!oldPosition.IsEqualTo(newPosition)) {
544 if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
545 oldCenter.Normalize(); // note weighted sum of normalized weight is not normalized
546 Rotation.first = oldCenter;
547 LOG(6, "DEBUG: Picking normalized oldCenter as Rotation.first " << oldCenter);
548 oldPosition.ProjectOntoPlane(Rotation.first);
549 newPosition.ProjectOntoPlane(Rotation.first);
550 LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
551 } else {
552 if (_bestmatching.size() == 2) {
553 // line situation
554 try {
555 Plane oldplane(oldPosition, oldCenter, newPosition);
556 Rotation.first = oldplane.getNormal();
557 LOG(6, "DEBUG: Plane is " << oldplane << " and normal is " << Rotation.first);
558 } catch (LinearDependenceException &e) {
559 LOG(6, "DEBUG: Vectors defining plane are linearly dependent.");
560 // oldPosition and newPosition are on a line, just flip when not equal
561 if (!oldPosition.IsEqualTo(newPosition)) {
562 Rotation.first.Zero();
563 Rotation.first.GetOneNormalVector(oldPosition);
564 LOG(6, "DEBUG: For flipping we use Rotation.first " << Rotation.first);
565 assert( Rotation.first.ScalarProduct(oldPosition) < std::numeric_limits<double>::epsilon()*1e4);
566 // Rotation.second = M_PI;
567 } else {
568 LOG(6, "DEBUG: oldPosition and newPosition are equivalent.");
569 }
570 }
571 } else {
572 // triangle situation
573 Plane oldplane(remainingold[0], remainingold[1], remainingold[2]);
574 Rotation.first = oldplane.getNormal();
575 LOG(6, "DEBUG: oldPlane is " << oldplane << " and normal is " << Rotation.first);
576 oldPosition.ProjectOntoPlane(Rotation.first);
577 LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
578 }
579 }
580 // construct reference vector to determine direction of rotation
581 const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
582 Rotation.second = sign * oldPosition.Angle(newPosition);
583 } else {
584 LOG(6, "DEBUG: oldPosition and newPosition are equivalent, hence no orientating rotation.");
585 }
586
587 return Rotation;
588}
589
590
591SphericalPointDistribution::Polygon_t
592SphericalPointDistribution::matchSphericalPointDistributions(
593 const SphericalPointDistribution::WeightedPolygon_t &_polygon,
594 SphericalPointDistribution::Polygon_t &_newpolygon
595 )
596{
597 SphericalPointDistribution::Polygon_t remainingpoints;
598 VectorArray_t remainingold;
599 for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
600 iter != _polygon.end(); ++iter)
601 remainingold.push_back(iter->first);
602 VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
603 LOG(2, "INFO: Matching old polygon " << _polygon
604 << " with new polygon " << _newpolygon);
605
606 if (_polygon.size() == _newpolygon.size()) {
607 // same number of points desired as are present? Do nothing
608 LOG(2, "INFO: There are no vacant points to return.");
609 return remainingpoints;
610 }
611
612 if (_polygon.size() > 0) {
613 IndexList_t bestmatching = findBestMatching(_polygon, _newpolygon);
614 LOG(2, "INFO: Best matching is " << bestmatching);
615
616 // determine rotation angles to align the two point distributions with
617 // respect to bestmatching:
618 // we use the center between the three first matching points
619 /// the first rotation brings these two centers to coincide
620 VectorArray_t rotated_newpolygon = remainingnew;
621 {
622 Rotation_t Rotation = findPlaneAligningRotation(
623 remainingold,
624 remainingnew,
625 bestmatching);
626 LOG(5, "DEBUG: Rotating coordinate system by " << Rotation.second
627 << " around axis " << Rotation.first);
628 Line Axis(zeroVec, Rotation.first);
629
630 // apply rotation angle to bring newCenter to oldCenter
631 for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
632 iter != rotated_newpolygon.end(); ++iter) {
633 Vector &current = *iter;
634 LOG(6, "DEBUG: Original point is " << current);
635 current = Axis.rotateVector(current, Rotation.second);
636 LOG(6, "DEBUG: Rotated point is " << current);
637 }
638
639#ifndef NDEBUG
640 // check: rotated "newCenter" should now equal oldCenter
641 {
642 Vector oldCenter;
643 Vector rotatednewCenter;
644 calculateOldAndNewCenters(
645 oldCenter, rotatednewCenter,
646 remainingold, rotated_newpolygon, bestmatching);
647 // NOTE: Center must not necessarily lie on the sphere with norm 1, hence, we
648 // have to normalize it just as before, as oldCenter and newCenter lengths may differ.
649 if ((!oldCenter.IsZero()) && (!rotatednewCenter.IsZero())) {
650 oldCenter.Normalize();
651 rotatednewCenter.Normalize();
652 LOG(4, "CHECK: rotatednewCenter is " << rotatednewCenter
653 << ", oldCenter is " << oldCenter);
654 ASSERT( (rotatednewCenter - oldCenter).NormSquared() < std::numeric_limits<double>::epsilon()*1e4,
655 "matchSphericalPointDistributions() - rotation does not work as expected by "
656 +toString((rotatednewCenter - oldCenter).NormSquared())+".");
657 }
658 }
659#endif
660 }
661 /// the second (orientation) rotation aligns the planes such that the
662 /// points themselves coincide
663 if (bestmatching.size() > 1) {
664 Rotation_t Rotation = findPointAligningRotation(
665 remainingold,
666 rotated_newpolygon,
667 bestmatching);
668
669 // construct RotationAxis and two points on its plane, defining the angle
670 Rotation.first.Normalize();
671 const Line RotationAxis(zeroVec, Rotation.first);
672
673 LOG(5, "DEBUG: Rotating around self is " << Rotation.second
674 << " around axis " << RotationAxis);
675
676#ifndef NDEBUG
677 // check: first bestmatching in rotated_newpolygon and remainingnew
678 // should now equal
679 {
680 const IndexList_t::const_iterator iter = bestmatching.begin();
681 Vector rotatednew = RotationAxis.rotateVector(
682 rotated_newpolygon[*iter],
683 Rotation.second);
684 LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
685 << " while old was " << remainingold[0]);
686 ASSERT( (rotatednew - remainingold[0]).Norm() < warn_amplitude,
687 "matchSphericalPointDistributions() - orientation rotation ends up off by more than "
688 +toString(warn_amplitude)+".");
689 }
690#endif
691
692 for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
693 iter != rotated_newpolygon.end(); ++iter) {
694 Vector &current = *iter;
695 LOG(6, "DEBUG: Original point is " << current);
696 current = RotationAxis.rotateVector(current, Rotation.second);
697 LOG(6, "DEBUG: Rotated point is " << current);
698 }
699 }
700
701 // remove all points in matching and return remaining ones
702 SphericalPointDistribution::Polygon_t remainingpoints =
703 removeMatchingPoints(rotated_newpolygon, bestmatching);
704 LOG(2, "INFO: Remaining points are " << remainingpoints);
705 return remainingpoints;
706 } else
707 return _newpolygon;
708}
709
710
Note: See TracBrowser for help on using the repository browser.