Changeset 4a44ed for src/Fragmentation/Exporters/SphericalPointDistribution.cpp

Timestamp:

Aug 20, 2014, 1:07:11 PM (11 years ago)

Author:

Frederik Heber <heber@…>

Children:

c56380

Parents:

d635829

git-author:

Frederik Heber <heber@…> (07/20/14 11:26:01)

git-committer:

Frederik Heber <heber@…> (08/20/14 13:07:11)

Message:

tempcommit: Fixes to calculateCenterOfMinimumDistance. Merge with 2056d8e9

• we also check whether all points lie on a great circle, which also results in non-uniqueness.

File:

• src/Fragmentation/Exporters/SphericalPointDistribution.cpp (modified) (5 diffs)

src/Fragmentation/Exporters/SphericalPointDistribution.cpp

@@ -144 +144 @@
       return center;
     // one position given: return it directly
-    if (_positions.size() == (size_t)1)
-      return _positions[0];
+    if (_indices.size() == (size_t)1)
+      return _positions[*_indices.begin()];
     // two positions on a line given: return closest point to (1.,0.,0.)
-    if (fabs(_positions[0].ScalarProduct(_positions[1]) + 1.)
-        < std::numeric_limits<double>::epsilon()*1e4) {
-      Vector candidate;
-      if (_positions[0].ScalarProduct(center) > _positions[1].ScalarProduct(center))
-        candidate = _positions[0];
-      else
-        candidate = _positions[1];
-      // non-uniqueness: all positions on great circle, normal to given line are valid
-      // so, we just pick one because returning a unique point is topmost priority
-      Vector normal;
-      normal.GetOneNormalVector(candidate);
-      Vector othernormal = candidate;
-      othernormal.VectorProduct(normal);
-      // now both normal and othernormal describe the plane containing the great circle
-      Plane greatcircle(normal, zeroVec, othernormal);
-      // check which axis is contained and pick the one not
-      if (greatcircle.isContained(center)) {
-        center = Vector(0.,1.,0.);
-        if (greatcircle.isContained(center))
-          center = Vector(0.,0.,1.);
-        // now we are done cause a plane cannot contain all three axis vectors
-      }
-      center = greatcircle.getClosestPoint(candidate);
-      // assure length of 1
-      center.Normalize();
+//    IndexList_t::const_iterator indexiter = _indices.begin();
+//    const unsigned int firstindex = *indexiter++;
+//    const unsigned int secondindex = *indexiter;
+//    if ( fabs(_positions[firstindex].ScalarProduct(_positions[secondindex]) + 1.)
+//        < std::numeric_limits<double>::epsilon()*1e4) {
+//      Vector candidate;
+//      if (_positions[firstindex].ScalarProduct(center) > _positions[secondindex].ScalarProduct(center))
+//        candidate = _positions[firstindex];
+//      else
+//        candidate = _positions[secondindex];
+//      // non-uniqueness: all positions on great circle, normal to given line are valid
+//      // so, we just pick one because returning a unique point is topmost priority
+//      Vector normal;
+//      normal.GetOneNormalVector(candidate);
+//      Vector othernormal = candidate;
+//      othernormal.VectorProduct(normal);
+//      // now both normal and othernormal describe the plane containing the great circle
+//      Plane greatcircle(normal, zeroVec, othernormal);
+//      // check which axis is contained and pick the one not
+//      if (greatcircle.isContained(center)) {
+//        center = Vector(0.,1.,0.);
+//        if (greatcircle.isContained(center))
+//          center = Vector(0.,0.,1.);
+//        // now we are done cause a plane cannot contain all three axis vectors
+//      }
+//      center = greatcircle.getClosestPoint(candidate);
+//      // assure length of 1
+//      center.Normalize();
+//
+//      return center;
+//    }
+    // given points lie on a great circle and go completely round it
+    // two or more positions on a great circle given: return closest point to (1.,0.,0.)
+    {
+      bool AllNormal = true;
+      Vector Normal;
+      {
+        IndexList_t::const_iterator indexiter = _indices.begin();
+        Normal = _positions[*indexiter++];
+        Normal.VectorProduct(_positions[*indexiter++]);
+        Normal.Normalize();
+        for (;(AllNormal) && (indexiter != _indices.end()); ++indexiter)
+          AllNormal &= _positions[*indexiter].IsNormalTo(Normal, 1e-8);
+      }
+      double AngleSum = 0.;
+      if (AllNormal) {
+        // check by angle sum whether points go round are cover just one half
+        IndexList_t::const_iterator indexiter = _indices.begin();
+        Vector CurrentVector = _positions[*indexiter++];
+        for(; indexiter != _indices.end(); ++indexiter) {
+          AngleSum += CurrentVector.Angle(_positions[*indexiter]);
+          CurrentVector = _positions[*indexiter];
+        }
+      }
+      if (AngleSum - M_PI > std::numeric_limits<double>::epsilon()*1e4) {
+//        Vector candidate;
+//        double oldSKP = -1.;
+//        for (IndexList_t::const_iterator iter = _indices.begin();
+//            iter != _indices.end(); ++iter) {
+//          const double newSKP = _positions[*iter].ScalarProduct(center);
+//          if (newSKP > oldSKP) {
+//            candidate = _positions[*iter];
+//            oldSKP = newSKP;
+//          }
+//        }
+        // non-uniqueness: all positions on great circle, normal to given line are valid
+        // so, we just pick one because returning a unique point is topmost priority
+//        Vector normal;
+//        normal.GetOneNormalVector(candidate);
+//        Vector othernormal = candidate;
+//        othernormal.VectorProduct(normal);
+//        // now both normal and othernormal describe the plane containing the great circle
+//        Plane greatcircle(normal, zeroVec, othernormal);
+        // check which axis is contained and pick the one not
+//        if (greatcircle.isContained(center)) {
+//          center = Vector(0.,1.,0.);
+//          if (greatcircle.isContained(center))
+//            center = Vector(0.,0.,1.);
+//          // now we are done cause a plane cannot contain all three axis vectors
+//        }
+//        center = greatcircle.getClosestPoint(candidate);
+//        center = greatcircle.getNormal();
+        center = Normal;
+        // assure length of 1
+        center.Normalize();
+
+        return center;
+      }
     }
   }
