Context Navigation

-              r6393ff
+              r2d8c4e
       _referencepositions, _currentpositions);
   if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
     LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
     oldCenter.Normalize();
     newCenter.Normalize();
     if (!oldCenter.IsEqualTo(newCenter)) {
       // calculate rotation axis and angle
       Rotation.first = oldCenter;
       Rotation.first.VectorProduct(newCenter);
       Rotation.second = oldCenter.Angle(newCenter); // /(M_PI/2.);
     } else {
       // no rotation required anymore
+    }
+  ASSERT( !oldCenter.IsZero() && !newCenter.IsZero(),
+      "findPlaneAligningRotation() - either old "+toString(oldCenter)
+      +" or new center "+toString(newCenter)+" are zero.");
+  LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
+  if (!oldCenter.IsEqualTo(newCenter)) {
+    // calculate rotation axis and angle
+    Rotation.first = oldCenter;
+    Rotation.first.VectorProduct(newCenter);
+    Rotation.first.Normalize();
+    // construct reference vector to determine direction of rotation
+    const double sign = determineSignOfRotation(newCenter, oldCenter, Rotation.first);
+    Rotation.second = sign * oldCenter.Angle(newCenter);
   } else {
+    LOG(4, "DEBUG: oldCenter is " << oldCenter << ", newCenter is " << newCenter);
+    if ((oldCenter.IsZero()) && (newCenter.IsZero())) {
+      // either oldCenter or newCenter (or both) is directly at origin
+      if (_currentpositions.indices.size() == 2) {
+        // line case
+        Vector oldPosition = _currentpositions.polygon[*_currentpositions.indices.begin()];
+        Vector newPosition = _referencepositions.polygon[0];
+        // check whether we need to rotate at all
+        if (!oldPosition.IsEqualTo(newPosition)) {
+          Rotation.first = oldPosition;
+          Rotation.first.VectorProduct(newPosition);
+          // orientation will fix the sign here eventually
+          Rotation.second = oldPosition.Angle(newPosition);
+        } else {
+          // no rotation required anymore
+        }
+      } else {
+        // triangle case
+        // both triangles/planes have same center, hence get axis by
+        // VectorProduct of Normals
+        Plane newplane(_referencepositions.polygon[0], _referencepositions.polygon[1], _referencepositions.polygon[2]);
+        VectorArray_t vectors;
+        for (IndexList_t::const_iterator iter = _currentpositions.indices.begin();
+            iter != _currentpositions.indices.end(); ++iter)
+          vectors.push_back(_currentpositions.polygon[*iter]);
+        Plane oldplane(vectors[0], vectors[1], vectors[2]);
+        Vector oldPosition = oldplane.getNormal();
+        Vector newPosition = newplane.getNormal();
+        // check whether we need to rotate at all
+        if (!oldPosition.IsEqualTo(newPosition)) {
+          Rotation.first = oldPosition;
+          Rotation.first.VectorProduct(newPosition);
+          Rotation.first.Normalize();
+          // construct reference vector to determine direction of rotation
+          const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
+          Rotation.second = sign * oldPosition.Angle(newPosition);
+          LOG(5, "DEBUG: Rotating plane normals by " << Rotation.second
+              << " around axis " << Rotation.first);
+        } else {
+          // else do nothing
+        }
+      }
+    } else {
+      // TODO: we can't do anything here, but this case needs to be dealt with when
+      // we have no ideal geometries anymore
+      if ((oldCenter-newCenter).Norm() > warn_amplitude)
+        ELOG(2, "oldCenter is " << oldCenter << ", yet newCenter is " << newCenter);
+#ifndef NDEBUG
+      // else they are considered close enough
+      dontcheck = true;
+#endif
+    }
+    // no rotation required anymore
+  }
 …
       << oldPosition.Norm() << ") and newPosition is " << newPosition << " length(: "
       << newPosition.Norm() << ")");
   if (!oldPosition.IsEqualTo(newPosition)) {
+    if ((!oldCenter.IsZero()) && (!newCenter.IsZero())) {
+      // we rotate at oldCenter and around the radial direction, which is again given
+      // by oldCenter.
+      Rotation.first = oldCenter;
+      Rotation.first.Normalize();  // note weighted sum of normalized weight is not normalized
+      LOG(6, "DEBUG: Using oldCenter " << oldCenter << " as rotation center and "
+          << Rotation.first << " as axis.");
+      oldPosition -= oldCenter;
+      newPosition -= oldCenter;
+      oldPosition = (oldPosition - oldPosition.Projection(Rotation.first));
+      newPosition = (newPosition - newPosition.Projection(Rotation.first));
+      LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
+    } else {
+      if (remainingnew.indices.size() == 2) {
+        // line situation
+        try {
+          Plane oldplane(oldPosition, oldCenter, newPosition);
+          Rotation.first = oldplane.getNormal();
+          LOG(6, "DEBUG: Plane is " << oldplane << " and normal is " << Rotation.first);
+        } catch (LinearDependenceException &e) {
+          LOG(6, "DEBUG: Vectors defining plane are linearly dependent.");
+          // oldPosition and newPosition are on a line, just flip when not equal
+          if (!oldPosition.IsEqualTo(newPosition)) {
+            Rotation.first.Zero();
+            Rotation.first.GetOneNormalVector(oldPosition);
+            LOG(6, "DEBUG: For flipping we use Rotation.first " << Rotation.first);
+            assert( Rotation.first.ScalarProduct(oldPosition) < std::numeric_limits<double>::epsilon()*1e4);
+  //              Rotation.second = M_PI;
+          } else {
+            LOG(6, "DEBUG: oldPosition and newPosition are equivalent.");
+          }
+        }
+      } else {
+        // triangle situation
+        Plane oldplane(remainingold.polygon[0], remainingold.polygon[1], remainingold.polygon[2]);
+        Rotation.first = oldplane.getNormal();
+        LOG(6, "DEBUG: oldPlane is " << oldplane << " and normal is " << Rotation.first);
+        oldPosition.ProjectOntoPlane(Rotation.first);
+        LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
+      }
+    }
+    // we rotate at oldCenter and around the radial direction, which is again given
+    // by oldCenter.
+    Rotation.first = oldCenter;
+    Rotation.first.Normalize();  // note weighted sum of normalized weight is not normalized
+    LOG(6, "DEBUG: Using oldCenter " << oldCenter << " as rotation center and "
+        << Rotation.first << " as axis.");
+    oldPosition -= oldCenter;
+    newPosition -= oldCenter;
+    oldPosition = (oldPosition - oldPosition.Projection(Rotation.first));
+    newPosition = (newPosition - newPosition.Projection(Rotation.first));
+    LOG(6, "DEBUG: Positions after projection are " << oldPosition << " and " << newPosition);
     // construct reference vector to determine direction of rotation
     const double sign = determineSignOfRotation(oldPosition, newPosition, Rotation.first);
 …
 #ifndef NDEBUG
       // check: rotated "newCenter" should now equal oldCenter
+      {
+      // we don't check in case of two points as these lie on a great circle
+      // and the center cannot stably be recalculated. We may reactivate this
+      // when we calculate centers only once
+      if (oldSet.indices.size() > 2) {
         Vector oldCenter;
         Vector rotatednewCenter;

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 2d8c4e for src/Fragmentation

Legend:

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Download in other formats: