1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkImplicitSelectionLoop.cxx
5 
6   Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7   All rights reserved.
8   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9 
10      This software is distributed WITHOUT ANY WARRANTY; without even
11      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12      PURPOSE.  See the above copyright notice for more information.
13 
14 =========================================================================*/
15 #include "vtkImplicitSelectionLoop.h"
16 
17 #include "vtkDoubleArray.h"
18 #include "vtkLine.h"
19 #include "vtkMath.h"
20 #include "vtkObjectFactory.h"
21 #include "vtkPlane.h"
22 #include "vtkPoints.h"
23 #include "vtkPolygon.h"
24 
25 vtkStandardNewMacro(vtkImplicitSelectionLoop);
26 vtkCxxSetObjectMacro(vtkImplicitSelectionLoop, Loop, vtkPoints);
27 
28 //------------------------------------------------------------------------------
29 // Instantiate object with no initial loop.
vtkImplicitSelectionLoop()30 vtkImplicitSelectionLoop::vtkImplicitSelectionLoop()
31 {
32   this->Loop = nullptr;
33   this->AutomaticNormalGeneration = 1;
34   this->Normal[0] = 0.0;
35   this->Normal[1] = 0.0;
36   this->Normal[2] = 1.0;
37   this->Polygon = vtkPolygon::New();
38 }
39 
40 //------------------------------------------------------------------------------
~vtkImplicitSelectionLoop()41 vtkImplicitSelectionLoop::~vtkImplicitSelectionLoop()
42 {
43   if (this->Loop)
44   {
45     this->Loop->Delete();
46   }
47   this->Polygon->Delete();
48   this->Polygon = nullptr;
49 }
50 
51 //------------------------------------------------------------------------------
52 #define VTK_DELTA 0.0001
53 // Generate plane equations only once to avoid a lot of extra work
Initialize()54 void vtkImplicitSelectionLoop::Initialize()
55 {
56   int i, numPts;
57   double x[3], xProj[3];
58 
59   numPts = this->Loop->GetNumberOfPoints();
60   this->Polygon->Points->SetDataTypeToDouble();
61   this->Polygon->Points->SetNumberOfPoints(numPts);
62 
63   if (this->AutomaticNormalGeneration)
64   {
65     // Make sure points define a loop with a normal
66     vtkPolygon::ComputeNormal(this->Loop, this->Normal);
67     if (this->Normal[0] == 0.0 && this->Normal[1] == 0.0 && this->Normal[2] == 0.0)
68     {
69       vtkErrorMacro(<< "Cannot determine inside/outside of loop");
70     }
71   }
72 
73   // Determine origin point by taking average
74   this->Origin[0] = this->Origin[1] = this->Origin[2] = 0.0;
75   for (i = 0; i < numPts; i++)
76   {
77     this->Loop->GetPoint(i, x);
78     this->Origin[0] += x[0];
79     this->Origin[1] += x[1];
80     this->Origin[2] += x[2];
81   }
82   this->Origin[0] /= numPts;
83   this->Origin[1] /= numPts;
84   this->Origin[2] /= numPts;
85 
86   // Project points onto plane generating new coordinates
87   for (i = 0; i < numPts; i++)
88   {
89     this->Loop->GetPoint(i, x);
90     vtkPlane::ProjectPoint(x, this->Origin, this->Normal, xProj);
91     this->Polygon->Points->SetPoint(i, xProj);
92   }
93 
94   this->Polygon->GetBounds(this->Bounds);
95   this->DeltaX = VTK_DELTA * (this->Bounds[1] - this->Bounds[0]);
96   this->DeltaY = VTK_DELTA * (this->Bounds[3] - this->Bounds[2]);
97   this->DeltaZ = VTK_DELTA * (this->Bounds[5] - this->Bounds[4]);
98   this->InitializationTime.Modified();
99 }
100 
101 //------------------------------------------------------------------------------
102 // Evaluate plane equations. Return smallest absolute value.
EvaluateFunction(double x[3])103 double vtkImplicitSelectionLoop::EvaluateFunction(double x[3])
104 {
105   int i, numPts;
106   double xProj[3];
107   double t, dist2, minDist2, closest[3];
108   int inside = 0;
109 
110   if (this->InitializationTime < this->GetMTime())
111   {
112     this->Initialize();
113   }
114   // Initialize may change the number of points
115   numPts = this->Polygon->Points->GetNumberOfPoints();
116 
117   // project point onto plane
118   vtkPlane::ProjectPoint(x, this->Origin, this->Normal, xProj);
119 
120   // determine whether it's in the selection loop and then evaluate point
121   // in polygon only if absolutely necessary.
122   if (xProj[0] >= this->Bounds[0] && xProj[0] <= this->Bounds[1] && xProj[1] >= this->Bounds[2] &&
123     xProj[1] <= this->Bounds[3] && xProj[2] >= this->Bounds[4] && xProj[2] <= this->Bounds[5] &&
124     vtkPolygon::PointInPolygon(xProj, numPts,
125       vtkArrayDownCast<vtkDoubleArray>(this->Polygon->Points->GetData())->GetPointer(0),
126       this->Bounds, this->Normal) == 1)
127   {
128     inside = 1;
129   }
130 
131   // determine distance to boundary
132   for (minDist2 = VTK_DOUBLE_MAX, i = 0; i < numPts; i++)
133   {
134     double p1[3], p2[3];
135     this->Polygon->Points->GetPoint(i, p1);
136     this->Polygon->Points->GetPoint((i + 1) % numPts, p2);
137     dist2 = vtkLine::DistanceToLine(xProj, p1, p2, t, closest);
138     if (dist2 < minDist2)
139     {
140       minDist2 = dist2;
141     }
142   }
143 
144   minDist2 = static_cast<double>(sqrt(minDist2));
145   return (inside ? -minDist2 : minDist2);
146 }
147 
148 //------------------------------------------------------------------------------
149 // Evaluate gradient of the implicit function. Use a numerical scheme: evaluate
150 // the function at four points (O,O+dx,O+dy,O+dz) and approximate the gradient.
151 // It's damn slow.
EvaluateGradient(double x[3],double n[3])152 void vtkImplicitSelectionLoop::EvaluateGradient(double x[3], double n[3])
153 {
154   double xp[3], yp[3], zp[3], g0, gx, gy, gz;
155   int i;
156 
157   g0 = this->EvaluateFunction(x); // side-effect is to compute DeltaX, Y, and Z
158 
159   for (i = 0; i < 3; i++)
160   {
161     xp[i] = yp[i] = zp[i] = x[i];
162   }
163   xp[0] += this->DeltaX;
164   yp[1] += this->DeltaY;
165   zp[2] += this->DeltaZ;
166 
167   gx = this->EvaluateFunction(xp);
168   gy = this->EvaluateFunction(yp);
169   gz = this->EvaluateFunction(zp);
170 
171   n[0] = (gx - g0) / this->DeltaX;
172   n[1] = (gy - g0) / this->DeltaY;
173   n[2] = (gz - g0) / this->DeltaZ;
174 }
175 
176 //------------------------------------------------------------------------------
GetMTime()177 vtkMTimeType vtkImplicitSelectionLoop::GetMTime()
178 {
179   vtkMTimeType mTime = this->vtkImplicitFunction::GetMTime();
180   vtkMTimeType time;
181 
182   if (this->Loop != nullptr)
183   {
184     time = this->Loop->GetMTime();
185     mTime = (time > mTime ? time : mTime);
186   }
187 
188   return mTime;
189 }
190 
191 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)192 void vtkImplicitSelectionLoop::PrintSelf(ostream& os, vtkIndent indent)
193 {
194   this->Superclass::PrintSelf(os, indent);
195 
196   if (this->Loop)
197   {
198     os << indent << "Loop of " << this->Loop->GetNumberOfPoints() << " points defined\n";
199   }
200   else
201   {
202     os << indent << "Loop not defined\n";
203   }
204 
205   os << indent
206      << "Automatic Normal Generation: " << (this->AutomaticNormalGeneration ? "On\n" : "Off\n");
207 
208   os << indent << "Normal: (" << this->Normal[0] << ", " << this->Normal[1] << ", "
209      << this->Normal[2] << ")\n";
210 }
211