text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    $RCSfile: vtkImageToPolyDataFilter.cxx,v $
  
    Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
    All rights reserved.
    See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
  
       This software is distributed WITHOUT ANY WARRANTY; without even
       the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
       PURPOSE.  See the above copyright notice for more information.
  
  =========================================================================*/
  #include "vtkImageToPolyDataFilter.h"
  
  #include "vtkAppendPolyData.h"
  #include "vtkCellArray.h"
  #include "vtkCellData.h"
  #include "vtkEdgeTable.h"
  #include "vtkImageData.h"
  #include "vtkInformation.h"
  #include "vtkInformationVector.h"
  #include "vtkLine.h"
  #include "vtkObjectFactory.h"
  #include "vtkPointData.h"
  #include "vtkPolyData.h"
  #include "vtkScalarsToColors.h"
  #include "vtkUnsignedCharArray.h"
  
  vtkCxxRevisionMacro(vtkImageToPolyDataFilter, "$Revision: 1.31 $");
  vtkStandardNewMacro(vtkImageToPolyDataFilter);
  
  vtkCxxSetObjectMacro(vtkImageToPolyDataFilter,LookupTable,vtkScalarsToColors);
  
  vtkImageToPolyDataFilter::vtkImageToPolyDataFilter()
  {
    this->OutputStyle = VTK_STYLE_POLYGONALIZE;
    this->ColorMode = VTK_COLOR_MODE_LINEAR_256;
    this->Smoothing = 1;
    this->NumberOfSmoothingIterations = 40;
    this->Decimation = 1;
    this->DecimationError = 1.5;
    this->Error = 100;
    this->SubImageSize = 250;
  
    this->Table = vtkUnsignedCharArray::New();
    this->LookupTable = NULL;
  }
  
  vtkImageToPolyDataFilter::~vtkImageToPolyDataFilter()
  {
    this->Table->Delete();
    if ( this->LookupTable )
      {
      this->LookupTable->Delete();
      }
  }
  
  // declare helper functions
  //
  
  int vtkImageToPolyDataFilter::RequestData(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector)
  {
    // get the info objects
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
  
    // get the input and ouptut
    vtkImageData *input = vtkImageData::SafeDownCast(
      inInfo->Get(vtkDataObject::DATA_OBJECT()));
    vtkPolyData *output = vtkPolyData::SafeDownCast(
      outInfo->Get(vtkDataObject::DATA_OBJECT()));
  
    vtkPolyData *tmpOutput;
    vtkPolyData *tmpInput;
    vtkAppendPolyData *append;
    vtkPolyData *appendOutput;
    vtkDataArray *inScalars = input->GetPointData()->GetScalars();
    vtkIdType numPixels=input->GetNumberOfPoints();
    int dims[3], numComp;
    double origin[3], spacing[3];
    vtkUnsignedCharArray *pixels;
    int type;
    int numPieces[2], extent[4];
    int i, j, newDims[3], totalPieces, pieceNum, abortExecute=0;
    double newOrigin[3];
  
    // Check input and initialize
    vtkDebugMacro(<<"Vectorizing image...");
    
    if ( inScalars == NULL || numPixels < 1 )
      {
      vtkDebugMacro(<<"Not enough input to create output");
      return 1;
      }
  
    append = vtkAppendPolyData::New();
    tmpOutput=vtkPolyData::New();
    tmpInput=vtkPolyData::New();    
    numComp=inScalars->GetNumberOfComponents();
    type=inScalars->GetDataType();
  
    appendOutput=append->GetOutput();
  
    input->GetDimensions(dims);
    input->GetOrigin(origin);
    input->GetSpacing(spacing);
    
    // Figure out how many pieces to break the image into (the image
    // might be too big to process). The filter does a series of appends
    // to join the pieces together.
    numPieces[0] = ((dims[0]-2) / this->SubImageSize) + 1;
    numPieces[1] = ((dims[1]-2) / this->SubImageSize) + 1;
    totalPieces = numPieces[0]*numPieces[1];
  
    appendOutput->Initialize(); //empty the output
    append->AddInput(tmpOutput); //output of piece
    append->AddInput(tmpInput); //output of previoius append
  
    // Loop over this many pieces
    for (pieceNum=j=0; j < numPieces[1] && !abortExecute; j++)
      {
      extent[2] = j*this->SubImageSize; //the y range
      extent[3] = (j+1)*this->SubImageSize;
      if ( extent[3] >= dims[1] )
        {
        extent[3] = dims[1] - 1;
        }
      
      for (i=0; i < numPieces[0] && !abortExecute; i++)
        {
        extent[0] = i*this->SubImageSize; //the x range
        extent[1] = (i+1)*this->SubImageSize;
        if ( extent[1] >= dims[0] )
          {
          extent[1] = dims[0] - 1;
          }
  
        vtkDebugMacro(<<"Processing #" << pieceNum);
        this->UpdateProgress ((double)pieceNum/totalPieces);
        if (this->GetAbortExecute())
          {
          abortExecute = 1;
          break;
          }
        pieceNum++;
  
        // Figure out characteristics of current sub-image
        newDims[0] = extent[1] - extent[0] + 1;
        newDims[1] = extent[3] - extent[2] + 1;
        newOrigin[0] = origin[0] + extent[0]*spacing[0];
        newOrigin[1] = origin[1] + extent[2]*spacing[1];
        newOrigin[2] = 0.0;
  
        // Create a quantized copy of the image based on the color table
        //
        pixels = this->QuantizeImage(inScalars, numComp, type, dims, extent);
        vtkDebugMacro(<<"Quantizing color...image size (" <<newDims[0] 
                      <<", " <<newDims[1] << ")");
  
        // Generate polygons according to mode setting
        //
        if ( this->OutputStyle == VTK_STYLE_PIXELIZE )
          {
          this->PixelizeImage(pixels, newDims, newOrigin, spacing, tmpOutput);
          }
        else if ( this->OutputStyle == VTK_STYLE_RUN_LENGTH )
          {
          this->RunLengthImage(pixels, newDims, newOrigin, spacing, tmpOutput);
          }
        else //VTK_STYLE_POLYGONALIZE
          {
          this->PolygonalizeImage(pixels, newDims, newOrigin, spacing,tmpOutput);
          }
  
        // Append pieces together
        //
        tmpInput->CopyStructure(appendOutput);
        tmpInput->GetPointData()->PassData(appendOutput->GetPointData());
        tmpInput->GetCellData()->PassData(appendOutput->GetCellData());
        append->Update();
  
        // Clean up this iteration
        //
        pixels->Delete();
        tmpInput->Initialize();
        tmpOutput->Initialize();
        } // for i pieces
      } // for j pieces
    
    // Create the final output contained in the append filter
    output->CopyStructure(appendOutput);
    output->GetPointData()->PassData(appendOutput->GetPointData());
    output->GetCellData()->PassData(appendOutput->GetCellData());
  
    append->Delete();
    tmpInput->Delete();
    tmpOutput->Delete();
  
    return 1;
  }
  
  void vtkImageToPolyDataFilter::PixelizeImage(vtkUnsignedCharArray *pixels, 
                                               int dims[3], double origin[3], 
                                               double spacing[3], 
                                               vtkPolyData *output)
  {
    int numPts, numCells, i, j, id;
    vtkIdType pts[4];
    vtkPoints *newPts;
    vtkCellArray *newPolys;
    double x[3];
    vtkUnsignedCharArray *polyColors;
    unsigned char *ptr, *colors=pixels->GetPointer(0);
    
    // create the points - see whether to create or append
    numPts = (dims[0]+1) * (dims[1]+1);
    newPts = vtkPoints::New();
    newPts->SetNumberOfPoints(numPts);
    
    x[2] = 0.0;
    for (id=0, j=0; j<=dims[1]; j++)
      {
      x[1] = origin[1] + j*spacing[1];
      for (i=0; i<=dims[0]; i++)
        {
        x[0] = origin[0] + i*spacing[0];
        newPts->SetPoint(id, x);
        id++;
        }
      }
    output->SetPoints(newPts);
    newPts->Delete();
  
    // create the cells and cell colors
    //
    numCells = dims[0] * dims[1];
    newPolys = vtkCellArray::New();
    newPolys->Allocate(newPolys->EstimateSize(numCells,4));
  
    polyColors = vtkUnsignedCharArray::New();
    polyColors->SetNumberOfValues(3*numCells); //for rgb
    polyColors->SetNumberOfComponents(3);
  
    // loop over all pixels, creating a quad per pixel.
    // Note: copying point data (pixel values) to cell data (quad colors).
    for (id=0, j=0; j<dims[1]; j++)
      {
      for (i=0; i<dims[0]; i++)
        {
        pts[0] = i + j*(dims[0]+1);
        pts[1] = pts[0] + 1;
        pts[2] = pts[1] + dims[0] + 1;
        pts[3] = pts[2] - 1;
        newPolys->InsertNextCell(4, pts);
        ptr = colors + 3*id;
        polyColors->SetValue(3*id, ptr[0]);
        polyColors->SetValue(3*id+1, ptr[1]);
        polyColors->SetValue(3*id+2, ptr[2]);
        id++;
        }
      }
  
    output->SetPolys(newPolys);
    newPolys->Delete();
  
    output->GetCellData()->SetScalars(polyColors);
    polyColors->Delete();
  }
  
  void vtkImageToPolyDataFilter::RunLengthImage(vtkUnsignedCharArray *pixels, 
                                               int dims[3], double origin[3], 
                                               double spacing[3], 
                                               vtkPolyData *output)
  {
    int i, j;
    vtkIdType pts[4], id;
    vtkPoints *newPts;
    vtkCellArray *newPolys;
    double x[3], minX, maxX, minY, maxY;
    vtkUnsignedCharArray *polyColors;
    unsigned char *ptr, *colors=pixels->GetPointer(0), *color;
    
    // Setup data
    newPts = vtkPoints::New();
    newPolys = vtkCellArray::New();
    newPolys->Allocate(newPolys->EstimateSize(dims[0]*dims[1]/10,4));
  
    polyColors = vtkUnsignedCharArray::New();
    polyColors->Allocate(3*dims[0]*dims[1]/10); //for rgb
    polyColors->SetNumberOfComponents(3);
  
    // Loop over row-by-row generating quad polygons
    x[2] = 0.0;
    for (j=0; j<dims[1]; j++)
      {
      if ( j == 0 )
        {
        minY = origin[1];
        maxY = origin[1] + 0.5*spacing[1];
        }
      else if ( j == (dims[1]-1) )
        {
        minY = origin[1] + j*spacing[1] - 0.5*spacing[1];
        maxY = origin[1] + j*spacing[1];
        }
      else
        {
        minY = origin[1] + j*spacing[1] - 0.5*spacing[1];
        maxY = origin[1] + j*spacing[1] + 0.5*spacing[1];
        }
      
      for ( i=0; i < dims[0]; )
        {
        if ( i == 0 )
          {
          minX = origin[0];
          }
        else
          {
          minX = origin[0] + i*spacing[0] - 0.5*spacing[0];
          }
        color = colors + 3*(i+j*dims[0]);
        ptr = color;
        while ( i < dims[0] )
          {
          ptr = colors + 3*(i+j*dims[0]);
          if ( ! this->IsSameColor(color,ptr) )
            {
            break;
            }
          else
            {
            i++;
            }
          }
  
        if ( i >= dims[0] )
          {
          maxX = origin[0] + (dims[0]-1)*spacing[0];
          }
        else
          {
          maxX = origin[0] + (i-1)*spacing[0] + 0.5*spacing[0];
          }
  
        // Create quad cell
        x[0] = minX;
        x[1] = minY;
        pts[0] = newPts->InsertNextPoint(x);
        x[0] = maxX;
        pts[1] = newPts->InsertNextPoint(x);
        x[1] = maxY;
        pts[2] = newPts->InsertNextPoint(x);
        x[0] = minX;
        pts[3] = newPts->InsertNextPoint(x);
        id = newPolys->InsertNextCell(4,pts);
        polyColors->InsertValue(3*id, color[0]);
        polyColors->InsertValue(3*id+1, color[1]);
        polyColors->InsertValue(3*id+2, color[2]);
        }
      }
    
    output->SetPoints(newPts);
    newPts->Delete();
    output->SetPolys(newPolys);
    newPolys->Delete();
  
    output->GetCellData()->SetScalars(polyColors);
    polyColors->Delete();
  }
  
  
  void vtkImageToPolyDataFilter::PolygonalizeImage(vtkUnsignedCharArray *pixels, 
                                 int dims[3], double origin[3], double spacing[3], 
                                 vtkPolyData *output)
  {
    int numPolys;
    int numPixels=dims[0]*dims[1];
  
    // Perform connected traversal on quantized points. This builds
    // the initial "polygons" in implicit form.
    //
    this->PolyColors = vtkUnsignedCharArray::New();
    this->PolyColors->SetNumberOfComponents(3);
    this->PolyColors->Allocate(5000);
  
    numPolys = this->ProcessImage(pixels, dims);
    vtkDebugMacro(<<"Visited regions..." << numPolys << " polygons");
  
    // Build edges around the boundary of the polygons. Also identify
    // junction points where 3 or 4 polygons meet.
    // 
    vtkPoints *points =  vtkPoints::New();
    points->Allocate(numPixels/2, numPixels/2);
  
    vtkUnsignedCharArray *pointDescr = vtkUnsignedCharArray::New();
    pointDescr->Allocate(numPixels/2, numPixels/2);
  
    vtkCellArray *edgeConn = vtkCellArray::New();
    edgeConn->Allocate(numPixels/2, numPixels/2);
    vtkPolyData *edges = vtkPolyData::New();
    edges->SetPoints(points);
    edges->SetLines(edgeConn);
    points->Delete();
    edgeConn->Delete();
  
    this->BuildEdges(pixels, dims, origin, spacing, pointDescr, edges);
    vtkDebugMacro(<<"Edges built...");
  
    // Now that we've got the edges, we have to build the "loops" around the
    // polygons that define the polygon explicitly.
    //
    vtkUnsignedCharArray *polyColors = vtkUnsignedCharArray::New();
    polyColors->SetNumberOfComponents(3);
    polyColors->SetNumberOfValues(numPolys*3);
  
    this->BuildPolygons(pointDescr, edges, numPolys, polyColors);
    this->PolyColors->Delete();
    delete [] this->Visited;
    vtkDebugMacro(<<"Constructed polygons...");
  
    // Smooth edge network. Some points are identified as fixed, others 
    // move using Laplacian smoothing.
    // 
    if ( this->Smoothing )
      {
      this->SmoothEdges(pointDescr, edges);
      vtkDebugMacro(<<"Edges smoothed...");
      }
  
    // Decimate edge network. There will be colinear vertices along edges.
    // These are eliminated.
    //
    if ( this->Decimation )
      {
      this->DecimateEdges(edges, pointDescr, this->DecimationError);
      }
  
    // Create output polydata. Each polyon is output with its edges.
    //
    this->GeneratePolygons(edges, numPolys, output, polyColors, pointDescr);
    vtkDebugMacro(<<"Output generated...");
  
    // clean up and get out 
    edges->Delete();
    polyColors->Delete();
    pointDescr->Delete();
  }
  
  // The following are private helper functions----------------------------------
  //
  vtkUnsignedCharArray *vtkImageToPolyDataFilter::QuantizeImage(
                            vtkDataArray *inScalars, int numComp, int type,
                            int dims[3], int extent[4])
  {
    int numPixels, i, j, idx, id;
    vtkUnsignedCharArray *pixels;
    unsigned char *ptr, *ptr2, *color, *outPixels;
    unsigned char *inPixels;
    
    // doing a portion of the image
    numPixels = (extent[1]-extent[0]+1) * (extent[3]-extent[2]+1);
    pixels = vtkUnsignedCharArray::New();
    pixels->SetNumberOfValues(3*numPixels);
    outPixels = pixels->GetPointer(0);
    
    // Figure out how to quantize
    //
    if ( this->ColorMode == VTK_COLOR_MODE_LINEAR_256 )
      {
      // Check scalar type
      if ( type != VTK_UNSIGNED_CHAR || numComp != 3 )
        {
        vtkErrorMacro(<<"Wrong input scalar type");
        return 0;
        }
      else
        {
        inPixels = static_cast<vtkUnsignedCharArray *>(inScalars)->GetPointer(0);
        }
  
      // Generate a color table used to quantize the points
      //
      if ( this->GetMTime() > this->TableMTime )
        {
        this->BuildTable(inPixels);
        }
  
      for (id=0, j=extent[2]; j <= extent[3]; j++)
        {
        for (i=extent[0]; i <= extent[1]; i++)
          {
          idx = i + j*dims[0];
          ptr = inPixels + 3*idx;
          ptr2 = outPixels + 3*id;
          color = this->GetColor(ptr);
          ptr2[0] = color[0];
          ptr2[1] = color[1];
          ptr2[2] = color[2];
          id++;
          }
        }
      }//using build in table
    
    else //using provided lookup table
      {
      if ( numComp != 1 || this->LookupTable == NULL )
        {
        vtkErrorMacro(<<"LUT mode requires single component scalar and LUT");
        return 0;
        }
  
      double s;
      for (id=0, j=extent[2]; j <= extent[3]; j++)
        {
        for (i=extent[0]; i <= extent[1]; i++)
          {
          idx = i + j*dims[0];
          s = inScalars->GetComponent(idx,0);
          color = this->LookupTable->MapValue(s);
          ptr2 = outPixels + 3*id;
          ptr2[0] = color[0];
          ptr2[1] = color[1];
          ptr2[2] = color[2];
          id++;
          }
        }
      }
  
    return pixels;
  }
  
  void vtkImageToPolyDataFilter::BuildTable(unsigned char *vtkNotUsed(inPixels))
  {
    int red, green, blue, idx=0; 
  
    this->Table->SetNumberOfValues(256*3);
  
    // use 3-3-2 bits for rgb
    for (blue=0; blue<256; blue+=64)
      {
      for (green=0; green<256; green+=32)
        {
        for (red=0; red<256; red+=32)
          {
          this->Table->SetValue(idx++, red);
          this->Table->SetValue(idx++, green);
          this->Table->SetValue(idx++, blue);
          }
        }
      }
  }
  
  int vtkImageToPolyDataFilter::FillInputPortInformation(int, vtkInformation *info)
  {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
    return 1;
  }
  
  void vtkImageToPolyDataFilter::PrintSelf(ostream& os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os,indent);
   
    os << indent << "Output Style: ";
    if ( this->OutputStyle == VTK_STYLE_PIXELIZE )
      {
      os << indent << "Pixelize\n";
      }
    else if ( this->OutputStyle == VTK_STYLE_RUN_LENGTH )
      {
      os << indent << "RunLength\n";
      }
    else // this->OutputStyle == VTK_STYLE_POLYGONALIZE
      {
      os << indent << "Polygonalize\n";
      }
  
    os << indent << "Color Mode: ";
    if ( this->ColorMode == VTK_STYLE_PIXELIZE )
      {
      os << indent << "LUT\n";
      }
    else // this->ColorMode == VTK_STYLE_POLYGONALIZE
      {
      os << indent << "Linear256\n";
      }
  
    os << indent << "Smoothing: " << (this->Smoothing ? "On\n" : "Off\n");
    os << indent << "Number of Smoothing Iterations: " 
       << this->NumberOfSmoothingIterations << "\n";
  
    os << indent << "Decimation: " << (this->Decimation ? "On\n" : "Off\n");
    os << indent << "Decimation Error: " 
       << (this->DecimationError ? "On\n" : "Off\n");
  
    os << indent << "Error: " << this->Error << "\n";
    os << indent << "Sub-Image Size: " << this->SubImageSize << "\n";
  
    if ( this->LookupTable )
      {
      os << indent << "LookupTable:\n";
      this->LookupTable->PrintSelf(os,indent.GetNextIndent());
      }
    else
      {
      os << indent << "LookupTable: (none)\n";
      }
  
  }
  
  
  //--------------------------private helper functions---------------------------
  // Determines whether two pixels are the same color
  int vtkImageToPolyDataFilter::IsSameColor(unsigned char *p1, unsigned char *p2)
  {
    int d2 = (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) + 
             (p1[2]-p2[2])*(p1[2]-p2[2]);
    
    return (d2 > this->Error ? 0 : 1);
  }
  
  unsigned char *vtkImageToPolyDataFilter::GetColor(unsigned char *rgb)
  {
    // round to nearest value
    int red = (rgb[0] + 16) / 32;
    red = (red > 7 ? 7 : red);
    int green =(rgb[1] + 16) / 32;
    green = (green > 7 ? 7 : green);
    int blue = (rgb[2] + 32) / 64;
    blue = (blue > 3 ? 3 : blue);
  
    return this->Table->GetPointer(3*(red + green*8 + blue*64));
  }
  
  void vtkImageToPolyDataFilter::GetIJ(int id, int &i, int &j, int dims[3])
  {
    i = id % dims[0];
    j = id / dims[0];
  }
  
  // Get the left-right-top-bottom neighboring pixels of a given pixel
  // The method has been modified to return right neighbor (mode==0);
  // or top neighbor (mode==1) or all neighbors (mode==2).
  int vtkImageToPolyDataFilter::GetNeighbors(unsigned char *ptr, int &i, int &j, 
                         int dims[2],unsigned char *neighbors[4], int mode)
  {
    int numNeis=0;
  
    if ( mode == 0 )
      {
      if ( (i+1) < dims[0] )
        {
        neighbors[numNeis++] = ptr + 3; //jump over rgb
        }
      if ( (i-1) >= 0 )
        {
        neighbors[numNeis++] = ptr - 3; //jump over rgb
        }
      }
    
    else if ( mode == 1 )
      {
      if ( (j+1) < dims[1] )
        {
        neighbors[numNeis++] = ptr + 3*dims[0];
        }
      }
    
    else
      {
      if ( (i+1) < dims[0] )
        {
        neighbors[numNeis++] = ptr + 3; //jump over rgb
        }
      if ( (i-1) >= 0 )
        {
        neighbors[numNeis++] = ptr - 3;
        }
  
      if ( (j+1) < dims[1] )
        {
        neighbors[numNeis++] = ptr + 3*dims[0];
        }
      if ( (j-1) >= 0 )
        {
        neighbors[numNeis++] = ptr - 3*dims[0];
        }
      }
  
    return numNeis;
    
  }
  
  
  // Marks connected regions with different colors.
  int vtkImageToPolyDataFilter::ProcessImage(vtkUnsignedCharArray *scalars, 
                                             int dims[2])
  {
    int numPixels = dims[0]*dims[1];
    vtkIdList *wave, *wave2, *tmpWave;
    int numIds, regionNumber, i, j, k, id, x, y, numNeighbors;
    unsigned char *neighbors[4], *ptr;
    unsigned char *pixels = scalars->GetPointer(0);
  
    // Collect groups of pixels together into similar colored regions. These
    // will be eventually grouped into polygons and/or lines.
    //
    // mark all pixels unvisited
    regionNumber = -1;
    this->Visited = new int [numPixels];
    memset (this->Visited, (int)-1, numPixels*sizeof(int));
    
    // set up the connected traversal
    wave = vtkIdList::New();
    wave->Allocate(static_cast<int>(numPixels/4.0), 
                   static_cast<int>(numPixels/4.0));
    wave2 = vtkIdList::New();
    wave2->Allocate(static_cast<int>(numPixels/4.0), 
                    static_cast<int>(numPixels/4.0));
  
    // visit connected pixels. Pixels are connected if they are topologically
    // adjacent and they have "equal" color values.
    for (i=0; i < numPixels; i++)
      {
      if ( this->Visited[i] == -1 )
        {//start a connected wave
        this->Visited[i] = ++regionNumber;
        ptr = pixels + 3*i;
        this->PolyColors->InsertValue(3*regionNumber, ptr[0]); //assign color
        this->PolyColors->InsertValue(3*regionNumber+1, ptr[1]);
        this->PolyColors->InsertValue(3*regionNumber+2, ptr[2]);
        wave->Reset(); wave2->Reset();
  
        // To prevent creating polygons with inner loops, we're going to start
        // the wave as a "vertical" stack of pixels, and then propogate the
        // wave horizontally only.
        wave->InsertId(0,i);
        this->GetIJ(i, x, y, dims);
        while ( (numNeighbors = this->GetNeighbors(ptr, x, y, dims, neighbors, 1)) )
          {
          id = (neighbors[0] - pixels) / 3;
          if ( this->Visited[id] == -1 && this->IsSameColor(ptr, neighbors[0]) )
            {
            this->Visited[id] = regionNumber;
            wave->InsertNextId(id);
            ptr = pixels + 3*id;
            this->GetIJ(id, x, y, dims);
            }
          else
            {
            break;
            }
          }
  
        // Okay, defined vertical wave, now propogate horizontally
        numIds = wave->GetNumberOfIds();
        while ( numIds > 0 )
          {
          for (j=0; j<numIds; j++) //propogate wave
            {
            id = wave->GetId(j);
            ptr = pixels + 3*id;
            this->GetIJ(id, x, y, dims);
            numNeighbors = this->GetNeighbors(ptr, x, y, dims, neighbors, 0);
            for (k=0; k<numNeighbors; k++)
              {
              id = (neighbors[k] - pixels) / 3;
              if ( this->Visited[id] == -1 && this->IsSameColor(ptr, neighbors[k]) )
                {
                this->Visited[id] = regionNumber;
                wave2->InsertNextId(id);
                }
              }//for each pixel neighbor
            }//for pixels left in wave
          numIds = wave2->GetNumberOfIds();
          tmpWave = wave;
          wave = wave2;
          wave2 = tmpWave;
          wave2->Reset();
          }//while still propogating
        }//if not, start wave
      }//for all pixels
    
    wave->Delete();
    wave2->Delete();
  
    return regionNumber+1;
  }
  
  // Create polygons and place into output
  void vtkImageToPolyDataFilter::GeneratePolygons(vtkPolyData *edges, 
                                 int vtkNotUsed(numPolys), vtkPolyData *output, 
                                 vtkUnsignedCharArray *polyColors,
                                 vtkUnsignedCharArray *pointDescr)
  {
    vtkCellArray *newPolys, *inPolys;
    int i, numPts;
    vtkIdType *pts = 0;
    vtkIdType npts = 0;
    
    // Copy the points via reference counting
    //
    output->SetPoints(edges->GetPoints());
  
    // Create the polygons - points may have been decimated so these
    // points have to be culled.
    //
    inPolys = edges->GetPolys();
    newPolys = vtkCellArray::New();
    newPolys->Allocate(inPolys->GetSize());
  
    for ( inPolys->InitTraversal(); inPolys->GetNextCell(npts,pts); )
      {
      newPolys->InsertNextCell(0);
      numPts = 0;
      for (i=0; i<npts; i++)
        {
        if ( pointDescr->GetValue(pts[i]) != 2 )
          {
          newPolys->InsertCellPoint(pts[i]);
          numPts++;
          }
        }
      newPolys->UpdateCellCount(numPts);
      }
  
    output->SetPolys(newPolys);
    newPolys->Delete();
    
    output->GetCellData()->SetScalars(polyColors);
  }
  
  // Uses clipping approach to build the polygon edges
  int vtkImageToPolyDataFilter::BuildEdges(vtkUnsignedCharArray *vtkNotUsed(pixels), 
                                           int dims[3], double origin[3],
                                           double spacing[3],
                                           vtkUnsignedCharArray *pointDescr, 
                                           vtkPolyData *edges)
  {
    double x[3];
    int i, j, edgeCount;
    vtkIdType ptId, p0, p1, p2, p3, startId, attrId, id[8], pts[4];
    vtkCellArray *edgeConn = edges->GetLines();
    vtkPoints *points = edges->GetPoints();
  
    // Build edges around perimeter of image. Note that the point ids
    // The first four points are the image corners and are inserted and 
    // marked so that they can't be moved during smoothing.
    points->InsertPoint(0, origin);
    pointDescr->InsertValue(0, 1);
  
    // Keep track of the polygons that use each edge as well as associated
    // intersection points on edge (if any)
    this->EdgeTable = vtkEdgeTable::New();
    this->EdgeTable->InitEdgeInsertion(dims[0]*dims[1],1);
  
    this->EdgeUseTable = vtkEdgeTable::New();
    this->EdgeUseTable->InitEdgeInsertion(dims[0]*dims[1],1);
  
    this->EdgeUses = vtkIntArray::New();
    this->EdgeUses->SetNumberOfComponents(2);
    this->EdgeUses->Allocate(4*dims[0]*dims[1],dims[0]*dims[1]);
  
    // Generate corner points of image
    x[0] = origin[0] + (dims[0]-1)*spacing[0];
    x[1] = origin[1];
    x[2] = 0.0;
    points->InsertPoint(1, x);
    pointDescr->InsertValue(1, 1);
  
    x[0] = origin[0] + (dims[0]-1)*spacing[0];
    x[1] = origin[1] + (dims[1]-1)*spacing[1];;
    x[2] = 0.0;
    points->InsertPoint(2, x);
    pointDescr->InsertValue(2, 1);
  
    x[0] = origin[0];
    x[1] = origin[1] + (dims[1]-1)*spacing[1];
    x[2] = 0.0;
    points->InsertPoint(3, x);
    pointDescr->InsertValue(3, 1);
    
    // Let's create perimeter edges - bottom x edge
    startId = 0;
    x[1] = origin[1];
    for (i=0; i<(dims[0]-1); i++)
      {
      p0 = i;
      p1 = i + 1;
      if ( this->Visited[p0] != this->Visited[p1] )
        {
        x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
        ptId = points->InsertNextPoint(x);
        this->EdgeTable->InsertEdge(p0,p1,ptId);
        pointDescr->InsertValue(ptId, 1); //can't be smoothed
        
        edgeConn->InsertNextCell(2);
        edgeConn->InsertCellPoint(startId);
        edgeConn->InsertCellPoint(ptId);
        attrId = this->EdgeUseTable->InsertEdge(startId,ptId);
        this->EdgeUses->InsertValue(2*attrId, this->Visited[p0]);
        this->EdgeUses->InsertValue(2*attrId+1, -1);
        startId = ptId;
        }
      }
    edgeConn->InsertNextCell(2); //finish off the edge
    edgeConn->InsertCellPoint(startId);
    edgeConn->InsertCellPoint(1);
    attrId = this->EdgeUseTable->InsertEdge(startId,1);
    this->EdgeUses->InsertValue(2*attrId, this->Visited[dims[0]-1]);
    this->EdgeUses->InsertValue(2*attrId+1, -1);
  
    // Let's create perimeter edges - top x edge
    startId = 3;
    x[1] = origin[1] + (dims[1]-1)*spacing[1];
    for (i=0; i<(dims[0]-1); i++)
      {
      p0 = i + dims[0]*(dims[1]-1);
      p1 = p0 + 1;
      if ( this->Visited[p0] != this->Visited[p1] )
        {
        x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
        ptId = points->InsertNextPoint(x);
        this->EdgeTable->InsertEdge(p0,p1,ptId);
        pointDescr->InsertValue(ptId, 1); //can't be smoothed
        
        edgeConn->InsertNextCell(2);
        edgeConn->InsertCellPoint(startId);
        edgeConn->InsertCellPoint(ptId);
        attrId = this->EdgeUseTable->InsertEdge(startId,ptId);
        this->EdgeUses->InsertValue(2*attrId, this->Visited[p0]);
        this->EdgeUses->InsertValue(2*attrId+1, -1);
        startId = ptId;
        }
      }
    edgeConn->InsertNextCell(2); //finish off the edge
    edgeConn->InsertCellPoint(startId);
    edgeConn->InsertCellPoint(2);
    attrId = this->EdgeUseTable->InsertEdge(startId,2);
    this->EdgeUses->InsertValue(2*attrId, this->Visited[dims[1]*dims[0]-1]);
    this->EdgeUses->InsertValue(2*attrId+1, -1);
  
    // Let's create perimeter edges - min y edge
    startId = 0;
    x[0] = origin[0];
    for (j=0; j<(dims[1]-1); j++)
      {
      p0 = j*dims[0];
      p1 = p0 + dims[0];
      if ( this->Visited[p0] != this->Visited[p1] )
        {
        x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
        ptId = points->InsertNextPoint(x);
        this->EdgeTable->InsertEdge(p0,p1,ptId);
        pointDescr->InsertValue(ptId, 1); //can't be smoothed
        
        edgeConn->InsertNextCell(2);
        edgeConn->InsertCellPoint(startId);
        edgeConn->InsertCellPoint(ptId);
        attrId = this->EdgeUseTable->InsertEdge(startId,ptId);
        this->EdgeUses->InsertValue(2*attrId, this->Visited[p0]);
        this->EdgeUses->InsertValue(2*attrId+1, -1);
        startId = ptId;
        }
      }
    edgeConn->InsertNextCell(2); //finish off the edge
    edgeConn->InsertCellPoint(startId);
    edgeConn->InsertCellPoint(3);
    attrId = this->EdgeUseTable->InsertEdge(startId,3);
    this->EdgeUses->InsertValue(2*attrId, this->Visited[(dims[1]-1)*dims[0]]);
    this->EdgeUses->InsertValue(2*attrId+1, -1);
  
    // Let's create perimeter edges - max y edge
    startId = 1;
    x[0] = origin[0] + (dims[0]-1)*spacing[0];
    for (j=0; j<(dims[1]-1); j++)
      {
      p0 = j*dims[0] + (dims[0]-1);
      p1 = p0 + dims[0];
      if ( this->Visited[p0] != this->Visited[p1] )
        {
        x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
        ptId = points->InsertNextPoint(x);
        this->EdgeTable->InsertEdge(p0,p1,ptId);
        pointDescr->InsertValue(ptId, 1); //can't be smoothed
        
        edgeConn->InsertNextCell(2);
        edgeConn->InsertCellPoint(startId);
        edgeConn->InsertCellPoint(ptId);
        attrId = this->EdgeUseTable->InsertEdge(startId,ptId);
        this->EdgeUses->InsertValue(2*attrId, this->Visited[p0]);
        this->EdgeUses->InsertValue(2*attrId+1, -1);
        startId = ptId;
        }
      }
    edgeConn->InsertNextCell(2); //finish off the edge
    edgeConn->InsertCellPoint(startId);
    edgeConn->InsertCellPoint(2);
    attrId = this->EdgeUseTable->InsertEdge(startId,2);
    this->EdgeUses->InsertValue(2*attrId, this->Visited[dims[1]*dims[0]-1]);
    this->EdgeUses->InsertValue(2*attrId+1, -1);
  
    // Loop over all edges generating intersection points and outer boundary
    // edge segments.
    //
    for (j=1; j<(dims[1]-1); j++) //loop over all x edges (except boundary)
      {
      x[1] = origin[1] + j*spacing[1];
      for (i=0; i<(dims[0]-1); i++)
        {
        p0 = i + j*dims[0];
        p1 = p0 + 1;
        if ( this->Visited[p0] != this->Visited[p1] )
          {
          x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
          ptId = points->InsertNextPoint(x);
          this->EdgeTable->InsertEdge(p0,p1,ptId);
          pointDescr->InsertValue(ptId, 0);
          }
        }
      }
    
    for (i=1; i<(dims[0]-1); i++) //loop over all y edges (except boundary)
      {
      x[0] = origin[0] + i*spacing[0];
      for (j=0; j<(dims[1]-1); j++)
        {
        p0 = i + j*dims[0];
        p1 = i + (j+1)*dims[0];
        if ( this->Visited[p0] != this->Visited[p1] )
          {
          x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
          ptId = points->InsertNextPoint(x);
          this->EdgeTable->InsertEdge(p0,p1,ptId);
          pointDescr->InsertValue(ptId, 0); //can be smoothed
          }
        }
      }
  
    // All intersection points are generated, now create edges. Use a clipping
    // approach to create line segments. Later we'll connect them into polylines.
    //
    for (j=0; j<(dims[1]-1); j++) //loop over all x edges
      {
      for (i=0; i<(dims[0]-1); i++)
        {
        edgeCount = 0;
        p0 = i   + j*dims[0];
        p1 = p0 + 1;
        p2 = i+1 + (j+1)*dims[0];
        p3 = p2 - 1;
   
        if ( (ptId=this->EdgeTable->IsEdge(p0,p1)) != -1 )
          {
          id[2*edgeCount] = p0; id[2*edgeCount+1] = p1;
          pts[edgeCount++] = ptId;
          }
        if ( (ptId=this->EdgeTable->IsEdge(p1,p2)) != -1 )
          {
          id[2*edgeCount] = p1; id[2*edgeCount+1] = p2;
          pts[edgeCount++] = ptId;
          }
        if ( (ptId=this->EdgeTable->IsEdge(p2,p3)) != -1 )
          {
          id[2*edgeCount] = p2; id[2*edgeCount+1] = p3;
          pts[edgeCount++] = ptId;
          }
        if ( (ptId=this->EdgeTable->IsEdge(p3,p0)) != -1 )
          {
          id[2*edgeCount] = p3; id[2*edgeCount+1] = p0;
          pts[edgeCount++] = ptId;
          }
  
        if ( edgeCount == 4 )
          {
          x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
          x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
          ptId = points->InsertNextPoint(x);
          pointDescr->InsertValue(ptId, 0); //intersection points are fixed
  
          edgeConn->InsertNextCell(2);
          edgeConn->InsertCellPoint(ptId);
          edgeConn->InsertCellPoint(pts[0]);
          attrId = this->EdgeUseTable->InsertEdge(ptId,pts[0]);
          this->EdgeUses->InsertValue(2*attrId, this->Visited[id[0]]);
          this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[1]]);
          
          edgeConn->InsertNextCell(2);
          edgeConn->InsertCellPoint(ptId);
          edgeConn->InsertCellPoint(pts[1]);
          attrId = this->EdgeUseTable->InsertEdge(ptId,pts[1]);
          this->EdgeUses->InsertValue(2*attrId, this->Visited[id[2]]);
          this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[3]]);
  
          edgeConn->InsertNextCell(2);
          edgeConn->InsertCellPoint(ptId);
          edgeConn->InsertCellPoint(pts[2]);
          attrId = this->EdgeUseTable->InsertEdge(ptId,pts[2]);
          this->EdgeUses->InsertValue(2*attrId, this->Visited[id[4]]);
          this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[5]]);
  
          edgeConn->InsertNextCell(2);
          edgeConn->InsertCellPoint(ptId);
          edgeConn->InsertCellPoint(pts[3]);
          attrId = this->EdgeUseTable->InsertEdge(ptId,pts[3]);
          this->EdgeUses->InsertValue(2*attrId, this->Visited[id[6]]);
          this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[7]]);
          }
  
        else if ( edgeCount == 3 )
          {
          x[0] = origin[0] + i*spacing[0] + 0.5*spacing[0];
          x[1] = origin[1] + j*spacing[1] + 0.5*spacing[1];
          ptId = points->InsertNextPoint(x);
          pointDescr->InsertValue(ptId, 0); //intersection points are fixed
  
          edgeConn->InsertNextCell(2);
          edgeConn->InsertCellPoint(ptId);
          edgeConn->InsertCellPoint(pts[0]);
          attrId = this->EdgeUseTable->InsertEdge(ptId, pts[0]);
          this->EdgeUses->InsertValue(2*attrId, this->Visited[id[0]]);
          this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[1]]);
  
          edgeConn->InsertNextCell(2);
          edgeConn->InsertCellPoint(ptId);
          edgeConn->InsertCellPoint(pts[1]);
          attrId = this->EdgeUseTable->InsertEdge(ptId, pts[1]);
          this->EdgeUses->InsertValue(2*attrId, this->Visited[id[2]]);
          this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[3]]);
  
          edgeConn->InsertNextCell(2);
          edgeConn->InsertCellPoint(ptId);
          edgeConn->InsertCellPoint(pts[2]);
          attrId = this->EdgeUseTable->InsertEdge(ptId, pts[2]);
          this->EdgeUses->InsertValue(2*attrId, this->Visited[id[4]]);
          this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[5]]);
          }
  
        else if ( edgeCount == 2 )
          {
          edgeConn->InsertNextCell(2);
          edgeConn->InsertCellPoint(pts[0]);
          edgeConn->InsertCellPoint(pts[1]);
          attrId = this->EdgeUseTable->InsertEdge(pts[0],pts[1]);
          this->EdgeUses->InsertValue(2*attrId, this->Visited[id[0]]);
          this->EdgeUses->InsertValue(2*attrId+1, this->Visited[id[1]]);
          }
  
        else if ( edgeCount == 1 )
          {
          vtkErrorMacro(<<"Bad mojo");
          return 0;
          }
        }
      }
    
    // Cleanup
    this->EdgeUseTable->Delete();
    this->EdgeTable->Delete();
  
    return 0;
  }
  
  void vtkImageToPolyDataFilter::BuildPolygons(vtkUnsignedCharArray *vtkNotUsed(pointDescr), 
                                               vtkPolyData *edges, int numPolys,
                                               vtkUnsignedCharArray *polyColors)
  {
    vtkPoints *points = edges->GetPoints();
    vtkIdType numPts = points->GetNumberOfPoints(), ptId;
    int i, j, k, *polyId, *polyId2, edgeId;
    vtkIdType *cells, *pts, *cells2, npts, cellId;
    int numPolyPts, p1, p2;
    unsigned short ncells, ncells2;
    unsigned char *polyVisited, *ptr;
    vtkCellArray *newPolys;
    
    // Make sure we can topological info
    edges->BuildLinks();
  
    // Mark all polygons as unvisited
    polyVisited = new unsigned char [numPolys];
    for (i=0; i<numPolys; i++)
      {
      polyVisited[i] = 0;
      }
  
    // Create connectivity array for polygon definition
    newPolys = vtkCellArray::New();
    newPolys->Allocate(newPolys->EstimateSize(numPolys,25));
  
    // Loop over all edge points tracking around each polygon
    for (ptId=0; ptId<numPts; ptId++)
      {
      edges->GetPointCells(ptId, ncells, cells);
      if (ncells < 2)
        {
        vtkErrorMacro(<<"Bad mojo");
        return;
        }
      //for each edge, walk around polygon (if not visited before)
      for (i=0; i<ncells; i++)
        {
        edgeId = cells[i];
        polyId = this->EdgeUses->GetPointer(2*edgeId);
        for (j=0; j<2; j++)
          {
          if ( polyId[j] != -1 && !polyVisited[polyId[j]] )
            {//build loop
            polyVisited[polyId[j]] = 1;
            numPolyPts = 1;
            cellId = newPolys->InsertNextCell(0); //will update count later
            newPolys->InsertCellPoint(ptId);
  
            // Update polygonal color
            ptr = this->PolyColors->GetPointer(3*polyId[j]);
            polyColors->SetValue(3*cellId, ptr[0]); //assign poly color
            polyColors->SetValue(3*cellId+1, ptr[1]);
            polyColors->SetValue(3*cellId+2, ptr[2]);
  
            p1 = ptId; p2 = -1;
            while ( 1 )
              {
              edges->GetCellPoints(edgeId, npts, pts);
              p2 = (pts[0] != p1 ? pts[0] : pts[1]);
              if (p2 == ptId) 
                {
                break;
                }
              
              newPolys->InsertCellPoint(p2);
              numPolyPts++;
              edges->GetPointCells(p2, ncells2, cells2);
              if (ncells < 2)
                {
                vtkErrorMacro(<<"Bad mojo");
                return;
                }
              for (k=0; k<ncells2; k++)
                {
                polyId2 = this->EdgeUses->GetPointer(2*cells2[k]);
                if ( cells2[k] != edgeId && 
                     (polyId2[0] == polyId[j] || polyId2[1] == polyId[j]) )
                  {
                  p1 = p2;
                  edgeId = cells2[k];
                  break;
                  }
                }
              }//while not completed loop
            newPolys->UpdateCellCount(numPolyPts);
            
            }//if polygon not yet visited
          }//for each use of edge by polygon (at most 2 polygons)
        }//for each edge connected to this point
      }//for all points in edge list
  
    edges->SetPolys(newPolys);
    newPolys->Delete();
    this->EdgeUses->Delete();
    delete [] polyVisited;
  }
  
  void vtkImageToPolyDataFilter::SmoothEdges(vtkUnsignedCharArray *pointDescr, 
                                             vtkPolyData *edges)
  {
    
    vtkPoints *points=edges->GetPoints();
    vtkIdType numPts=points->GetNumberOfPoints(), ptId;
    int i, iterNum;
    int connId;
    double x[3], xconn[3], xave[3], factor;
    unsigned short int ncells;
    vtkIdType *cells, *pts, npts;
  
  
    // For each smoothing operation, loop over points. Points that can be 
    // smoothed are moved in the direction of the average of their neighbor
    // points.
    for ( iterNum=0; iterNum<this->NumberOfSmoothingIterations; iterNum++ )
      {
      if ( (iterNum % 2) ) //alternate smoothing direction
        {
        factor = -0.331;
        }
      else
        {
        factor = 0.330;
        }
      
      for (ptId=0; ptId<numPts; ptId++)
        {
        if ( pointDescr->GetValue(ptId) == 0 ) //can smooth
          {
          points->GetPoint(ptId, x);
          edges->GetPointCells(ptId, ncells, cells);
          xave[0] = xave[1] = xave[2] = 0.0;
          for (i=0; i<ncells; i++)
            {
            edges->GetCellPoints(cells[i], npts, pts);
            if (pts[0] != ptId)
              {
              connId = pts[0];
              }
            else if (npts > 1)
              {
              connId = pts[1];
              }
            else
              {
              vtkErrorMacro("Bad cell in smoothing operation");
              connId = pts[0];
              }
            points->GetPoint(connId, xconn);
            xave[0] += xconn[0]; xave[1] += xconn[1]; xave[2] += xconn[2];
            }
          if ( ncells > 0 )
            {
            xave[0] /= ncells; xave[1] /= ncells; xave[2] /= ncells;
            x[0] = x[0] + factor * (xave[0] - x[0]);
            x[1] = x[1] + factor * (xave[1] - x[1]);
            x[2] = x[2] + factor * (xave[2] - x[2]);
            points->SetPoint(ptId, x);
            }
          
          }//if smoothable point
        }//for all points
      }//for all smoothing operations
  }
  
  
  // Remove points that are nearly co-linear to reduce the total point count
  //
  void vtkImageToPolyDataFilter::DecimateEdges(vtkPolyData *edges, 
                                               vtkUnsignedCharArray *pointDescr,
                                               double tol2)
  {
    vtkPoints *points=edges->GetPoints();
    vtkIdType numPts=points->GetNumberOfPoints(), ptId, prevId, nextId;
    vtkIdType npts;
    double x[3], xPrev[3], xNext[3];
    unsigned short int ncells;
    vtkIdType *cells, *pts;
  
    // Loop over all points, finding those that are connected to just two
    // edges. If the point is colinear to the previous and next edge point,
    // then mark it as deleted.
    for (ptId=0; ptId<numPts; ptId++)
      {
      if ( pointDescr->GetValue(ptId) == 0 )
        {
        points->GetPoint(ptId, x);
        edges->GetPointCells(ptId, ncells, cells);
        if ( ncells == 2 )
          {
          edges->GetCellPoints(cells[0], npts, pts);
          prevId = (pts[0] != ptId ? pts[0] : pts[1]);
          points->GetPoint(prevId, xPrev);
  
          edges->GetCellPoints(cells[1], npts, pts);
          nextId = (pts[0] != ptId ? pts[0] : pts[1]);
          points->GetPoint(nextId, xNext);
  
          if ( vtkLine::DistanceToLine(x, xPrev, xNext) <= tol2 )
            {
            pointDescr->SetValue(ptId, 2); //mark deleted
            }
          }
        } //if manifold
      } //for all points
  }
  