| 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 |  * Shape.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 "CodePatterns/Assert.hpp"
 | 
|---|
| 23 | #include "LinearAlgebra/Vector.hpp"
 | 
|---|
| 24 | 
 | 
|---|
| 25 | #include "Shapes/Shape.hpp"
 | 
|---|
| 26 | #include "Shapes/Shape_impl.hpp"
 | 
|---|
| 27 | #include "Shapes/ShapeExceptions.hpp"
 | 
|---|
| 28 | 
 | 
|---|
| 29 | #include <string>
 | 
|---|
| 30 | 
 | 
|---|
| 31 | 
 | 
|---|
| 32 | Shape::Shape(const Shape& src) :
 | 
|---|
| 33 |   impl(src.getImpl())
 | 
|---|
| 34 | {}
 | 
|---|
| 35 | 
 | 
|---|
| 36 | Shape::~Shape(){}
 | 
|---|
| 37 | 
 | 
|---|
| 38 | bool Shape::isInside(const Vector &point) const{
 | 
|---|
| 39 |   return impl->isInside(point);
 | 
|---|
| 40 | }
 | 
|---|
| 41 | 
 | 
|---|
| 42 | bool Shape::isOnSurface(const Vector &point) const{
 | 
|---|
| 43 |   return impl->isOnSurface(point);
 | 
|---|
| 44 | }
 | 
|---|
| 45 | 
 | 
|---|
| 46 | Vector Shape::getNormal(const Vector &point) const throw (NotOnSurfaceException){
 | 
|---|
| 47 |   return impl->getNormal(point);
 | 
|---|
| 48 | }
 | 
|---|
| 49 | 
 | 
|---|
| 50 | LineSegmentSet Shape::getLineIntersections(const Line &line){
 | 
|---|
| 51 |   return impl->getLineIntersections(line);
 | 
|---|
| 52 | }
 | 
|---|
| 53 | 
 | 
|---|
| 54 | std::vector<Vector> Shape::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 55 |   return impl->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 56 | }
 | 
|---|
| 57 | 
 | 
|---|
| 58 | Shape::Shape(Shape::impl_ptr _impl) :
 | 
|---|
| 59 |     impl(_impl)
 | 
|---|
| 60 | {}
 | 
|---|
| 61 | 
 | 
|---|
| 62 | Shape &Shape::operator=(const Shape& rhs){
 | 
|---|
| 63 |   if(&rhs!=this){
 | 
|---|
| 64 |     impl=rhs.getImpl();
 | 
|---|
| 65 |   }
 | 
|---|
| 66 |   return *this;
 | 
|---|
| 67 | }
 | 
|---|
| 68 | 
 | 
|---|
| 69 | std::string Shape::toString() const{
 | 
|---|
| 70 |   return impl->toString();
 | 
|---|
| 71 | }
 | 
|---|
| 72 | 
 | 
|---|
| 73 | Shape::impl_ptr Shape::getImpl() const{
 | 
|---|
| 74 |   return impl;
 | 
|---|
| 75 | }
 | 
|---|
| 76 | 
 | 
|---|
| 77 | // allows arbitrary friendship, but only if implementation is known
 | 
|---|
| 78 | Shape::impl_ptr getShapeImpl(const Shape &shape){
 | 
|---|
| 79 |   return shape.getImpl();
 | 
|---|
| 80 | }
 | 
|---|
| 81 | 
 | 
|---|
| 82 | /***************************** Some simple Shapes ***************************/
 | 
|---|
| 83 | 
 | 
|---|
| 84 | Shape Everywhere(){
 | 
|---|
| 85 |   static Shape::impl_ptr impl = Shape::impl_ptr(new Everywhere_impl());
 | 
|---|
| 86 |   return Shape(impl);
 | 
|---|
| 87 | }
 | 
|---|
| 88 | 
 | 
|---|
| 89 | Shape Nowhere(){
 | 
|---|
| 90 |   static Shape::impl_ptr impl = Shape::impl_ptr(new Nowhere_impl());
 | 
|---|
| 91 |   return Shape(impl);
 | 
|---|
| 92 | }
 | 
|---|
| 93 | 
 | 
|---|
| 94 | /****************************** Operators ***********************************/
 | 
|---|
| 95 | 
 | 
|---|
| 96 | // AND
 | 
|---|
| 97 | 
 | 
|---|
| 98 | AndShape_impl::AndShape_impl(const Shape::impl_ptr &_lhs, const Shape::impl_ptr &_rhs) :
 | 
