src/INTERP_KERNEL/Geometric2D/EdgeLin.cxx

   1 //  Copyright (C) 2007-2008  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.
   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 email : webmaster.salome@opencascade.com
  18 //
  19 #include "EdgeLin.hxx"
  20 #include "Node.hxx"
  21 #include "InterpKernelException.hxx"
  22
  23 using namespace std;
  24 using namespace INTERP_KERNEL;
  25
  26 namespace INTERP_KERNEL
  27 {
  28   extern const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
  29 }
  30
  31 SegSegIntersector::SegSegIntersector(const EdgeLin& e1, const EdgeLin& e2):SameTypeEdgeIntersector(e1,e2)
  32 {
  33   _matrix[0]=(*(e2.getStartNode()))[0]-(*(e2.getEndNode()))[0];
  34   _matrix[1]=(*(e1.getEndNode()))[0]-(*(e1.getStartNode()))[0];
  35   _matrix[2]=(*(e2.getStartNode()))[1]-(*(e2.getEndNode()))[1];
  36   _matrix[3]=(*(e1.getEndNode()))[1]-(*(e1.getStartNode()))[1];
  37   _col[0]=_matrix[3]*(*(e1.getStartNode()))[0]-_matrix[1]*(*(e1.getStartNode()))[1];
  38   _col[1]=-_matrix[2]*(*(e2.getStartNode()))[0]+_matrix[0]*(*(e2.getStartNode()))[1];
  39   //Little trick to avoid problems if 'e1' and 'e2' are colinears and along Ox or Oy axes.
  40   if(fabs(_matrix[3])>fabs(_matrix[1]))
  41     _ind=0;
  42   else
  43     _ind=1;
  44 }
  45
  46 /*!
  47  * Must be called when 'this' and 'other' have been detected to be at least colinear. Typically they are overlapped.
  48  * Must be called after call of areOverlappedOrOnlyColinears.
  49  */
  50 bool SegSegIntersector::haveTheySameDirection() const
  51 {
  52   return (_matrix[3]*_matrix[1]+_matrix[2]*_matrix[0])>0.;
  53   //return (_matrix[_ind?1:0]>0. && _matrix[_ind?3:2]>0.) || (_matrix[_ind?1:0]<0. && _matrix[_ind?3:2]<0.);
  54 }
  55
  56 /*!
  57  * Precondition start and end must be so that there predecessor was in the same direction than 'e1'
  58  */
  59 void SegSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
  60 {
  61   getCurveAbscisse(start,whereStart,commonNode);
  62   getCurveAbscisse(end,whereEnd,commonNode);
  63 }
  64
  65 void SegSegIntersector::getCurveAbscisse(Node *node, TypeOfLocInEdge& where, MergePoints& commonNode) const
  66 {
  67   bool obvious;
  68   obviousCaseForCurvAbscisse(node,where,commonNode,obvious);
  69   if(obvious)
  70     return ;
  71   double ret=((*node)[!_ind]-(*_e1.getStartNode())[!_ind])/((*_e1.getEndNode())[!_ind]-(*_e1.getStartNode())[!_ind]);
  72   if(ret>0. && ret <1.)
  73     where=INSIDE;
  74   else if(ret<0.)
  75     where=OUT_BEFORE;
  76   else
  77     where=OUT_AFTER;
  78 }
  79
  80 /*!
  81  * areColinears method should be called before with a returned colinearity equal to false to avoid bad news.
  82  */
  83 std::list< IntersectElement > SegSegIntersector::getIntersectionsCharacteristicVal() const
  84 {
  85   std::list< IntersectElement > ret;
  86   double x=_matrix[0]*_col[0]+_matrix[1]*_col[1];
  87   double y=_matrix[2]*_col[0]+_matrix[3]*_col[1];
  88   //Only one intersect point possible
  89   Node *node=new Node(x,y);
  90   node->declareOn();
  91   bool i_1S=_e1.getStartNode()->isEqual(*node);
  92   bool i_1E=_e1.getEndNode()->isEqual(*node);
  93   bool i_2S=_e2.getStartNode()->isEqual(*node);
  94   bool i_2E=_e2.getEndNode()->isEqual(*node);
  95   ret.push_back(IntersectElement(_e1.getCharactValue(*node),
  96                                  _e2.getCharactValue(*node),
  97                                  i_1S,i_1E,i_2S,i_2E,node,_e1,_e2,keepOrder()));
  98   return ret;
  99 }
 100
 101 /*!
 102  * retrieves if segs are colinears.
 103  * WARNING !!! Contrary to areOverlappedOrOnlyColinears method, this method use an
 104  * another precision to detect colinearity !
 105  */
 106 bool SegSegIntersector::areColinears() const
 107 {
 108   double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
 109   return fabs(determinant)<QUADRATIC_PLANAR::_arc_detection_precision;
 110 }
 111
 112 /*!
 113  * Should be called \b once ! non const method.
 114  * \param whereToFind specifies the box where final seek should be done. Essentially it is used for caracteristic reason.
 115  * \param colinearity returns if regarding QUADRATIC_PLANAR::_precision ; e1 and e2 are colinears
 116  *                    If true 'this' is modified ! So this method be called once above all if true is returned for this parameter.
 117  * \param areOverlapped if colinearity if true, this parameter looks if e1 and e2 are overlapped.
 118  */
 119 void SegSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& colinearity, bool& areOverlapped)
 120 {
 121   double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
 122   if(fabs(determinant)>2.*QUADRATIC_PLANAR::_precision)//2*_precision due to max of offset on _start and _end
 123     {
 124       colinearity=false; areOverlapped=false;
 125       _matrix[0]/=determinant; _matrix[1]/=determinant; _matrix[2]/=determinant; _matrix[3]/=determinant;
 126     }
 127   else
 128     {
 129       colinearity=true;
 130       //retrieving initial matrix
 131       double tmp=_matrix[0]; _matrix[0]=_matrix[3]; _matrix[3]=tmp;
 132       _matrix[1]=-_matrix[1]; _matrix[2]=-_matrix[2];
 133       //
 134       double deno=sqrt(_matrix[0]*_matrix[0]+_matrix[1]*_matrix[1]);
 135       double x=(*(_e1.getStartNode()))[0]-(*(_e2.getStartNode()))[0];
 136       double y=(*(_e1.getStartNode()))[1]-(*(_e2.getStartNode()))[1];
 137       areOverlapped=fabs((_matrix[1]*y+_matrix[0]*x)/deno)<QUADRATIC_PLANAR::_precision;
 138     }
 139 }
 140
 141 EdgeLin::EdgeLin(std::istream& lineInXfig)
 142 {
 143   char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
 144   lineInXfig.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
 145   _start=new Node(lineInXfig);
 146   _end=new Node(lineInXfig);
 147   updateBounds();
 148 }
 149
 150 EdgeLin::EdgeLin(Node *start, Node *end, bool direction):Edge(start,end,direction)
 151 {
 152   updateBounds();
 153 }
 154
 155 EdgeLin::EdgeLin(double sX, double sY, double eX, double eY):Edge(sX,sY,eX,eY)
 156 {
 157   updateBounds();
 158 }
 159
 160 EdgeLin::~EdgeLin()
 161 {
 162 }
 163
 164 /*!
 165  * Characteristic for edges is relative position btw 0.;1.
 166  */
 167 bool EdgeLin::isIn(double characterVal) const
 168 {
 169   return characterVal>0. && characterVal<1.;
 170 }
 171
 172 Node *EdgeLin::buildRepresentantOfMySelf() const
 173 {
 174   return new Node(((*(_start))[0]+(*(_end))[0])/2.,((*(_start))[1]+(*(_end))[1])/2.);
 175 }
 176
 177 double EdgeLin::getCharactValue(const Node& node) const
 178 {
 179   return getCharactValueEng(node);
 180 }
 181
 182 double EdgeLin::getDistanceToPoint(const double *pt) const
 183 {
 184   double loc=getCharactValueEng(pt);
 185   if(loc>0. && loc<1.)
 186     {
 187       double tmp[2];
 188       tmp[0]=(*_start)[0]*(1-loc)+loc*(*_end)[0];
 189       tmp[1]=(*_start)[1]*(1-loc)+loc*(*_end)[1];
 190       return Node::distanceBtw2Pt(pt,tmp);
 191     }
 192   else
 193     {
 194       double dist1=Node::distanceBtw2Pt(*_start,pt);
 195       double dist2=Node::distanceBtw2Pt(*_end,pt);
 196       return std::min(dist1,dist2);
 197     }
 198 }
 199
 200 bool EdgeLin::isNodeLyingOn(const double *coordOfNode) const
 201 {
 202   double dBase=sqrt(_start->distanceWithSq(*_end));
 203   double d1=Node::distanceBtw2Pt(*_start,coordOfNode);
 204   d1+=Node::distanceBtw2Pt(*_end,coordOfNode);
 205   return Node::areDoubleEquals(dBase,d1);
 206 }
 207
 208 void EdgeLin::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
 209 {
 210   stream << "2 1 0 1 ";
 211   fillXfigStreamForLoc(stream);
 212   stream << " 7 50 -1 -1 0.000 0 0 -1 0 0 2" << endl;
 213   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
 214   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
 215   stream << endl;
 216 }
 217
 218 void EdgeLin::update(Node *m)
 219 {
 220   updateBounds();
 221 }
 222
 223 double EdgeLin::getNormSq() const
 224 {
 225   return _start->distanceWithSq(*_end);
 226 }
 227
 228 /*!
 229  * This methods computes :
 230  * \f[
 231  * \int_{Current Edge} -ydx
 232  * \f]
 233  */
 234 double EdgeLin::getAreaOfZone() const
 235 {
 236   return ((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
 237 }
 238
 239 void EdgeLin::getBarycenter(double *bary) const
 240 {
 241   bary[0]=((*_start)[0]+(*_end)[0])/2.;
 242   bary[1]=((*_start)[1]+(*_end)[1])/2.;
 243 }
 244
 245 /*!
 246  * \f[
 247  * bary[0]=\int_{Current Edge} -yxdx
 248  * \f]
 249  * \f[
 250  * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
 251  * \f]
 252  * To compute these 2 expressions in this class we have :
 253  * \f[
 254  * y=y_{1}+\frac{y_{2}-y_{1}}{x_{2}-x_{1}}(x-x_{1})
 255  * \f]
 256  */
 257 void EdgeLin::getBarycenterOfZone(double *bary) const
 258 {
 259   double x1=(*_start)[0];
 260   double y1=(*_start)[1];
 261   double x2=(*_end)[0];
 262   double y2=(*_end)[1];
 263   bary[0]=(x1-x2)*(y1*(2.*x1+x2)+y2*(2.*x2+x1))/6.;
 264   //bary[0]+=(y1-y2)*(x2*x2/3.-(x1*x2+x1*x1)/6.)+y1*(x1*x1-x2*x2)/2.;
 265   //bary[0]+=(y1-y2)*((x2*x2+x1*x2+x1*x1)/3.-(x2+x1)*x1/2.)+y1*(x1*x1-x2*x2)/2.;
 266   bary[1]=(x1-x2)*(y1*(y1+y2)+y2*y2)/6.;
 267 }
 268
 269 double EdgeLin::getCurveLength() const
 270 {
 271   double x=(*_start)[0]-(*_end)[0];
 272   double y=(*_start)[1]-(*_end)[1];
 273   return sqrt(x*x+y*y);
 274 }
 275
 276 Edge *EdgeLin::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
 277 {
 278   return new EdgeLin(start,end,direction);
 279 }
 280
 281 /*!
 282  * No precision should be introduced here. Just think as if precision was perfect.
 283  */
 284 void EdgeLin::updateBounds()
 285 {
 286   _bounds.setValues(std::min((*_start)[0],(*_end)[0]),std::max((*_start)[0],(*_end)[0]),std::min((*_start)[1],(*_end)[1]),std::max((*_start)[1],(*_end)[1]));
 287 }
 288
 289 double EdgeLin::getCharactValueEng(const double *node) const
 290 {
 291   double car1_1x=node[0]-(*(_start))[0]; double car1_2x=(*(_end))[0]-(*(_start))[0];
 292   double car1_1y=node[1]-(*(_start))[1]; double car1_2y=(*(_end))[1]-(*(_start))[1];
 293   return (car1_1x*car1_2x+car1_1y*car1_2y)/(car1_2x*car1_2x+car1_2y*car1_2y);
 294 }