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

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkElevationFilter.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 "vtkBoundaryRegionsFilter.h"
#include "vtkCellData.h"
#include "vtkDataSet.h"
#include "vtkIntArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkSmartPointer.h"
#include "vtkObjectFactory.h"
#include "vtkPolyDataNormals.h"


vtkCxxRevisionMacro(vtkBoundaryRegionsFilter, "$Revision: 0.0 $");
vtkStandardNewMacro(vtkBoundaryRegionsFilter);

//----------------------------------------------------------------------------
vtkBoundaryRegionsFilter::vtkBoundaryRegionsFilter()
{
  this->normals = vtkPolyDataNormals::New();
}

//----------------------------------------------------------------------------
vtkBoundaryRegionsFilter::~vtkBoundaryRegionsFilter()
{
  if (this->normals)
	{
	  this->normals->Delete();
	}
}

//----------------------------------------------------------------------------
void vtkBoundaryRegionsFilter::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
}

//----------------------------------------------------------------------------
int vtkBoundaryRegionsFilter::RequestData(vtkInformation*,
                                    vtkInformationVector** inputVector,
                                    vtkInformationVector* outputVector)
{

  // Definition of some variables used by the algorithm
  double  cellBounds[6], xb, yb, zb;
  double *cellNormals, *bounds;
  double  one = 0.99;
  int     reg;

  // Get the input and output data objects.
  vtkPolyData* input = vtkPolyData::GetData(inputVector[0]);
  vtkPolyData* output = vtkPolyData::GetData(outputVector);

  // Check the number of cells of the input
  vtkIdType numCells = input->GetNumberOfCells();
  if(numCells < 1)
    {
    vtkDebugMacro("No input!");
    return 1;
    }

  // Allocate space for the scalar per cell data
  vtkSmartPointer<vtkIntArray> newScalars = vtkSmartPointer<vtkIntArray>::New();
  newScalars->SetNumberOfTuples(numCells);

  // Computing polygon normals to be used during the filtering process
  vtkPolyData* inputCopy = input->NewInstance();
  inputCopy->ShallowCopy(input);
  this->normals->SetInput(inputCopy);
  this->normals->ComputeCellNormalsOn();
  this->normals->ComputePointNormalsOff();
  inputCopy->Delete();
  this->normals->Update();

  // Getting the geometry bounds
  bounds = this->normals->GetOutput()->GetBounds();

  // Support progress and abort.
  vtkIdType tenth = (numCells >= 10? numCells/10 : 1);
  double numPtsInv = 1.0/numCells;
  int abort = 0;

  // Notify that the filter is doing something
  vtkDebugMacro("Generating elevation scalars!");

  for(vtkIdType i=0; i < numCells && !abort; ++i)
    {

    // Periodically update progress and check for an abort request.
    if(i % tenth == 0)
      {
      this->UpdateProgress((i+1)*numPtsInv);
      abort = this->GetAbortExecute();
      }

      cellNormals = this->normals->GetOutput()->GetCellData()->GetNormals()->GetTuple(i);

	  this->normals->GetOutput()->GetCellBounds(i,cellBounds);

      // Computing the baricenter of the bounding box of the element
	  xb = cellBounds[0] + (cellBounds[1]-cellBounds[0])*0.5;
      yb = cellBounds[2] + (cellBounds[3]-cellBounds[2])*0.5;
	  zb = cellBounds[4] + (cellBounds[5]-cellBounds[4])*0.5;

	  // Filtering process according to regions
	  if      ( cellNormals[0] >= +one && xb >= bounds[1] ) // Front (region 0)
		  reg = 0;
      else if ( cellNormals[0] <= -one && xb <= bounds[0] ) // Back (region 1)
		  reg = 1;
      else if ( cellNormals[1] >= +one && yb >= bounds[3] ) // Right (region 2)
		  reg = 2;
      else if ( cellNormals[1] <= -one && yb <= bounds[2] ) // Left (region 3)
		  reg = 3;
      else if ( cellNormals[2] >= +one && zb >= bounds[5] ) // Top (region 4)
		  reg = 4;
      else if ( cellNormals[2] <= -one && zb <= bounds[4] ) // Bottom (region 5)
		  reg = 5;
      else  // Body (region 6)
		  reg = 6;

		  newScalars->InsertValue(i,reg);
    }

  // Copy all the input geometry and data to the output.
  output->CopyStructure(input);
  output->GetPointData()->PassData(input->GetPointData());
  output->GetCellData()->PassData(input->GetCellData());

  // Add the new scalars array to the output.
  newScalars->SetName("Regions");
  output->GetCellData()->AddArray(newScalars);
  output->GetCellData()->SetActiveScalars("Regions");

  return 1;
}