source: src/Shapes/BaseShapes.cpp@ d6203a

/*
 * 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.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include "Shapes/BaseShapes.hpp"
#include "Shapes/BaseShapes_impl.hpp"
#include "Shapes/ShapeExceptions.hpp"
#include "Shapes/ShapeOps.hpp"

#include "Helpers/defs.hpp"

#include "CodePatterns/Assert.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/Line.hpp"
#include "LinearAlgebra/Plane.hpp"
#include "LinearAlgebra/LineSegment.hpp"
#include "LinearAlgebra/LineSegmentSet.hpp"

#include <cmath>
#include <algorithm>

bool Sphere_impl::isInside(const Vector &point) const{
  return point.NormSquared()<=1.;
}

bool Sphere_impl::isOnSurface(const Vector &point) const{
  return fabs(point.NormSquared()-1.)<MYEPSILON;
}

Vector Sphere_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException){
  if(!isOnSurface(point)){
    throw NotOnSurfaceException() << ShapeVector(&point);
  }
  return point;
}

Vector Sphere_impl::getCenter() const
{
  return Vector(0.,0.,0.);
}

double Sphere_impl::getRadius() const
{
  return 1.;
}

double Sphere_impl::getVolume() const
{
        return (4./3.)*M_PI; // 4/3 pi r^3
}

double Sphere_impl::getSurfaceArea() const
{
        return 2.*M_PI; // 2 pi r^2
}


LineSegmentSet Sphere_impl::getLineIntersections(const Line &line) const{
  LineSegmentSet res(line);
  std::vector<Vector> intersections = line.getSphereIntersections();
  if(intersections.size()==2){
    res.insert(LineSegment(intersections[0],intersections[1]));
  }
  return res;
}

std::string Sphere_impl::toString() const{
  return "Sphere()";
}

enum ShapeType Sphere_impl::getType() const
{
        return SphereType;
}

/**
 * algorithm taken from http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere
 * \param N number of points on surface
 */
std::vector<Vector> Sphere_impl::getHomogeneousPointsOnSurface(const size_t N) const
{
  std::vector<Vector> PointsOnSurface;
  if (true) {
    // Exactly N points but not symmetric.

    // This formula is derived by finding a curve on the sphere that spirals down from
    // the north pole to the south pole keeping a constant distance between consecutive turns.
    // The curve is then parametrized by arch length and evaluated in constant intervals.
    double a = sqrt(N) * 2;
    for (int i=0; i<N; i++){
      double t0 = ((double)i + 0.5) / (double)N;
      double t = (sqrt(t0) - sqrt(1.0 - t0) + 1.0) / 2.0 * M_PI;
      Vector point;
      point.Zero();
      point[0] = sin(t) * sin(t * a);
      point[1] = sin(t) * cos(t * a);
      point[2] = cos(t);
      PointsOnSurface.push_back(point);
    }
    ASSERT(PointsOnSurface.size() == N,
        "Sphere_impl::getHomogeneousPointsOnSurface() did not create "
        +::toString(N)+" but "+::toString(PointsOnSurface.size())+" points.");
  } else {
    // Symmetric but only approximately N points.
    double a=4*M_PI/N;
    double d= sqrt(a);
    int Mtheta=int(M_PI/d);
    double dtheta=M_PI/Mtheta;
    double dphi=a/dtheta;
    for (int m=0; m<Mtheta; m++)
    {
      double theta=M_PI*(m+0.5)/Mtheta;
      int Mphi=int(2*M_PI*sin(theta)/dphi);
      for (int n=0; n<Mphi;n++)
      {
        double phi= 2*M_PI*n/Mphi;
        Vector point;
        point.Zero();
        point[0]=sin(theta)*cos(phi);
        point[1]=sin(theta)*sin(phi);
        point[2]=cos(theta);
        PointsOnSurface.push_back(point);
      }
    }
  }
  return PointsOnSurface;
}

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

Shape Sphere(){
  Shape::impl_ptr impl = Shape::impl_ptr(new Sphere_impl());
  return Shape(impl);
}

Shape Sphere(const Vector &center,double radius){
  return translate(resize(Sphere(),radius),center);
}

