source: src/Shapes/ShapeOps.cpp@ 987145

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. 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/>.
 */

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

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

//#include "CodePatterns/MemDebug.hpp"

#include <algorithm>
#include <boost/bind.hpp>
#include <cmath>

#include "Shapes/ShapeExceptions.hpp"
#include "Shapes/ShapeOps.hpp"
#include "Shapes/ShapeOps_impl.hpp"

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

/*************** Base case ***********************/

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

ShapeOpsBase_impl::~ShapeOpsBase_impl(){}

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

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

Vector ShapeOpsBase_impl::getNormal(const Vector &point) const throw (NotOnSurfaceException){
  Vector helper = translateIn(point);
  if(!arg->isOnSurface(helper)){
    throw NotOnSurfaceException() << ShapeVector(&helper);
  }
  Vector res = translateOutNormal(arg->getNormal(helper));
  res.Normalize();
  return res;
}

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

double ShapeOpsBase_impl::getRadius() const
{
  return translateOutPos(Vector(arg->getRadius(), 0., 0.)).Norm();
}


LineSegmentSet ShapeOpsBase_impl::getLineIntersections(const Line &line) const{
  Line newLine(translateIn(line.getOrigin()),translateIn(line.getDirection()));
  LineSegmentSet res(line);
  LineSegmentSet helper = getArg()->getLineIntersections(newLine);
  for(LineSegmentSet::iterator iter = helper.begin();iter!=helper.end();++iter){
    LinePoint lpBegin = iter->getBegin();
    LinePoint lpEnd = iter->getBegin();
    // translate both linepoints
    lpBegin = lpBegin.isNegInfinity()?
                line.negEndpoint():
                line.getLinePoint(translateOutPos(lpBegin.getPoint()));
    lpEnd = lpEnd.isPosInfinity()?
              line.posEndpoint():
              line.getLinePoint(translateOutPos(lpEnd.getPoint()));
    res.insert(LineSegment(lpBegin,lpEnd));
  }
  return res;
}

enum ShapeType ShapeOpsBase_impl::getType() const {
        return getArg()->getType();
}

std::vector<Vector> ShapeOpsBase_impl::getHomogeneousPointsOnSurface(const size_t N) const {
  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
  std::transform(PointsOnSurface.begin(), PointsOnSurface.end(), PointsOnSurface.begin(),
      boost::bind(&ShapeOpsBase_impl::translateOutPos, this, _1) );
  return PointsOnSurface;
}

std::vector<Vector> ShapeOpsBase_impl::getHomogeneousPointsInVolume(const size_t N) const {
  std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsInVolume(N);
  std::transform(PointsOnSurface.begin(), PointsOnSurface.end(), PointsOnSurface.begin(),
      boost::bind(&ShapeOpsBase_impl::translateOutPos, this, _1) );
  return PointsOnSurface;
}

Shape::impl_ptr ShapeOpsBase_impl::getArg() const{
  return arg;
}

/********************* Resize ********************/

Resize_impl::Resize_impl(const Shape::impl_ptr &_arg,double _size) :
  ShapeOpsBase_impl(_arg), size(_size)
{
  ASSERT(size>0,"Cannot resize a Shape to size zero or below");
}

Resize_impl::~Resize_impl(){}

double Resize_impl::getVolume() const
{
        return getArg()->getVolume() * size * size * size;
}

double Resize_impl::getSurfaceArea() const
{
        return getArg()->getSurfaceArea() * size * size;
}


bool Resize_impl::isInside(const Vector& point) const{
  return getArg()->isInside((1./size) * point);
}

Vector Resize_impl::translateIn(const Vector& point) const{
  return (1./size) * point;
}

Vector Resize_impl::translateOutPos(const Vector& point) const{
  return size * point;
}

Vector Resize_impl::translateOutNormal(const Vector& point) const{
  return point;
}

std::string Resize_impl::toString() const{
  std::stringstream sstr;
  sstr << "resize(" << getArg()->toString() << "," << size << ")";
  return sstr.str();
}

Shape resize(const Shape &arg,double size){
  Shape::impl_ptr impl = Shape::impl_ptr(new Resize_impl(getShapeImpl(arg),size));
  return Shape(impl);
}

/*************************** translate *******************/

Translate_impl::Translate_impl(const Shape::impl_ptr &_arg, const Vector &_offset) :
  ShapeOpsBase_impl(_arg),offset(_offset)
{}

Translate_impl::~Translate_impl(){}

bool Translate_impl::isInside(const Vector& point) const{
  return getArg()->isInside(point-offset);
}