|---|
| 99 |   lhs(_lhs),rhs(_rhs)
 | 
|---|
| 100 | {}
 | 
|---|
| 101 | 
 | 
|---|
| 102 | AndShape_impl::~AndShape_impl(){}
 | 
|---|
| 103 | 
 | 
|---|
| 104 | bool AndShape_impl::isInside(const Vector &point){
 | 
|---|
| 105 |   return lhs->isInside(point) && rhs->isInside(point);
 | 
|---|
| 106 | }
 | 
|---|
| 107 | 
 | 
|---|
| 108 | bool AndShape_impl::isOnSurface(const Vector &point){
 | 
|---|
| 109 |   // check the number of surfaces that this point is on
 | 
|---|
| 110 |   int surfaces =0;
 | 
|---|
| 111 |   surfaces += lhs->isOnSurface(point);
 | 
|---|
| 112 |   surfaces += rhs->isOnSurface(point);
 | 
|---|
| 113 | 
 | 
|---|
| 114 |   switch(surfaces){
 | 
|---|
| 115 |     case 0:
 | 
|---|
| 116 |       return false;
 | 
|---|
| 117 |       // no break necessary
 | 
|---|
| 118 |     case 1:
 | 
|---|
| 119 |       // if it is inside for the object where it does not lie on
 | 
|---|
| 120 |       // the surface the whole point lies inside
 | 
|---|
| 121 |       return (lhs->isOnSurface(point) && rhs->isInside(point)) ||
 | 
|---|
| 122 |              (rhs->isOnSurface(point) && lhs->isInside(point));
 | 
|---|
| 123 |       // no break necessary
 | 
|---|
| 124 |     case 2:
 | 
|---|
| 125 |       {
 | 
|---|
| 126 |         // it lies on both Shapes... could be an edge or an inner point
 | 
|---|
| 127 |         // test the direction of the normals
 | 
|---|
| 128 |         Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
 | 
|---|
| 129 |         // if the directions are opposite we lie on the inside
 | 
|---|
| 130 |         return !direction.IsZero();
 | 
|---|
| 131 |       }
 | 
|---|
| 132 |       // no break necessary
 | 
|---|
| 133 |     default:
 | 
|---|
| 134 |       // if this happens there is something wrong
 | 
|---|
| 135 |       ASSERT(0,"Default case should have never been used");
 | 
|---|
| 136 |   }
 | 
|---|
| 137 |   return false; // never reached
 | 
|---|
| 138 | }
 | 
|---|
| 139 | 
 | 
|---|
| 140 | Vector AndShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
 | 
|---|
| 141 |   Vector res;
 | 
|---|
| 142 |   if(!isOnSurface(point)){
 | 
|---|
| 143 |     throw NotOnSurfaceException() << ShapeVector(&point);
 | 
|---|
| 144 |   }
 | 
|---|
| 145 |   res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
 | 
|---|
| 146 |   res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
 | 
|---|
| 147 |   res.Normalize();
 | 
|---|
| 148 |   return res;
 | 
|---|
| 149 | }
 | 
|---|
| 150 | 
 | 
|---|
| 151 | LineSegmentSet AndShape_impl::getLineIntersections(const Line &line){
 | 
|---|
| 152 |   return intersect(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
 | 
|---|
| 153 | }
 | 
|---|
| 154 | 
 | 
|---|
| 155 | string AndShape_impl::toString(){
 | 
|---|
| 156 |   return string("(") + lhs->toString() + string("&&") + rhs->toString() + string(")");
 | 
|---|
| 157 | }
 | 
|---|
| 158 | 
 | 
|---|
| 159 | std::vector<Vector> AndShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 160 |   std::vector<Vector> PointsOnSurface_lhs = lhs->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 161 |   std::vector<Vector> PointsOnSurface_rhs = rhs->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 162 |   std::vector<Vector> PointsOnSurface;
 | 
|---|
| 163 | 
 | 
|---|
| 164 |   for (std::vector<Vector>::const_iterator iter = PointsOnSurface_lhs.begin(); iter != PointsOnSurface_lhs.end(); ++iter) {
 | 
|---|
| 165 |     if (rhs->isInside(*iter))
 | 
|---|
| 166 |       PointsOnSurface.push_back(*iter);
 | 
|---|
| 167 |   }
 | 
|---|
| 168 |   for (std::vector<Vector>::const_iterator iter = PointsOnSurface_rhs.begin(); iter != PointsOnSurface_rhs.end(); ++iter) {
 | 
|---|
| 169 |     if (lhs->isInside(*iter))
 | 
|---|
| 170 |       PointsOnSurface.push_back(*iter);
 | 
|---|
| 171 |   }
 | 
|---|
| 172 | 
 | 
|---|
| 173 |   return PointsOnSurface;
 | 
|---|
| 174 | }
 | 
