source: src/Shapes/Shape.cpp@ e32fa6

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * Shape.cpp
 *
 *  Created on: Jun 18, 2010
 *      Author: crueger
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "CodePatterns/Assert.hpp"
#include "LinearAlgebra/Vector.hpp"

#include "Shapes/Shape.hpp"
#include "Shapes/Shape_impl.hpp"
#include "Shapes/ShapeExceptions.hpp"
#include "Shapes/ShapeType.hpp"

#include "Tesselation/ApproximateShapeArea.hpp"
#include "Tesselation/ApproximateShapeVolume.hpp"

#include <algorithm>
#include <limits>
#include <string>


Shape::Shape(const Shape& src) :
  impl(src.getImpl())
{}

Shape::~Shape(){}

bool Shape::isInside(const Vector &point) const{
  return impl->isInside(point);
}

bool Shape::isOnSurface(const Vector &point) const{
  return impl->isOnSurface(point);
}

Vector Shape::getNormal(const Vector &point) const throw (NotOnSurfaceException){
  return impl->getNormal(point);
}

Vector Shape::getCenter() const{
  return impl->getCenter();
}

double Shape::getRadius() const{
  return impl->getRadius();
}

/** Returns the volume of the Shape.
 *
 * If the underlying implementation does not have a working implementation,
 * i.e. returns -1., then we use an approximate method to calculate the
 * volume via a mesh of grid points and checking for isInside (basically
 * a Monte-Carlo integration of the volume).
 *
 * \return volume of the shape
 */
double Shape::getVolume() const
{
        const double volume = impl->getVolume();
        if (volume  != -1.) {
                return volume;
        } else {
                ApproximateShapeVolume Approximator(*this);
                return Approximator();
        }
}

/** Returns the surface area of the Shape.
 *
 * If the underlying implementation does not have a working implementation,
 * i.e. returns -1., then we use the working filling of the shapes surface
 * with points and subsequent tesselation and obtaining the approximate
 * surface area therefrom.
 *
 * @return surface area of the Shape
 */
double Shape::getSurfaceArea() const
{
        const double surfacearea = impl->getSurfaceArea();
        if (surfacearea != -1.) {
                return surfacearea;
        } else {
                ApproximateShapeArea Approximator(*this);
                return Approximator();
        }
}

LineSegmentSet Shape::getLineIntersections(const Line &line) const{
  return impl->getLineIntersections(line);
}

std::vector<Vector> Shape::getHomogeneousPointsOnSurface(const size_t N) const {
  return impl->getHomogeneousPointsOnSurface(N);
}

std::vector<Vector> Shape::getHomogeneousPointsInVolume(const size_t N) const {
        return impl->getHomogeneousPointsInVolume(N);
}

Shape::Shape(Shape::impl_ptr _impl) :
    impl(_impl)
{}

Shape &Shape::operator=(const Shape& rhs){
  if(&rhs!=this){
    impl=rhs.getImpl();
  }
  return *this;
}

bool Shape::operator==(const Shape &rhs) const{
        return (this->getType() == rhs.getType());
}

std::string Shape::toString() const{
  return impl->toString();
}

enum ShapeType Shape::getType() const{
        return impl->getType();
}

Shape::impl_ptr Shape::getImpl() const{
  return impl;
}

// allows arbitrary friendship, but only if implementation is known
Shape::impl_ptr getShapeImpl(const Shape &shape){
  return shape.getImpl();
}

/***************************** Some simple Shapes ***************************/

Shape Everywhere(){
  static Shape::impl_ptr impl = Shape::impl_ptr(new Everywhere_impl());
  return Shape(impl);
}

Shape Nowhere(){
  static Shape::impl_ptr impl = Shape::impl_ptr(new Nowhere_impl());
  return Shape(impl);
}

/****************************** Operators ***********************************/

// AND

AndShape_impl::AndShape_impl(const Shape::impl_ptr &_lhs, const Shape::impl_ptr &_rhs) :
  lhs(_lhs),rhs(_rhs)
{}

AndShape_impl::~AndShape_impl(){}

bool AndShape_impl::isInside(const Vector &point) const{
  return lhs->isInside(point) && rhs->isInside(point);
}

bool AndShape_impl::isOnSurface(const Vector &point) const{
  // check the number of surfaces that this point is on
  int surfaces =0;
  surfaces += lhs->isOnSurface(point);
  surfaces += rhs->isOnSurface(point);

  switch(surfaces){
    case 0:
      return false;
      // no break necessary
    case 1:
      // if it is inside for the object where it does not lie on
      // the surface the whole point lies inside
      return (lhs->isOnSurface(point) && rhs->isInside(point)) ||
             (rhs->isOnSurface(point) && lhs->isInside(point));
      // no break necessary
    case 2:
      {
        // it lies on both Shapes... could be an edge or an inner point
        // test the direction of the normals
        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
        // if the directions are opposite we lie on the inside
        return !direction.IsZero();
      }
      // no break necessary
    default:
      // if this happens there is something wrong
      ASSERT(0,"Default case should have never been used");
  }
  return false; // never reached
}

