src/IntPatch/IntPatch_PointLine.cxx

   1 // Created on: 2015-02-18
   2 // Created by: Nikolai BUKHALOV
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2015 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <IntPatch_PointLine.hxx>
  18
  19 #include <Adaptor3d_Surface.hxx>
  20 #include <IntSurf_PntOn2S.hxx>
  21 #include <Precision.hxx>
  22
  23 IMPLEMENT_STANDARD_RTTIEXT(IntPatch_PointLine,IntPatch_Line)
  24
  25 IntPatch_PointLine::IntPatch_PointLine (const Standard_Boolean Tang,
  26                                         const IntSurf_TypeTrans Trans1,
  27                                         const IntSurf_TypeTrans Trans2) :
  28   IntPatch_Line(Tang, Trans1, Trans2)
  29 {}
  30
  31 IntPatch_PointLine::IntPatch_PointLine (const Standard_Boolean Tang,
  32                                         const IntSurf_Situation Situ1,
  33                                         const IntSurf_Situation Situ2) :
  34   IntPatch_Line(Tang, Situ1, Situ2)
  35 {}
  36
  37 IntPatch_PointLine::IntPatch_PointLine (const Standard_Boolean Tang) :
  38   IntPatch_Line(Tang)
  39 {}
  40
  41 //=======================================================================
  42 //function : CurvatureRadiusOfIntersLine
  43 //purpose  :
  44 //        ATTENTION!!!
  45 //            Returns negative value if computation is not possible
  46 //=======================================================================
  47 Standard_Real IntPatch_PointLine::
  48                 CurvatureRadiusOfIntersLine(const Handle(Adaptor3d_Surface)& theS1,
  49                                             const Handle(Adaptor3d_Surface)& theS2,
  50                                             const IntSurf_PntOn2S& theUVPoint)
  51 {
  52   const Standard_Real aSmallValue = 1.0/Precision::Infinite();
  53   const Standard_Real aSqSmallValue = aSmallValue*aSmallValue;
  54
  55   Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
  56   theUVPoint.Parameters(aU1, aV1, aU2, aV2);
  57   gp_Pnt aPt;
  58   gp_Vec aDU1, aDV1, aDUU1, aDUV1, aDVV1;
  59   gp_Vec aDU2, aDV2, aDUU2, aDUV2, aDVV2;
  60
  61   theS1->D2(aU1, aV1, aPt, aDU1, aDV1, aDUU1, aDVV1, aDUV1);
  62   theS2->D2(aU2, aV2, aPt, aDU2, aDV2, aDUU2, aDVV2, aDUV2);
  63
  64 #if 0
  65   //The code in this block contains TEST CASES for
  66   //this algorithm only. It is stupedly to create OCCT-test for
  67   //the method, which will be changed possibly never.
  68   //However, if we do something in this method we can check its
  69   //functionality easily. For that:
  70   //  1. Initialize aTestID variable by the correct value;
  71   //  2. Compile this test code fragment.
  72
  73   int aTestID = 0;
  74   Standard_Real anExpectedSqRad = -1.0;
  75   switch(aTestID)
  76   {
  77   case 1:
  78     //Intersection between two spherical surfaces: O1(0.0, 0.0, 0.0), R1 = 3
  79     //and O2(5.0, 0.0, 0.0), R2 = 5.0.
  80     //Considered point has coordinates: (0.9, 0.0, 0.3*sqrt(91.0)).
  81
  82     aDU1.SetCoord(0.00000000000000000, 0.90000000000000002, 0.00000000000000000);
  83     aDV1.SetCoord(-2.8618176042508372, 0.00000000000000000, 0.90000000000000002);
  84     aDUU1.SetCoord(-0.90000000000000002, 0.00000000000000000, 0.00000000000000000);
  85     aDUV1.SetCoord(0.00000000000000000, -2.8618176042508372, 0.00000000000000000);
  86     aDVV1.SetCoord(-0.90000000000000002, 0.00000000000000000, -2.8618176042508372);
  87     aDU2.SetCoord(0.00000000000000000, -4.0999999999999996, 0.00000000000000000);
  88     aDV2.SetCoord(-2.8618176042508372, 0.00000000000000000, -4.0999999999999996);
  89     aDUU2.SetCoord(4.0999999999999996, 0.00000000000000000, 0.00000000000000000);
  90     aDUV2.SetCoord(0.00000000000000000, -2.8618176042508372, 0.00000000000000000);
  91     aDVV2.SetCoord(4.0999999999999996, 0.00000000000000000, -2.8618176042508372);
  92     anExpectedSqRad = 819.0/100.0;
  93     break;
  94   case 2:
  95     //Intersection between spherical surfaces: O1(0.0, 0.0, 0.0), R1 = 10
  96     //and the plane 3*x+4*y+z=26.
  97     //Considered point has coordinates: (-1.68, 5.76, 8.0).
  98
  99     aDU1.SetCoord(-5.76, -1.68, 0.0);
 100     aDV1.SetCoord(2.24, -7.68, 6.0);
 101     aDUU1.SetCoord(1.68, -5.76, 0.0);
 102     aDUV1.SetCoord(7.68, 2.24, 0.0);
 103     aDVV1.SetCoord(1.68, -5.76, -8.0);
 104     aDU2.SetCoord(1.0, 0.0, -3.0);
 105     aDV2.SetCoord(0.0, 1.0, -4.0);
 106     aDUU2.SetCoord(0.0, 0.0, 0.0);
 107     aDUV2.SetCoord(0.0, 0.0, 0.0);
 108     aDVV2.SetCoord(0.0, 0.0, 0.0);
 109     anExpectedSqRad = 74.0;
 110     break;
 111   default:
 112     aTestID = 0;
 113     break;
 114   }
 115 #endif
 116
 117   const gp_Vec aN1(aDU1.Crossed(aDV1)), aN2(aDU2.Crossed(aDV2));
 118   //Tangent vector to the intersection curve
 119   const gp_Vec aCTan(aN1.Crossed(aN2));
 120   const Standard_Real aSqMagnFDer = aCTan.SquareMagnitude();
 121
 122   if (aSqMagnFDer < 1.0e-8)
 123   {
 124     // Use 1.0e-4 (instead of aSmallValue) to provide
 125     // stable computation between different platforms.
 126     // See test bugs modalg_7 bug29807_sc01
 127     return -1.0;
 128   }
 129
 130   Standard_Real aDuS1 = 0.0, aDvS1 = 0.0, aDuS2 = 0.0, aDvS2 = 1.0;
 131
 132   {
 133     //This algorithm is described in NonSingularProcessing() function
 134     //in ApproxInt_ImpPrmSvSurfaces.gxx file
 135     Standard_Real aSqNMagn = aN1.SquareMagnitude();
 136     gp_Vec aTgU(aCTan.Crossed(aDU1)), aTgV(aCTan.Crossed(aDV1));
 137     Standard_Real aDeltaU = aTgV.SquareMagnitude()/aSqNMagn;
 138     Standard_Real aDeltaV = aTgU.SquareMagnitude()/aSqNMagn;
 139
 140     aDuS1 = Sign(sqrt(aDeltaU), aTgV.Dot(aN1));
 141     aDvS1 = -Sign(sqrt(aDeltaV), aTgU.Dot(aN1));
 142
 143     aSqNMagn = aN2.SquareMagnitude();
 144     aTgU.SetXYZ(aCTan.Crossed(aDU2).XYZ());
 145     aTgV.SetXYZ(aCTan.Crossed(aDV2).XYZ());
 146     aDeltaU = aTgV.SquareMagnitude()/aSqNMagn;
 147     aDeltaV = aTgU.SquareMagnitude()/aSqNMagn;
 148
 149     aDuS2 = Sign(sqrt(aDeltaU), aTgV.Dot(aN2));
 150     aDvS2 = -Sign(sqrt(aDeltaV), aTgU.Dot(aN2));
 151   }
 152
 153   //According to "Marching along surface/surface intersection curves
 154   //with an adaptive step length"
 155   //by Tz.E.Stoyagov
 156   //(http://www.sciencedirect.com/science/article/pii/016783969290046R)
 157   //we obtain the system:
 158   //            {A*a+B*b=F1
 159   //            {B*a+C*b=F2
 160   //where a and b should be found.
 161   //After that, 2nd derivative of the intersection curve can be computed as
 162   //            r''(t)=a*aN1+b*aN2.
 163
 164   const Standard_Real aA = aN1.Dot(aN1), aB = aN1.Dot(aN2), aC = aN2.Dot(aN2);
 165   const Standard_Real aDetSyst = aB*aB - aA*aC;
 166
 167   if(Abs(aDetSyst) < aSmallValue)
 168   {
 169     //Undetermined system solution
 170     return -1.0;
 171   }
 172
 173   const Standard_Real aF1 = aDuS1*aDuS1*aDUU1.Dot(aN1) +
 174                             2.0*aDuS1*aDvS1*aDUV1.Dot(aN1) +
 175                             aDvS1*aDvS1*aDVV1.Dot(aN1);
 176   const Standard_Real aF2 = aDuS2*aDuS2*aDUU2.Dot(aN2) +
 177                             2.0*aDuS2*aDvS2*aDUV2.Dot(aN2) +
 178                             aDvS2*aDvS2*aDVV2.Dot(aN2);
 179
 180   //Principal normal to the intersection curve
 181   const gp_Vec aCNorm((aF1*aC-aF2*aB)/aDetSyst*aN1 + (aA*aF2-aF1*aB)/aDetSyst*aN2);
 182   const Standard_Real aSqMagnSDer = aCNorm.CrossSquareMagnitude(aCTan);
 183
 184   if(aSqMagnSDer < aSqSmallValue)
 185   {//Intersection curve has null curvature in observed point
 186     return Precision::Infinite();
 187   }
 188
 189   //square of curvature radius
 190   const Standard_Real aFactSqRad = aSqMagnFDer*aSqMagnFDer*aSqMagnFDer/aSqMagnSDer;
 191
 192 #if 0
 193   if(aTestID)
 194   {
 195     if(Abs(aFactSqRad - anExpectedSqRad) < Precision::Confusion())
 196     {
 197       printf("OK: Curvature radius is equal to expected (%5.10g)", anExpectedSqRad);
 198     }
 199     else
 200     {
 201       printf("Error: Curvature radius is not equal to expected: %5.10g != %5.10g",
 202               aFactSqRad, anExpectedSqRad);
 203     }
 204   }
 205 #endif
 206
 207   return sqrt(aFactSqRad);
 208 }