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

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


Copyright (c) 1993-2000 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 REGENTS 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 "vtkVoxelContoursToSurfaceFilter.h"
#include "vtkStructuredPoints.h"
#include "vtkContourFilter.h"
#include "vtkAppendPolyData.h"
#include "vtkObjectFactory.h"



//------------------------------------------------------------------------------
vtkVoxelContoursToSurfaceFilter* vtkVoxelContoursToSurfaceFilter::New()
{
  // First try to create the object from the vtkObjectFactory
  vtkObject* ret = vtkObjectFactory::CreateInstance("vtkVoxelContoursToSurfaceFilter");
  if(ret)
    {
    return (vtkVoxelContoursToSurfaceFilter*)ret;
    }
  // If the factory was unable to create the object, then create it here.
  return new vtkVoxelContoursToSurfaceFilter;
}




vtkVoxelContoursToSurfaceFilter::vtkVoxelContoursToSurfaceFilter()
{
  this->MemoryLimitInBytes = 10000000;
  this->Spacing[0]         = 1.0;
  this->Spacing[1]         = 1.0;
  this->Spacing[2]         = 1.0;
  this->LineList           = new float[4*1000];
  this->LineListLength     = 0;
  this->LineListSize       = 1000;
  this->SortedXList        = NULL;
  this->SortedYList        = NULL;
  this->WorkingList        = NULL;
  this->IntersectionList   = NULL;
  this->SortedListSize     = 0;
}

vtkVoxelContoursToSurfaceFilter::~vtkVoxelContoursToSurfaceFilter()
{
  if ( this->LineList )
    {
    delete [] this->LineList;
    }
  if ( this->SortedXList )
    {
    delete [] this->SortedXList;
    }
  if ( this->SortedYList )
    {
    delete [] this->SortedYList;
    }
  if ( this->WorkingList )
    {
    delete [] this->WorkingList;
    }
  if ( this->IntersectionList )
    {
    delete [] this->IntersectionList;
    }
} 

void vtkVoxelContoursToSurfaceFilter::AddLineToLineList( float x1, float y1,
						   float x2, float y2 )
{
  // Do we need to increase the size of our list?
  if ( this->LineListLength >= this->LineListSize )
    {
    // Double the space we had before
    float *newList = new float[this->LineListSize*4*2];
    memcpy( newList, this->LineList, 4*this->LineListSize*sizeof(float) );
    delete [] this->LineList;
    this->LineList = newList;
    this->LineListSize *= 2;
    }

  // Now we are sure we have space - add the line
  this->LineList[4*this->LineListLength + 0] = x1;
  this->LineList[4*this->LineListLength + 1] = y1;
  this->LineList[4*this->LineListLength + 2] = x2;
  this->LineList[4*this->LineListLength + 3] = y2;
  this->LineListLength++;
}

