/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkImageToPolyDataFilter.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


Copyright (c) 1993-2001 Ken Martin, Will Schroeder, Bill Lorensen 
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.

 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

 * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
   of any contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

 * Modified source versions must be plainly marked as such, and must not be
   misrepresented as being the original software.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

=========================================================================*/
#include "vtkImageToPolyDataFilter.h"
#include "vtkUnsignedCharArray.h"
#include "vtkLine.h"
#include "vtkEdgeTable.h"
#include "vtkAppendPolyData.h"

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
//

void vtkImageToPolyDataFilter::Execute()
{
  vtkStructuredPoints *input=this->GetInput();
  vtkPolyData *output=this->GetOutput();
  vtkPolyData *tmpOutput=vtkPolyData::New();
  vtkPolyData *tmpInput=vtkPolyData::New();
  vtkAppendPolyData *append = vtkAppendPolyData::New();
  vtkPolyData *appendOutput;
  vtkScalars *inScalars = input->GetPointData()->GetScalars();
  int numPixels=input->GetNumberOfPoints();
  int dims[3], numComp=inScalars->GetData()->GetNumberOfComponents();
  float origin[3], spacing[3];
  vtkUnsignedCharArray *pixels;
  int type=inScalars->GetData()->GetDataType();
  int numPieces[2], extent[4];
  int i, j, newDims[3], totalPieces, pieceNum, abortExecute=0;
  float newOrigin[3];

  // Check input and initialize
  vtkDebugMacro(<<"Vectorizing image...");
  
  if ( inScalars == NULL || numPixels < 1 )
    {
    vtkErrorMacro(<<"Not enough input to create output");
    return;
    }
  
  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 ((float)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();
}

void vtkImageToPolyDataFilter::PixelizeImage(vtkUnsignedCharArray *pixels, 
                                             int dims[3], float origin[3], 
                                             float spacing[3], 
                                             vtkPolyData *output)
{
  int numPts, numCells, pts[4], i, j, id;
  vtkPoints *newPts;
  vtkCellArray *newPolys;
  float x[3];
  vtkUnsignedCharArray *polyColors;
  unsigned char *ptr, *colors=pixels->GetPointer(0);
  vtkScalars *polyScalars;
  
  // 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);
  polyScalars = vtkScalars::New();
  polyScalars->SetData(polyColors);

  // 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(polyScalars);
  polyScalars->Delete();
  polyColors->Delete();
}

void vtkImageToPolyDataFilter::RunLengthImage(vtkUnsignedCharArray *pixels, 
                                             int dims[3], float origin[3], 
                                             float spacing[3], 
                                             vtkPolyData *output)
{
  int pts[4], i, j, id;
  vtkPoints *newPts;
  vtkCellArray *newPolys;
  float 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();

  vtkScalars *polyScalars = vtkScalars::New();
  polyScalars->SetData(polyColors);
  output->GetCellData()->SetScalars(polyScalars);
  polyScalars->Delete();
  polyColors->Delete();
}


void vtkImageToPolyDataFilter::PolygonalizeImage(vtkUnsignedCharArray *pixels, 
                               int dims[3], float origin[3], float 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(
                          vtkScalars *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 = ((vtkUnsignedCharArray *)inScalars->GetData())->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;
      }

    float 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->GetScalar(idx);
        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);
        }
      }
    }
}

void vtkImageToPolyDataFilter::PrintSelf(ostream& os, vtkIndent indent)
{
  vtkStructuredPointsToPolyDataFilter::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(numPixels/4.0, numPixels/4.0);
  wave2 = vtkIdList::New();
  wave2->Allocate(numPixels/4.0, 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, npts, *pts, numPts;
  
  // 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();
  
  vtkScalars *polyScalars = vtkScalars::New();
  polyScalars->SetData(polyColors);
  output->GetCellData()->SetScalars(polyScalars);
  polyScalars->Delete();
}

// Uses clipping approach to build the polygon edges
int vtkImageToPolyDataFilter::BuildEdges(vtkUnsignedCharArray *vtkNotUsed(pixels), 
                                         int dims[3], float origin[3],
                                         float spacing[3],
                                         vtkUnsignedCharArray *pointDescr, 
                                         vtkPolyData *edges)
{
  float x[3];
  int i, j, ptId, edgeCount, p0, p1, p2, p3;
  vtkCellArray *edgeConn = edges->GetLines();
  int pts[4], startId, attrId, id[8];
  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();
  int numPts = points->GetNumberOfPoints();
  int i, j, k, ptId, cellId, *cells, *cells2, *polyId, *polyId2, edgeId;
  int numPolyPts, *pts, npts, 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(npts); //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();
  int numPts=points->GetNumberOfPoints();
  int i, ptId, iterNum;
  int npts, *pts, connId;
  float x[3], *xconn, xave[3], factor;
  unsigned short int ncells;
  int *cells;

  // 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);
          connId = (pts[0] != ptId ? pts[0] : pts[1]);
          xconn = points->GetPoint(connId);
          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,
                                             float tol2)
{
  vtkPoints *points=edges->GetPoints();
  int numPts=points->GetNumberOfPoints();
  int ptId, prevId, nextId;
  int npts, *pts;
  float *x, *xPrev, *xNext;
  unsigned short int ncells;
  int *cells;

  // 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 )
      {
      x = points->GetPoint(ptId);
      edges->GetPointCells(ptId, ncells, cells);
      if ( ncells == 2 )
        {
        edges->GetCellPoints(cells[0], npts, pts);
        prevId = (pts[0] != ptId ? pts[0] : pts[1]);
        xPrev = points->GetPoint(prevId);

        edges->GetCellPoints(cells[1], npts, pts);
        nextId = (pts[0] != ptId ? pts[0] : pts[1]);
        xNext = points->GetPoint(nextId);

        if ( vtkLine::DistanceToLine(x, xPrev, xNext) <= tol2 )
          {
          pointDescr->SetValue(ptId, 2); //mark deleted
          }
        }
      } //if manifold
    } //for all points
}


