source: src/Tesselation/triangleintersectionlist.cpp@ d94e7e

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * triangleintersectionlist.cpp
 *
 * This files encapsulates the TriangleIntersectionList class where all matters related to points in space being how close to
 * a tesselated surface are answered, i.e. distance, is inside, ...
 *
 *  Created on: Mar 1, 2010
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include <boost/scoped_ptr.hpp>

#include "triangleintersectionlist.hpp"

#include "tesselation.hpp"
#include "BoundaryMaps.hpp"
#include "BoundaryPointSet.hpp"
#include "BoundaryLineSet.hpp"
#include "BoundaryTriangleSet.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Helpers/defs.hpp"
#include "LinearAlgebra/Vector.hpp"

/** Constructor for class TriangleIntersectionList.
 *
 */
TriangleIntersectionList::TriangleIntersectionList(const Vector &x, const Tesselation *surface, const LinkedCell_deprecated* LC) :
  Point(x),
  Tess(surface),
  Vicinity(LC)
{
  //Info FunctionInfo(__func__);
  GatherIntersectionsWithTriangles();
};

/** Destructor for class TriangleIntersectionList.
 * Free's all allocated intersection vectors
 */
TriangleIntersectionList::~TriangleIntersectionList()
{
  //Info FunctionInfo(__func__);
  for (TriangleVectorMap::iterator Runner = IntersectionList.begin(); Runner != IntersectionList.end(); Runner++)
    delete((*Runner).second);
};

/** Returns the smallest distance between TriangleIntersectionList::Point and the surface given by
 * TriangleIntersectionList::Tess.
 * \return distance >=0, -1 if no specific distance could be found
 */
double TriangleIntersectionList::GetSmallestDistance() const
{
  //Info FunctionInfo(__func__);
  FillDistanceList();
  if (!DistanceList.empty())
    return DistanceList.begin()->first;
  else
    // return reference, if none can be found
    return -1.;
};

/** Returns the vector of closest intersection between TriangleIntersectionList::Point and the surface
 * TriangleIntersectionList::Tess.
 * \return Intersection vector or TriangleIntersectionList::Point if none could be found
 */
Vector TriangleIntersectionList::GetClosestIntersection() const
{
  //Info FunctionInfo(__func__);
  TriangleVectorMap::const_iterator runner;
  runner = GetIteratortoSmallestDistance();
  if (runner != IntersectionList.end()) {
    return *((*runner).second);
  }
  // return reference, if none can be found
  return Point;
};

/** Returns the triangle closest to TriangleIntersectionList::Point and the surface
 * TriangleIntersectionList::Tess.
 * \return pointer to triangle or NULL if none could be found
 */
BoundaryTriangleSet * TriangleIntersectionList::GetClosestTriangle() const
{
  //Info FunctionInfo(__func__);
  TriangleVectorMap::const_iterator runner;
  runner = GetIteratortoSmallestDistance();
  if (runner != IntersectionList.end()) {
    return ((*runner).first);
  }
  // return reference, if none can be found
  return NULL;
};

/** Checks whether reference TriangleIntersectionList::Point of this class is inside or outside.
 * \return true - TriangleIntersectionList::Point is inside, false - is outside
 */
bool TriangleIntersectionList::IsInside() const
{
  //Info FunctionInfo(__func__);
  TriangleVectorMap::const_iterator runner = GetIteratortoSmallestDistance();
  if (runner != IntersectionList.end()) {
    // if we have found one, check Scalarproduct between the vector
    Vector TestVector = (Point) - (*(*runner).second);
    LOG(1, "INFO: Distance vector between point " << Point
        << " and triangle intersection at " << (*(*runner).second) << " is "
        << TestVector << ".");
    if (fabs(TestVector.NormSquared()) < MYEPSILON)  {//
      LOG(1, "ACCEPT: Point is on the intersected triangle.");
      return true;
    }
    const double sign = (*runner).first->NormalVector.ScalarProduct(TestVector);
    LOG(1, "INFO: Checking sign of SKP between Distance vector and triangle's NormalVector "
        << (*runner).first->NormalVector << ": " << sign << ".");
    if (sign < 0) {
      LOG(1, "ACCEPT: Point lies on inner side of intersected triangle.");
      return true;
    } else {
      LOG(1, "REJECT: Point lies on outer side of intersected triangle.");
      return false;
    }
  }
  // so far away from surface that there are no triangles in list
  return false;
};

