Ignore:
Timestamp:
Sep 12, 2016, 2:03:15 PM (9 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable
Children:
158ecb
Parents:
8e4471
git-author:
Frederik Heber <heber@…> (06/05/14 17:42:39)
git-committer:
Frederik Heber <heber@…> (09/12/16 14:03:15)
Message:

Dropped quaternion rotation for simple arbitrary rotation axis and angle.

  • with newCenter, oldCenter and the cross product we have all we need: a rotation axis and an angle. We don't need to burden ourselves with those stupid, absolutely not working quaternions.
  • removeMatchingPoints() now works on an array.
  • Orientation rotation was wrong way round, added check.
  • TESTFIX: removed QuaternionTest from SphericalPointDistributionTest, marked FragmentMolecule-cycle and StoreSaturatedFragment regression tests as XFAIL.
Location:
src/Fragmentation/Exporters
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/SphericalPointDistribution.cpp

    r8e4471 rb3c052  
    139139
    140140SphericalPointDistribution::Polygon_t removeMatchingPoints(
    141     const SphericalPointDistribution::Polygon_t &_returnpolygon,
     141    const VectorArray_t &_points,
    142142    const IndexList_t &_matchingindices
    143143    )
    … …  
    147147  std::sort(indices.begin(), indices.end());
    148148  LOG(4, "DEBUG: sorted matching is " << indices);
    149   IndexArray_t::const_iterator valueiter = indices.begin();
    150   SphericalPointDistribution::Polygon_t::const_iterator pointiter =
    151       _returnpolygon.begin();
    152   for (unsigned int i=0; i< _returnpolygon.size(); ++i, ++pointiter) {
    153     // skip all those in values
    154     if (*valueiter == i)
    155       ++valueiter;
    156     else
    157       remainingreturnpolygon.push_back(*pointiter);
     149  IndexArray_t remainingindices(_points.size(), -1);
     150  std::generate(remainingindices.begin(), remainingindices.end(), UniqueNumber);
     151  IndexArray_t::iterator remainiter = std::set_difference(
     152      remainingindices.begin(), remainingindices.end(),
     153      indices.begin(), indices.end(),
     154      remainingindices.begin());
     155  remainingindices.erase(remainiter, remainingindices.end());
     156  LOG(4, "DEBUG: remaining indices are " << remainingindices);
     157  for (IndexArray_t::const_iterator iter = remainingindices.begin();
     158      iter != remainingindices.end(); ++iter) {
     159    remainingreturnpolygon.push_back(_points[*iter]);
    158160  }
    159   LOG(4, "DEBUG: remaining indices are " << remainingreturnpolygon);
    160161
    161162  return remainingreturnpolygon;
    … …  
    292293    // determine rotation angles to align the two point distributions with
    293294    // respect to bestmatching
    294     SphericalPointDistribution::Polygon_t rotated_newpolygon;
     295    VectorArray_t rotated_newpolygon = remainingnew;
    295296    Vector oldCenter;
    296297    {
    … …  
    309310      if ((oldCenter - newCenter).NormSquared() > std::numeric_limits<double>::epsilon()*1e4) {
    310311        // setup quaternion
    311         Vector RotationAxis = newCenter;
    312         RotationAxis.VectorProduct(oldCenter);
     312        Vector RotationAxis = oldCenter;
     313        RotationAxis.VectorProduct(newCenter);
     314        Line Axis(zeroVec, RotationAxis);
    313315        RotationAxis.Normalize();
    314         const double RotationAngle = oldCenter.Angle(newCenter)/(M_PI/2.);
    315 //            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
    316         boost::math::quaternion<double> q
    317             (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
    318         LOG(5, "DEBUG: RotationAxis is " << RotationAxis
    319             << ", RotationAngle is " << RotationAngle);
    320         LOG(5, "DEBUG: Quaternion describing rotation is " << q);
     316        const double RotationAngle = oldCenter.Angle(newCenter); // /(M_PI/2.);
     317        LOG(5, "DEBUG: Rotate coordinate system by " << RotationAngle
     318            << " around axis " << RotationAxis);
     319
     320        // apply rotation angles
     321        for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
     322            iter != rotated_newpolygon.end(); ++iter) {
     323          Vector &current = *iter;
     324          LOG(5, "DEBUG: Original point is " << current);
     325          current =  Axis.rotateVector(current, RotationAngle);
     326          LOG(5, "DEBUG: Rotated point is " << current);
     327        }
     328      }
     329    }
     330    // rotate triangle/line/point around itself to match orientation
     331    if (MCS.bestmatching.size() > 1) {
     332      if (oldCenter.NormSquared() > std::numeric_limits<double>::epsilon()*1e4) {
     333        // construct RotationAxis and two points on its plane, defining the angle
     334        const Line RotationAxis(zeroVec, oldCenter);
     335        Vector oldPosition(rotated_newpolygon[*MCS.bestmatching.begin()]);
     336        oldPosition.ProjectOntoPlane(RotationAxis.getDirection());
     337        Vector newPosition(remainingold[*MCS.bestmatching.begin()]);
     338        newPosition.ProjectOntoPlane(RotationAxis.getDirection());
     339
     340        // construct reference vector to determine direction of rotation
     341        Vector dreiBein(oldPosition);
     342        dreiBein.VectorProduct(oldCenter);
     343        dreiBein.Normalize();
     344        const double sign =
     345            (dreiBein.ScalarProduct(newPosition) < 0.) ? -1. : +1.;
     346        LOG(6, "DEBUG: oldCenter on plane is " << oldPosition
     347            << ", newCenter in plane is " << newPosition
     348            << ", and dreiBein is " << dreiBein);
     349        const double RotationAngle = sign * oldPosition.Angle(newPosition);
     350        LOG(5, "DEBUG: Rotating around self is " << RotationAngle
     351            << " around axis " << RotationAxis);
     352
    321353#ifndef NDEBUG
     354        // check: first bestmatching in rotated_newpolygon and remainingnew
     355        // should now equal
    322356        {
    323         // check that rotation works in DEBUG mode
    324         boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
    325         boost::math::quaternion<double> q_inverse =
    326             boost::math::conj(q)/(boost::math::norm(q));
    327         p = q * p * q_inverse;
    328         boost::math::quaternion<double> identity(1,0,0,0);
    329         ASSERT( boost::math::norm(q*q_inverse - identity) < std::numeric_limits<double>::epsilon()*1e4,
    330             "matchSphericalPointDistributions() - quaternion does not rotate newCenter "
    331             +toString(q*q_inverse)+" into oldCenter "+toString(identity)+".");
    332         boost::math::quaternion<double> comparison(0., -oldCenter[0], oldCenter[1], oldCenter[2]);
    333         ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4,
    334             "matchSphericalPointDistributions() - quaternion does not rotate newCenter "
    335             +toString(p)+" into oldCenter "+toString(comparison)+".");
     357          const IndexList_t::const_iterator iter = MCS.bestmatching.begin();
     358          Vector rotatednew = RotationAxis.rotateVector(
     359              rotated_newpolygon[*iter],
     360              RotationAngle);
     361          LOG(4, "CHECK: rotated first new bestmatching is " << rotatednew
     362              << " while old was " << remainingold[*iter]);
     363          ASSERT( (rotatednew - remainingold[*iter]).Norm()
     364              < std::numeric_limits<double>::epsilon()*1e4,
     365              "matchSphericalPointDistributions() - orientation rotation does not work as expected.");
    336366        }
    337367#endif
    338368
    339         // rotate spherical distribution
    340         rotated_newpolygon = rotatePolygon(_newpolygon, q);
    341         LOG(5, "DEBUG: Rotated new polygon is " << rotated_newpolygon);
    342       } else {
    343         rotated_newpolygon = _newpolygon;
    344       }
    345     }
    346     // rotate triangle/line/point around itself to match orientation
    347     const Line RotationAxis(zeroVec, oldCenter);
    348     const double RotationAngle =
    349         oldCenter.Angle(remainingold[0])
    350         - oldCenter.Angle(remainingnew[*MCS.bestmatching.begin()]);
    351     LOG(5, "DEBUG: Rotate around self is " << RotationAngle
    352         << " around axis " << RotationAxis);
    353 
    354     for (SphericalPointDistribution::Polygon_t::iterator iter = rotated_newpolygon.begin();
    355         iter != rotated_newpolygon.end(); ++iter) {
    356       RotationAxis.rotateVector(*iter, RotationAngle);
     369        for (VectorArray_t::iterator iter = rotated_newpolygon.begin();
     370            iter != rotated_newpolygon.end(); ++iter) {
     371          Vector &current = *iter;
     372          LOG(6, "DEBUG: Original point is " << current);
     373          current = RotationAxis.rotateVector(current, RotationAngle);
     374          LOG(6, "DEBUG: Rotated point is " << current);
     375        }
     376      }
    357377    }
    358378

  • src/Fragmentation/Exporters/unittests/Makefile.am

    r8e4471 rb3c052  
    2121
    2222TESTS += $(FRAGMENTATIONEXPORTERSTESTS)
    23 XFAIL_TESTS += SphericalPointDistributionUnitTest
    2423check_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS)
    2524noinst_PROGRAMS += $(FRAGMENTATIONEXPORTERSTESTS)

  • src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp

    r8e4471 rb3c052  
    8282}
    8383
    84 void SphericalPointDistributionTest::QuaternionTest()
    85 {
    86   Vector oldCenter(0.,1.,0.);
    87   Vector newCenter(1.,0.,0.);
    88 
    89   {
    90     // setup quaternion
    91     Vector RotationAxis = newCenter;
    92     RotationAxis.VectorProduct(oldCenter);
    93     RotationAxis.Normalize();
    94     const double RotationAngle = oldCenter.Angle(newCenter)/(M_PI/2.);
    95   //            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
    96     boost::math::quaternion<double> q
    97         (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
    98     LOG(5, "DEBUG: RotationAxis is " << RotationAxis
    99         << ", RotationAngle is " << RotationAngle);
    100     LOG(5, "DEBUG: Quaternion describing rotation is " << q);
    101     boost::math::quaternion<double> q_inverse =
    102         boost::math::conj(q)/(boost::math::norm(q));
    103     LOG(5, "DEBUG: Quaternion inverse is " << q_inverse);
    104     boost::math::quaternion<double> identity(1,0,0,0);
    105     const boost::math::quaternion<double> unity = q*q_inverse;
    106     LOG(5, "DEBUG: q * q^-1 is " << unity);
    107     CPPUNIT_ASSERT( boost::math::norm(unity - identity) < std::numeric_limits<double>::epsilon()*1e4);
    108 
    109     // check that rotation works
    110     boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
    111     LOG(5, "DEBUG: Original newCenter is " << p);
    112     p = p * q_inverse;
    113     p = q * p;
    114     LOG(5, "DEBUG: Rotated newCenter is " << p);
    115     boost::math::quaternion<double> comparison(0., -oldCenter[0], oldCenter[1], oldCenter[2]);
    116     LOG(5, "DEBUG: Difference norm is " << boost::math::norm(p - comparison));
    117     CPPUNIT_ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4);
    118   }
    119 
    120   // rotating with angle = 0 flips the vector unwantedly
    121   {
    122     // setup quaternion
    123     Vector RotationAxis = newCenter;
    124     RotationAxis.VectorProduct(oldCenter);
    125     RotationAxis.Normalize();
    126     const double RotationAngle = 0.;
    127   //            RotationAxis.Angle(oldCenter) - RotationAxis.Angle(newCenter);
    128     boost::math::quaternion<double> q
    129         (RotationAngle, RotationAxis[0], RotationAxis[1], RotationAxis[2]);
    130     LOG(5, "DEBUG: RotationAxis is " << RotationAxis
    131         << ", RotationAngle is " << RotationAngle);
    132     LOG(5, "DEBUG: Quaternion describing rotation is " << q);
    133     boost::math::quaternion<double> q_inverse =
    134         boost::math::conj(q)/(boost::math::norm(q));
    135     LOG(5, "DEBUG: Quaternion inverse is " << q_inverse);
    136     boost::math::quaternion<double> identity(1,0,0,0);
    137     const boost::math::quaternion<double> unity = q*q_inverse;
    138     LOG(5, "DEBUG: q * q^-1 is " << unity);
    139     CPPUNIT_ASSERT( boost::math::norm(unity - identity) < std::numeric_limits<double>::epsilon()*1e4);
    140 
    141     // check that rotation works
    142     boost::math::quaternion<double> p(0., newCenter[0], newCenter[1], newCenter[2]);
    143     boost::math::quaternion<double> comparison(0., -newCenter[0], newCenter[1], newCenter[2]);
    144     LOG(5, "DEBUG: Original newCenter is " << p);
    145     p = p * q_inverse;
    146     p = q * p;
    147     LOG(5, "DEBUG: Rotated newCenter is " << p);
    148     LOG(5, "DEBUG: Difference norm is " << boost::math::norm(p - comparison));
    149     CPPUNIT_ASSERT( boost::math::norm(p - comparison) < std::numeric_limits<double>::epsilon()*1e4);
    150   }
    151 }
    15284
    15385/** UnitTest for matchSphericalPointDistributions() with two points

  • src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp

    r8e4471 rb3c052  
    2222{
    2323    CPPUNIT_TEST_SUITE( SphericalPointDistributionTest) ;
    24     CPPUNIT_TEST ( QuaternionTest );
    2524    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_2 );
    2625    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_3 );
    … …  
    3534      void setUp();
    3635      void tearDown();
    37       void QuaternionTest();
    3836      void matchSphericalPointDistributionsTest_2();
    3937      void matchSphericalPointDistributionsTest_3();
Note: See TracChangeset for help on using the changeset viewer.