@@ -179 +242 @@
   if (!center.IsZero()) {
     center.Normalize();
-    LOG(4, "DEBUG: Starting with geographical midpoint of " << _positions << " under indices "
+    LOG(5, "DEBUG: Starting with geographical midpoint of " << _positions << " under indices "
         << _indices << " is " << center);
   } else {
     // any point is good actually
     center = _positions[0];
-    LOG(4, "DEBUG: Starting with first position " << center);
+    LOG(5, "DEBUG: Starting with first position " << center);
   }
 
   // calculate initial MinimumDistance
   double MinimumDistance = calculateMinimumDistance(center, _positions, _indices);
-  LOG(5, "DEBUG: MinimumDistance to this center is " << MinimumDistance);
+  LOG(6, "DEBUG: MinimumDistance to this center is " << MinimumDistance);
 
   // check all present points whether they may serve as a better center
@@ -199 +262 @@
       MinimumDistance = candidateDistance;
       center = centerCandidate;
-      LOG(5, "DEBUG: new MinimumDistance to current test point " << MinimumDistance
+      LOG(6, "DEBUG: new MinimumDistance to current test point " << MinimumDistance
           << " is " << center);
     }
   }
-  LOG(5, "DEBUG: new MinimumDistance to center " << center << " is " << MinimumDistance);
+  LOG(6, "DEBUG: new MinimumDistance to center " << center << " is " << MinimumDistance);
 
   // now iterate over TestDistance
@@ -231 +294 @@
         center = centerCandidate;
         updatedflag = true;
-        LOG(5, "DEBUG: new MinimumDistance to test point at " << *angleiter/M_PI*180.
+        LOG(7, "DEBUG: new MinimumDistance to test point at " << *angleiter/M_PI*180.
             << "° is " << MinimumDistance);
       }
@@ -411 +474 @@
       errors.second += gap*gap;
     }
+    // there is at least one distance, we've checked that before
+    errors.second *= 1./(double)firstdistances.size();
   } else {
     // check whether we have any zero centers: Combining points on new sphere may result