/** Puts all triangles in vicinity (by \a *LC) into a hash map with Intersection vectors.
 * Private function for constructor of TriangleIntersectionList.
 */
void TriangleIntersectionList::GatherIntersectionsWithTriangles()
{
  //Info FunctionInfo(__func__);

  // get closest points
  boost::scoped_ptr< DistanceToPointMap > points(Tess->FindClosestBoundaryPointsToVector(Point,Vicinity));
  if (points == NULL) {
    ELOG(1, "There is no nearest point: too far away from the surface.");
    return;
  }

  // gather all triangles in *LC-neighbourhood
  TriangleSet triangles;
  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++)
    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++)
      for (TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
        triangles.insert(TriangleRunner->second);

  Vector *Intersection = NULL;
  for (TriangleSet::const_iterator TriangleRunner = triangles.begin(); TriangleRunner != triangles.end(); TriangleRunner++) {
    // get intersection with triangle plane
    Intersection = new Vector;
    (*TriangleRunner)->GetClosestPointInsideTriangle(Point, *Intersection);
    //LOG(1, "Intersection between " << Point << " and " << **TriangleRunner << " is at " << *Intersection << ".");
    IntersectionList.insert( pair<BoundaryTriangleSet *, Vector * > (*TriangleRunner, Intersection) );
  }
};

/** Fills the private distance list
 */
void TriangleIntersectionList::FillDistanceList() const
{
  //Info FunctionInfo(__func__);
  if (DistanceList.empty())
    for (TriangleVectorMap::const_iterator runner = IntersectionList.begin(); runner != IntersectionList.end(); runner++) {
      // when the closest point is on the edge of a triangle (and hence
      // we find two closes triangles due to it having an adjacent one)
      // we should make sure that both triangles end up in the same entry
      // in the distance multimap. Hence, we round to 6 digit precision.
      const double distance =
          1e-6*floor(Point.distance(*(*runner).second)*1e+6);
      DistanceList.insert( pair<double, BoundaryTriangleSet *> (distance, (*runner).first) );
    }
  //for (DistanceTriangleMap::const_iterator runner = DistanceList.begin(); runner != DistanceList.end(); runner++)
  //  LOG(1, (*runner).first << " away from " << *(*runner).second);
};


/** Returns the iterator in VectorTriangleMap to the smallest distance.
 * This is a private helper function as the iterator is then evaluated in some other functions.
 * \return iterator or TriangleIntersectionList::IntersectionList.end() if none could be found
 */
TriangleVectorMap::const_iterator TriangleIntersectionList::GetIteratortoSmallestDistance() const
{
  //Info FunctionInfo(__func__);
  FillDistanceList();

  // do we have any intersections?
  TriangleVectorMap::const_iterator runner = IntersectionList.end();
  if (!DistanceList.empty()) {
    // get closest triangle(s)
    std::pair< DistanceTriangleMap::const_iterator, DistanceTriangleMap::const_iterator >  DistanceRange =
        DistanceList.equal_range(DistanceList.begin()->first);
    {
      size_t count = 0;
      for (DistanceTriangleMap::const_iterator iter = DistanceRange.first;
          iter != DistanceRange.second; ++iter) {
        LOG(3, "DEBUG: " << *(iter->second) << " is in the list with " << iter->first << ".");
        ++count;
      }
      LOG(1, "INFO: There are " << count << " possible triangles at the smallest distance.");
    }
    // if there is more than one, check all of their normal vectors
    // return the one with positive SKP if present because if the point
    // would be truely inside, all of the closest triangles would have to show
    // their inner side
    LOG(1, "INFO: Looking for first SKP greater than zero ...");
    for (DistanceTriangleMap::const_iterator iter = DistanceRange.first;
        iter != DistanceRange.second; ++iter) {
      // find the entry by the triangle key in IntersectionList to get the Intersection point
      runner = IntersectionList.find(iter->second);
      Vector TestVector = (Point) - (*(*runner).second);
      // construct SKP
      const double sign = (*iter).second->NormalVector.ScalarProduct(TestVector);
      LOG(1, "INFO: Checking SKP of closest triangle " << *(iter->second) << " = " << sign << ".");
      // if positive return
      if (sign > 0)
        break;
    }
    // if there's just one or all are negative, runner is not set yet
    if (runner == IntersectionList.end())
      runner = IntersectionList.find(DistanceList.begin()->second);
  }
  return runner;
};