void vtkVoxelContoursToSurfaceFilter::SortLineList()
{
  int     i, j;
  float   tmp[4];
  float   tmpval;


  // Make sure we have enough space in our sorted list
  if ( this->SortedListSize < this->LineListLength )
    {
    if ( this->SortedXList )
      {
      delete [] this->SortedXList;
      }
    if ( this->SortedYList )
      {
      delete [] this->SortedYList;
      }
    if ( this->WorkingList )
      {
      delete [] this->WorkingList;
      }
    if ( this->IntersectionList )
      {
      delete [] this->IntersectionList;
      }

    this->SortedXList = new float[4*this->LineListLength];
    this->SortedYList = new float[4*this->LineListLength];
    this->SortedListSize = this->LineListLength;

    // Create the space we'll need for our working list of indices
    // The is the list of lines that we are currently considering
    // for intersections. Lines move in then out of the list as
    // we pass the first then the second endpoint. This will be
    // used during the CastXLines and CastYLines methods, and is 
    // the same size as the number of lines.
    this->WorkingList = new int[this->LineListLength];    

    // Create the space we'll need for the intersection list
    // There can't be more intersections than there are possible
    // lines. Although it is highly doubtful we'll actually use
    // all this space, it isn't much and it makes the code simpler
    // not to have to worry about exceeding the bounds. This will be
    // used during the CastXLines and CastYLines methods, and is 
    // the same size as the number of lines.
    this->IntersectionList = new float[this->LineListLength];    
    }

  // Copy the lines into the lists
  memcpy( this->SortedXList, this->LineList, 
	  4*this->LineListLength*sizeof(float) );
  memcpy( this->SortedYList, this->LineList, 
	  4*this->LineListLength*sizeof(float) );

  // Now sort on x and y
  // Use a simple bubble sort - will improve if necessary
  for ( i = 0; i < this->LineListLength; i++ )
    {
    // swap x entry if necessary to keep min x the first endpoint
    if ( this->SortedXList[4*i + 0] > this->SortedXList[4*i + 2] )
      {
      tmpval = this->SortedXList[4*i];
      this->SortedXList[4*i] = this->SortedXList[4*i + 2];
      this->SortedXList[4*i + 2] = tmpval;
      tmpval = this->SortedXList[4*i + 1];
      this->SortedXList[4*i + 1] = this->SortedXList[4*i + 3];
      this->SortedXList[4*i + 3] = tmpval;
      }

    // swap y entry if necessary to keep min y the first endpoint
    if ( this->SortedYList[4*i + 1] > this->SortedYList[4*i + 3] )
      {
      tmpval = this->SortedYList[4*i];
      this->SortedYList[4*i] = this->SortedYList[4*i + 2];
      this->SortedYList[4*i + 2] = tmpval;
      tmpval = this->SortedYList[4*i + 1];
      this->SortedYList[4*i + 1] = this->SortedYList[4*i + 3];
      this->SortedYList[4*i + 3] = tmpval;
      }

    // Sort x list
    for ( j = i; j > 0; j-- )
      {
      if ( this->SortedXList[j*4] < this->SortedXList[(j-1)*4] )
	{
	memcpy( tmp, this->SortedXList + j*4, 4*sizeof(float) );
	memcpy( this->SortedXList + j*4, 
		this->SortedXList + (j-1)*4, 4*sizeof(float) );
	memcpy( this->SortedXList + (j-1)*4, tmp, 4*sizeof(float) );
	}
      else
	{
	break;
	}
      }

    // Sort y list
    for ( j = i; j > 0; j-- )
      {
      if ( this->SortedYList[j*4+1] < this->SortedYList[(j-1)*4+1] )
	{
	memcpy( tmp, this->SortedYList + j*4, 4*sizeof(float) );
	memcpy( this->SortedYList + j*4, this->SortedYList + (j-1)*4, 
		4*sizeof(float) );
	memcpy( this->SortedYList + (j-1)*4, tmp, 4*sizeof(float) );
	}
      else
	{
	break;
	}
      }
    }
}


