Changeset 1f591b for molecuilder/src/tesselationhelpers.cpp

Timestamp:

Apr 13, 2010, 1:22:42 PM (16 years ago)

Author:

Tillmann Crueger <crueger@…>

Children:

e7ea64

Parents:

0f55b2

Message:

Prepared interface of Vector Class for transition to VectorComposites

File:

• molecuilder/src/tesselationhelpers.cpp (modified) (19 diffs)

molecuilder/src/tesselationhelpers.cpp

@@ -88 +88 @@
   center->at(2) =  0.5 * m14/ m11;
 
-  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
-    eLog() << Verbose(1) << "The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
+  if (fabs(a.Distance(*center) - RADIUS) > MYEPSILON)
+    eLog() << Verbose(1) << "The given center is further way by " << fabs(a.Distance(*center) - RADIUS) << " from a than RADIUS." << endl;
 
   gsl_matrix_free(A);
@@ -121 +121 @@
   Vector OtherCenter;
   Center->Zero();
-  helper.CopyVector(&a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(&b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(&c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
+  helper = sin(2.*alpha) * a;
+  (*Center) += helper;
+  helper = sin(2.*beta) * b;
+  (*Center) += helper;
+  helper = sin(2.*gamma) * c;
+  (*Center) += helper;
   //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
   Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
-  NewUmkreismittelpunkt->CopyVector(Center);
+  (*NewUmkreismittelpunkt) = (*Center);
   Log() << Verbose(1) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
   // Here we calculated center of circumscribing circle, using barycentric coordinates
   Log() << Verbose(1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
 
-  TempNormal.CopyVector(&a);
-  TempNormal.SubtractVector(&b);
-  helper.CopyVector(&a);
-  helper.SubtractVector(&c);
-  TempNormal.VectorProduct(&helper);
+  TempNormal = a - b;
+  helper = a - c;
+  TempNormal.VectorProduct(helper);
   if (fabs(HalfplaneIndicator) < MYEPSILON)
     {
-      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(AlternativeDirection) >0 && AlternativeIndicator <0))
+      if ((TempNormal.ScalarProduct(*AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(*AlternativeDirection) >0 && AlternativeIndicator <0))
         {
-          TempNormal.Scale(-1);
+          TempNormal *= -1;
         }
     }
   else
     {
-      if (((TempNormal.ScalarProduct(Direction)<0) && (HalfplaneIndicator >0)) || ((TempNormal.ScalarProduct(Direction)>0) && (HalfplaneIndicator<0)))
+      if (((TempNormal.ScalarProduct(*Direction)<0) && (HalfplaneIndicator >0)) || ((TempNormal.ScalarProduct(*Direction)>0) && (HalfplaneIndicator<0)))
         {
-          TempNormal.Scale(-1);
+          TempNormal *= -1;
         }
     }
@@ -163 +158 @@
   Log() << Verbose(1) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
 
-  Center->AddVector(&TempNormal);
+  (*Center) += TempNormal;
   Log() << Verbose(1) << "Center of sphere of circumference is " << *Center << ".\n";
   GetSphere(&OtherCenter, a, b, c, RADIUS);
@@ -181 +176 @@
   Vector helper;
   double alpha, beta, gamma;
-  Vector SideA, SideB, SideC;
-  SideA.CopyVector(b);
-  SideA.SubtractVector(&c);
-  SideB.CopyVector(c);
-  SideB.SubtractVector(&a);
-  SideC.CopyVector(a);
-  SideC.SubtractVector(&b);
-  alpha = M_PI - SideB.Angle(&SideC);
-  beta = M_PI - SideC.Angle(&SideA);
-  gamma = M_PI - SideA.Angle(&SideB);
+  Vector SideA = b - c;
+  Vector SideB = c - a;
+  Vector SideC = a - b;
+  alpha = M_PI - SideB.Angle(SideC);
+  beta = M_PI - SideC.Angle(SideA);
+  gamma = M_PI - SideA.Angle(SideB);
   //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
   if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) {
@@ -197 +188 @@
 
   Center->Zero();
-  helper.CopyVector(a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
+  helper = sin(2.*alpha) * a;
+  (*Center) += helper;
+  helper = sin(2.*beta) * b;
+  (*Center) += helper;
+  helper = sin(2.*gamma) * c;
+  (*Center) += helper;
   Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
 };
@@ -227 +215 @@
   Vector helper;
   double radius, alpha;
-  Vector RelativeOldSphereCenter;
-  Vector RelativeNewSphereCenter;
-
-  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
-  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
-  RelativeNewSphereCenter.CopyVector(&NewSphereCenter);
-  RelativeNewSphereCenter.SubtractVector(&CircleCenter);
-  helper.CopyVector(&RelativeNewSphereCenter);
+
+  Vector RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
+  Vector RelativeNewSphereCenter = NewSphereCenter - CircleCenter;
+  helper = RelativeNewSphereCenter;
   // test whether new center is on the parameter circle's plane
-  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-    eLog() << Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
-    helper.ProjectOntoPlane(&CirclePlaneNormal);
+  if (fabs(helper.ScalarProduct(CirclePlaneNormal)) > HULLEPSILON) {
+    eLog() << Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(CirclePlaneNormal))  << "!" << endl;
+    helper.ProjectOntoPlane(CirclePlaneNormal);
   }
   radius = helper.NormSquared();
@@ -244 +228 @@
   if (fabs(radius - CircleRadius) > HULLEPSILON)
     eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-  alpha = helper.Angle(&RelativeOldSphereCenter);
+  alpha = helper.Angle(RelativeOldSphereCenter);
   // make the angle unique by checking the halfplanes/search direction
-  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
+  if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
     alpha = 2.*M_PI - alpha;
   Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl;
-  radius = helper.Distance(&RelativeOldSphereCenter);
-  helper.ProjectOntoPlane(&NormalVector);
+  radius = helper.Distance(RelativeOldSphereCenter);
+  helper.ProjectOntoPlane(NormalVector);
   // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
   if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
@@ -280 +264 @@
   struct Intersection *I = (struct Intersection *)params;
   Vector intersection;
-  Vector SideA,SideB,HeightA, HeightB;
   for (int i=0;i<NDIM;i++)
     intersection[i] = gsl_vector_get(x, i);
 
-  SideA.CopyVector(&(I->x1));
-  SideA.SubtractVector(&I->x2);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&I->x1);
-  HeightA.ProjectOntoPlane(&SideA);
-
-  SideB.CopyVector(&I->x3);
-  SideB.SubtractVector(&I->x4);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&I->x3);
-  HeightB.ProjectOntoPlane(&SideB);
-
-  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
+  Vector SideA = I->x1 -I->x2 ;
+  Vector HeightA = intersection - I->x1;
+  HeightA.ProjectOntoPlane(SideA);
+
+  Vector SideB = I->x3 - I->x4;
+  Vector HeightB = intersection - I->x3;
+  HeightB.ProjectOntoPlane(SideB);
+
+  retval = HeightA.ScalarProduct(HeightA) + HeightB.ScalarProduct(HeightB);
   //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
 
