| 1 | /*
 | 
|---|
| 2 |  * Project: MoleCuilder
 | 
|---|
| 3 |  * Description: creates and alters molecular systems
 | 
|---|
| 4 |  * Copyright (C)  2010 University of Bonn. All rights reserved.
 | 
|---|
| 5 |  * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 | 
|---|
| 6 |  */
 | 
|---|
| 7 | 
 | 
|---|
| 8 | /*
 | 
|---|
| 9 |  * ShapeOps.cpp
 | 
|---|
| 10 |  *
 | 
|---|
| 11 |  *  Created on: Jun 18, 2010
 | 
|---|
| 12 |  *      Author: crueger
 | 
|---|
| 13 |  */
 | 
|---|
| 14 | 
 | 
|---|
| 15 | // include config.h
 | 
|---|
| 16 | #ifdef HAVE_CONFIG_H
 | 
|---|
| 17 | #include <config.h>
 | 
|---|
| 18 | #endif
 | 
|---|
| 19 | 
 | 
|---|
| 20 | #include "CodePatterns/MemDebug.hpp"
 | 
|---|
| 21 | 
 | 
|---|
| 22 | #include "Shapes/ShapeExceptions.hpp"
 | 
|---|
| 23 | #include "Shapes/ShapeOps.hpp"
 | 
|---|
| 24 | #include "Shapes/ShapeOps_impl.hpp"
 | 
|---|
| 25 | 
 | 
|---|
| 26 | #include "LinearAlgebra/Vector.hpp"
 | 
|---|
| 27 | #include "CodePatterns/Assert.hpp"
 | 
|---|
| 28 | 
 | 
|---|
| 29 | /*************** Base case ***********************/
 | 
|---|
| 30 | 
 | 
|---|
| 31 | ShapeOpsBase_impl::ShapeOpsBase_impl(const Shape::impl_ptr &_arg) :
 | 
|---|
| 32 |   arg(_arg){}
 | 
|---|
| 33 | 
 | 
|---|
| 34 | ShapeOpsBase_impl::~ShapeOpsBase_impl(){}
 | 
|---|
| 35 | 
 | 
|---|
| 36 | bool ShapeOpsBase_impl::isInside(const Vector &point){
 | 
|---|
| 37 |   return arg->isInside(translateIn(point));
 | 
|---|
| 38 | }
 | 
|---|
| 39 | 
 | 
|---|
| 40 | bool ShapeOpsBase_impl::isOnSurface(const Vector &point){
 | 
|---|
| 41 |   return arg->isOnSurface(translateIn(point));
 | 
|---|
| 42 | }
 | 
|---|
| 43 | 
 | 
|---|
| 44 | Vector ShapeOpsBase_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
 | 
|---|
| 45 |   Vector helper = translateIn(point);
 | 
|---|
| 46 |   if(!arg->isOnSurface(helper)){
 | 
|---|
| 47 |     throw NotOnSurfaceException() << ShapeVector(&helper);
 | 
|---|
| 48 |   }
 | 
|---|
| 49 |   Vector res = translateOutNormal(arg->getNormal(helper));
 | 
|---|
| 50 |   res.Normalize();
 | 
|---|
| 51 |   return res;
 | 
|---|
| 52 | }
 | 
|---|
| 53 | 
 | 
|---|
| 54 | LineSegmentSet ShapeOpsBase_impl::getLineIntersections(const Line &line){
 | 
|---|
| 55 |   Line newLine(translateIn(line.getOrigin()),translateIn(line.getDirection()));
 | 
|---|
| 56 |   LineSegmentSet res(line);
 | 
|---|
| 57 |   LineSegmentSet helper = getArg()->getLineIntersections(newLine);
 | 
|---|
| 58 |   for(LineSegmentSet::iterator iter = helper.begin();iter!=helper.end();++iter){
 | 
|---|
| 59 |     LinePoint lpBegin = iter->getBegin();
 | 
|---|
| 60 |     LinePoint lpEnd = iter->getBegin();
 | 
|---|
| 61 |     // translate both linepoints
 | 
|---|
| 62 |     lpBegin = lpBegin.isNegInfinity()?
 | 
|---|
| 63 |                 line.negEndpoint():
 | 
|---|
| 64 |                 line.getLinePoint(translateOutPos(lpBegin.getPoint()));
 | 
|---|
| 65 |     lpEnd = lpEnd.isPosInfinity()?
 | 
|---|
| 66 |               line.posEndpoint():
 | 
|---|
| 67 |               line.getLinePoint(translateOutPos(lpEnd.getPoint()));
 | 
|---|
| 68 |     res.insert(LineSegment(lpBegin,lpEnd));
 | 
|---|
| 69 |   }
 | 
|---|
| 70 |   return res;
 | 
|---|
| 71 | }
 | 