void vtkVoxelContoursToSurfaceFilter::CastLines( float *slicePtr, 
						 float gridOrigin[3], 
						 int gridSize[3],
						 int type )
{
  float   axis1, axis2;
  float   d1, d2;
  int     index;
  int     i, j;
  float   tmp;
  float   *line;
  float   *currSlicePtr;
  int     currSlice;
  int     currentIntersection;
  float   sign;
  float   *sortedList;
  float   low1, low2, high1, high2;
  int     increment1, increment2;
  int     offset1, offset2, offset3, offset4;

  // this is the x direction
  if ( type == 0 )
    {
    low1 = gridOrigin[0];
    high1 = gridOrigin[0] + (float)gridSize[0];
    low2 = gridOrigin[1];
    high2 = gridOrigin[1] + (float)gridSize[1];
    increment1 = gridSize[0];
    increment2 = 1;
    sortedList = this->SortedXList;
    offset1 = 0;
    offset2 = 2;
    offset3 = 1;
    offset4 = 3;
    }
  // This is the y direction
  else
    {
    low1 = gridOrigin[1];
    high1 = gridOrigin[1] + (float)gridSize[1];
    low2 = gridOrigin[0];
    high2 = gridOrigin[0] + (float)gridSize[0];
    increment1 = 1;
    increment2 = gridSize[0];
    sortedList = this->SortedYList;
    offset1 = 1;
    offset2 = 3;
    offset3 = 0;
    offset4 = 2;
    }

  // Initialize the working list to nothing. We will start
  // looking at index = 0 for the next line to add to the 
  // working list
  this->WorkingListLength = 0;
  index = 0;
  
  // Loop through the x or y lines
  for ( axis1 = low1, currSlice = 0; axis1 < high1; axis1 += 1.0, currSlice++ )
    {
    // Initialize the intersection list to nothing
    this->IntersectionListLength = 0;

    // Add lines to the working list if necessary
    for ( ; index < this->LineListLength; index++ )
      {
      if ( sortedList[4*index + offset1] < axis1 )
	{
	this->WorkingList[this->WorkingListLength] = index;
	this->WorkingListLength++;
	}
      else
	{
	break;
	}
      }

    // Do the intersections, removing lines from the
    // working list if necessary
    for ( i = 0; i < this->WorkingListLength; i++ )
      {
      line = sortedList + 4*this->WorkingList[i];

      // Yes, it intersects, add it to the intersection list
      if ( line[offset1] < axis1 && line[offset2] > axis1 )
	{
	// Compute the intersection distance
	// For x lines this is y = y1 + (y2 - y1)*((x - x1)/(x2 - x1))
	// For y lines this is x = x1 + (x2 - x1)*((y - y1)/(y2 - y1))
	this->IntersectionList[this->IntersectionListLength] = 
	  line[offset3] + (line[offset4] - line[offset3]) *
	  ((axis1 - line[offset1]) / (line[offset2] - line[offset1] ));

	// Make sure this distance is sorted
	for ( j = this->IntersectionListLength; j > 0; j-- )
	  {
	  if ( this->IntersectionList[j] < this->IntersectionList[j-1] )
	    {
	    tmp = this->IntersectionList[j];
	    this->IntersectionList[j] = this->IntersectionList[j-1];
	    this->IntersectionList[j-1] = tmp;
	    }
	  else
	    {
	    break;
	    }
	  }
	this->IntersectionListLength++;
	}
      // No, it doesn't intersect, remove it from the working list
      else
	{
	for ( j = i; j < (this->WorkingListLength-1); j++ )
	  {
	  this->WorkingList[j] = this->WorkingList[j+1];
	  }
	this->WorkingListLength--;
	i--;
	}
      }

    // Now we have all the intersections for the x or y line, in sorted
    // order. Use them to fill in distances (as long as there are
    // any)
    if ( this->IntersectionListLength )
      {
      currSlicePtr = slicePtr + currSlice*increment2;
      currentIntersection = 0;
      // We are starting outside which has a negative distance
      sign = -1.0;
      for ( axis2 = low2; axis2 < high2; axis2 += 1.0 )
	{
	while( currentIntersection < this->IntersectionListLength &&
	       this->IntersectionList[currentIntersection] < axis2 )

	  {
	  currentIntersection++;

	  // Each time we cross a line we are moving across an 
	  // inside/outside boundary
	  sign *= -1.0;
	  }
	// We are now positioned at an x or y value between currentIntersection
	// and currentIntersection - 1 (except at boundaries where we are
	// before intersection 0 or after the last intersection)

	if ( currentIntersection == 0 )
	  {
	  d1 = axis2 - this->IntersectionList[currentIntersection];
	  *currSlicePtr = (*currSlicePtr > d1 )?(*currSlicePtr):(d1);
	  }
	else if ( currentIntersection == this->IntersectionListLength )
	  {
	  d1 = this->IntersectionList[currentIntersection-1] - axis2;
	  *currSlicePtr = (*currSlicePtr > d1 )?(*currSlicePtr):(d1);
	  }
	else
	  {
	  d1 = axis2 - this->IntersectionList[currentIntersection-1];
	  d2 = this->IntersectionList[currentIntersection] - axis2;
	  d1 = ( d1 < d2 )?(d1):(d2);
	  if ( type == 0 )
	    {
	    *currSlicePtr = sign*d1;
	    }
	  else 
	    {
	    *currSlicePtr = 
	      (sign*(*currSlicePtr) < d1 )?(*currSlicePtr):(sign*d1);
	    }
	  }

	currSlicePtr += increment1;
	}
      }
    }
}