Vector AndShape_impl::getNormal(const Vector &point) const throw (NotOnSurfaceException){
  Vector res;
  if(!isOnSurface(point)){
    throw NotOnSurfaceException() << ShapeVector(&point);
  }
  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
  res.Normalize();
  return res;
}

Vector AndShape_impl::getCenter() const
{
  // calculate closest position on sphere surface to other center ..
  const Vector rhsDistance = rhs->getCenter() + rhs->getRadius()*((lhs->getCenter() - rhs->getCenter()).getNormalized());
  const Vector lhsDistance = lhs->getCenter() + lhs->getRadius()*((rhs->getCenter() - lhs->getCenter()).getNormalized());
  // .. and then it's right in between those two
  return 0.5*(rhsDistance + lhsDistance);
}

double AndShape_impl::getRadius() const
{
  const double distance = (lhs->getCenter() - rhs->getCenter()).Norm();
  const double minradii = std::min(lhs->getRadius(), rhs->getRadius());
  // if no intersection
  if (distance > (lhs->getRadius() + rhs->getRadius()))
    return 0.;
  else // if intersection it can only be the smaller one
    return minradii;
}

double AndShape_impl::getVolume() const
{
        // TODO
        return -1.;
}

double AndShape_impl::getSurfaceArea() const
{
        // TODO
        return -1.;
}

LineSegmentSet AndShape_impl::getLineIntersections(const Line &line) const{
  return intersect(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
}

std::string AndShape_impl::toString() const{
  return std::string("(") + lhs->toString() + std::string("&&") + rhs->toString() + std::string(")");
}

enum ShapeType AndShape_impl::getType() const{
        return CombinedType;
}

std::vector<Vector> AndShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
  std::vector<Vector> PointsOnSurface_lhs = lhs->getHomogeneousPointsOnSurface(N);
  std::vector<Vector> PointsOnSurface_rhs = rhs->getHomogeneousPointsOnSurface(N);
  std::vector<Vector> PointsOnSurface;

  for (std::vector<Vector>::const_iterator iter = PointsOnSurface_lhs.begin(); iter != PointsOnSurface_lhs.end(); ++iter) {
    if (rhs->isInside(*iter))
      PointsOnSurface.push_back(*iter);
  }
  for (std::vector<Vector>::const_iterator iter = PointsOnSurface_rhs.begin(); iter != PointsOnSurface_rhs.end(); ++iter) {
    if (lhs->isInside(*iter))
      PointsOnSurface.push_back(*iter);
  }

  return PointsOnSurface;
}

std::vector<Vector> AndShape_impl::getHomogeneousPointsInVolume(const size_t N) const {
        ASSERT(0,
                        "AndShape_impl::getHomogeneousPointsInVolume() - not implemented.");
        return std::vector<Vector>();
}


