/*# ---------------------------------------------------------------------
#
# Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
#                        pour la Sant�)
# Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
# Previous Authors : Laurent Guigues, Jean-Pierre Roux
# CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
#
#  This software is governed by the CeCILL-B license under French law and
#  abiding by the rules of distribution of free software. You can  use,
#  modify and/ or redistribute the software under the terms of the CeCILL-B
#  license as circulated by CEA, CNRS and INRIA at the following URL
#  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
#  or in the file LICENSE.txt.
#
#  As a counterpart to the access to the source code and  rights to copy,
#  modify and redistribute granted by the license, users are provided only
#  with a limited warranty  and the software's author,  the holder of the
#  economic rights,  and the successive licensors  have only  limited
#  liability.
#
#  The fact that you are presently reading this means that you have had
#  knowledge of the CeCILL-B license and that you accept its terms.
# ------------------------------------------------------------------------ */

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

  Program:   wxMaracas
  Module:    $RCSfile: vtkImagePolyDataSeedConnectivity.cxx,v $
  Language:  C++
  Date:      $Date: 2012/11/15 14:15:18 $
  Version:   $Revision: 1.2 $

  Copyright: (c) 2002, 2003
  License:
  
     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 "vtkImagePolyDataSeedConnectivity.h"

#include <vtkPolyData.h>
#include <vtkImageData.h>
#include <vtkMath.h>
#include <vtkPlane.h>
#include <vtkSphere.h>
#include <vtkImplicitBoolean.h>
#include <vtkImplicitFunctionToImageStencil.h>
#include <vtkImageStencil.h>
#include <vtkObjectFactory.h>
#include <vtkImageSeedConnectivity.h>
#include <vtkImageIterator.h>
#include <vtkImageGaussianSmooth.h>
#include <vtkMarchingCubes.h>
#include <vtkPolyDataConnectivityFilter.h>
#include <vtkCleanPolyData.h>
#include <vtkTriangleFilter.h>
#include <vtkPointData.h>
#include <vtkImageContinuousDilate3D.h>

vtkCxxRevisionMacro(vtkImagePolyDataSeedConnectivity, "$Revision: 1.2 $");
vtkStandardNewMacro(vtkImagePolyDataSeedConnectivity);
vtkCxxSetObjectMacro(vtkImagePolyDataSeedConnectivity, Axis,vtkPolyData);

//----------------------------------------------------------------------------
vtkImagePolyDataSeedConnectivity::vtkImagePolyDataSeedConnectivity()
{
  this->Axis = NULL;
  this->ClipImageData = vtkImageData::New();
  this->OuterMold = vtkPolyData::New();
  this->ThresholdRatio = 0.45;
}
//----------------------------------------------------------------------------
vtkImagePolyDataSeedConnectivity::~vtkImagePolyDataSeedConnectivity()
{
  if(this->Axis)
  {
    this->Axis->Delete();
  }
  this->ClipImageData->Delete();
  this->ClipImageData = NULL;
  this->OuterMold->Delete();
  this->OuterMold = NULL;
}
 /**
  * Utilisation de static pour accel�rer l'execution ?
  * cf: vtkDividingCubes.cxx
  */

