1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkTextureMapToPlane.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 "vtkTextureMapToPlane.h"
16 
17 #include "vtkCellData.h"
18 #include "vtkDataSet.h"
19 #include "vtkFloatArray.h"
20 #include "vtkInformation.h"
21 #include "vtkInformationVector.h"
22 #include "vtkMath.h"
23 #include "vtkObjectFactory.h"
24 #include "vtkPointData.h"
25 
26 vtkStandardNewMacro(vtkTextureMapToPlane);
27 
28 // Construct with s,t range=(0,1) and automatic plane generation turned on.
vtkTextureMapToPlane()29 vtkTextureMapToPlane::vtkTextureMapToPlane()
30 {
31   // all zero - indicates that using normal is preferred and automatic is off
32   this->Origin[0] = this->Origin[1] = this->Origin[2] = 0.0;
33   this->Point1[0] = this->Point1[1] = this->Point1[2] = 0.0;
34   this->Point2[0] = this->Point2[1] = this->Point2[2] = 0.0;
35 
36   this->Normal[0] = 0.0;
37   this->Normal[1] = 0.0;
38   this->Normal[2] = 1.0;
39 
40   this->SRange[0] = 0.0;
41   this->SRange[1] = 1.0;
42 
43   this->TRange[0] = 0.0;
44   this->TRange[1] = 1.0;
45 
46   this->AutomaticPlaneGeneration = 1;
47 }
48 
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)49 int vtkTextureMapToPlane::RequestData(vtkInformation* vtkNotUsed(request),
50   vtkInformationVector** inputVector, vtkInformationVector* outputVector)
51 {
52   // get the info objects
53   vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
54   vtkInformation* outInfo = outputVector->GetInformationObject(0);
55 
56   // get the input and output
57   vtkDataSet* input = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
58   vtkDataSet* output = vtkDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
59 
60   double tcoords[2];
61   vtkIdType numPts;
62   vtkFloatArray* newTCoords;
63   vtkIdType i;
64   int j;
65   double proj, minProj, axis[3], sAxis[3], tAxis[3];
66   int dir = 0;
67   double s, t, sSf, tSf, p[3];
68   int abort = 0;
69   vtkIdType progressInterval;
70 
71   vtkDebugMacro(<< "Generating texture coordinates!");
72 
73   // First, copy the input to the output as a starting point
74   output->CopyStructure(input);
75 
76   if ((numPts = input->GetNumberOfPoints()) < 3 && this->AutomaticPlaneGeneration)
77   {
78     vtkErrorMacro(<< "Not enough points for automatic plane mapping\n");
79     return 1;
80   }
81 
82   //  Allocate texture data
83   //
84   newTCoords = vtkFloatArray::New();
85   newTCoords->SetName("Texture Coordinates");
86   newTCoords->SetNumberOfComponents(2);
87   newTCoords->SetNumberOfTuples(numPts);
88   progressInterval = numPts / 20 + 1;
89 
90   //  Compute least squares plane if on automatic mode; otherwise use
91   //  normal specified or plane specified
92   //
93   if (this->AutomaticPlaneGeneration &&
94     (this->Origin[0] == 0.0 && this->Origin[1] == 0.0 && this->Origin[2] == 0.0 &&
95       this->Point1[0] == 0.0 && this->Point1[1] == 0.0 && this->Point1[2] == 0.0))
96   {
97     this->ComputeNormal(output);
98 
99     vtkMath::Normalize(this->Normal);
100 
101     //  Now project each point onto plane generating s,t texture coordinates
102     //
103     //  Create local s-t coordinate system.  Need to find the two axes on
104     //  the plane and encompassing all the points.  Hence use the bounding
105     //  box as a reference.
106     //
107     for (minProj = 1.0, i = 0; i < 3; i++)
108     {
109       axis[0] = axis[1] = axis[2] = 0.0;
110       axis[i] = 1.0;
111       if ((proj = fabs(vtkMath::Dot(this->Normal, axis))) < minProj)
112       {
113         minProj = proj;
114         dir = i;
115       }
116     }
117     axis[0] = axis[1] = axis[2] = 0.0;
118     axis[dir] = 1.0;
119 
120     vtkMath::Cross(this->Normal, axis, tAxis);
121     vtkMath::Normalize(tAxis);
122 
123     vtkMath::Cross(tAxis, this->Normal, sAxis);
124 
125     //  Construct projection matrices
126     //
127     //  Arrange s-t axes so that parametric location of points will fall
128     //  between s_range and t_range.  Simplest to do by projecting maximum
129     //  corner of bounding box unto plane and backing out scale factors.
130     //
131     const double* bounds = output->GetBounds();
132     for (i = 0; i < 3; i++)
133     {
134       axis[i] = bounds[2 * i + 1] - bounds[2 * i];
135     }
136 
137     s = vtkMath::Dot(sAxis, axis);
138     t = vtkMath::Dot(tAxis, axis);
139 
140     sSf = (this->SRange[1] - this->SRange[0]) / s;
141     tSf = (this->TRange[1] - this->TRange[0]) / t;
142 
143     //  Now can loop over all points, computing parametric coordinates.
144     //
145     for (i = 0; i < numPts && !abort; i++)
146     {
147       if (!(i % progressInterval))
148       {
149         this->UpdateProgress((double)i / numPts);
150         abort = this->GetAbortExecute();
151       }
152 
153       output->GetPoint(i, p);
154       for (j = 0; j < 3; j++)
155       {
156         axis[j] = p[j] - bounds[2 * j];
157       }
158 
159       tcoords[0] = this->SRange[0] + vtkMath::Dot(sAxis, axis) * sSf;
160       tcoords[1] = this->TRange[0] + vtkMath::Dot(tAxis, axis) * tSf;
161 
162       newTCoords->SetTuple(i, tcoords);
163     }
164   } // compute plane and/or parametric range
165 
166   else // use the axes specified
167   {
168     double num, sDenom, tDenom;
169 
170     for (i = 0; i < 3; i++) // compute axes
171     {
172       sAxis[i] = this->Point1[i] - this->Origin[i];
173       tAxis[i] = this->Point2[i] - this->Origin[i];
174     }
175 
176     sDenom = vtkMath::Dot(sAxis, sAxis);
177     tDenom = vtkMath::Dot(tAxis, tAxis);
178 
179     if (sDenom == 0.0 || tDenom == 0.0)
180     {
181       vtkErrorMacro("Bad plane definition");
182       sDenom = tDenom = 1.0;
183     }
184 
185     // compute s-t coordinates
186     for (i = 0; i < numPts && !abort; i++)
187     {
188       if (!(i % progressInterval))
189       {
190         this->UpdateProgress((double)i / numPts);
191         abort = this->GetAbortExecute();
192       }
193       output->GetPoint(i, p);
194       for (j = 0; j < 3; j++)
195       {
196         axis[j] = p[j] - this->Origin[j];
197       }
198 
199       // s-coordinate
200       num = sAxis[0] * axis[0] + sAxis[1] * axis[1] + sAxis[2] * axis[2];
201       tcoords[0] = num / sDenom;
202 
203       // t-coordinate
204       num = tAxis[0] * axis[0] + tAxis[1] * axis[1] + tAxis[2] * axis[2];
205       tcoords[1] = num / tDenom;
206 
207       newTCoords->SetTuple(i, tcoords);
208     }
209   }
210 
211   // Update ourselves
212   //
213   output->GetPointData()->CopyTCoordsOff();
214   output->GetPointData()->PassData(input->GetPointData());
215   output->GetCellData()->PassData(input->GetCellData());
216 
217   output->GetPointData()->SetTCoords(newTCoords);
218   newTCoords->Delete();
219 
220   return 1;
221 }
222 
223 #define VTK_TOLERANCE 1.0e-03
224 
ComputeNormal(vtkDataSet * output)225 void vtkTextureMapToPlane::ComputeNormal(vtkDataSet* output)
226 {
227   vtkIdType numPts = output->GetNumberOfPoints();
228   double m[9], v[3], x[3];
229   vtkIdType ptId;
230   int dir = 0, i;
231   double length, w, *c1, *c2, *c3, det;
232 
233   //  First thing to do is to get an initial normal and point to define
234   //  the plane.  Then, use this information to construct better
235   //  matrices.  If problem occurs, then the point and plane becomes the
236   //  fallback value.
237   //
238   //  Get minimum width of bounding box.
239   const double* bounds = output->GetBounds();
240   length = output->GetLength();
241 
242   for (w = length, i = 0; i < 3; i++)
243   {
244     this->Normal[i] = 0.0;
245     if ((bounds[2 * i + 1] - bounds[2 * i]) < w)
246     {
247       dir = i;
248       w = bounds[2 * i + 1] - bounds[2 * i];
249     }
250   }
251 
252   //  If the bounds is perpendicular to one of the axes, then can
253   //  quickly compute normal.
254   //
255   this->Normal[dir] = 1.0;
256   if (w <= (length * VTK_TOLERANCE))
257   {
258     return;
259   }
260 
261   //  Need to compute least squares approximation.  Depending on major
262   //  normal direction (dir), construct matrices appropriately.
263   //
264   //  Compute 3x3 least squares matrix
265   v[0] = v[1] = v[2] = 0.0;
266   for (i = 0; i < 9; i++)
267   {
268     m[i] = 0.0;
269   }
270 
271   for (ptId = 0; ptId < numPts; ptId++)
272   {
273     output->GetPoint(ptId, x);
274 
275     v[0] += x[0] * x[2];
276     v[1] += x[1] * x[2];
277     v[2] += x[2];
278 
279     m[0] += x[0] * x[0];
280     m[1] += x[0] * x[1];
281     m[2] += x[0];
282 
283     m[3] += x[0] * x[1];
284     m[4] += x[1] * x[1];
285     m[5] += x[1];
286 
287     m[6] += x[0];
288     m[7] += x[1];
289   }
290   m[8] = numPts;
291 
292   //  Solve linear system using Kramers rule
293   //
294   c1 = m;
295   c2 = m + 3;
296   c3 = m + 6;
297   if ((det = vtkMath::Determinant3x3(c1, c2, c3)) <= VTK_TOLERANCE)
298   {
299     return;
300   }
301 
302   this->Normal[0] = vtkMath::Determinant3x3(v, c2, c3) / det;
303   this->Normal[1] = vtkMath::Determinant3x3(c1, v, c3) / det;
304   this->Normal[2] = -1.0; // because of the formulation
305 }
306 
PrintSelf(ostream & os,vtkIndent indent)307 void vtkTextureMapToPlane::PrintSelf(ostream& os, vtkIndent indent)
308 {
309   this->Superclass::PrintSelf(os, indent);
310 
311   os << indent << "Origin: (" << this->Origin[0] << ", " << this->Origin[1] << ", "
312      << this->Origin[2] << " )\n";
313 
314   os << indent << "Axis Point 1: (" << this->Point1[0] << ", " << this->Point1[1] << ", "
315      << this->Point1[2] << " )\n";
316 
317   os << indent << "Axis Point 2: (" << this->Point2[0] << ", " << this->Point2[1] << ", "
318      << this->Point2[2] << " )\n";
319 
320   os << indent << "S Range: (" << this->SRange[0] << ", " << this->SRange[1] << ")\n";
321   os << indent << "T Range: (" << this->TRange[0] << ", " << this->TRange[1] << ")\n";
322   os << indent
323      << "Automatic Normal Generation: " << (this->AutomaticPlaneGeneration ? "On\n" : "Off\n");
324   os << indent << "Normal: (" << this->Normal[0] << ", " << this->Normal[1] << ", "
325      << this->Normal[2] << ")\n";
326 }
327