void vtkVoxelContoursToSurfaceFilter::PushDistances( float *volumePtr,
						     int gridSize[3],
						     int chunkSize )
{
  int    i, j, k;
  float  *vptr;

  // Push distances along x (both ways) and y (both ways) on each slice
  for ( k = 0; k < chunkSize; k++ )
    {
    // Do the x rows
    for ( j = 0; j < gridSize[1]; j++ )
      {
      vptr = volumePtr + k*gridSize[0]*gridSize[1] + j*gridSize[0];
      vptr++;

      // first one way
      for ( i = 1; i < gridSize[0]; i++ )
	{
	if ( *vptr > 0 &&  *(vptr-1) + 1 < *(vptr) )
	  {
	  *vptr = *(vptr-1) + 1;
	  }
	else if ( *vptr < 0 && *(vptr-1) - 1 > *(vptr) )
	  {
	  *vptr = *(vptr-1) - 1;
	  }
	vptr++;
	}

      vptr -= 2;
      i    -= 2;

      // then the other
      for ( ; i >= 0; i-- )
	{
	if ( *vptr > 0 &&  *(vptr+1) + 1 < *(vptr) )
	  {
	  *vptr = *(vptr+1) + 1;
	  }
	else if ( *vptr < 0 && *(vptr+1) - 1 > *(vptr) )
	  {
	  *vptr = *(vptr+1) - 1;
	  }
	}

      }


    // Do the y columns
    for ( i = 0; i < gridSize[0]; i++ )
      {
      vptr = volumePtr + k*gridSize[0]*gridSize[1] + i;

      vptr+=gridSize[0];

      // first one way
      for ( j = 1; j < gridSize[1]; j++ )
	{
	if ( *vptr > 0 &&  *(vptr-gridSize[0]) + 1 < *(vptr) )
	  {
	  *vptr = *(vptr-gridSize[0]) + 1;
	  }
	else if ( *vptr < 0 && *(vptr-gridSize[0]) - 1 > *(vptr) )
	  {
	  *vptr = *(vptr-gridSize[0]) - 1;
	  }
	vptr += gridSize[0];
	}

      vptr -= 2*gridSize[0];
      j    -= 2;

      // then the other
      for ( ; j >= 0; j-- )
	{
	if ( *vptr > 0 &&  *(vptr+gridSize[0]) + 1 < *(vptr) )
	  {
	  *vptr = *(vptr+gridSize[0]) + 1;
	  }
	else if ( *vptr < 0 && *(vptr+gridSize[0]) - 1 > *(vptr) )
	  {
	  *vptr = *(vptr+gridSize[0]) - 1;
	  }
	}

      }
    }
} 