@@ -321 +300 @@
 
   struct Intersection par;
-    par.x1.CopyVector(&point1);
-    par.x2.CopyVector(&point2);
-    par.x3.CopyVector(&point3);
-    par.x4.CopyVector(&point4);
+    par.x1 = point1;
+    par.x2 = point2;
+    par.x3 = point3;
+    par.x4 = point4;
 
     const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
@@ -370 +349 @@
 
     // check whether intersection is in between or not
-  Vector intersection, SideA, SideB, HeightA, HeightB;
+  Vector intersection;
   double t1, t2;
   for (int i = 0; i < NDIM; i++) {
@@ -376 +355 @@
   }
 
-  SideA.CopyVector(&par.x2);
-  SideA.SubtractVector(&par.x1);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&par.x1);
-
-  t1 = HeightA.ScalarProduct(&SideA)/SideA.ScalarProduct(&SideA);
-
-  SideB.CopyVector(&par.x4);
-  SideB.SubtractVector(&par.x3);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&par.x3);
-
-  t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
+  Vector SideA = par.x2 - par.x1;
+  Vector HeightA = intersection - par.x1;
+
+  t1 = HeightA.ScalarProduct(SideA)/SideA.ScalarProduct(SideA);
+
+  Vector SideB = par.x4 - par.x3;
+  Vector HeightB = intersection - par.x3;
+
+  t2 = HeightB.ScalarProduct(SideB)/SideB.ScalarProduct(SideB);
 
   Log() << Verbose(1) << "Intersection " << intersection << " is at "