//Macros definitions
#define PSC_MIN(x,y)      ((x)<(y) ? (x) : (y))
#define PSC_MAX(x,y)      ((x)>(y) ? (x) : (y))
//----------------------------------------------------------------------------
void vtkImagePolyDataSeedConnectivity::Execute()
{
  vtkImageData *input = this->GetInput();
  vtkPolyData *output = this->GetOutput();
  int intervalle = 2;
  int dim[3];

  vtkDebugMacro(<< "Executing polydata seed connectivity...");

  if (input == NULL)
    {
    vtkErrorMacro(<<"Input is NULL");
    return;
    }

  // just deal with volumes
  if ( input->GetDataDimension() != 3 )
    {
    vtkErrorMacro("Bad input: only treats 3D structured point datasets");
    return;
    }

  input->GetDimensions(dim);

  //Allocate a temporary image that will be filled in b&w
  //this is the base of the algothm, from a graysccale input + an axis
  //we produce a b&w image the the marching cubes can contour:

  vtkImageData *imagedata = vtkImageData::New();
  imagedata->SetOrigin( input->GetOrigin() );
  imagedata->SetSpacing( input->GetSpacing() );
  imagedata->SetExtent( input->GetExtent() );
  imagedata->SetScalarTypeToUnsignedChar (); //for connectityseed + reduce mem usage
  imagedata->AllocateScalars();

  /*
  Result: piece of crap !!
  There was some leakage from another vessel that was near enough from the heart.
  */

  vtkImageSeedConnectivity *seed = vtkImageSeedConnectivity::New();
  seed->SetInput( imagedata );
  seed->SetInputConnectValue( 255 ); //If 0 -> inverse video, 255 -> rien !

  // Before starting the loop we should make sure the image has been clipped
  // with the polydata:
  this->ClipImageWithAxis();

  vtkPoints *points = this->Axis->GetPoints();
  vtkDataArray* scalars = this->Axis->GetPointData()->GetScalars();
  
  int numPts = points->GetNumberOfPoints();
  double p1[3], p2[3], p3[3], tuple[2];
  double radius, signalvalue;
  int region[6];
  int i_p2; //allow saving of p2 index
  unsigned short* inSI, *inSIEnd;
  unsigned char*  outSI,*outSIEnd;
    
  for(int i=0; i<numPts; i+=intervalle)
  {
    i_p2 = i+intervalle/2;
    points->GetPoint(i, p1);
    if(numPts <= i+2*intervalle)
    {
      i_p2 = i+intervalle;
      points->GetPoint(numPts - 1, p3);
      i = numPts;	//no more iteration
    }
    else
    {
      points->GetPoint(i+intervalle, p3);
    }
    points->GetPoint(i_p2, p2);
    
    //Now we can add p2 point to the see list:
    //seed->AddSeed( 3, p2); //too bad :(
    seed->AddSeed( (int)p2[0], (int)p2[1], (int)p2[2] );

    //this parameter seems to be ok for most image
    //furthermore it doesn't make any change to make it bigger.
    scalars->GetTuple(i_p2, tuple);
    radius		= tuple[0]; 
	signalvalue = tuple[1];
    radius *= 4;		//multiply for a factor 4 for aneurism case

    region[0] = (int)PSC_MAX(p2[0] - radius, 0);
    region[2] = (int)PSC_MAX(p2[1] - radius, 0);
    region[4] = (int)PSC_MAX(p2[2] - radius, 0);

    region[1] = (int)PSC_MIN(p2[0] + radius, dim[0] - 1);
    region[3] = (int)PSC_MIN(p2[1] + radius, dim[1] - 1);
    region[5] = (int)PSC_MIN(p2[2] + radius, dim[2] - 1);

    // multiply signal by a factor this allow user to accuratley find the best
    // range
    signalvalue *= ThresholdRatio;  

      /**
       * Instead of working on the real image, we might take advantage of
       * class Stencil made by dgobbi
       *
       * Steps:
       *  - Create a fake disk (black) of radius ~ Optimale sphere and normal
       *    orientation = last point first eigen vector
       * See: Disk-Stencil.py for more info
       */
      //Solution: this is too much a trouble: rather do it once on the entry (=first  axis point)

      vtkImageIterator<unsigned short> inIt(this->ClipImageData, region);//image_woheart
      vtkImageIterator<unsigned char> outIt(imagedata, region);
      // vtkImageProgressIterator should be nicer in conjonction with a wxGauge

      // Loop through ouput pixels
      while (!inIt.IsAtEnd())
      {
        inSI = inIt.BeginSpan(); inSIEnd = inIt.EndSpan();
        outSI = outIt.BeginSpan(); outSIEnd = outIt.EndSpan();

        // Pixel operation
        while (inSI != inSIEnd)
          *outSI++ = (signalvalue <= *inSI++) ? 255 : 0;

        inIt.NextSpan(); outIt.NextSpan();
      }
   }
   
  //NEW: introduced by Maciek/Marcela: use a vtkImageGaussian to get rid 
  //of the crenelage/steps introduced by MarchingCubes
  
  vtkImageGaussianSmooth *gauss = vtkImageGaussianSmooth::New();
  gauss->SetInput( seed->GetOutput() );
  gauss->SetDimensionality(3);
  //gauss->SetStandardDeviation(1.0);
  //gauss->SetRadiusFactor(3.);

  vtkMarchingCubes *contour = vtkMarchingCubes::New();
  contour->SetInput ( gauss->GetOutput() );
  contour->SetNumberOfContours( 1 );
  contour->SetValue(0, 128);
  contour->ComputeGradientsOn ();
  contour->ComputeScalarsOn ();

  vtkCleanPolyData *clean = vtkCleanPolyData::New();
  clean->SetInput ( contour->GetOutput() );

  vtkTriangleFilter *tf = vtkTriangleFilter::New();
  tf->SetInput( clean->GetOutput() );
  tf->Update();

  //Need to do a deep copy as we'll do a ShallowCopy afterwards
  output->DeepCopy( tf->GetOutput() );

  //output->Squeeze();  //usefull ?

  // Now we do a dilatation of seed->GetOutput() in order to have inner (=output)
  // and outer mold (=OuterMold):
 /**
  * Internal function that take an input of a b&w or grayscale image (should be 
  * unsigned short/ unsigned char and apply a dilatation (see math morpho for
  * info).
  * NB: dilation o gaussian = gaussian o dilatation 
  */

  //First apply a morpho math: dilatation of 2mm:
  vtkImageContinuousDilate3D *dilate = vtkImageContinuousDilate3D ::New();
  dilate->SetInput ( seed->GetOutput() );
  //TODO: determine real kernel based on scale = cm/pixel
  //See: http://public.kitware.com/pipermail/vtkusers/2000-September/004210.html
  //[Spacing is x and y is the distance between the samples. The units of measure
  //are your choice. Often medical images use millimeters. Let's say your medical
  //data resolution is 512 x 512 and your study has a 250 mm field of view. Then
  //the spacing in x and y would be 250/512 or .488 mm.]
  //therefore we have here:

  int kernelsize[3];	
  double spacing[3];
  seed->GetOutput()->GetSpacing( spacing );  //beware input should be Update()
  for(int j=0; j<3; j++)
  {
    //2mm thickness, we assume spacing units is in mm !!
    kernelsize[j] = 3; //(int) (2. / spacing[i]);
    //std::cout << "kernel=" << kernelsize[i] << ",et:" << (int)(2./0.3334) << ",";
    //beware if spacing=0.3334 -> 2/0.3334 = 5 !!!! bad bad +/-0.5 to see
  }
  dilate->SetKernelSize (kernelsize[0], kernelsize[1], kernelsize[2]);

  //VTK filters rules !!!
  gauss->SetInput( dilate->GetOutput());
  tf->Update();
  // shallow copy result into output.
  // Could be re-written, as we know there is no vert neither strips...
  OuterMold->ShallowCopy( tf->GetOutput() );
  
  //OuterMold->Squeeze();  //usefull ?

  //Delete temp stuff:
  imagedata		-> Delete();
  seed			-> Delete();
  gauss			-> Delete();
  contour		-> Delete();
  clean			-> Delete();
  tf			-> Delete();
  dilate		-> Delete();
  
  //TODO: if this filter is used many times when working on different image size
  //I guess we need to find something brighter for this->ClipImageData
}
//end of macros use
#undef PSC_MAX
#undef PSC_MIN
//----------------------------------------------------------------------------
//We need to add a function that from the raw data is able to clip with
//the begining and ending of polydata to produce sharp edge.
void vtkImagePolyDataSeedConnectivity::ClipImageWithAxis()
{
  //no input need we'll work with member data:
  //Now we only need the region growing part:
 /**
  * Instead of working directly on leo we might remove the heart by a  trick:
  * If we could cut with a disk at place where the last(first) axis point is
  * we could avoid the flood fill algo to spread in the heart...smart  !
  */

 /**
  * First we construct a virtual disk (=half sphere for now)
  * See: Disk-Stencil.py
  */
  double normal[3];

 /**
  * What we do here:
  * Get first point (p2) and second point (p1) : assume it exists !
  * Then calculate vector p2->p1
  * normalize it
  * use it to define intersection plane between sphere and plane
  */
  this->Axis->Update();
  vtkPoints *points = this->Axis->GetPoints();
  int numPts = points->GetNumberOfPoints();
  double p1[3];
  double p2[3];
  
  points->GetPoint( 0, p2 );
  points->GetPoint( 1, p1 );
  
  normal[0] = p1[0] - p2[0];
  normal[1] = p1[1] - p2[1];
  normal[2] = p1[2] - p2[2];
  vtkMath::Normalize( normal );

  vtkPlane *plane1 = vtkPlane::New();
  plane1->SetOrigin( p2[ 0 ], p2[ 1 ], p2[ 2 ] );
  plane1->SetNormal( normal[0], normal[1], normal[2]);

  vtkDataArray* scalars = this->Axis->GetPointData()->GetScalars();
  double tuple[2];
  scalars->GetTuple(0, tuple);
  double radius = tuple[0];

   vtkSphere *sphere1 = vtkSphere::New();
   sphere1->SetCenter( p2[ 0 ], p2[ 1 ], p2[ 2 ] );
   sphere1->SetRadius( radius + 1);

   vtkImplicitBoolean *disk1 = vtkImplicitBoolean::New();
   disk1->SetOperationTypeToIntersection();
   disk1->AddFunction(sphere1);
   disk1->AddFunction(plane1);

 /**
  * What we do here:
  * Get last point (p2) and before last point (p1) : assume it exists !
  * Then calculate vector p2->p1
  * normalize it
  * use it to define intersection plane between sphere and plane
  */

  //ax->getControlPoint(numPts - 1, p2, p2+3);
  points->GetPoint( numPts - 1, p2);
  //ax->getControlPoint(numPts - 2, p1, p1+3); //assume it exists
  points->GetPoint( numPts - 2, p1);
  normal[0] = p1[0] - p2[0];
  normal[1] = p1[1] - p2[1];
  normal[2] = p1[2] - p2[2];
  vtkMath::Normalize( normal );

  vtkPlane *plane2 = vtkPlane::New();
  plane2->SetOrigin( p2[ 0 ], p2[ 1 ], p2[ 2 ] );
  plane2->SetNormal( normal[0], normal[1], normal[2]);

  vtkSphere *sphere2 = vtkSphere::New();
  sphere2->SetCenter( p2[ 0 ], p2[ 1 ], p2[ 2 ] );
  sphere2->SetRadius( radius + 1 );

  vtkImplicitBoolean *disk2 = vtkImplicitBoolean::New();
  disk2->SetOperationTypeToIntersection();
  disk2->AddFunction(sphere2);
  disk2->AddFunction(plane2);

  vtkImplicitBoolean *disk = vtkImplicitBoolean::New();
  disk->SetOperationTypeToUnion();
  disk->AddFunction(disk1);
  disk->AddFunction(disk2);

 /**
  * Now use dgobbi class: Stencil stuff to remove the 'bad' heart part
  */

  vtkImplicitFunctionToImageStencil *functionToStencil = vtkImplicitFunctionToImageStencil::New ();
  functionToStencil->SetInput ( disk );

  vtkImageStencil *stencil = vtkImageStencil::New();
  stencil->SetInput( this->GetInput() );
  stencil->SetStencil( functionToStencil->GetOutput() );
  stencil->ReverseStencilOn();
  stencil->SetBackgroundValue (0);
  //don't forget to call ->Update before using an iterator !!!
  stencil->Update();
  
  this->ClipImageData->ShallowCopy( stencil->GetOutput() );
  
 /**
  * Should I do my own vtkCappedCylinder ??
  * solution: no , cause ... I can't remember why... :p
  */
  
  //Delete stuff:
  plane1			-> Delete();
  sphere1			-> Delete();
  disk1				-> Delete();
  plane2			-> Delete();
  sphere2			-> Delete();
  disk2				-> Delete();
  disk				-> Delete();
  functionToStencil	-> Delete();
  stencil			-> Delete();
}
//----------------------------------------------------------------------------
void vtkImagePolyDataSeedConnectivity::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
  if ( this->Axis )
    {
    os << indent << "Axis of " << this->Axis->GetNumberOfPoints()
       << "points defined\n";
    }
  else
    {
    os << indent << "Axis not defined\n";
    }
}