1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkImageToPolyDataFilter.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 "vtkImageToPolyDataFilter.h"
16
17 #include "vtkAppendPolyData.h"
18 #include "vtkCellArray.h"
19 #include "vtkCellData.h"
20 #include "vtkEdgeTable.h"
21 #include "vtkImageData.h"
22 #include "vtkInformation.h"
23 #include "vtkInformationVector.h"
24 #include "vtkLine.h"
25 #include "vtkObjectFactory.h"
26 #include "vtkPointData.h"
27 #include "vtkPolyData.h"
28 #include "vtkScalarsToColors.h"
29 #include "vtkUnsignedCharArray.h"
30
31 vtkStandardNewMacro(vtkImageToPolyDataFilter);
32
33 vtkCxxSetObjectMacro(vtkImageToPolyDataFilter,LookupTable,vtkScalarsToColors);
34
vtkImageToPolyDataFilter()35 vtkImageToPolyDataFilter::vtkImageToPolyDataFilter()
36 {
37 this->OutputStyle = VTK_STYLE_POLYGONALIZE;
38 this->ColorMode = VTK_COLOR_MODE_LINEAR_256;
39 this->Smoothing = 1;
40 this->NumberOfSmoothingIterations = 40;
41 this->Decimation = 1;
42 this->DecimationError = 1.5;
43 this->Error = 100;
44 this->SubImageSize = 250;
45
46 this->Table = vtkUnsignedCharArray::New();
47 this->LookupTable = NULL;
48 }
49
~vtkImageToPolyDataFilter()50 vtkImageToPolyDataFilter::~vtkImageToPolyDataFilter()
51 {
52 this->Table->Delete();
53 if ( this->LookupTable )
54 {
55 this->LookupTable->Delete();
56 }
57 }
58
59 // declare helper functions
60 //
61
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)62 int vtkImageToPolyDataFilter::RequestData(
63 vtkInformation *vtkNotUsed(request),
64 vtkInformationVector **inputVector,
65 vtkInformationVector *outputVector)
66 {
67 // get the info objects
68 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
69 vtkInformation *outInfo = outputVector->GetInformationObject(0);
70
71 // get the input and output
72 vtkImageData *input = vtkImageData::SafeDownCast(
73 inInfo->Get(vtkDataObject::DATA_OBJECT()));
74 vtkPolyData *output = vtkPolyData::SafeDownCast(
75 outInfo->Get(vtkDataObject::DATA_OBJECT()));
76
77 vtkPolyData *tmpOutput;
78 vtkPolyData *tmpInput;
79 vtkAppendPolyData *append;
80 vtkPolyData *appendOutput;
81 vtkDataArray *inScalars = input->GetPointData()->GetScalars();
82 vtkIdType numPixels=input->GetNumberOfPoints();
83 int dims[3], numComp;
84 double origin[3], spacing[3];
85 vtkUnsignedCharArray *pixels;
86 int type;
87 int numPieces[2], extent[4];
88 int i, j, newDims[3], totalPieces, pieceNum, abortExecute=0;
89 double newOrigin[3];
90
91 // Check input and initialize
92 vtkDebugMacro(<<"Vectorizing image...");
93
94 if ( inScalars == NULL || numPixels < 1 )
95 {
96 vtkDebugMacro(<<"Not enough input to create output");
97 return 1;
98 }
99
100 append = vtkAppendPolyData::New();
101 tmpOutput=vtkPolyData::New();
102 tmpInput=vtkPolyData::New();
103 numComp=inScalars->GetNumberOfComponents();
104 type=inScalars->GetDataType();
105
106 appendOutput=append->GetOutput();
107
108 input->GetDimensions(dims);
109 input->GetOrigin(origin);
110 input->GetSpacing(spacing);
111
112 // Figure out how many pieces to break the image into (the image
113 // might be too big to process). The filter does a series of appends
114 // to join the pieces together.
115 numPieces[0] = ((dims[0]-2) / this->SubImageSize) + 1;
116 numPieces[1] = ((dims[1]-2) / this->SubImageSize) + 1;
117 totalPieces = numPieces[0]*numPieces[1];
118
119 appendOutput->Initialize(); //empty the output
120 append->AddInputData(tmpOutput); //output of piece
121 append->AddInputData(tmpInput); //output of previoius append
122
123 // Loop over this many pieces
124 for (pieceNum=j=0; j < numPieces[1] && !abortExecute; j++)
125 {
126 extent[2] = j*this->SubImageSize; //the y range
127 extent[3] = (j+1)*this->SubImageSize;
128 if ( extent[3] >= dims[1] )
129 {
130 extent[3] = dims[1] - 1;
131 }
132
133 for (i=0; i < numPieces[0] && !abortExecute; i++)
134 {
135 extent[0] = i*this->SubImageSize; //the x range
136 extent[1] = (i+1)*this->SubImageSize;
137 if ( extent[1] >= dims[0] )
138 {
139 extent[1] = dims[0] - 1;
140 }
141
142 vtkDebugMacro(<<"Processing #" << pieceNum);
143 this->UpdateProgress ((double)pieceNum/totalPieces);
144 if (this->GetAbortExecute())
145 {
146 abortExecute = 1;
147 break;
148 }
149 pieceNum++;
150
151 // Figure out characteristics of current sub-image
152 newDims[0] = extent[1] - extent[0] + 1;
153 newDims[1] = extent[3] - extent[2] + 1;
154 newOrigin[0] = origin[0] + extent[0]*spacing[0];
155 newOrigin[1] = origin[1] + extent[2]*spacing[1];
156 newOrigin[2] = 0.0;
157
158 // Create a quantized copy of the image based on the color table
159 //
160 pixels = this->QuantizeImage(inScalars, numComp, type, dims, extent);
161 vtkDebugMacro(<<"Quantizing color...image size (" <<newDims[0]
162 <<", " <<newDims[1] << ")");
163
164 // Generate polygons according to mode setting
165 //
166 if ( this->OutputStyle == VTK_STYLE_PIXELIZE )
167 {
168 this->PixelizeImage(pixels, newDims, newOrigin, spacing, tmpOutput);
169 }
170 else if ( this->OutputStyle == VTK_STYLE_RUN_LENGTH )
171 {
172 this->RunLengthImage(pixels, newDims, newOrigin, spacing, tmpOutput);
173 }
174 else //VTK_STYLE_POLYGONALIZE
175 {
176 this->PolygonalizeImage(pixels, newDims, newOrigin, spacing,tmpOutput);
177 }
178
179 // Append pieces together
180 //
181 tmpInput->CopyStructure(appendOutput);
182 tmpInput->GetPointData()->PassData(appendOutput->GetPointData());
183 tmpInput->GetCellData()->PassData(appendOutput->GetCellData());
184 append->Update();
185
186 // Clean up this iteration
187 //
188 pixels->Delete();
189 tmpInput->Initialize();
190 tmpOutput->Initialize();
191 } // for i pieces
192 } // for j pieces
193
194 // Create the final output contained in the append filter
195 output->CopyStructure(appendOutput);
196 output->GetPointData()->PassData(appendOutput->GetPointData());
197 output->GetCellData()->PassData(appendOutput->GetCellData());
198
199 append->Delete();
200 tmpInput->Delete();
201 tmpOutput->Delete();
202
203 return 1;
204 }
205
PixelizeImage(vtkUnsignedCharArray * pixels,int dims[3],double origin[3],double spacing[3],vtkPolyData * output)206 void vtkImageToPolyDataFilter::PixelizeImage(vtkUnsignedCharArray *pixels,
207 int dims[3], double origin[3],
208 double spacing[3],
209 vtkPolyData *output)
210 {
211 int numPts, numCells, i, j, id;
212 vtkIdType pts[4];
213 vtkPoints *newPts;
214 vtkCellArray *newPolys;
215 double x[3];
216 vtkUnsignedCharArray *polyColors;
217 unsigned char *ptr, *colors=pixels->GetPointer(0);
218
219 // create the points - see whether to create or append
220 numPts = (dims[0]+1) * (dims[1]+1);
221 newPts = vtkPoints::New();
222 newPts->SetNumberOfPoints(numPts);
223
224 x[2] = 0.0;
225 for (id=0, j=0; j<=dims[1]; j++)
226 {
227 x[1] = origin[1] + j*spacing[1];
228 for (i=0; i<=dims[0]; i++)
229 {
230 x[0] = origin[0] + i*spacing[0];
231 newPts->SetPoint(id, x);
232 id++;
233 }
234 }
235 output->SetPoints(newPts);
236 newPts->Delete();
237
238 // create the cells and cell colors
239 //
240 numCells = dims[0] * dims[1];
241 newPolys = vtkCellArray::New();
242 newPolys->Allocate(newPolys->EstimateSize(numCells,4));
243
244 polyColors = vtkUnsignedCharArray::New();
245 polyColors->SetNumberOfValues(3*numCells); //for rgb
246 polyColors->SetNumberOfComponents(3);
247
248 // loop over all pixels, creating a quad per pixel.
249 // Note: copying point data (pixel values) to cell data (quad colors).
250 for (id=0, j=0; j<dims[1]; j++)
251 {
252 for (i=0; i<dims[0]; i++)
253 {
254 pts[0] = i + j*(dims[0]+1);
255 pts[1] = pts[0] + 1;
256 pts[2] = pts[1] + dims[0] + 1;
257 pts[3] = pts[2] - 1;
258 newPolys->InsertNextCell(4, pts);
259 ptr = colors + 3*id;
260 polyColors->SetValue(3*id, ptr[0]);
261 polyColors->SetValue(3*id+1, ptr[1]);
262 polyColors->SetValue(3*id+2, ptr[2]);
263 id++;
264 }
265 }
266
267 output->SetPolys(newPolys);
268 newPolys->Delete();
269
270 output->GetCellData()->SetScalars(polyColors);
271 polyColors->Delete();
272 }
273
RunLengthImage(vtkUnsignedCharArray * pixels,int dims[3],double origin[3],double spacing[3],vtkPolyData * output)274 void vtkImageToPolyDataFilter::RunLengthImage(vtkUnsignedCharArray *pixels,
275 int dims[3], double origin[3],
276 double spacing[3],
277 vtkPolyData *output)
278 {
279 int i, j;
280 vtkIdType pts[4], id;
281 vtkPoints *newPts;
282 vtkCellArray *newPolys;
283 double x[3], minX, maxX, minY, maxY;
284 vtkUnsignedCharArray *polyColors;
285 unsigned char *colors=pixels->GetPointer(0), *color;
286
287 // Setup data
288 newPts = vtkPoints::New();
289 newPolys = vtkCellArray::New();
290 newPolys->Allocate(newPolys->EstimateSize(dims[0]*dims[1]/10,4));
291
292 polyColors = vtkUnsignedCharArray::New();
293 polyColors->Allocate(3*dims[0]*dims[1]/10); //for rgb
294 polyColors->SetNumberOfComponents(3);
295
296 // Loop over row-by-row generating quad polygons
297 x[2] = 0.0;
298 for (j=0; j<dims[1]; j++)
299 {
300 if ( j == 0 )
301 {
302 minY = origin[1];
303 maxY = origin[1] + 0.5*spacing[1];
304 }
305 else if ( j == (dims[1]-1) )
306 {
307 minY = origin[1] + j*spacing[1] - 0.5*spacing[1];
308 maxY = origin[1] + j*spacing[1];
309 }
310 else
311 {
312 minY = origin[1] + j*spacing[1] - 0.5*spacing[1];
313 maxY = origin[1] + j*spacing[1] + 0.5*spacing[1];
314 }
315
316 for ( i=0; i < dims[0]; )
317 {
318 if ( i == 0 )
319 {
320 minX = origin[0];
321 }
322 else
323 {
324 minX = origin[0] + i*spacing[0] - 0.5*spacing[0];
325 }
326 color = colors + 3*(i+j*dims[0]);
327 while ( i < dims[0] )
328 {
329 unsigned char *ptr = colors + 3*(i+j*dims[0]);
330 if ( ! this->IsSameColor(color,ptr) )
331 {
332 break;
333 }
334 else
335 {
336 i++;
337 }
338 }
339
340 if ( i >= dims[0] )
341 {
342 maxX = origin[0] + (dims[0]-1)*spacing[0];
343 }
344 else
345 {
346 maxX = origin[0] + (i-1)*spacing[0] + 0.5*spacing[0];
347 }
348
349 // Create quad cell
350 x[0] = minX;
351 x[1] = minY;
352 pts[0] = newPts->InsertNextPoint(x);
353 x[0] = maxX;
354 pts[1] = newPts->InsertNextPoint(x);
355 x[1] = maxY;
356 pts[2] = newPts->InsertNextPoint(x);
357 x[0] = minX;
358 pts[3] = newPts->InsertNextPoint(x);
359 id = newPolys->InsertNextCell(4,pts);
360 polyColors->InsertValue(3*id, color[0]);
361 polyColors->InsertValue(3*id+1, color[1]);
362 polyColors->InsertValue(3*id+2, color[2]);
363 }
364 }
365
366 output->SetPoints(newPts);
367 newPts->Delete();
368 output->SetPolys(newPolys);
369 newPolys->Delete();
370
371 output->GetCellData()->SetScalars(polyColors);
372 polyColors->Delete();
373 }
374
375
PolygonalizeImage(vtkUnsignedCharArray * pixels,int dims[3],double origin[3],double spacing[3],vtkPolyData * output)376 void vtkImageToPolyDataFilter::PolygonalizeImage(vtkUnsignedCharArray *pixels,
377 int dims[3], double origin[3], double spacing[3],
378 vtkPolyData *output)
379 {
380 int numPolys;
381 int numPixels=dims[0]*dims[1];
382
383 // Perform connected traversal on quantized points. This builds
384 // the initial "polygons" in implicit form.
385 //
386 this->PolyColors = vtkUnsignedCharArray::New();
387 this->PolyColors->SetNumberOfComponents(3);
388 this->PolyColors->Allocate(5000);
389
390 numPolys = this->ProcessImage(pixels, dims);
391 vtkDebugMacro(<<"Visited regions..." << numPolys << " polygons");
392
393 // Build edges around the boundary of the polygons. Also identify
394 // junction points where 3 or 4 polygons meet.
395 //
396 vtkPoints *points = vtkPoints::New();
397 points->Allocate(numPixels/2, numPixels/2);
398
399 vtkUnsignedCharArray *pointDescr = vtkUnsignedCharArray::New();
400 pointDescr->Allocate(numPixels/2, numPixels/2);
401
402 vtkCellArray *edgeConn = vtkCellArray::New();
403 edgeConn->Allocate(numPixels/2, numPixels/2);
404 vtkPolyData *edges = vtkPolyData::New();
405 edges->SetPoints(points);
406 edges->SetLines(edgeConn);
407 points->Delete();
408 edgeConn->Delete();
409
410 this->BuildEdges(pixels, dims, origin, spacing, pointDescr, edges);
411 vtkDebugMacro(<<"Edges built...");
412
413 // Now that we've got the edges, we have to build the "loops" around the
414 // polygons that define the polygon explicitly.
415 //
416 vtkUnsignedCharArray *polyColors = vtkUnsignedCharArray::New();
417 polyColors->SetNumberOfComponents(3);
418 polyColors->SetNumberOfValues(numPolys*3);
419
420 this->BuildPolygons(pointDescr, edges, numPolys, polyColors);
421 this->PolyColors->Delete();
422 delete [] this->Visited;
423 vtkDebugMacro(<<"Constructed polygons...");
424
425 // Smooth edge network. Some points are identified as fixed, others
426 // move using Laplacian smoothing.
427 //
428 if ( this->Smoothing )
429 {
430 this->SmoothEdges(pointDescr, edges);
431 vtkDebugMacro(<<"Edges smoothed...");
432 }
433
434 // Decimate edge network. There will be colinear vertices along edges.
435 // These are eliminated.
436 //
437 if ( this->Decimation )
438 {
439 this->DecimateEdges(edges, pointDescr, this->DecimationError);
440 }
441
442 // Create output polydata. Each polyon is output with its edges.
443 //
444 this->GeneratePolygons(edges, numPolys, output, polyColors, pointDescr);
445 vtkDebugMacro(<<"Output generated...");
446
447 // clean up and get out
448 edges->Delete();
449 polyColors->Delete();
450 pointDescr->Delete();
451 }
452
453 // The following are private helper functions----------------------------------
454 //
QuantizeImage(vtkDataArray * inScalars,int numComp,int type,int dims[3],int extent[4])455 vtkUnsignedCharArray *vtkImageToPolyDataFilter::QuantizeImage(
456 vtkDataArray *inScalars, int numComp, int type,
457 int dims[3], int extent[4])
458 {
459 int numPixels, i, j, idx, id;
460 vtkUnsignedCharArray *pixels;
461 unsigned char *ptr, *ptr2, *outPixels;
462 unsigned char *inPixels;
463
464 // doing a portion of the image
465 numPixels = (extent[1]-extent[0]+1) * (extent[3]-extent[2]+1);
466 pixels = vtkUnsignedCharArray::New();
467 pixels->SetNumberOfValues(3*numPixels);
468 outPixels = pixels->GetPointer(0);
469
470 // Figure out how to quantize
471 //
472 if ( this->ColorMode == VTK_COLOR_MODE_LINEAR_256 )
473 {
474 // Check scalar type
475 if ( type != VTK_UNSIGNED_CHAR || numComp != 3 )
476 {
477 vtkErrorMacro(<<"Wrong input scalar type");
478 return 0;
479 }
480 else
481 {
482 inPixels = static_cast<vtkUnsignedCharArray *>(inScalars)->GetPointer(0);
483 }
484
485 // Generate a color table used to quantize the points
486 //
487 if ( this->GetMTime() > this->TableMTime )
488 {
489 this->BuildTable(inPixels);
490 }
491
492 for (id=0, j=extent[2]; j <= extent[3]; j++)
493 {
494 for (i=extent[0]; i <= extent[1]; i++)
495 {
496 idx = i + j*dims[0];
497 ptr = inPixels + 3*idx;
498 ptr2 = outPixels + 3*id;
499 const unsigned char *color = this->GetColor(ptr);
500 ptr2[0] = color[0];
501 ptr2[1] = color[1];
502 ptr2[2] = color[2];
503 id++;
504 }
505 }
506 }//using build in table
507
508 else //using provided lookup table
509 {
510 if ( numComp != 1 || this->LookupTable == NULL )
511 {
512 vtkErrorMacro(<<"LUT mode requires single component scalar and LUT");
513 return 0;
514 }
515
516 double s;
517 for (id=0, j=extent[2]; j <= extent[3]; j++)
518 {
519 for (i=extent[0]; i <= extent[1]; i++)
520 {
521 idx = i + j*dims[0];
522 s = inScalars->GetComponent(idx,0);
523 const unsigned char *color = this->LookupTable->MapValue(s);
524 ptr2 = outPixels + 3*id;
525 ptr2[0] = color[0];
526 ptr2[1] = color[1];
527 ptr2[2] = color[2];
528 id++;
529 }
530 }
531 }
532
533 return pixels;
534 }
535
BuildTable(unsigned char * vtkNotUsed (inPixels))536 void vtkImageToPolyDataFilter::BuildTable(unsigned char *vtkNotUsed(inPixels))
537 {
538 int red, green, blue, idx=0;
539
540 this->Table->SetNumberOfValues(256*3);
541
542 // use 3-3-2 bits for rgb
543 for (blue=0; blue<256; blue+=64)
544 {
545 for (green=0; green<256; green+=32)
546 {
547 for (red=0; red<256; red+=32)
548 {
549 this->Table->SetValue(idx++, red);
550 this->Table->SetValue(idx++, green);
551 this->Table->SetValue(idx++, blue);
552 }
553 }
554 }
555 }
556
FillInputPortInformation(int,vtkInformation * info)557 int vtkImageToPolyDataFilter::FillInputPortInformation(int, vtkInformation *info)
558 {
559 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
560 return 1;
561 }
562
PrintSelf(ostream & os,vtkIndent indent)563 void vtkImageToPolyDataFilter::PrintSelf(ostream& os, vtkIndent indent)
564 {
565 this->Superclass::PrintSelf(os,indent);
566
567 os << indent << "Output Style: ";
568 if ( this->OutputStyle == VTK_STYLE_PIXELIZE )
569 {
570 os << indent << "Pixelize\n";
571 }
572 else if ( this->OutputStyle == VTK_STYLE_RUN_LENGTH )
573 {
574 os << indent << "RunLength\n";
575 }
576 else // this->OutputStyle == VTK_STYLE_POLYGONALIZE
577 {
578 os << indent << "Polygonalize\n";
579 }
580
581 os << indent << "Color Mode: ";
582 if ( this->ColorMode == VTK_STYLE_PIXELIZE )
583 {
584 os << indent << "LUT\n";
585 }
586 else // this->ColorMode == VTK_STYLE_POLYGONALIZE
587 {
588 os << indent << "Linear256\n";
589 }
590
591 os << indent << "Smoothing: " << (this->Smoothing ? "On\n" : "Off\n");
592 os << indent << "Number of Smoothing Iterations: "
593 << this->NumberOfSmoothingIterations << "\n";
594
595 os << indent << "Decimation: " << (this->Decimation ? "On\n" : "Off\n");
596 os << indent << "Decimation Error: "
597 << (this->DecimationError ? "On\n" : "Off\n");
598
599 os << indent << "Error: " << this->Error << "\n";
600 os << indent << "Sub-Image Size: " << this->SubImageSize << "\n";
601
602 if ( this->LookupTable )
603 {
604 os << indent << "LookupTable:\n";
605 this->LookupTable->PrintSelf(os,indent.GetNextIndent());
606 }
607 else
608 {
609 os << indent << "LookupTable: (none)\n";
610 }
611
612 }
613
614
615 //--------------------------private helper functions---------------------------
616 // Determines whether two pixels are the same color
IsSameColor(unsigned char * p1,unsigned char * p2)617 int vtkImageToPolyDataFilter::IsSameColor(unsigned char *p1, unsigned char *p2)
618 {
619 int d2 = (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) +
620 (p1[2]-p2[2])*(p1[2]-p2[2]);
621
622 return (d2 > this->Error ? 0 : 1);
623 }
624
GetColor(unsigned char * rgb)625 unsigned char *vtkImageToPolyDataFilter::GetColor(unsigned char *rgb)
626 {
627 // round to nearest value
628 int red = (rgb[0] + 16) / 32;
629 red = (red > 7 ? 7 : red);
630 int green =(rgb[1] + 16) / 32;
631 green = (green > 7 ? 7 : green);
632 int blue = (rgb[2] + 32) / 64;
633 blue = (blue > 3 ? 3 : blue);
634
635 return this->Table->GetPointer(3*(red + green*8 + blue*64));
636 }
637
GetIJ(int id,int & i,int & j,int dims[3])638 void vtkImageToPolyDataFilter::GetIJ(int id, int &i, int &j, int dims[3])
639 {
640 i = id % dims[0];
641 j = id / dims[0];
642 }
643
644 // Get the left-right-top-bottom neighboring pixels of a given pixel
645 // The method has been modified to return right neighbor (mode==0);
646 // or top neighbor (mode==1) or all neighbors (mode==2).
GetNeighbors(unsigned char * ptr,int & i,int & j,int dims[2],unsigned char * neighbors[4],int mode)647 int vtkImageToPolyDataFilter::GetNeighbors(unsigned char *ptr, int &i, int &j,
648 int dims[2],unsigned char *neighbors[4], int mode)
649 {
650 int numNeis=0;
651
652 if ( mode == 0 )
653 {
654 if ( (i+1) < dims[0] )
655 {
656 neighbors[numNeis++] = ptr + 3; //jump over rgb
657 }
658 if ( (i-1) >= 0 )
659 {
660 neighbors[numNeis++] = ptr - 3; //jump over rgb
661 }
662 }
663
664 else if ( mode == 1 )
665 {
666 if ( (j+1) < dims[1] )
667 {
668 neighbors[numNeis++] = ptr + 3*dims[0];
669 }
670 }
671
672 else
673 {
674 if ( (i+1) < dims[0] )
675 {
676 neighbors[numNeis++] = ptr + 3; //jump over rgb
677 }
678 if ( (i-1) >= 0 )
679 {
680 neighbors[numNeis++] = ptr - 3;
681 }
682
683 if ( (j+1) < dims[1] )
684 {
685 neighbors[numNeis++] = ptr + 3*dims[0];
686 }
687 if ( (j-1) >= 0 )
688 {
689 neighbors[numNeis++] = ptr - 3*dims[0];
690 }
691 }
692
693 return numNeis;
694
695 }
696
697
698 // Marks connected regions with different colors.
ProcessImage(vtkUnsignedCharArray * scalars,int dims[2])699 int vtkImageToPolyDataFilter::ProcessImage(vtkUnsignedCharArray *scalars,
700 int dims[2])
701 {
702 int numPixels = dims[0]*dims[1];
703 vtkIdList *wave, *wave2, *tmpWave;
704 int numIds, regionNumber, i, j, k, id, x, y, numNeighbors;
705 unsigned char *neighbors[4], *ptr;
706 unsigned char *pixels = scalars->GetPointer(0);
707
708 // Collect groups of pixels together into similar colored regions. These
709 // will be eventually grouped into polygons and/or lines.
710 //
711 // mark all pixels unvisited
712 regionNumber = -1;
713 this->Visited = new int [numPixels];
714 memset (this->Visited, (int)-1, numPixels*sizeof(int));
715
716 // set up the connected traversal
717 wave = vtkIdList::New();
718 wave->Allocate(static_cast<int>(numPixels/4.0),
719 static_cast<int>(numPixels/4.0));
720 wave2 = vtkIdList::New();
721 wave2->Allocate(static_cast<int>(numPixels/4.0),
722 static_cast<int>(numPixels/4.0));
723
724 // visit connected pixels. Pixels are connected if they are topologically
725 // adjacent and they have "equal" color values.
726 for (i=0; i < numPixels; i++)
727 {
728 if ( this->Visited[i] == -1 )
729 {//start a connected wave
730 this->Visited[i] = ++regionNumber;
731 ptr = pixels + 3*i;
732 this->PolyColors->InsertValue(3*regionNumber, ptr[0]); //assign color
733 this->PolyColors->InsertValue(3*regionNumber+1, ptr[1]);
734 this->PolyColors->InsertValue(3*regionNumber+2, ptr[2]);
735 wave->Reset(); wave2->Reset();
736
737 // To prevent creating polygons with inner loops, we're going to start
738 // the wave as a "vertical" stack of pixels, and then propagate the
739 // wave horizontally only.
740 wave->InsertId(0,i);
741 this->GetIJ(i, x, y, dims);
742 while ( (numNeighbors = this->GetNeighbors(ptr, x, y, dims, neighbors, 1)) )
743 {
744 id = (neighbors[0] - pixels) / 3;
745 if ( this->Visited[id] == -1 && this->IsSameColor(ptr, neighbors[0]) )
746 {
747 this->Visited[id] = regionNumber;
748 wave->InsertNextId(id);
749 ptr = pixels + 3*id;
750 this->GetIJ(id, x, y, dims);
751 }
752 else
753 {
754 break;
755 }
756 }
757
758 // Okay, defined vertical wave, now propagate horizontally
759 numIds = wave->GetNumberOfIds();
760 while ( numIds > 0 )
761 {
762 for (j=0; j<numIds; j++) //propagate wave
763 {
764 id = wave->GetId(j);
765 ptr = pixels + 3*id;
766 this->GetIJ(id, x, y, dims);
767 numNeighbors = this->GetNeighbors(ptr, x, y, dims, neighbors, 0);
768 for (k=0; k<numNeighbors; k++)
769 {
770 id = (neighbors[k] - pixels) / 3;
771 if ( this->Visited[id] == -1 && this->IsSameColor(ptr, neighbors[k]) )
772 {
773 this->Visited[id] = regionNumber;
774 wave2->InsertNextId(id);
775 }
776 }//for each pixel neighbor
777 }//for pixels left in wave
778 numIds = wave2->GetNumberOfIds();
779 tmpWave = wave;
780 wave = wave2;
781 wave2 = tmpWave;
782 wave2->Reset();
783 }//while still propogating
784 }//if not, start wave
785 }//for all pixels
786
787 wave->Delete();
788 wave2->Delete();
789
790 return regionNumber+1;
791 }
792
793 // Create polygons and place into output
GeneratePolygons(vtkPolyData * edges,int vtkNotUsed (numPolys),vtkPolyData * output,vtkUnsignedCharArray * polyColors,vtkUnsignedCharArray * pointDescr)794 void vtkImageToPolyDataFilter::GeneratePolygons(vtkPolyData *edges,
795 int vtkNotUsed(numPolys), vtkPolyData *output,
796 vtkUnsignedCharArray *polyColors,
797 vtkUnsignedCharArray *pointDescr)
798 {
799 vtkCellArray *newPolys, *inPolys;
800 int i, numPts;
801 vtkIdType *pts = 0;
802 vtkIdType npts = 0;
803
804 // Copy the points via reference counting
805 //
806 output->SetPoints(edges->GetPoints());
807
808 // Create the polygons - points may have been decimated so these
809 // points have to be culled.
810 //
811 inPolys = edges->GetPolys();
812 newPolys = vtkCellArray::New();
813 newPolys->Allocate(inPolys->GetSize());
814
815 for ( inPolys->InitTraversal(); inPolys->GetNextCell(npts,pts); )
816 {
817 newPolys->InsertNextCell(0);
818 numPts = 0;
819 for (i=0; i<npts; i++)
820 {
821 if ( pointDescr->GetValue(pts[i]) != 2 )
822 {
823 newPolys->InsertCellPoint(pts[i]);
824 numPts++;
825 }
826 }
827 newPolys->UpdateCellCount(numPts);
828 }
829
830 output->SetPolys(newPolys);
831 newPolys->Delete();
832
833 output->GetCellData()->SetScalars(polyColors);
834 }
835
836 // Uses clipping approach to build the polygon edges
BuildEdges(vtkUnsignedCharArray * vtkNotUsed (pixels),int dims[3],double origin[3],double spacing[3],vtkUnsignedCharArray * pointDescr,vtkPolyData * edges)837 int vtkImageToPolyDataFilter::BuildEdges(vtkUnsignedCharArray *vtkNotUsed(pixels),
838 int dims[3], double origin[3],
839 double spacing[3],
840 vtkUnsignedCharArray *pointDescr,
841 vtkPolyData *edges)
842 {
843 double x[3];
844 int i, j, edgeCount;
845 vtkIdType ptId, p0, p1, p2, p3, startId, attrId, id[8], pts[4];
846 vtkCellArray *edgeConn = edges->GetLines();
847 vtkPoints *points = edges->GetPoints();
848
849 // Build edges around perimeter of image. Note that the point ids
850 // The first four points are the image corners and are inserted and
851 // marked so that they can't be moved during smoothing.
852 points->InsertPoint(0, origin);
853 pointDescr->InsertValue(0, 1);
854
855 // Keep track of the polygons that use each edge as well as associated
856 // intersection points on edge (if any)
857 this->EdgeTable = vtkEdgeTable::New();
858 this->EdgeTable->InitEdgeInsertion(dims[0]*dims[1],1);
859
860 this->EdgeUseTable = vtkEdgeTable::New();
861 this->EdgeUseTable->InitEdgeInsertion(dims[0]*dims[1],1);
862
863 this->EdgeUses = vtkIntArray::New();
864 this->EdgeUses->SetNumberOfComponents(2);
865 this->EdgeUses->Allocate(4*dims[0]*dims[1],dims[0]*dims[1]);
866
867 // Generate corner points of image
868 x[0] = origin[0] + (dims[0]-1)*spacing[0];
869 x[1] = origin[1];
870 x[2] = 0.0;
871 points->InsertPoint(1, x);
872 pointDescr->InsertValue(1, 1);
873
874 x[0] = origin[0] + (dims[0]-1)*spacing[0];
875 x[1] = origin[1] + (dims[1]-1)*spacing[1];;
876 x[2] = 0.0;
877 points->InsertPoint(2, x);
878 pointDescr->InsertValue(2, 1);
879
880 x[0] = origin[0];
881 x[1] = origin[1] + (dims[1]-1)*spacing[1];
882 x[2] = 0.0;
883 points->InsertPoint(3, x);
884 pointDescr->InsertValue(3, 1);
885
886 // Let's create perimeter edges - bottom x edge
887 startId = 0;
888 x[1] = origin[1];
889 for (i=0; i<(dims[0]-1); i++)
890 {
891 p0 = i;
892 p1 = i + 1;
893 if ( this->Visited[p0] != this->Visited[p1] )
894 {
895 x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
896 ptId = points->InsertNextPoint(x);
897 this->EdgeTable->InsertEdge(p0,p1,ptId);
898 pointDescr->InsertValue(ptId, 1); //can't be smoothed
899
900 edgeConn->InsertNextCell(2);
901 edgeConn->InsertCellPoint(startId);
902 edgeConn->InsertCellPoint(ptId);
903 attrId = this->EdgeUseTable->InsertEdge(startId,ptId);
904 this->EdgeUses->InsertValue(2*attrId, this->Visited[p0]);
905 this->EdgeUses->InsertValue(2*attrId+1, -1);
906 startId = ptId;
907 }
908 }
909 edgeConn->InsertNextCell(2); //finish off the edge
910 edgeConn->InsertCellPoint(startId);
911 edgeConn->InsertCellPoint(1);
912 attrId = this->EdgeUseTable->InsertEdge(startId,1);
913 this->EdgeUses->InsertValue(2*attrId, this->Visited[dims[0]-1]);
914 this->EdgeUses->InsertValue(2*attrId+1, -1);
915
916 // Let's create perimeter edges - top x edge
917 startId = 3;
918 x[1] = origin[1] + (dims[1]-1)*spacing[1];
919 for (i=0; i<(dims[0]-1); i++)
920 {
921 p0 = i + dims[0]*(dims[1]-1);
922 p1 = p0 + 1;
923 if ( this->Visited[p0] != this->Visited[p1] )
924 {
925 x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
926 ptId = points->InsertNextPoint(x);
927 this->EdgeTable->InsertEdge(p0,p1,ptId);
928 pointDescr->InsertValue(ptId, 1); //can't be smoothed
929
930 edgeConn->InsertNextCell(2);
931 edgeConn->InsertCellPoint(startId);
932 edgeConn->InsertCellPoint(ptId);
933 attrId = this->EdgeUseTable->InsertEdge(startId,ptId);
934 this->EdgeUses->InsertValue(2*attrId, this->Visited[p0]);
935 this->EdgeUses->InsertValue(2*attrId+1, -1);
936 startId = ptId;
937 }
938 }
939 edgeConn->InsertNextCell(2); //finish off the edge
940 edgeConn->InsertCellPoint(startId);
941 edgeConn->InsertCellPoint(2);
942 attrId = this->EdgeUseTable->InsertEdge(startId,2);
943 this->EdgeUses->InsertValue(2*attrId, this->Visited[dims[1]*dims[0]-1]);
944 this->EdgeUses->InsertValue(2*attrId+1, -1);
945
946 // Let's create perimeter edges - min y edge
947 startId = 0;
948 x[0] = origin[0];
949 for (j=0; j<(dims[1]-1); j++)
950 {
951 p0 = j*dims[0];
952 p1 = p0 + dims[0];
953 if ( this->Visited[p0] != this->Visited[p1] )
954 {
955 x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
956 ptId = points->InsertNextPoint(x);
957 this->EdgeTable->InsertEdge(p0,p1,ptId);
958 pointDescr->InsertValue(ptId, 1); //can't be smoothed
959
960 edgeConn->InsertNextCell(2);
961 edgeConn->InsertCellPoint(startId);
962 edgeConn->InsertCellPoint(ptId);
963 attrId = this->EdgeUseTable->InsertEdge(startId,ptId);
964 this->EdgeUses->InsertValue(2*attrId, this->Visited[p0]);
965 this->EdgeUses->InsertValue(2*attrId+1, -1);
966 startId = ptId;
967 }
968 }
969 edgeConn->InsertNextCell(2); //finish off the edge
970 edgeConn->InsertCellPoint(startId);
971 edgeConn->InsertCellPoint(3);
972 attrId = this->EdgeUseTable->InsertEdge(startId,3);
973 this->EdgeUses->InsertValue(2*attrId, this->Visited[(dims[1]-1)*dims[0]]);
974 this->EdgeUses->InsertValue(2*attrId+1, -1);
975
976 // Let's create perimeter edges - max y edge
977 startId = 1;
978 x[0] = origin[0] + (dims[0]-1)*spacing[0];
979 for (j=0; j<(dims[1]-1); j++)
980 {
981 p0 = j*dims[0] + (dims[0]-1);
982 p1 = p0 + dims[0];
983 if ( this->Visited[p0] != this->Visited[p1] )
984 {
985 x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
986 ptId = points->InsertNextPoint(x);
987 this->EdgeTable->InsertEdge(p0,p1,ptId);
988 pointDescr->InsertValue(ptId, 1); //can't be smoothed
989
990 edgeConn->InsertNextCell(2);
991 edgeConn->InsertCellPoint(startId);
992 edgeConn->InsertCellPoint(ptId);
993 attrId = this->EdgeUseTable->InsertEdge(startId,ptId);
994 this->EdgeUses->InsertValue(2*attrId, this->Visited[p0]);
995 this->EdgeUses->InsertValue(2*attrId+1, -1);
996 startId = ptId;
997 }
998 }
999 edgeConn->InsertNextCell(2); //finish off the edge
1000 edgeConn->InsertCellPoint(startId);
1001 edgeConn->InsertCellPoint(2);
1002 attrId = this->EdgeUseTable->InsertEdge(startId,2);
1003 this->EdgeUses->InsertValue(2*attrId, this->Visited[dims[1]*dims[0]-1]);
1004 this->EdgeUses->InsertValue(2*attrId+1, -1);
1005
1006 // Loop over all edges generating intersection points and outer boundary
1007 // edge segments.
1008 //
1009 for (j=1; j<(dims[1]-1); j++) //loop over all x edges (except boundary)
1010 {
1011 x[1] = origin[1] + j*spacing[1];
1012 for (i=0; i<(dims[0]-1); i++)
1013 {
1014 p0 = i + j*dims[0];
1015 p1 = p0 + 1;
1016 if ( this->Visited[p0] != this->Visited[p1] )
1017 {
1018 x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
1019 ptId = points->InsertNextPoint(x);
1020 this->EdgeTable->InsertEdge(p0,p1,ptId);
1021 pointDescr->InsertValue(ptId, 0);
1022 }
1023 }
1024 }
1025
1026 for (i=1; i<(dims[0]-1); i++) //loop over all y edges (except boundary)
1027 {
1028 x[0] = origin[0] + i*spacing[0];
1029 for (j=0; j<(dims[1]-1); j++)
1030 {
1031 p0 = i + j*dims[0];
1032 p1 = i + (j+1)*dims[0];
1033 if ( this->Visited[p0] != this->Visited[p1] )
1034 {
1035 x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
1036 ptId = points->InsertNextPoint(x);
1037 this->EdgeTable->InsertEdge(p0,p1,ptId);
1038 pointDescr->InsertValue(ptId, 0); //can be smoothed
1039 }
1040 }
1041 }
1042
1043 // All intersection points are generated, now create edges. Use a clipping
1044 // approach to create line segments. Later we'll connect them into polylines.
1045 //
1046 for (j=0; j<(dims[1]-1); j++) //loop over all x edges
1047 {
1048 for (i=0; i<(dims[0]-1); i++)
1049 {
1050 edgeCount = 0;
1051 p0 = i + j*dims[0];
1052 p1 = p0 + 1;
1053 p2 = i+1 + (j+1)*dims[0];
1054 p3 = p2 - 1;
1055
1056 if ( (ptId=this->EdgeTable->IsEdge(p0,p1)) != -1 )
1057 {
1058 id[2*edgeCount] = p0; id[2*edgeCount+1] = p1;
1059 pts[edgeCount++] = ptId;
1060 }
1061 if ( (ptId=this->EdgeTable->IsEdge(p1,p2)) != -1 )
1062 {
1063 id[2*edgeCount] = p1; id[2*edgeCount+1] = p2;
1064 pts[edgeCount++] = ptId;
1065 }
1066 if ( (ptId=this->EdgeTable->IsEdge(p2,p3)) != -1 )
1067 {
1068 id[2*edgeCount] = p2; id[2*edgeCount+1] = p3;
1069 pts[edgeCount++] = ptId;
1070 }
1071 if ( (ptId=this->EdgeTable->IsEdge(p3,p0)) != -1 )
1072 {
1073 id[2*edgeCount] = p3; id[2*edgeCount+1] = p0;
1074 pts[edgeCount++] = ptId;
1075 }
1076
1077 if ( edgeCount == 4 )
1078 {
1079 x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
1080 x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
1081 ptId = points->InsertNextPoint(x);
1082 pointDescr->InsertValue(ptId, 0); //intersection points are fixed
1083
1084 edgeConn->InsertNextCell(2);
1085 edgeConn->InsertCellPoint(ptId);
1086 edgeConn->InsertCellPoint(pts[0]);
1087 attrId = this->EdgeUseTable->InsertEdge(ptId,pts[0]);
1088 this->EdgeUses->InsertValue(2*attrId, this->Visited[id[0]]);
1089 this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[1]]);
1090
1091 edgeConn->InsertNextCell(2);
1092 edgeConn->InsertCellPoint(ptId);
1093 edgeConn->InsertCellPoint(pts[1]);
1094 attrId = this->EdgeUseTable->InsertEdge(ptId,pts[1]);
1095 this->EdgeUses->InsertValue(2*attrId, this->Visited[id[2]]);
1096 this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[3]]);
1097
1098 edgeConn->InsertNextCell(2);
1099 edgeConn->InsertCellPoint(ptId);
1100 edgeConn->InsertCellPoint(pts[2]);
1101 attrId = this->EdgeUseTable->InsertEdge(ptId,pts[2]);
1102 this->EdgeUses->InsertValue(2*attrId, this->Visited[id[4]]);
1103 this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[5]]);
1104
1105 edgeConn->InsertNextCell(2);
1106 edgeConn->InsertCellPoint(ptId);
1107 edgeConn->InsertCellPoint(pts[3]);
1108 attrId = this->EdgeUseTable->InsertEdge(ptId,pts[3]);
1109 this->EdgeUses->InsertValue(2*attrId, this->Visited[id[6]]);
1110 this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[7]]);
1111 }
1112
1113 else if ( edgeCount == 3 )
1114 {
1115 x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
1116 x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
1117 ptId = points->InsertNextPoint(x);
1118 pointDescr->InsertValue(ptId, 0); //intersection points are fixed
1119
1120 edgeConn->InsertNextCell(2);
1121 edgeConn->InsertCellPoint(ptId);
1122 edgeConn->InsertCellPoint(pts[0]);
1123 attrId = this->EdgeUseTable->InsertEdge(ptId, pts[0]);
1124 this->EdgeUses->InsertValue(2*attrId, this->Visited[id[0]]);
1125 this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[1]]);
1126
1127 edgeConn->InsertNextCell(2);
1128 edgeConn->InsertCellPoint(ptId);
1129 edgeConn->InsertCellPoint(pts[1]);
1130 attrId = this->EdgeUseTable->InsertEdge(ptId, pts[1]);
1131 this->EdgeUses->InsertValue(2*attrId, this->Visited[id[2]]);
1132 this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[3]]);
1133
1134 edgeConn->InsertNextCell(2);
1135 edgeConn->InsertCellPoint(ptId);
1136 edgeConn->InsertCellPoint(pts[2]);
1137 attrId = this->EdgeUseTable->InsertEdge(ptId, pts[2]);
1138 this->EdgeUses->InsertValue(2*attrId, this->Visited[id[4]]);
1139 this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[5]]);
1140 }
1141
1142 else if ( edgeCount == 2 )
1143 {
1144 edgeConn->InsertNextCell(2);
1145 edgeConn->InsertCellPoint(pts[0]);
1146 edgeConn->InsertCellPoint(pts[1]);
1147 attrId = this->EdgeUseTable->InsertEdge(pts[0],pts[1]);
1148 this->EdgeUses->InsertValue(2*attrId, this->Visited[id[0]]);
1149 this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[1]]);
1150 }
1151
1152 else if ( edgeCount == 1 )
1153 {
1154 vtkErrorMacro(<<"Bad mojo");
1155 return 0;
1156 }
1157 }
1158 }
1159
1160 // Cleanup
1161 this->EdgeUseTable->Delete();
1162 this->EdgeTable->Delete();
1163
1164 return 0;
1165 }
1166
BuildPolygons(vtkUnsignedCharArray * vtkNotUsed (pointDescr),vtkPolyData * edges,int numPolys,vtkUnsignedCharArray * polyColors)1167 void vtkImageToPolyDataFilter::BuildPolygons(vtkUnsignedCharArray *vtkNotUsed(pointDescr),
1168 vtkPolyData *edges, int numPolys,
1169 vtkUnsignedCharArray *polyColors)
1170 {
1171 vtkPoints *points = edges->GetPoints();
1172 vtkIdType numPts = points->GetNumberOfPoints(), ptId;
1173 int i, j, k, *polyId, *polyId2, edgeId;
1174 vtkIdType *cells, *pts, *cells2, npts, cellId;
1175 int numPolyPts, p1, p2;
1176 unsigned short ncells, ncells2;
1177 unsigned char *polyVisited, *ptr;
1178 vtkCellArray *newPolys;
1179
1180 // Make sure we can topological info
1181 edges->BuildLinks();
1182
1183 // Mark all polygons as unvisited
1184 polyVisited = new unsigned char [numPolys];
1185 for (i=0; i<numPolys; i++)
1186 {
1187 polyVisited[i] = 0;
1188 }
1189
1190 // Create connectivity array for polygon definition
1191 newPolys = vtkCellArray::New();
1192 newPolys->Allocate(newPolys->EstimateSize(numPolys,25));
1193
1194 // Loop over all edge points tracking around each polygon
1195 for (ptId=0; ptId<numPts; ptId++)
1196 {
1197 edges->GetPointCells(ptId, ncells, cells);
1198 if (ncells < 2)
1199 {
1200 vtkErrorMacro(<<"Bad mojo");
1201 return;
1202 }
1203 //for each edge, walk around polygon (if not visited before)
1204 for (i=0; i<ncells; i++)
1205 {
1206 edgeId = cells[i];
1207 polyId = this->EdgeUses->GetPointer(2*edgeId);
1208 for (j=0; j<2; j++)
1209 {
1210 if ( polyId[j] != -1 && !polyVisited[polyId[j]] )
1211 {//build loop
1212 polyVisited[polyId[j]] = 1;
1213 numPolyPts = 1;
1214 cellId = newPolys->InsertNextCell(0); //will update count later
1215 newPolys->InsertCellPoint(ptId);
1216
1217 // Update polygonal color
1218 ptr = this->PolyColors->GetPointer(3*polyId[j]);
1219 polyColors->SetValue(3*cellId, ptr[0]); //assign poly color
1220 polyColors->SetValue(3*cellId+1, ptr[1]);
1221 polyColors->SetValue(3*cellId+2, ptr[2]);
1222
1223 p1 = ptId;
1224 while ( true )
1225 {
1226 edges->GetCellPoints(edgeId, npts, pts);
1227 p2 = (pts[0] != p1 ? pts[0] : pts[1]);
1228 if (p2 == ptId)
1229 {
1230 break;
1231 }
1232
1233 newPolys->InsertCellPoint(p2);
1234 numPolyPts++;
1235 edges->GetPointCells(p2, ncells2, cells2);
1236 if (ncells < 2)
1237 {
1238 vtkErrorMacro(<<"Bad mojo");
1239 return;
1240 }
1241 for (k=0; k<ncells2; k++)
1242 {
1243 polyId2 = this->EdgeUses->GetPointer(2*cells2[k]);
1244 if ( cells2[k] != edgeId &&
1245 (polyId2[0] == polyId[j] || polyId2[1] == polyId[j]) )
1246 {
1247 p1 = p2;
1248 edgeId = cells2[k];
1249 break;
1250 }
1251 }
1252 }//while not completed loop
1253 newPolys->UpdateCellCount(numPolyPts);
1254
1255 }//if polygon not yet visited
1256 }//for each use of edge by polygon (at most 2 polygons)
1257 }//for each edge connected to this point
1258 }//for all points in edge list
1259
1260 edges->SetPolys(newPolys);
1261 newPolys->Delete();
1262 this->EdgeUses->Delete();
1263 delete [] polyVisited;
1264 }
1265
SmoothEdges(vtkUnsignedCharArray * pointDescr,vtkPolyData * edges)1266 void vtkImageToPolyDataFilter::SmoothEdges(vtkUnsignedCharArray *pointDescr,
1267 vtkPolyData *edges)
1268 {
1269
1270 vtkPoints *points=edges->GetPoints();
1271 vtkIdType numPts=points->GetNumberOfPoints(), ptId;
1272 int i, iterNum;
1273 int connId;
1274 double x[3], xconn[3], xave[3], factor;
1275 unsigned short int ncells;
1276 vtkIdType *cells, *pts, npts;
1277
1278
1279 // For each smoothing operation, loop over points. Points that can be
1280 // smoothed are moved in the direction of the average of their neighbor
1281 // points.
1282 for ( iterNum=0; iterNum<this->NumberOfSmoothingIterations; iterNum++ )
1283 {
1284 if ( (iterNum % 2) ) //alternate smoothing direction
1285 {
1286 factor = -0.331;
1287 }
1288 else
1289 {
1290 factor = 0.330;
1291 }
1292
1293 for (ptId=0; ptId<numPts; ptId++)
1294 {
1295 if ( pointDescr->GetValue(ptId) == 0 ) //can smooth
1296 {
1297 points->GetPoint(ptId, x);
1298 edges->GetPointCells(ptId, ncells, cells);
1299 xave[0] = xave[1] = xave[2] = 0.0;
1300 for (i=0; i<ncells; i++)
1301 {
1302 edges->GetCellPoints(cells[i], npts, pts);
1303 if (pts[0] != ptId)
1304 {
1305 connId = pts[0];
1306 }
1307 else if (npts > 1)
1308 {
1309 connId = pts[1];
1310 }
1311 else
1312 {
1313 vtkErrorMacro("Bad cell in smoothing operation");
1314 connId = pts[0];
1315 }
1316 points->GetPoint(connId, xconn);
1317 xave[0] += xconn[0]; xave[1] += xconn[1]; xave[2] += xconn[2];
1318 }
1319 if ( ncells > 0 )
1320 {
1321 xave[0] /= ncells; xave[1] /= ncells; xave[2] /= ncells;
1322 x[0] = x[0] + factor * (xave[0] - x[0]);
1323 x[1] = x[1] + factor * (xave[1] - x[1]);
1324 x[2] = x[2] + factor * (xave[2] - x[2]);
1325 points->SetPoint(ptId, x);
1326 }
1327
1328 }//if smoothable point
1329 }//for all points
1330 }//for all smoothing operations
1331 }
1332
1333
1334 // Remove points that are nearly co-linear to reduce the total point count
1335 //
DecimateEdges(vtkPolyData * edges,vtkUnsignedCharArray * pointDescr,double tol2)1336 void vtkImageToPolyDataFilter::DecimateEdges(vtkPolyData *edges,
1337 vtkUnsignedCharArray *pointDescr,
1338 double tol2)
1339 {
1340 vtkPoints *points=edges->GetPoints();
1341 vtkIdType numPts=points->GetNumberOfPoints(), ptId, prevId, nextId;
1342 vtkIdType npts;
1343 double x[3], xPrev[3], xNext[3];
1344 unsigned short int ncells;
1345 vtkIdType *cells, *pts;
1346
1347 // Loop over all points, finding those that are connected to just two
1348 // edges. If the point is colinear to the previous and next edge point,
1349 // then mark it as deleted.
1350 for (ptId=0; ptId<numPts; ptId++)
1351 {
1352 if ( pointDescr->GetValue(ptId) == 0 )
1353 {
1354 points->GetPoint(ptId, x);
1355 edges->GetPointCells(ptId, ncells, cells);
1356 if ( ncells == 2 )
1357 {
1358 edges->GetCellPoints(cells[0], npts, pts);
1359 prevId = (pts[0] != ptId ? pts[0] : pts[1]);
1360 points->GetPoint(prevId, xPrev);
1361
1362 edges->GetCellPoints(cells[1], npts, pts);
1363 nextId = (pts[0] != ptId ? pts[0] : pts[1]);
1364 points->GetPoint(nextId, xNext);
1365
1366 if ( vtkLine::DistanceToLine(x, xPrev, xNext) <= tol2 )
1367 {
1368 pointDescr->SetValue(ptId, 2); //mark deleted
1369 }
1370 }
1371 } //if manifold
1372 } //for all points
1373 }
1374