Ignore:
Timestamp:
Jan 8, 2010, 2:04:22 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Children:
22b47e
Parents:
eda56a
Message:

Attempt to fix the tesselation::IsInnerPoint().

We try the IsInnerPoint() as follows:

  1. Find nearest BoundaryPoints - working
  2. Find Closest BoundaryLine's - working
  3. Find closest Triangle that is well aligned (wrt to NormalVector and Distance) - unsure whether correctly working
  4. Check whether alignment is on boundary or inside/outside - working
  5. If on boundary, we check whether it's inside of triangle by intersecting with boundary lines - not working

Hence, we code a wrapper for GSL routines, to - finally - allow for solution of linear system of equations.

Signed-off-by: Frederik Heber <heber@…>

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/tesselation.cpp

    reda56a ra1acc5  
    388388};
    389389
    390 /** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
     390/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
    391391 * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
    392392 * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
    … …  
    411411  }
    412412
     413  Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl;
     414  Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl;
     415  Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
     416
    413417  // Calculate cross point between one baseline and the line from the third endpoint to intersection
    414418  int i=0;
    415419  do {
    416420    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
     421      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
    417422      helper.CopyVector(endpoints[(i+1)%3]->node->node);
    418423      helper.SubtractVector(endpoints[i%3]->node->node);
     424      break;
    419425    } else
    420426      i++;
    421     if (i>2)
    422       break;
    423   } while (CrossPoint.NormSquared() < MYEPSILON);
     427  } while (i<3);
    424428  if (i==3) {
    425429    eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
    426430  }
    427   CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
     431  Log() << Verbose(1) << "INFO: Crosspoint is " << CrossPoint << "." << endl;
    428432
    429433  // check whether intersection is inside or not by comparing length of intersection and length of cross point
    … …  
    29872991 * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
    29882992 */
    2989 DistanceMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
     2993DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
    29902994{
    29912995  Info FunctionInfo(__func__);
    … …  
    30043008  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    30053009
    3006   DistanceMap * points = new DistanceMap;
     3010  DistanceToPointMap * points = new DistanceToPointMap;
    30073011  LC->GetNeighbourBounds(Nlower, Nupper);
    30083012  //Log() << Verbose(1) << endl;
    … …  
    30153019          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    30163020            FindPoint = PointsOnBoundary.find((*Runner)->nr);
    3017             if (FindPoint != PointsOnBoundary.end())
    3018               points->insert(DistancePair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
     3021            if (FindPoint != PointsOnBoundary.end()) {
     3022              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
     3023              Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
     3024            }
    30193025          }
    30203026        } else {
    … …  
    30423048
    30433049  // get closest points
    3044   DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
     3050  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
    30453051  if (points == NULL) {
    30463052    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
    … …  
    30553061  Vector Center;
    30563062  Vector BaseLine;
    3057   for (DistanceMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
     3063  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
    30583064    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
    30593065      // calculate closest point on line to desired point
    … …  
    31093115
    31103116        // get closest points
    3111         DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
     3117        DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
    31123118  if (points == NULL) {
    31133119    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
    … …  
    31183124  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
    31193125  BoundaryLineSet *ClosestLine = NULL;
    3120   double MinDistance = -1.;
    3121   Vector helper;
     3126  double MinDistance = 1e+16;
     3127  Vector BaseLineIntersection;
    31223128  Vector Center;
    31233129  Vector BaseLine;
    3124   for (DistanceMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
     3130  Vector BaseLineCenter;
     3131  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
    31253132    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
    3126       for (TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++) {
     3133
     3134      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
     3135      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
     3136      const double lengthBase = BaseLine.NormSquared();
     3137
     3138      BaseLineIntersection.CopyVector(x);
     3139      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
     3140      const double lengthEndA = BaseLineIntersection.NormSquared();
     3141
     3142      BaseLineIntersection.CopyVector(x);
     3143      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
     3144      const double lengthEndB = BaseLineIntersection.NormSquared();
     3145
     3146      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
     3147        if ((lengthEndA < MinDistance) && (lengthEndA < lengthEndB)) {
     3148          ClosestLine = LineRunner->second;
     3149          MinDistance = lengthEndA;
     3150          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEndA << "." << endl;
     3151        } else if (lengthEndB < MinDistance) {
     3152          ClosestLine = LineRunner->second;
     3153          MinDistance = lengthEndB;
     3154          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is closer with " << lengthEndB << "." << endl;
     3155        } else {
     3156          Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
     3157        }
     3158      } else { // intersection is closer, calculate
    31273159        // calculate closest point on line to desired point
    3128         helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
    3129         helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
    3130         helper.Scale(0.5);
    3131         Center.CopyVector(x);
    3132         Center.SubtractVector(&helper);
    3133         BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
    3134         BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
     3160        BaseLineIntersection.CopyVector(x);
     3161        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
     3162        Center.CopyVector(&BaseLineIntersection);
    31353163        Center.ProjectOntoPlane(&BaseLine);
    3136         const double distance = Center.NormSquared();
     3164        BaseLineIntersection.SubtractVector(&Center);
     3165        const double distance = BaseLineIntersection.NormSquared();
     3166        if (Center.NormSquared() > BaseLine.NormSquared()) {
     3167          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
     3168        }
    31373169        if ((ClosestLine == NULL) || (distance < MinDistance)) {
    3138           // additionally calculate intersection on line (whether it's on the line section or not)
    3139           helper.CopyVector(x);
    3140           helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
    3141           helper.SubtractVector(&Center);
    3142           const double lengthA = helper.ScalarProduct(&BaseLine);
    3143           helper.CopyVector(x);
    3144           helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
    3145           helper.SubtractVector(&Center);
    3146           const double lengthB = helper.ScalarProduct(&BaseLine);
    3147           if (lengthB*lengthA < 0) {  // if have different sign
    3148             ClosestLine = LineRunner->second;
    3149             MinDistance = distance;
    3150             Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
    3151           } else {
    3152             Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
    3153           }
     3170          ClosestLine = LineRunner->second;
     3171          MinDistance = distance;
     3172          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
    31543173        } else {
    3155           Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
     3174          Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
    31563175        }
    31573176      }
    … …  
    31903209    return NULL;
    31913210
    3192   // reject all which are not closest
    3193   double MinDistance = -1.;
     3211  // go through all and pick the one with the best alignment to x
     3212  double MinAlignment = 2.*M_PI;
    31943213  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
    31953214    (*Runner)->GetCenter(&Center);
    31963215    helper.CopyVector(x);
    31973216    helper.SubtractVector(&Center);
    3198     helper.ProjectOntoPlane(&(*Runner)->NormalVector);
    3199     const double distance = helper.NormSquared();
    3200     if (fabs(distance - MinDistance) < MYEPSILON) {  // degenerate case, keep both
    3201       candidates.push_back(*Runner);
    3202     } else if (distance < MinDistance) {
    3203       candidates.clear();
    3204       candidates.push_back(*Runner);
    3205       MinDistance = distance;
     3217    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
     3218    if (Alignment < MinAlignment) {
     3219      result = *Runner;
     3220      MinAlignment = Alignment;
     3221      Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl;
     3222    } else {
     3223      Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl;
    32063224    }
    32073225  }
    32083226  delete(triangles);
    3209   if (!candidates.empty()) {
    3210     if (candidates.size() == 1) { // there is no degenerate case
    3211       result = candidates.front();
    3212       Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
    3213     } else {
    3214       result = candidates.front();
    3215       result->GetCenter(&Center);
    3216       Center.SubtractVector(x);
    3217       Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
    3218       if (Center.ScalarProduct(&result->NormalVector) < 0) {
    3219         result = candidates.back();
    3220         Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
    3221         if (Center.ScalarProduct(&result->NormalVector) < 0) {
    3222           eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
    3223         }
    3224       }
    3225     }
    3226   } else {
    3227     result = NULL;
    3228     Log() << Verbose(0) << "No closest triangle found." << endl;
    3229   }
     3227
    32303228  return result;
    32313229};
    32323230
    32333231/** Checks whether the provided Vector is within the tesselation structure.
     3232 * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
     3233 * towards or away from the given \a &Point.
    32343234 *
    32353235 * @param point of which to check the position
    … …  
    32453245  Vector Center;
    32463246  Vector helper;
     3247  Vector DistanceToCenter;
     3248  Vector Intersection;
    32473249
    32483250  if (result == NULL) {// is boundary point or only point in point cloud?
    … …  
    32553257  result->GetCenter(&Center);
    32563258  Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
    3257   Center.SubtractVector(&Point);
     3259  DistanceToCenter.CopyVector(&Center);
     3260  DistanceToCenter.SubtractVector(&Point);
    32583261  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << Center << "." << endl;
    3259   if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
     3262
     3263  // check whether we are on boundary
     3264  if (fabs(DistanceToCenter.ScalarProduct(&result->NormalVector)) < MYEPSILON) {
     3265    // calculate whether inside of triangle
     3266    DistanceToCenter.ProjectOntoPlane(&result->NormalVector);
     3267    DistanceToCenter.AddVector(&Center);
     3268    Center.CopyVector(&DistanceToCenter);
     3269    Center.SubtractVector(&result->NormalVector); // points towards MolCenter
     3270    DistanceToCenter.AddVector(&result->NormalVector); // points outside
     3271    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
     3272    if (result->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
     3273      Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
     3274      return true;
     3275    } else {
     3276      Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
     3277      return false;
     3278    }
     3279  }
     3280
     3281  // then check direction to boundary
     3282  if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
    32603283    Log() << Verbose(1) << Point << " is an inner point." << endl;
    32613284    return true;
    32623285  } else {
    3263     for (int i=0;i<3;i++) {
    3264       helper.CopyVector(result->endpoints[i]->node->node);
    3265       helper.SubtractVector(&Point);
    3266       if (Center.NormSquared() > helper.NormSquared())
    3267         Center.CopyVector(&helper);
    3268     }
    3269     if (Center.NormSquared() < epsilon*epsilon) {
    3270       Log() << Verbose(1) << Point << " is inner point, being sufficiently close (wrt " << epsilon << ") to boundary." << endl;
    3271       return true;
    3272     }
    32733286    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
    32743287    return false;
    32753288  }
    3276 }
    3277 
    3278 /** Checks whether the provided TesselPoint is within the tesselation structure.
    3279  *
    3280  * @param *Point of which to check the position
    3281  * @param *LC Linked Cell structure
    3282  * @param epsilon Distance of \a &Point to Center of triangle is tested greater against this value (standard: -MYEPSILON)
    3283  *
    3284  * @return true if the point is inside the tesselation structure, false otherwise
    3285  */
    3286 bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon) const
    3287 {
    3288         Info FunctionInfo(__func__);
    3289   return IsInnerPoint(*(Point->node), LC, epsilon);
    32903289}
    32913290
    … …  
    34603459  }
    34613460
     3461  // check whether there's something to do
     3462  if (SetOfNeighbours->size() < 3) {
     3463    for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
     3464      connectedCircle->push_back(*TesselRunner);
     3465    return connectedCircle;
     3466  }
     3467
     3468  Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
    34623469  // calculate central point
    3463   for (TesselPointSet::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
    3464     center.AddVector((*TesselRunner)->node);
     3470
     3471  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
     3472  TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
     3473  TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
     3474  TesselB++;
     3475  TesselC++;
     3476  TesselC++;
     3477  int counter = 0;
     3478  while (TesselC != SetOfNeighbours->end()) {
     3479    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
     3480    Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
     3481    counter++;
     3482    TesselA++;
     3483    TesselB++;
     3484    TesselC++;
     3485    PlaneNormal.AddVector(&helper);
     3486  }
    34653487  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
    3466   //  << "; scale factor " << 1.0/connectedPoints.size();
    3467   center.Scale(1.0/SetOfNeighbours->size());
    3468   Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
    3469 
    3470   // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
    3471   PlaneNormal.CopyVector(Point->node);
    3472   PlaneNormal.SubtractVector(&center);
    3473   PlaneNormal.Normalize();
     3488  //  << "; scale factor " << counter;
     3489  PlaneNormal.Scale(1.0/(double)counter);
     3490//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
     3491//
     3492//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
     3493//  PlaneNormal.CopyVector(Point->node);
     3494//  PlaneNormal.SubtractVector(&center);
     3495//  PlaneNormal.Normalize();
    34743496  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
    34753497
    … …  
    35043526    helper.ProjectOntoPlane(&PlaneNormal);
    35053527    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     3528    if (angle > M_PI) // the correction is of no use here (and not desired)
     3529      angle = 2.*M_PI - angle;
    35063530    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
    35073531    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
Note: See TracChangeset for help on using the changeset viewer.