// Append data sets into single unstructured grid
void vtkVoxelContoursToSurfaceFilter::Execute()
{
  vtkPolyData          *output = this->GetOutput();
  vtkPolyData          *input  = this->GetInput();
  vtkCellArray         *inputPolys = input->GetPolys();
  int                  gridSize[3];
  float                gridOrigin[3];
  float                contourBounds[6];
  int                  chunkSize;
  int                  currentSlice, lastSlice, currentIndex;
  int                  i, j;
  int                  numberOfInputCells;
  int                  currentInputCellIndex;
  int                  npts, *pts;
  float                point1[3], point2[3];
  float                currentZ;
  vtkStructuredPoints  *volume;
  float                *volumePtr, *slicePtr;
  vtkContourFilter     *contourFilter;
  vtkPolyData          *contourOutput;
  vtkAppendPolyData    *appendFilter;

  vtkDebugMacro(<<"Creating surfaces from contours");

  // Get the bounds of the input contours
  input->GetBounds( contourBounds );

  // From the bounds, compute the grid size, and origin
  
  // The origin of the grid should be (-0.5, -0.5, 0.0) away from the 
  // lower bounds of the contours. This is because we want the grid
  // to lie halfway between integer endpoint locations of the line
  // segments on each plane. Also, we want an extra plane on each end
  // for capping
  gridOrigin[0] = contourBounds[0] - 0.5;
  gridOrigin[1] = contourBounds[2] - 0.5;
  gridOrigin[2] = contourBounds[4] - 1.0;

  // The difference between the bounds, plus one to account a 
  // sample on the first and last location, plus one to account 
  // for the larger grid size ( the 0.5 unit border ) On Z, we 
  // want to sample exactly on the contours so we don't need to 
  // add the extra 1, but we have added two extra planes so we
  // need another 2.
  gridSize[0] = (int) (contourBounds[1] - contourBounds[0] + 2);
  gridSize[1] = (int) (contourBounds[3] - contourBounds[2] + 2);
  gridSize[2] = (int) (contourBounds[5] - contourBounds[4] + 3);

  // How many slices in a chunk? This will later be decremented 
  // by one to account for the fact that the last slice in the 
  // previous chuck is copied to the first slice in the next chunk. 
  // Stay within memory limit. There are 4 bytes per float.
  chunkSize = this->MemoryLimitInBytes / ( gridSize[0] * gridSize[1] * 4 );
  if ( chunkSize > gridSize[2] )
    {
    chunkSize = gridSize[2];
    }

  currentSlice          = 0;
  currentZ              = contourBounds[4] - 1.0;
  currentIndex          = 0;
  lastSlice             = gridSize[2] - 1;
  numberOfInputCells    = inputPolys->GetNumberOfCells();
  currentInputCellIndex = 0;

  volume = vtkStructuredPoints::New();
  volume->SetDimensions( gridSize[0], gridSize[1], chunkSize );
  volume->SetSpacing( this->Spacing );
  volume->SetScalarType( VTK_FLOAT );
  volume->AllocateScalars();
  volumePtr = 
    (float *)(volume->GetPointData()->GetScalars()->GetVoidPointer(0));


  contourFilter = vtkContourFilter::New();
  contourFilter->SetInput( volume );
  contourFilter->SetNumberOfContours(1);
  contourFilter->SetValue( 0, 0.0 );

  appendFilter = vtkAppendPolyData::New();

  inputPolys->InitTraversal();
  inputPolys->GetNextCell( npts, pts );
 
  while ( currentSlice <= lastSlice )
    {
    // Make sure the origin of the volume is in the right
    // place so that the appended polydata all matches up
    // nicely.
    volume->SetOrigin( gridOrigin[0], gridOrigin[1],
		       gridOrigin[2] + 
		       this->Spacing[2] * (currentSlice - (currentSlice!=0)) );

    for ( i = currentIndex; i < chunkSize; i++ )
      {
      slicePtr = volumePtr + i * gridSize[0] * gridSize[1];

      // Clear out the slice memory - set it all to a large negative
      // value indicating no surfaces are nearby, and we assume we
      // are outside of any surface
      for ( j = 0; j < gridSize[0] * gridSize[1]; j++ )
	{
	*(slicePtr+j) = -9.99e10;
	}

      // If we are past the end, don't do anything
      if ( currentSlice > lastSlice )
	{
	continue;
	}

      this->LineListLength = 0;

      // Read in the lines for the contours on this slice
      while ( currentInputCellIndex < numberOfInputCells )
	{	  
	// Check if we are still on the right z slice
	input->GetPoint( pts[0], point1 );
	if ( point1[2] != currentZ )
	  {
	  break;
	  }

	// This contour is on the right z slice - add the lines
	// to our list
	for ( j = 0; j < npts; j++ )
	  {
	  input->GetPoint( pts[j], point1 );
	  input->GetPoint( pts[(j+1)%npts], point2 );
	  this->AddLineToLineList( point1[0], point1[1],
				   point2[0], point2[1] );
	  }

	inputPolys->GetNextCell( npts, pts );
	currentInputCellIndex++;
	}

      // Sort the contours in x and y
      this->SortLineList();

      // Cast lines in x and y filling in distance
      this->CastLines( slicePtr, gridOrigin, gridSize, 0 );
      this->CastLines( slicePtr, gridOrigin, gridSize, 1 );

      // Move on to the next slice
      currentSlice++;
      currentIndex++;
      currentZ += 1.0;
      }

    this->PushDistances( volumePtr, gridSize, chunkSize );

    // Update the contour filter and grab the output
    // Make a new output for it, then grab the output and
    // add it to the append filter, then delete the output
    // which is ok since it was registered by the appendFilter
    contourOutput = vtkPolyData::New();
    contourFilter->SetOutput( contourOutput );
    contourFilter->Update();
    appendFilter->AddInput( contourFilter->GetOutput() );
    contourFilter->SetOutput( NULL );
    contourOutput->Delete();


    if ( currentSlice <= lastSlice )
      {
      // Copy last slice to first slice
      memcpy( volumePtr, volumePtr + (chunkSize-1)*gridSize[0]*gridSize[1],
	      sizeof(float) * gridSize[0] * gridSize[1] );

      // reset currentIndex to 1
      currentIndex = 1;
      }
    }

  appendFilter->Update();

  // Grab the appended data as the output to this filter
  output->SetPoints( appendFilter->GetOutput()->GetPoints() );
  output->SetVerts(  appendFilter->GetOutput()->GetVerts() );
  output->SetLines(  appendFilter->GetOutput()->GetLines() );
  output->SetPolys(  appendFilter->GetOutput()->GetPolys() );
  output->SetStrips( appendFilter->GetOutput()->GetStrips() );
  output->GetPointData()->PassData(appendFilter->GetOutput()->GetPointData());

  contourFilter->Delete();
  appendFilter->Delete();
  volume->Delete();
}


void vtkVoxelContoursToSurfaceFilter::PrintSelf(ostream& os, vtkIndent indent)
{
  vtkPolyDataToPolyDataFilter::PrintSelf(os,indent);

  os << indent << "Memory Limit (in bytes): " << 
    this->MemoryLimitInBytes << endl;

  os << indent << "Spacing: " << this->Spacing[0] << " " <<
    this->Spacing[1] << " " << this->Spacing[2] << endl;
}