|---|
| 175 | 
 | 
|---|
| 176 | 
 | 
|---|
| 177 | Shape operator&&(const Shape &lhs,const Shape &rhs){
 | 
|---|
| 178 |   Shape::impl_ptr newImpl = Shape::impl_ptr(new AndShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
 | 
|---|
| 179 |   return Shape(newImpl);
 | 
|---|
| 180 | }
 | 
|---|
| 181 | 
 | 
|---|
| 182 | // OR
 | 
|---|
| 183 | 
 | 
|---|
| 184 | OrShape_impl::OrShape_impl(const Shape::impl_ptr &_lhs, const Shape::impl_ptr &_rhs) :
 | 
|---|
| 185 |   lhs(_lhs),rhs(_rhs)
 | 
|---|
| 186 | {}
 | 
|---|
| 187 | 
 | 
|---|
| 188 | OrShape_impl::~OrShape_impl(){}
 | 
|---|
| 189 | 
 | 
|---|
| 190 | bool OrShape_impl::isInside(const Vector &point){
 | 
|---|
| 191 |   return rhs->isInside(point) || lhs->isInside(point);
 | 
|---|
| 192 | }
 | 
|---|
| 193 | 
 | 
|---|
| 194 | bool OrShape_impl::isOnSurface(const Vector &point){
 | 
|---|
| 195 |   // check the number of surfaces that this point is on
 | 
|---|
| 196 |   int surfaces =0;
 | 
|---|
| 197 |   surfaces += lhs->isOnSurface(point);
 | 
|---|
| 198 |   surfaces += rhs->isOnSurface(point);
 | 
|---|
| 199 | 
 | 
|---|
| 200 |   switch(surfaces){
 | 
|---|
| 201 |     case 0:
 | 
|---|
| 202 |       return false;
 | 
|---|
| 203 |       // no break necessary
 | 
|---|
| 204 |     case 1:
 | 
|---|
| 205 |       // if it is inside for the object where it does not lie on
 | 
|---|
| 206 |       // the surface the whole point lies inside
 | 
|---|
| 207 |       return (lhs->isOnSurface(point) && !rhs->isInside(point)) ||
 | 
|---|
| 208 |              (rhs->isOnSurface(point) && !lhs->isInside(point));
 | 
|---|
| 209 |       // no break necessary
 | 
|---|
| 210 |     case 2:
 | 
|---|
| 211 |       {
 | 
|---|
| 212 |         // it lies on both Shapes... could be an edge or an inner point
 | 
|---|
| 213 |         // test the direction of the normals
 | 
|---|
| 214 |         Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
 | 
|---|
| 215 |         // if the directions are opposite we lie on the inside
 | 
|---|
| 216 |         return !direction.IsZero();
 | 
|---|
| 217 |       }
 | 
|---|
| 218 |       // no break necessary
 | 
|---|
| 219 |     default:
 | 
|---|
| 220 |       // if this happens there is something wrong
 | 
|---|
| 221 |       ASSERT(0,"Default case should have never been used");
 | 
|---|
| 222 |   }
 | 
|---|
| 223 |   return false; // never reached
 | 
|---|
| 224 | }
 | 
|---|
| 225 | 
 | 
|---|
| 226 | Vector OrShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
 | 
|---|
| 227 |   Vector res;
 | 
|---|
| 228 |   if(!isOnSurface(point)){
 | 
|---|
| 229 |     throw NotOnSurfaceException() << ShapeVector(&point);
 | 
|---|
| 230 |   }
 | 
|---|
| 231 |   res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
 | 
|---|
| 232 |   res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
 | 
|---|
| 233 |   res.Normalize();
 | 
|---|
| 234 |   return res;
 | 
|---|
| 235 | }
 | 
|---|
| 236 | 
 | 
