src/GeomAPI/GeomAPI_Shell.cpp

   1 // Copyright (C) 2018-20xx  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or
  18 // email : webmaster.salome@opencascade.com<mailto:webmaster.salome@opencascade.com>
  19 //
  20
  21 #include "GeomAPI_Shell.h"
  22
  23 #include "GeomAPI_Ax3.h"
  24 #include "GeomAPI_Box.h"
  25 #include "GeomAPI_Cone.h"
  26 #include "GeomAPI_Cylinder.h"
  27 #include "GeomAPI_Face.h"
  28 #include "GeomAPI_Pnt.h"
  29 #include "GeomAPI_Sphere.h"
  30 #include "GeomAPI_Torus.h"
  31 #include "GeomAPI_Wire.h"
  32 #include "GeomAPI_XYZ.h"
  33
  34 #include <BRep_Builder.hxx>
  35 #include <Precision.hxx>
  36 #include <TopExp_Explorer.hxx>
  37 #include <TopoDS.hxx>
  38
  39 #include <map>
  40
  41 //=================================================================================================
  42 GeomAPI_Shell::GeomAPI_Shell()
  43 {
  44   TopoDS_Shell* aShell = new TopoDS_Shell();
  45
  46   BRep_Builder aBuilder;
  47   aBuilder.MakeShell(*aShell);
  48
  49   this->setImpl(aShell);
  50 }
  51
  52 //=================================================================================================
  53 GeomAPI_Shell::GeomAPI_Shell(const std::shared_ptr<GeomAPI_Shape>& theShape)
  54 {
  55   if (!theShape->isNull() && theShape->isShell()) {
  56     setImpl(new TopoDS_Shape(theShape->impl<TopoDS_Shape>()));
  57   }
  58 }
  59
  60 //=================================================================================================
  61 std::shared_ptr<GeomAPI_Sphere> GeomAPI_Shell::getSphere() const
  62 {
  63   bool isSphere = true;
  64   bool isFirstFace = true;
  65   GeomSpherePtr aSphere;
  66
  67   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_FACE); anExp.More(); anExp.Next()) {
  68     GeomFacePtr aFace(new GeomAPI_Face);
  69     aFace->setImpl(new TopoDS_Shape(anExp.Current()));
  70
  71     GeomSpherePtr aCurSphere = aFace->getSphere();
  72     if (!aCurSphere) {
  73       isSphere = false;
  74       break;
  75     }
  76
  77     if (isFirstFace) {
  78       aSphere = aCurSphere;
  79       isFirstFace = false;
  80     }
  81     else if (aSphere->center()->distance(aCurSphere->center()) >= Precision::Confusion() ||
  82              Abs(aSphere->radius() - aCurSphere->radius()) >= Precision::Confusion()) {
  83       isSphere = false;
  84       break;
  85     }
  86   }
  87
  88   return isSphere ? aSphere : GeomSpherePtr();
  89 }
  90
  91 //=================================================================================================
  92 std::shared_ptr<GeomAPI_Cylinder> GeomAPI_Shell::getCylinder() const
  93 {
  94   bool isCylinder = true;
  95   bool isFirstFace = true;
  96
  97   GeomPointPtr aLocation;
  98   GeomDirPtr anAxis;
  99   double aRadius;
 100   double aHeight;
 101
 102   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_FACE); anExp.More(); anExp.Next()) {
 103     GeomFacePtr aFace(new GeomAPI_Face);
 104     aFace->setImpl(new TopoDS_Shape(anExp.Current()));
 105
 106     GeomCylinderPtr aCurCyl = aFace->getCylinder();
 107     if (!aCurCyl) {
 108       isCylinder = false;
 109       break;
 110     }
 111
 112     if (isFirstFace) {
 113       aLocation = aCurCyl->location();
 114       anAxis = aCurCyl->axis();
 115       aRadius = aCurCyl->radius();
 116       aHeight = aCurCyl->height();
 117       isFirstFace = false;
 118     }
 119     else {
 120       // compare radii
 121       if (Abs(aRadius - aCurCyl->radius()) >= Precision::Confusion() ||
 122         // check directions are collinear
 123         !anAxis->isParallel(aCurCyl->axis()) ||
 124         // check current center is on the main axis
 125         anAxis->xyz()->cross(aLocation->xyz()->decreased(aCurCyl->location()->xyz())
 126         )->squareModulus() >= Precision::SquareConfusion()) {
 127         isCylinder = false;
 128         break;
 129       }
 130
 131       double aMinHeight = 0.0;
 132       double aMaxHeight = aHeight;
 133
 134       std::shared_ptr<GeomAPI_XYZ> aCurCylLoc = aCurCyl->location()->xyz();
 135       std::shared_ptr<GeomAPI_XYZ> aCurCylLocHeight =
 136         aCurCylLoc->added(aCurCyl->axis()->xyz()->multiplied(aCurCyl->height()));
 137
 138       std::shared_ptr<GeomAPI_XYZ> anAxisXYZ = anAxis->xyz();
 139
 140       double aDist = anAxisXYZ->dot(aCurCylLoc->decreased(aLocation->xyz()));
 141       if (aDist < aMinHeight)
 142         aMinHeight = aDist;
 143       else if (aDist > aMaxHeight)
 144         aMaxHeight = aDist;
 145
 146       aDist = anAxisXYZ->dot(aCurCylLocHeight->decreased(aLocation->xyz()));
 147       if (aDist < aMinHeight)
 148         aMinHeight = aDist;
 149       else if (aDist > aMaxHeight)
 150         aMaxHeight = aDist;
 151
 152       if (aMinHeight < 0.0) {
 153         // move location of full cylinder
 154         aLocation->setX(aLocation->x() + aMinHeight * anAxis->x());
 155         aLocation->setY(aLocation->y() + aMinHeight * anAxis->y());
 156         aLocation->setZ(aLocation->z() + aMinHeight * anAxis->z());
 157       }
 158
 159       aHeight = aMaxHeight - aMinHeight;
 160     }
 161   }
 162
 163   GeomCylinderPtr aCylinder;
 164   if (isCylinder)
 165     aCylinder = GeomCylinderPtr(new GeomAPI_Cylinder(aLocation, anAxis, aRadius, aHeight));
 166   return aCylinder;
 167 }
 168
 169 //=================================================================================================
 170 std::shared_ptr<GeomAPI_Cone> GeomAPI_Shell::getCone() const
 171 {
 172   bool isCone = true;
 173   bool isFirstFace = true;
 174
 175   GeomPointPtr anApex;
 176   GeomDirPtr anAxis;
 177   double aSemiAngle, aCosSemiAngle;
 178   double aHeight1, aHeight2;
 179
 180   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_FACE); anExp.More(); anExp.Next()) {
 181     GeomFacePtr aFace(new GeomAPI_Face);
 182     aFace->setImpl(new TopoDS_Shape(anExp.Current()));
 183
 184     GeomConePtr aCurCone = aFace->getCone();
 185     if (!aCurCone) {
 186       isCone = false;
 187       break;
 188     }
 189
 190     if (isFirstFace) {
 191       anApex = aCurCone->apex();
 192       anAxis = aCurCone->axis();
 193       aSemiAngle = aCurCone->semiAngle();
 194       aCosSemiAngle = Cos(aSemiAngle);
 195       aHeight1 = aCurCone->radius1() * aCosSemiAngle;
 196       aHeight2 = aCurCone->radius2() * aCosSemiAngle;
 197       isFirstFace = false;
 198     }
 199     else {
 200       // check equal locations
 201       if (anApex->distance(aCurCone->apex()) >= Precision::Confusion() ||
 202       // check equal angles
 203           Abs(aSemiAngle - aCurCone->semiAngle()) >= Precision::Confusion() ||
 204       // check directions are collinear
 205           !anAxis->isParallel(aCurCone->axis())) {
 206         isCone = false;
 207         break;
 208       }
 209
 210       double aSign = anAxis->dot(aCurCone->axis());
 211       double aCurSemiAngle = aCurCone->semiAngle();
 212       double aCosCurSemiAngle = Cos(aSemiAngle);
 213
 214       double aH = aCurCone->radius1() * aCosCurSemiAngle * aSign;
 215       if (aH < aHeight1)
 216         aHeight1 = aH;
 217       else if (aH > aHeight2)
 218         aHeight2 = aH;
 219
 220       aH = aCurCone->radius2() * aCosCurSemiAngle * aSign;
 221       if (aH < aHeight1)
 222         aHeight1 = aH;
 223       else if (aH > aHeight2)
 224         aHeight2 = aH;
 225     }
 226   }
 227
 228   GeomConePtr aCone;
 229   if (isCone) {
 230     GeomPointPtr aLocation(new GeomAPI_Pnt(
 231       anApex->xyz()->added(anAxis->xyz()->multiplied(aHeight1))));
 232     double aRadius1 = aHeight1 * Tan(aSemiAngle);
 233     double aRadius2 = aHeight2 * Tan(aSemiAngle);
 234
 235     aCone = GeomConePtr(new GeomAPI_Cone(aLocation, anAxis, aSemiAngle, aRadius1, aRadius2));
 236   }
 237   return aCone;
 238 }
 239
 240 //=================================================================================================
 241 std::shared_ptr<GeomAPI_Torus> GeomAPI_Shell::getTorus() const
 242 {
 243   bool isTorus = true;
 244   bool isFirstFace = true;
 245   GeomTorusPtr aTorus;
 246
 247   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_FACE); anExp.More(); anExp.Next()) {
 248     GeomFacePtr aFace(new GeomAPI_Face);
 249     aFace->setImpl(new TopoDS_Shape(anExp.Current()));
 250
 251     GeomTorusPtr aCurTorus = aFace->getTorus();
 252     if (!aCurTorus) {
 253       isTorus = false;
 254       break;
 255     }
 256
 257     if (isFirstFace) {
 258       aTorus = aCurTorus;
 259       isFirstFace = false;
 260     }
 261     else {
 262       // compare radii
 263       if (Abs(aTorus->majorRadius() - aCurTorus->majorRadius()) >= Precision::Confusion() ||
 264           Abs(aTorus->minorRadius() - aCurTorus->minorRadius()) >= Precision::Confusion() ||
 265       // check equal centers
 266           aTorus->center()->distance(aCurTorus->center()) >= Precision::SquareConfusion() ||
 267       // check directions are collinear
 268           !aTorus->direction()->isParallel(aCurTorus->direction())) {
 269         isTorus = false;
 270         break;
 271       }
 272     }
 273   }
 274
 275   return isTorus ? aTorus : GeomTorusPtr();
 276 }
 277
 278 //=================================================================================================
 279 std::shared_ptr<GeomAPI_Box> GeomAPI_Shell::getParallelepiped() const
 280 {
 281   struct Plane
 282   {
 283     std::shared_ptr<GeomAPI_Ax3> myAxes;
 284     double myWidth;
 285     double myDepth;
 286     double myHeight;
 287   } aPlanes[6];
 288   std::map<int, int> aParallelPlanes;
 289
 290   int aNbPlanes = 0;
 291   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_WIRE);
 292        anExp.More() && aNbPlanes < 6; anExp.Next()) {
 293     GeomWirePtr aWire(new GeomAPI_Wire);
 294     aWire->setImpl(new TopoDS_Shape(anExp.Current()));
 295
 296     std::list<GeomPointPtr> aCorners;
 297     if (aWire->isRectangle(aCorners)) {
 298       // convert rectangle to plane with dimensions
 299       GeomPointPtr anOrigin = aCorners.front();
 300       aCorners.pop_front();
 301
 302       GeomPointPtr aFront = aCorners.front();
 303       GeomPointPtr aBack = aCorners.back();
 304
 305       aPlanes[aNbPlanes].myWidth = aBack->distance(anOrigin);
 306       aPlanes[aNbPlanes].myDepth = aFront->distance(anOrigin);
 307       aPlanes[aNbPlanes].myHeight = Precision::Infinite();
 308
 309       GeomDirPtr aDX(new GeomAPI_Dir(aBack->x() - anOrigin->x(),
 310                                      aBack->y() - anOrigin->y(),
 311                                      aBack->z() - anOrigin->z()));
 312       GeomDirPtr aDY(new GeomAPI_Dir(aFront->x() - anOrigin->x(),
 313                                      aFront->y() - anOrigin->y(),
 314                                      aFront->z() - anOrigin->z()));
 315       GeomDirPtr aDZ(new GeomAPI_Dir(aDX->cross(aDY)));
 316       aPlanes[aNbPlanes].myAxes =
 317           std::shared_ptr<GeomAPI_Ax3>(new GeomAPI_Ax3(anOrigin, aDX, aDZ));
 318
 319       // find parallel plane
 320       for (int i = 0; i < aNbPlanes; ++i) {
 321         double aDot = aPlanes[i].myAxes->normal()->dot(aDZ);
 322         if (Abs(aDot + 1.0) < Precision::Angular()) {
 323           if (aParallelPlanes.find(i) == aParallelPlanes.end())
 324             aParallelPlanes[i] = aNbPlanes;
 325           else
 326             break; // parallel planes already exist
 327         }
 328       }
 329
 330       ++aNbPlanes;
 331
 332     } else
 333       break;
 334   }
 335
 336   if (aNbPlanes != 6 || aParallelPlanes.size() != 3) // not a parallelepiped
 337     return GeomBoxPtr();
 338
 339   // calculate heights for planes computed by rectangles
 340   for (std::map<int, int>::iterator it = aParallelPlanes.begin();
 341        it != aParallelPlanes.end(); ++it) {
 342     GeomDirPtr aNormal = aPlanes[it->first].myAxes->normal();
 343     GeomPointPtr anOrigin = aPlanes[it->first].myAxes->origin();
 344     GeomPointPtr aNeighbor = aPlanes[it->second].myAxes->origin();
 345
 346     aPlanes[it->first].myHeight =
 347     aPlanes[it->second].myHeight =
 348         aNormal->xyz()->dot( aNeighbor->xyz()->decreased(anOrigin->xyz()) );
 349   }
 350
 351   // check if the box is oriented in the main axes
 352   int anIndex = 0;
 353   for (int i = 0; i < 6; ++i) {
 354     if (Abs(aPlanes[i].myAxes->dirX()->x() - 1.) < Precision::Angular() &&
 355         Abs(aPlanes[i].myAxes->normal()->z() - 1.) < Precision::Angular()) {
 356       anIndex = i;
 357       break;
 358     }
 359   }
 360
 361   // construct a box
 362   GeomBoxPtr aBox(new GeomAPI_Box(aPlanes[anIndex].myAxes,
 363                                   aPlanes[anIndex].myWidth,
 364                                   aPlanes[anIndex].myDepth,
 365                                   aPlanes[anIndex].myHeight));
 366   return aBox;
 367 }