|---|
| 72 | 
 | 
|---|
| 73 | std::vector<Vector> ShapeOpsBase_impl::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 74 |   return getArg()->getHomogeneousPointsOnSurface(N);;
 | 
|---|
| 75 | }
 | 
|---|
| 76 | 
 | 
|---|
| 77 | Shape::impl_ptr ShapeOpsBase_impl::getArg() const{
 | 
|---|
| 78 |   return arg;
 | 
|---|
| 79 | }
 | 
|---|
| 80 | 
 | 
|---|
| 81 | /********************* Resize ********************/
 | 
|---|
| 82 | 
 | 
|---|
| 83 | Resize_impl::Resize_impl(const Shape::impl_ptr &_arg,double _size) :
 | 
|---|
| 84 |   ShapeOpsBase_impl(_arg), size(_size)
 | 
|---|
| 85 | {
 | 
|---|
| 86 |   ASSERT(size>0,"Cannot resize a Shape to size zero or below");
 | 
|---|
| 87 | }
 | 
|---|
| 88 | 
 | 
|---|
| 89 | Resize_impl::~Resize_impl(){}
 | 
|---|
| 90 | 
 | 
|---|
| 91 | bool Resize_impl::isInside(const Vector& point){
 | 
|---|
| 92 |   return getArg()->isInside((1/size) * point);
 | 
|---|
| 93 | }
 | 
|---|
| 94 | 
 | 
|---|
| 95 | Vector Resize_impl::translateIn(const Vector& point){
 | 
|---|
| 96 |   return (1/size) * point;
 | 
|---|
| 97 | }
 | 
|---|
| 98 | 
 | 
|---|
| 99 | Vector Resize_impl::translateOutPos(const Vector& point){
 | 
|---|
| 100 |   return size * point;
 | 
|---|
| 101 | }
 | 
|---|
| 102 | 
 | 
|---|
| 103 | Vector Resize_impl::translateOutNormal(const Vector& point){
 | 
|---|
| 104 |   return point;
 | 
|---|
| 105 | }
 | 
|---|
| 106 | 
 | 
|---|
| 107 | string Resize_impl::toString(){
 | 
|---|
| 108 |   stringstream sstr;
 | 
|---|
| 109 |   sstr << "resize(" << getArg()->toString() << "," << size << ")";
 | 
|---|
| 110 |   return sstr.str();
 | 
|---|
| 111 | }
 | 
|---|
| 112 | 
 | 
|---|
| 113 | std::vector<Vector> Resize_impl::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 114 |   std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 115 |   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
 | 
|---|
| 116 |     *iter *= size;
 | 
|---|
| 117 |   }
 | 
|---|
| 118 |   return PointsOnSurface;
 | 
|---|
| 119 | }
 | 
|---|
| 120 | 
 | 
|---|
| 121 | 
 | 
|---|
| 122 | Shape resize(const Shape &arg,double size){
 | 
|---|
| 123 |   Shape::impl_ptr impl = Shape::impl_ptr(new Resize_impl(getShapeImpl(arg),size));
 | 
|---|
| 124 |   return Shape(impl);
 | 
|---|
| 125 | }
 | 
|---|
| 126 | 
 | 
|---|
| 127 | /*************************** translate *******************/
 | 
|---|
| 128 | 
 | 
|---|
| 129 | Translate_impl::Translate_impl(const Shape::impl_ptr &_arg, const Vector &_offset) :
 | 
|---|
| 130 |   ShapeOpsBase_impl(_arg),offset(_offset)
 | 
|---|
| 131 | {}
 | 
|---|
| 132 | 
 | 
|---|
| 133 | Translate_impl::~Translate_impl(){}
 | 
|---|
| 134 | 
 | 
|---|
| 135 | bool Translate_impl::isInside(const Vector& point){
 | 
|---|
| 136 |   return getArg()->isInside(point-offset);
 | 
|---|
| 137 | }
 | 
|---|
| 138 | 
 | 