Vector Translate_impl::getCenter() const
{
  return getArg()->getCenter()+offset;
}

double Translate_impl::getRadius() const
{
  return getArg()->getRadius();
}

double Translate_impl::getVolume() const
{
        return getArg()->getVolume();
}

double Translate_impl::getSurfaceArea() const
{
        return getArg()->getSurfaceArea();
}

Vector Translate_impl::translateIn(const Vector& point) const{
  return point-offset;
}

Vector Translate_impl::translateOutPos(const Vector& point) const{
  return point+offset;
}

Vector Translate_impl::translateOutNormal(const Vector& point) const{
  return point;
}

std::string Translate_impl::toString() const{
  std::stringstream sstr;
  sstr << "translate(" << getArg()->toString() << "," << offset << ")";
  return sstr.str();
}

Shape translate(const Shape &arg, const Vector &offset){
  Shape::impl_ptr impl = Shape::impl_ptr(new Translate_impl(getShapeImpl(arg),offset));
  return Shape(impl);
}

/*********************** stretch ******************/

Stretch_impl::Stretch_impl(const Shape::impl_ptr &_arg, const Vector &_factors) :
  ShapeOpsBase_impl(_arg),factors(_factors)
{
  for(int i = NDIM;i--;){
    ASSERT(factors[i]>0.,"cannot stretch a shape by a negative amount");
    reciFactors[i] = 1./factors[i];
  }
}

Stretch_impl::~Stretch_impl(){}

double Stretch_impl::getVolume() const
{
        // TODO
        return getArg()->getVolume() * factors[0] * factors[1]  * factors[2];
}

double Stretch_impl::getSurfaceArea() const
{
        // TODO
        return getArg()->getSurfaceArea() * pow(factors[0] * factors[1]  * factors[2], 2./3.);
}

bool Stretch_impl::isInside(const Vector& point) const{
  Vector helper=point;
  helper.ScaleAll(reciFactors);
  return getArg()->isInside(helper);
}

Vector Stretch_impl::translateIn(const Vector& point) const{
  Vector helper=point;
  helper.ScaleAll(reciFactors);
  return helper;
}

Vector Stretch_impl::translateOutPos(const Vector& point) const{
  Vector helper=point;
  helper.ScaleAll(factors);
  return helper;
}

Vector Stretch_impl::translateOutNormal(const Vector& point) const{
  Vector helper=point;
  // the normalFactors are derived from appearances of the factors
  // with in the vectorproduct
  Vector normalFactors;
  normalFactors[0]=factors[1]*factors[2];
  normalFactors[1]=factors[0]*factors[2];
  normalFactors[2]=factors[0]*factors[1];
  helper.ScaleAll(normalFactors);
  return helper;
}

std::string Stretch_impl::toString() const{
  std::stringstream sstr;
  sstr << "stretch(" << getArg()->toString() << "," << factors << ")";
  return sstr.str();
}

Shape stretch(const Shape &arg, const Vector &factors){
  Shape::impl_ptr impl = Shape::impl_ptr(new Stretch_impl(getShapeImpl(arg),factors));
  return Shape(impl);
}

/************************* transform *****************/

Transform_impl::Transform_impl(const Shape::impl_ptr &_arg, const RealSpaceMatrix &_transformation) :
  ShapeOpsBase_impl(_arg),transformation(_transformation)
{
  transformationInv = transformation.invert();
}

Transform_impl::~Transform_impl(){}

double Transform_impl::getVolume() const
{
        return getArg()->getVolume() * transformation.determinant();
}

double Transform_impl::getSurfaceArea() const
{
        return getArg()->getSurfaceArea() * pow(transformation.determinant(), 2./3.);
}

bool Transform_impl::isInside(const Vector& point) const{
  return getArg()->isInside(transformationInv * point);
}

Vector Transform_impl::translateIn(const Vector& point) const{
  return transformationInv * point;
}

Vector Transform_impl::translateOutPos(const Vector& point) const{
  return transformation * point;
}

Vector Transform_impl::translateOutNormal(const Vector& point) const
{
  RealSpaceMatrix mat = transformation.invert().transpose();
  return mat * point;
}

std::string Transform_impl::toString() const{
  std::stringstream sstr;
  sstr << "transform(" << getArg()->toString() << "," << transformation << ")";
  return sstr.str();
}

Shape transform(const Shape &arg, const RealSpaceMatrix &transformation){
  Shape::impl_ptr impl = Shape::impl_ptr(new Transform_impl(getShapeImpl(arg),transformation));
  return Shape(impl);
}