source: src/Shapes/BaseShapes.cpp@ c521e1

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 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){
  return point.NormSquared()<=1;
}

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

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

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

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

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

    double PI=3.14159265358979323846;
    double a=4*PI/N;
    double d= sqrt(a);
    int Mtheta=int(PI/d);
    double dtheta=PI/Mtheta;
    double dphi=a/dtheta;
    for (int m=0; m<Mtheta; m++)
    {
          double theta=PI*(m+0.5)/Mtheta;
          int Mphi=int(2*PI*sin(theta)/dphi);
          for (int n=0; n<Mphi;n++)
          {
                  double phi= 2*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);
          }
    }
  /*const double dlength = M_PI*(3.-sqrt(5.));
  double length = 0;
  const double dz = 2.0/N;
  double z = 1. - dz/2.;
  Vector point;
  for (size_t ka = 0; ka<N; ka++){
    const double r = sqrt(1.-z*z);
    point.Zero();
    point[0] = cos(length)*r;
    point[1] = sin(length)*r;
    point[2] = z;
    PointsOnSurface.push_back(point);
    z = z - dz;
    length = length + dlength;*/

  //ASSERT(PointsOnSurface.size() == N, "Sphere_impl::getHomogeneousPointsOnSurface() did not create enough points.");
  return PointsOnSurface;
}


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){
  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){
  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) 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;
}

LineSegmentSet Cuboid_impl::getLineIntersections(const Line &line){
  LineSegmentSet res(line);
  // get the intersection on each of the six faces
  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;
}

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

/**
 * \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;
}

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] = min(corner1[i],corner2[i]);
    sortedC2[i] = 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);
}