|---|
| 139 | Vector Translate_impl::translateIn(const Vector& point){
 | 
|---|
| 140 |   return point-offset;
 | 
|---|
| 141 | }
 | 
|---|
| 142 | 
 | 
|---|
| 143 | Vector Translate_impl::translateOutPos(const Vector& point){
 | 
|---|
| 144 |   return point+offset;
 | 
|---|
| 145 | }
 | 
|---|
| 146 | 
 | 
|---|
| 147 | Vector Translate_impl::translateOutNormal(const Vector& point){
 | 
|---|
| 148 |   return point;
 | 
|---|
| 149 | }
 | 
|---|
| 150 | 
 | 
|---|
| 151 | string Translate_impl::toString(){
 | 
|---|
| 152 |   stringstream sstr;
 | 
|---|
| 153 |   sstr << "translate(" << getArg()->toString() << "," << offset << ")";
 | 
|---|
| 154 |   return sstr.str();
 | 
|---|
| 155 | }
 | 
|---|
| 156 | 
 | 
|---|
| 157 | std::vector<Vector> Translate_impl::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 158 |   std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 159 |   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
 | 
|---|
| 160 |     *iter += offset;
 | 
|---|
| 161 |   }
 | 
|---|
| 162 |   return PointsOnSurface;
 | 
|---|
| 163 | }
 | 
|---|
| 164 | 
 | 
|---|
| 165 | Shape translate(const Shape &arg, const Vector &offset){
 | 
|---|
| 166 |   Shape::impl_ptr impl = Shape::impl_ptr(new Translate_impl(getShapeImpl(arg),offset));
 | 
|---|
| 167 |   return Shape(impl);
 | 
|---|
| 168 | }
 | 
|---|
| 169 | 
 | 
|---|
| 170 | /*********************** stretch ******************/
 | 
|---|
| 171 | 
 | 
|---|
| 172 | Stretch_impl::Stretch_impl(const Shape::impl_ptr &_arg, const Vector &_factors) :
 | 
|---|
| 173 |   ShapeOpsBase_impl(_arg),factors(_factors)
 | 
|---|
| 174 | {
 | 
|---|
| 175 |   ASSERT(factors[0]>0,"cannot stretch a shape by a negative amount");
 | 
|---|
| 176 |   ASSERT(factors[1]>0,"cannot stretch a shape by a negative amount");
 | 
|---|
| 177 |   ASSERT(factors[2]>0,"cannot stretch a shape by a negative amount");
 | 
|---|
| 178 |   for(int i = NDIM;i--;){
 | 
|---|
| 179 |     reciFactors[i] = 1/factors[i];
 | 
|---|
| 180 |   }
 | 
|---|
| 181 | }
 | 
|---|
| 182 | 
 | 
|---|
| 183 | Stretch_impl::~Stretch_impl(){}
 | 
|---|
| 184 | 
 | 
|---|
| 185 | bool Stretch_impl::isInside(const Vector& point){
 | 
|---|
| 186 |   Vector helper=point;
 | 
|---|
| 187 |   helper.ScaleAll(reciFactors);
 | 
|---|
| 188 |   return getArg()->isInside(helper);
 | 
|---|
| 189 | }
 | 
|---|
| 190 | 
 | 
|---|
| 191 | Vector Stretch_impl::translateIn(const Vector& point){
 | 
|---|
| 192 |   Vector helper=point;
 | 
|---|
| 193 |   helper.ScaleAll(reciFactors);
 | 
|---|
| 194 |   return helper;
 | 
|---|
| 195 | }
 | 
|---|
| 196 | 
 | 
|---|
| 197 | Vector Stretch_impl::translateOutPos(const Vector& point){
 | 
|---|
| 198 |   Vector helper=point;
 | 
|---|
| 199 |   helper.ScaleAll(factors);
 | 
|---|
| 200 |   return helper;
 | 
|---|
| 201 | }
 | 
|---|
| 202 | 
 | 
|---|
| 203 | Vector Stretch_impl::translateOutNormal(const Vector& point){
 | 
|---|
| 204 |   Vector helper=point;
 | 
|---|
| 205 |   // the normalFactors are derived from appearances of the factors
 | 
|---|
| 206 |   // with in the vectorproduct
 | 
|---|
| 207 |   Vector normalFactors;
 | 
|---|
| 208 |   normalFactors[0]=factors[1]*factors[2];
 | 
|---|
| 209 |   normalFactors[1]=factors[0]*factors[2];
 | 
|---|
| 210 |   normalFactors[2]=factors[0]*factors[1];
 | 
|---|
| 211 |   helper.ScaleAll(normalFactors);
 | 
|---|
| 212 |   return helper;
 | 
|---|
| 213 | }
 | 