Shape Ellipsoid(const Vector &center, const Vector &radius){
  return translate(stretch(Sphere(),radius),center);
}

bool Cuboid_impl::isInside(const Vector &point) const{
  return (point[0]>=0 && point[0]<=1) && (point[1]>=0 && point[1]<=1) && (point[2]>=0 && point[2]<=1);
}

bool Cuboid_impl::isOnSurface(const Vector &point) const{
  bool retVal = isInside(point);
  // test all borders of the cuboid
  // double fabs
  retVal = retVal &&
           (((fabs(point[0]-1.)  < MYEPSILON) || (fabs(point[0])  < MYEPSILON)) ||
            ((fabs(point[1]-1.)  < MYEPSILON) || (fabs(point[1])  < MYEPSILON)) ||
            ((fabs(point[2]-1.)  < MYEPSILON) || (fabs(point[2])  < MYEPSILON)));
  return retVal;
}

Vector Cuboid_impl::getNormal(const Vector &point) const throw(NotOnSurfaceException){
  if(!isOnSurface(point)){
    throw NotOnSurfaceException() << ShapeVector(&point);
  }
  Vector res;
  // figure out on which sides the Vector lies (maximum 3, when it is in a corner)
  for(int i=NDIM;i--;){
    if(fabs(fabs(point[i])-1)<MYEPSILON){
      // add the scaled (-1/+1) Vector to the set of surface vectors
      res[i] = point[i];
    }
  }
  ASSERT(res.NormSquared()>=1 && res.NormSquared()<=3,"To many or to few sides found for this Vector");

  res.Normalize();
  return res;
}


Vector Cuboid_impl::getCenter() const
{
  return Vector(0.5,0.5,0.5);
}

double Cuboid_impl::getRadius() const
{
  return .5;
}

double Cuboid_impl::getVolume() const
{
        return 1.; // l^3
}

double Cuboid_impl::getSurfaceArea() const
{
        return 6.;      // 6 * l^2
}

LineSegmentSet Cuboid_impl::getLineIntersections(const Line &line) const{
  LineSegmentSet res(line);
  // get the intersection on each of the six faces
  std::vector<Vector> intersections;
  intersections.resize(2);
  int c=0;
  int x[2]={-1,+1};
  for(int i=NDIM;i--;){
    for(int p=0;p<2;++p){
      if(c==2) goto end; // I know this sucks, but breaking two loops is stupid
      Vector base;
      base[i]=x[p];
      // base now points to the surface and is normal to it at the same time
      Plane p(base,base);
      Vector intersection = p.GetIntersection(line);
      if(isInside(intersection)){
        // if we have a point on the edge it might already be contained
        if(c==1 && intersections[0]==intersection)
          continue;
        intersections[c++]=intersection;
      }
    }
  }
  end:
  if(c==2){
    res.insert(LineSegment(intersections[0],intersections[1]));
  }
  return res;
}

std::string Cuboid_impl::toString() const{
  return "Cuboid()";
}

enum ShapeType Cuboid_impl::getType() const
{
        return CuboidType;
}

/**
 * \param N number of points on surface
 */
std::vector<Vector> Cuboid_impl::getHomogeneousPointsOnSurface(const size_t N) const {
  std::vector<Vector> PointsOnSurface;
  ASSERT(false, "Cuboid_impl::getHomogeneousPointsOnSurface() not implemented yet");
  return PointsOnSurface;
}

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

Shape Cuboid(){
  Shape::impl_ptr impl = Shape::impl_ptr(new Cuboid_impl());
  return Shape(impl);
}

Shape Cuboid(const Vector &corner1, const Vector &corner2){
  // make sure the two edges are upper left front and lower right back
  Vector sortedC1;
  Vector sortedC2;
  for(int i=NDIM;i--;){
    sortedC1[i] = std::min(corner1[i],corner2[i]);
    sortedC2[i] = std::max(corner1[i],corner2[i]);
    ASSERT(corner1[i]!=corner2[i],"Given points for cuboid edges did not define a valid space");
  }
  // get the middle point
  Vector middle = (1./2.)*(sortedC1+sortedC2);
  Vector factors = sortedC2-middle;
  return translate(stretch(Cuboid(),factors),middle);
}