Shape operator&&(const Shape &lhs,const Shape &rhs){
  Shape::impl_ptr newImpl = Shape::impl_ptr(new AndShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
  return Shape(newImpl);
}

// OR

OrShape_impl::OrShape_impl(const Shape::impl_ptr &_lhs, const Shape::impl_ptr &_rhs) :
  lhs(_lhs),rhs(_rhs)
{}

OrShape_impl::~OrShape_impl(){}

bool OrShape_impl::isInside(const Vector &point) const{
  return rhs->isInside(point) || lhs->isInside(point);
}

bool OrShape_impl::isOnSurface(const Vector &point) const{
  // check the number of surfaces that this point is on
  int surfaces =0;
  surfaces += lhs->isOnSurface(point);
  surfaces += rhs->isOnSurface(point);

  switch(surfaces){
    case 0:
      return false;
      // no break necessary
    case 1:
      // if it is inside for the object where it does not lie on
      // the surface the whole point lies inside
      return (lhs->isOnSurface(point) && !rhs->isInside(point)) ||
             (rhs->isOnSurface(point) && !lhs->isInside(point));
      // no break necessary
    case 2:
      {
        // it lies on both Shapes... could be an edge or an inner point
        // test the direction of the normals
        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
        // if the directions are opposite we lie on the inside
        return !direction.IsZero();
      }
      // no break necessary
    default:
      // if this happens there is something wrong
      ASSERT(0,"Default case should have never been used");
  }
  return false; // never reached
}

Vector OrShape_impl::getNormal(const Vector &point) const throw (NotOnSurfaceException){
  Vector res;
  if(!isOnSurface(point)){
    throw NotOnSurfaceException() << ShapeVector(&point);
  }
  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
  res.Normalize();
  return res;
}

Vector OrShape_impl::getCenter() const
{
  // calculate furthest position on sphere surface to other center ..
  const Vector rhsDistance = rhs->getCenter() + rhs->getRadius()*((rhs->getCenter() - lhs->getCenter()).getNormalized());
  const Vector lhsDistance = lhs->getCenter() + lhs->getRadius()*((lhs->getCenter() - rhs->getCenter()).getNormalized());
  // .. and then it's right in between those two
  return .5*(rhsDistance + lhsDistance);
}

double OrShape_impl::getRadius() const
{
  const Vector rhsDistance = rhs->getCenter() + rhs->getRadius()*((rhs->getCenter() - lhs->getCenter()).getNormalized());
  const Vector lhsDistance = lhs->getCenter() + lhs->getRadius()*((lhs->getCenter() - rhs->getCenter()).getNormalized());
  return .5*(rhsDistance - lhsDistance).Norm();
}

double OrShape_impl::getVolume() const
{
        // TODO
        return -1.;
}

double OrShape_impl::getSurfaceArea() const
{
        // TODO
        return -1.;
}

LineSegmentSet OrShape_impl::getLineIntersections(const Line &line) const{
  return merge(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
}

std::string OrShape_impl::toString() const{
  return std::string("(") + lhs->toString() + std::string("||") + rhs->toString() + std::string(")");
}

enum ShapeType OrShape_impl::getType() const{
        return CombinedType;
}

std::vector<Vector> OrShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
  std::vector<Vector> PointsOnSurface_lhs = lhs->getHomogeneousPointsOnSurface(N);
  std::vector<Vector> PointsOnSurface_rhs = rhs->getHomogeneousPointsOnSurface(N);
  std::vector<Vector> PointsOnSurface;

  for (std::vector<Vector>::const_iterator iter = PointsOnSurface_lhs.begin(); iter != PointsOnSurface_lhs.end(); ++iter) {
    if (!rhs->isInside(*iter))
      PointsOnSurface.push_back(*iter);
  }
  for (std::vector<Vector>::const_iterator iter = PointsOnSurface_rhs.begin(); iter != PointsOnSurface_rhs.end(); ++iter) {
    if (!lhs->isInside(*iter))
      PointsOnSurface.push_back(*iter);
  }

  return PointsOnSurface;
}

std::vector<Vector> OrShape_impl::getHomogeneousPointsInVolume(const size_t N) const {
        ASSERT(0,
                        "OrShape_impl::getHomogeneousPointsInVolume() - not implemented.");
        return std::vector<Vector>();
}

Shape operator||(const Shape &lhs,const Shape &rhs){
  Shape::impl_ptr newImpl = Shape::impl_ptr(new OrShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
  return Shape(newImpl);
}

// NOT

NotShape_impl::NotShape_impl(const Shape::impl_ptr &_arg) :
  arg(_arg)
{}

NotShape_impl::~NotShape_impl(){}

bool NotShape_impl::isInside(const Vector &point) const{
  return !arg->isInside(point);
}

bool NotShape_impl::isOnSurface(const Vector &point) const{
  return arg->isOnSurface(point);
}

Vector NotShape_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException){
  return -1.*arg->getNormal(point);
}

Vector NotShape_impl::getCenter() const
{
  return arg->getCenter();
}

double NotShape_impl::getRadius() const
{
  return std::numeric_limits<double>::infinity();
}

double NotShape_impl::getVolume() const
{
        // TODO
        return -1.; //-arg->getVolume();
}

double NotShape_impl::getSurfaceArea() const
{
        // TODO
        return -1.; // -arg->getSurfaceArea();
}

LineSegmentSet NotShape_impl::getLineIntersections(const Line &line) const{
  return invert(arg->getLineIntersections(line));
}

std::string NotShape_impl::toString() const{
  return std::string("!") + arg->toString();
}

enum ShapeType NotShape_impl::getType() const{
        return CombinedType;
}

std::vector<Vector> NotShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
  // surfaces are the same, only normal direction is different
  return arg->getHomogeneousPointsOnSurface(N);
}

std::vector<Vector> NotShape_impl::getHomogeneousPointsInVolume(const size_t N) const {
        ASSERT(0,
                        "NotShape_impl::getHomogeneousPointsInVolume() - not implemented.");
        return std::vector<Vector>();
}

Shape operator!(const Shape &arg){
  Shape::impl_ptr newImpl = Shape::impl_ptr(new NotShape_impl(getShapeImpl(arg)));
  return Shape(newImpl);
}

/**************** global operations *********************************/
std::ostream &operator<<(std::ostream &ost,const Shape &shape){
  ost << shape.toString();
  return ost;
}