|---|
| 214 | 
 | 
|---|
| 215 | string Stretch_impl::toString(){
 | 
|---|
| 216 |   stringstream sstr;
 | 
|---|
| 217 |   sstr << "stretch(" << getArg()->toString() << "," << factors << ")";
 | 
|---|
| 218 |   return sstr.str();
 | 
|---|
| 219 | }
 | 
|---|
| 220 | 
 | 
|---|
| 221 | std::vector<Vector> Stretch_impl::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 222 |   std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 223 |   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
 | 
|---|
| 224 |     (*iter).ScaleAll(reciFactors);
 | 
|---|
| 225 |   }
 | 
|---|
| 226 |   return PointsOnSurface;
 | 
|---|
| 227 | }
 | 
|---|
| 228 | 
 | 
|---|
| 229 | Shape stretch(const Shape &arg, const Vector &factors){
 | 
|---|
| 230 |   Shape::impl_ptr impl = Shape::impl_ptr(new Stretch_impl(getShapeImpl(arg),factors));
 | 
|---|
| 231 |   return Shape(impl);
 | 
|---|
| 232 | }
 | 
|---|
| 233 | 
 | 
|---|
| 234 | /************************* transform *****************/
 | 
|---|
| 235 | 
 | 
|---|
| 236 | Transform_impl::Transform_impl(const Shape::impl_ptr &_arg, const RealSpaceMatrix &_transformation) :
 | 
|---|
| 237 |   ShapeOpsBase_impl(_arg),transformation(_transformation)
 | 
|---|
| 238 | {
 | 
|---|
| 239 |   transformationInv = transformation.invert();
 | 
|---|
| 240 | }
 | 
|---|
| 241 | 
 | 
|---|
| 242 | Transform_impl::~Transform_impl(){}
 | 
|---|
| 243 | 
 | 
|---|
| 244 | bool Transform_impl::isInside(const Vector& point){
 | 
|---|
| 245 |   return getArg()->isInside(transformationInv * point);
 | 
|---|
| 246 | }
 | 
|---|
| 247 | 
 | 
|---|
| 248 | Vector Transform_impl::translateIn(const Vector& point){
 | 
|---|
| 249 |   return transformationInv * point;
 | 
|---|
| 250 | }
 | 
|---|
| 251 | 
 | 
|---|
| 252 | Vector Transform_impl::translateOutPos(const Vector& point){
 | 
|---|
| 253 |   return transformation * point;
 | 
|---|
| 254 | }
 | 
|---|
| 255 | 
 | 
|---|
| 256 | Vector Transform_impl::translateOutNormal(const Vector& point){
 | 
|---|
| 257 |   RealSpaceMatrix mat = transformation.invert().transpose();
 | 
|---|
| 258 |   return mat * point;
 | 
|---|
| 259 | }
 | 
|---|
| 260 | 
 | 
|---|
| 261 | string Transform_impl::toString(){
 | 
|---|
| 262 |   stringstream sstr;
 | 
|---|
| 263 |   sstr << "transform(" << getArg()->toString() << "," << transformation << ")";
 | 
|---|
| 264 |   return sstr.str();
 | 
|---|
| 265 | }
 | 
|---|
| 266 | 
 | 
|---|
| 267 | std::vector<Vector> Transform_impl::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 268 |   std::vector<Vector> PointsOnSurface = getArg()->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 269 |   for(std::vector<Vector>::iterator iter = PointsOnSurface.begin(); iter != PointsOnSurface.end(); ++iter) {
 | 
|---|
| 270 |     *iter = transformation * (*iter);
 | 
|---|
| 271 |   }
 | 
|---|
| 272 |   return PointsOnSurface;
 | 
|---|
| 273 | }
 | 
|---|
| 274 | 
 | 
|---|
| 275 | Shape transform(const Shape &arg, const RealSpaceMatrix &transformation){
 | 
|---|
| 276 |   Shape::impl_ptr impl = Shape::impl_ptr(new Transform_impl(getShapeImpl(arg),transformation));
 | 
|---|
| 277 |   return Shape(impl);
 | 
|---|
| 278 | }
 | 
|---|