|---|
| 237 | LineSegmentSet OrShape_impl::getLineIntersections(const Line &line){
 | 
|---|
| 238 |   return merge(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
 | 
|---|
| 239 | }
 | 
|---|
| 240 | 
 | 
|---|
| 241 | string OrShape_impl::toString(){
 | 
|---|
| 242 |   return string("(") + lhs->toString() + string("||") + rhs->toString() + string(")");
 | 
|---|
| 243 | }
 | 
|---|
| 244 | 
 | 
|---|
| 245 | std::vector<Vector> OrShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 246 |   std::vector<Vector> PointsOnSurface_lhs = lhs->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 247 |   std::vector<Vector> PointsOnSurface_rhs = rhs->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 248 |   std::vector<Vector> PointsOnSurface;
 | 
|---|
| 249 | 
 | 
|---|
| 250 |   for (std::vector<Vector>::const_iterator iter = PointsOnSurface_lhs.begin(); iter != PointsOnSurface_lhs.end(); ++iter) {
 | 
|---|
| 251 |     if (!rhs->isInside(*iter))
 | 
|---|
| 252 |       PointsOnSurface.push_back(*iter);
 | 
|---|
| 253 |   }
 | 
|---|
| 254 |   for (std::vector<Vector>::const_iterator iter = PointsOnSurface_rhs.begin(); iter != PointsOnSurface_rhs.end(); ++iter) {
 | 
|---|
| 255 |     if (!lhs->isInside(*iter))
 | 
|---|
| 256 |       PointsOnSurface.push_back(*iter);
 | 
|---|
| 257 |   }
 | 
|---|
| 258 | 
 | 
|---|
| 259 |   return PointsOnSurface;
 | 
|---|
| 260 | }
 | 
|---|
| 261 | 
 | 
|---|
| 262 | Shape operator||(const Shape &lhs,const Shape &rhs){
 | 
|---|
| 263 |   Shape::impl_ptr newImpl = Shape::impl_ptr(new OrShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
 | 
|---|
| 264 |   return Shape(newImpl);
 | 
|---|
| 265 | }
 | 
|---|
| 266 | 
 | 
|---|
| 267 | // NOT
 | 
|---|
| 268 | 
 | 
|---|
| 269 | NotShape_impl::NotShape_impl(const Shape::impl_ptr &_arg) :
 | 
|---|
| 270 |   arg(_arg)
 | 
|---|
| 271 | {}
 | 
|---|
| 272 | 
 | 
|---|
| 273 | NotShape_impl::~NotShape_impl(){}
 | 
|---|
| 274 | 
 | 
|---|
| 275 | bool NotShape_impl::isInside(const Vector &point){
 | 
|---|
| 276 |   return !arg->isInside(point);
 | 
|---|
| 277 | }
 | 
|---|
| 278 | 
 | 
|---|
| 279 | bool NotShape_impl::isOnSurface(const Vector &point){
 | 
|---|
| 280 |   return arg->isOnSurface(point);
 | 
|---|
| 281 | }
 | 
|---|
| 282 | 
 | 
|---|
| 283 | Vector NotShape_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
 | 
|---|
| 284 |   return -1.*arg->getNormal(point);
 | 
|---|
| 285 | }
 | 
|---|
| 286 | 
 | 
|---|
| 287 | LineSegmentSet NotShape_impl::getLineIntersections(const Line &line){
 | 
|---|
| 288 |   return invert(arg->getLineIntersections(line));
 | 
|---|
| 289 | }
 | 
|---|
| 290 | 
 | 
|---|
| 291 | string NotShape_impl::toString(){
 | 
|---|
| 292 |   return string("!") + arg->toString();
 | 
|---|
| 293 | }
 | 
|---|
| 294 | 
 | 
|---|
| 295 | std::vector<Vector> NotShape_impl::getHomogeneousPointsOnSurface(const size_t N) const {
 | 
|---|
| 296 |   // surfaces are the same, only normal direction is different
 | 
|---|
| 297 |   return arg->getHomogeneousPointsOnSurface(N);
 | 
|---|
| 298 | }
 | 
|---|
| 299 | 
 | 
|---|
| 300 | Shape operator!(const Shape &arg){
 | 
|---|
| 301 |   Shape::impl_ptr newImpl = Shape::impl_ptr(new NotShape_impl(getShapeImpl(arg)));
 | 
|---|
| 302 |   return Shape(newImpl);
 | 
|---|
| 303 | }
 | 
|---|
| 304 | 
 | 
|---|
| 305 | /**************** global operations *********************************/
 | 
|---|
| 306 | ostream &operator<<(ostream &ost,const Shape &shape){
 | 
|---|
| 307 |   ost << shape.toString();
 | 
|---|
| 308 |   return ost;
 | 
|---|
| 309 | }
 | 
|---|