@@ -428 +403 @@
   if (point.IsZero())
     return M_PI;
-  double phi = point.Angle(&reference);
-  if (OrthogonalVector.ScalarProduct(&point) > 0) {
+  double phi = point.Angle(reference);
+  if (OrthogonalVector.ScalarProduct(point) > 0) {
     phi = 2.*M_PI - phi;
   }
@@ -456 +431 @@
   TetraederVector[2].CopyVector(c);
   for (int j=0;j<3;j++)
-    TetraederVector[j].SubtractVector(&d);
-  Point.CopyVector(&TetraederVector[0]);
-  Point.VectorProduct(&TetraederVector[1]);
-  volume = 1./6. * fabs(Point.ScalarProduct(&TetraederVector[2]));
+    TetraederVector[j].SubtractVector(d);
+  Point.CopyVector(TetraederVector[0]);
+  Point.VectorProduct(TetraederVector[1]);
+  volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));
   return volume;
 };
@@ -583 +558 @@
         if (List != NULL) {
           for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
+            helper = (*Point) - (*(*Runner)->node);
+            double currentNorm = helper.Norm();
             if (currentNorm < distance) {
               // remember second point
@@ -638 +612 @@
         if (List != NULL) {
           for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
+            helper = (*Point) - (*(*Runner)->node);
             double currentNorm = helper.NormSquared();
             if (currentNorm < distance) {
@@ -678 +651 @@
         Info FunctionInfo(__func__);
   // construct the plane of the two baselines (i.e. take both their directional vectors)
-  Vector Normal;
-  Vector Baseline, OtherBaseline;
-  Baseline.CopyVector(Base->endpoints[1]->node->node);
-  Baseline.SubtractVector(Base->endpoints[0]->node->node);
-  OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
-  OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
-  Normal.CopyVector(&Baseline);
-  Normal.VectorProduct(&OtherBaseline);
+  Vector Baseline = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
+  Vector OtherBaseline = (*OtherBase->endpoints[1]->node->node) - (*OtherBase->endpoints[0]->node->node);
+  Vector Normal = Baseline;
+  Normal.VectorProduct(OtherBaseline);
   Normal.Normalize();
   Log() << Verbose(1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
 
   // project one offset point of OtherBase onto this plane (and add plane offset vector)
-  Vector NewOffset;
-  NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
-  NewOffset.SubtractVector(Base->endpoints[0]->node->node);
-  NewOffset.ProjectOntoPlane(&Normal);
-  NewOffset.AddVector(Base->endpoints[0]->node->node);
-  Vector NewDirection;
-  NewDirection.CopyVector(&NewOffset);
-  NewDirection.AddVector(&OtherBaseline);
+  Vector NewOffset = (*OtherBase->endpoints[0]->node->node) - (*Base->endpoints[0]->node->node);
+  NewOffset.ProjectOntoPlane(Normal);
+  NewOffset += (*Base->endpoints[0]->node->node);
+  Vector NewDirection = NewOffset + OtherBaseline;
 
   // calculate the intersection between this projected baseline and Base
@@ -704 +669 @@
                                                    *(Base->endpoints[1]->node->node),
                                                      NewOffset, NewDirection);
-  Normal.CopyVector(Intersection);
-  Normal.SubtractVector(Base->endpoints[0]->node->node);
-  Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
+  Normal = (*Intersection) - (*Base->endpoints[0]->node->node);
+  Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl;
 
   return Intersection;
@@ -724 +688 @@
     return -1;
   }
-  distance = x->DistanceToPlane(&triangle->NormalVector, triangle->endpoints[0]->node->node);
+  distance = x->DistanceToPlane(triangle->NormalVector, *triangle->endpoints[0]->node->node);
   return distance;
 };
@@ -787 +751 @@
     Vector *center = cloud->GetCenter();
     // make the circumsphere's center absolute again
-    helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
-    helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
-    helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
-    helper.Scale(1./3.);
-    helper.SubtractVector(center);
+    Vector helper = (1./3.) * ((*Tess->LastTriangle->endpoints[0]->node->node) +
+                               (*Tess->LastTriangle->endpoints[1]->node->node) +
+                               (*Tess->LastTriangle->endpoints[2]->node->node));
+    helper -= (*center);
     // and add to file plus translucency object
     *rasterfile << "# current virtual sphere\n";