// marAxisCT.cpp: implementation of the marAxisCT class.
//
//////////////////////////////////////////////////////////////////////


#include <vtkImageCast.h>
#include <vtkImageGaussianSmooth.h>
#include <vtkCleanPolyData.h>
#include <vtkPolyDataConnectivityFilter.h>
#include <vtkStripper.h>
#include <vtkCell.h>
#include <vtkIdList.h>
#include <vtkMatrix4x4.h>
#include <vtkTransform.h>
#include <vtkPlaneSource.h>
#include <vtkContourFilter.h>
#include <marLine.h>
#include <marUtils.h>
#include <vtkImageContinuousErode3D.h>
#include <vtkImageThreshold.h>
#include <vtkImageSobel2D.h>
#include <vtkCellArray.h>
#include <vtkImageLogic.h>


#include "marAxisCT.h"
#include "marGdcmDicom.h"




marAxisCT::marAxisCT(  ) : marAxis(  )                             
{
}

// ----------------------------------------------------------------------------
marContour* marAxisCT::getContour( int point , kVolume* vol, int index ) {
	if (quantContours[point] == NULL || quantContours[point]->getContour(index)->getContour() == NULL)
	{
		createContours(point,vol);
	}

	return quantContours[point]->getContour(index)->getContour();
} 

// ----------------------------------------------------------------------------
vtkPoints* marAxisCT::get3Dcontour( int point , kVolume* vol, int index ) {
	if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get3DContour() == NULL)
	{
		create3Dcontours(point,vol); 
	}

	marAxisContours* cont = quantContours[point];

	return cont->getContour(index)->get3DContour();
}

// ----------------------------------------------------------------------------
vtkPolyData* marAxisCT::get2Dcontour( int point , kVolume* vol, int index ) {
	if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DContour() == NULL)
	{
		create2Dcontours(point,vol); 
	}

	marAxisContours* cont = quantContours[point];

	

	return cont->getContour(index)->get2DContour();
}

// ----------------------------------------------------------------------------
vtkPoints* marAxisCT::get2DDiameterMin( int point , kVolume* vol, int index ) {
//	if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DDiameterMin() == NULL)
//	{
		create2DDiametersMin(point,vol); 
//	}

	marAxisContours* cont = quantContours[point];

	return cont->getContour(index)->get2DDiameterMin();
}

// ----------------------------------------------------------------------------
vtkPoints* marAxisCT::get2DDiameterMax( int point , kVolume* vol, int index ) {
//	if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DDiameterMax() == NULL)
//	{
		create2DDiametersMax(point,vol); 
//	}

	marAxisContours* cont = quantContours[point];

	return cont->getContour(index)->get2DDiameterMax();
}

// ----------------------------------------------------------------------------
int marAxisCT::getSignal(int point, int index, kVolume* vol) {
	if (quantContours[point] == NULL )
	{
		createSignals(point,vol); 
	}

	return (int) ( quantContours[point]->getContour(index)->getSignal() );
}

// ----------------------------------------------------------------------------
void marAxisCT::createContours(	int point, kVolume* vol ) {

	int j,id;
	int numberOfPoints,numberOfCells;
	vtkCell* cell;
	vtkIdList* pids;
	double p[ 3 ];
	vtkPolyData* vtkPolyData_2DContour;
	double x,y;

	int sizeIma				=	getParameters( )->getSizeIma( );
	double dimIma			=	getParameters( )->getDimIma( );

	if (quantContours.size() < point - 1 )
	{
		create2Dcontours(point,vol);
	}

	marAxisContours* axisCont = quantContours[point];
	
	for (int i = 0; i < axisCont->getSize(); i++)
	{
		marContourVO* aContourVO = axisCont->getContour(i);
		aContourVO->setContour(new marContour( point, getParameters( ) ));
		vtkPolyData_2DContour = aContourVO->get2DContour();
		numberOfCells			=	vtkPolyData_2DContour->GetNumberOfCells( );
		cell					=	vtkPolyData_2DContour->GetCell( 0 );
		pids					=	cell->GetPointIds( );
		numberOfPoints			=	pids->GetNumberOfIds( );
		for( j = 0; j < numberOfPoints; j++ ) 
		{
			id = pids->GetId( j );
			vtkPolyData_2DContour->GetPoint( id, p);
			x=p[0]-64.0;
			y=p[1]-64.0;
			x=x * dimIma / ( float ) sizeIma;
			y=y * dimIma / ( float ) sizeIma;
			aContourVO->getContour()->addContourPoint( x , y );
		}

		aContourVO->getContour()->do_spline();
		aContourVO->getContour()->calculateVariables();
	}
	
}

// ----------------------------------------------------------------------------
void marAxisCT::create3Dcontours(	int point	, kVolume* vol ) {
	
	if (quantContours[point] == NULL )
	{
		createContours(point,vol);
	}

	marAxisContours* axisCont = quantContours[point];

	for (int i = 0; i < axisCont->getSize(); i++)
	{
		vtkPoints *_vtkPoints;
		double *o;
		double *n;
		vtkMatrix4x4* mat;
		vtkTransform* trans;
		double nn[3];
		double c[3],p1[3],p2[3];
		double nA,nB,nC;

		_vtkPoints						= vtkPoints::New();
		o								= getPoints(point); //PENDIENTE  _points[i];
		n								= getNormal( i );
		mat								= vtkMatrix4x4::New();
		trans							= vtkTransform::New();


		nn[0]=n[0];  nn[1]=n[1];  nn[2]=n[2];
		nC=sqrt( nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2] );
		nn[0]=nn[0]/nC;  nn[1]=nn[1]/nC;  nn[2]=nn[2]/nC;

		vtkPlaneSource* pSource = vtkPlaneSource::New( );
		pSource->SetOrigin( 0, 0	, 0				);
		pSource->SetPoint1( 1, 0	, 0				);
		pSource->SetPoint2( 0, 0	, 1.0	);
		pSource->SetCenter( 0,0,0 );
		pSource->SetNormal( nn[ 0 ], nn[ 1 ], nn[ 2 ] );
		pSource->Update( );
		pSource->Update( );
		pSource->GetOrigin(c);
		pSource->GetPoint1(p1);
		pSource->GetPoint2(p2);
		pSource->Delete();
		p1[0]=p1[0]-c[0];  p1[1]=p1[1]-c[1];  p1[2]=p1[2]-c[2];
		p2[0]=p2[0]-c[0];  p2[1]=p2[1]-c[1];  p2[2]=p2[2]-c[2];
		nA=sqrt( p1[0]*p1[0] + p1[1]*p1[1] + p1[2]*p1[2] );
		nB=sqrt( p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2] );
		p1[0]=p1[0]/nA;  p1[1]=p1[1]/nA;  p1[2]=p1[2]/nA;
		p2[0]=p2[0]/nB;  p2[1]=p2[1]/nB;  p2[2]=p2[2]/nB;

		mat->SetElement (0,0, nn[0]);
		mat->SetElement (1,0, nn[1]);
		mat->SetElement (2,0, nn[2]);
		mat->SetElement (3,0, 0);
		mat->SetElement (0,1, p2[0]);
		mat->SetElement (1,1, p2[1]);
		mat->SetElement (2,1, p2[2]);
		mat->SetElement (3,1, 0);
		mat->SetElement (0,2, p1[0]);
		mat->SetElement (1,2, p1[1]);
		mat->SetElement (2,2, p1[2]);
		mat->SetElement (3,2, 0);
		mat->SetElement (0,3, 0);
		mat->SetElement (1,3, 0);
		mat->SetElement (2,3, 0);
		mat->SetElement (3,3, 1);

		double deter=mat->Determinant(mat);
		trans->Identity();
		trans->SetMatrix(mat);
		float ang;
		ang=-90;
		trans->RotateWXYZ  ( ang,0,1,0); 
		trans->Update();

		vtkMatrix4x4 *m=vtkMatrix4x4::New();
		trans->GetMatrix  (  m  );
    
		int j,numberOfPoints;

		marContour* contour2D	= getContour(point,vol,i);
		marContourVO* aContourVO = axisCont->getContour(i);
		numberOfPoints			= contour2D->getNumberOfPoints( );
		double fpA[3];
		double pA[3],ppp[3];

		for( j = 0; j <= numberOfPoints; j++ ) {
			contour2D->getPoint( fpA,j % numberOfPoints);
			pA[0]=fpA[0];
			pA[1]=fpA[1];
			pA[2]=0;

			ppp[0] = m->GetElement(0,0)*pA[0] + m->GetElement(0,1)*pA[1] + m->GetElement(0,2)*pA[2] + m->GetElement(0,3);
			ppp[1] = m->GetElement(1,0)*pA[0] + m->GetElement(1,1)*pA[1] + m->GetElement(1,2)*pA[2] + m->GetElement(1,3);
			ppp[2] = m->GetElement(2,0)*pA[0] + m->GetElement(2,1)*pA[1] + m->GetElement(2,2)*pA[2] + m->GetElement(2,3);

			ppp[0]=ppp[0]+o[0]-c[0]; ppp[1]=ppp[1]+o[1]-c[1]; ppp[2]=ppp[2]+o[2]-c[2];  
			_vtkPoints->InsertNextPoint( (float)ppp[0],(float)ppp[1],(float)ppp[2] );

		
		}

	
		m->Delete();
		mat->Delete();
		trans->Delete();

		if (aContourVO->get3DContour()!=NULL){ aContourVO->get3DContour()->Delete(); }
			aContourVO->set3DContour(_vtkPoints);
	
	}


	
	
}

// ----------------------------------------------------------------------------
void marAxisCT::create2Dcontours(	int point	, kVolume* vol ) {
	std::vector <marIsocontour *> contours;

	

	generatePoints(point,vol,&contours);
	
	vtkImageData* imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
	marAxisContours* axisCont = new marAxisContours();

	for (int i = 0; i < contours.size(); i++)
	{
		marContourVO* contourVO = new marContourVO();
		marIsocontour* iso = contours[i];
	
		contourVO->setIsocontour(vtkContourFilter::New());
		contourVO->getIsocontour()->SetInput( imagedata );
		contourVO->getIsocontour()->SetNumberOfContours( 1 );
		contourVO->getIsocontour()->SetValue( 0, iso->getMaxIntensity() );	
		contourVO->getIsocontour()->Update( );
		
		contourVO->setIsocontourCpd(vtkCleanPolyData::New());
		contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
		contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
		contourVO->getIsocontourCpd()->Update( );

		double cgx, cgy;

		iso->getCG(&cgx,&cgy);
		
		contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
		contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
		
		contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
		contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
		contourVO->getIsocontourDcf()->Update( );
		
		contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
		contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
		contourVO->getIsocontourCpd2()->Update();

		contourVO->setIsocontourStripped(vtkStripper::New( ));
		contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
		contourVO->getIsocontourStripped()->Update();
	
		contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
		contourVO->get2DContour()->Update();
		contourVO->setSignal(iso->getMaxIntensity());
		
		if (iso->getType() == marAxisContours::LUMEN)
		{
			contourVO->setType(marAxisContours::WALL);
		}
		else
		{
			contourVO->setType(iso->getType());
		}
		
		
		
		axisCont->addContour(contourVO);

		
		delete iso;
		contours[i] = NULL;
	}

	

	quantContours[point] = axisCont;
	


	updateLumenPercentage(point,vol);
	
	updateCalcPercentage(point,vol);
//	adjustWall(point,vol);
	if (calibration && point < startIndex)
	{
		this->adjustContour(point,vol);	
	}
	else
	{
		adjustWall(point,vol);
	}
	contours.clear();
	imagedata->Delete();


//	marIsocontour* cont = parsePolyDataToMarIsocontour(quantContours[point]->getContour(0));

	/*FILE *archivo;
	char nomArchivo[30], temp[3];

	strcpy(nomArchivo,"./points");
	itoa(point,temp,10);
	strcat(nomArchivo,temp);
	strcat(nomArchivo,".csv");
	
	archivo = fopen(nomArchivo,"w");
	
	for (int h = 0; h < cont->getSize(); h++)
	{
		double x, y;
		cont->getPoint(h,&x,&y);
		fprintf(archivo,"%f;%f\n",x,y);
	}

	fclose(archivo);*/
}

// ----------------------------------------------------------------------------
void marAxisCT::create2DDiametersMin(int point	, kVolume* vol ) {
	
if (quantContours[point] == NULL )
	{
		createContours(point,vol);
	}

	marAxisContours* axisCont = quantContours[point];

	for (int i = 0; i < axisCont->getSize(); i++)
	{
		double p1[2],p2[2];
		marContour *marcontour=getContour(point,vol,i);
		marContourVO* aContourVO = axisCont->getContour(i);

		marcontour->getMinimumLine(p1,p2);
		vtkPoints *_vtkPoints = vtkPoints::New();
		
		int sizeIma = getParameters( )->getSizeIma( );
		double dimIma = getParameters( )->getDimIma( );
		double factor=( float ) sizeIma / dimIma;
		
		p1[0]=p1[0]*factor+sizeIma/2;
		p1[1]=p1[1]*factor+sizeIma/2;
		p2[0]=p2[0]*factor+sizeIma/2;
		p2[1]=p2[1]*factor+sizeIma/2;
		_vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
		_vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );

		if ( aContourVO->get2DDiameterMin()!=NULL ) { aContourVO->get2DDiameterMin()->Delete(); }
			aContourVO->set2DDiameterMin(_vtkPoints);

	}

}

// ----------------------------------------------------------------------------
void marAxisCT::create2DDiametersMax(int point	, kVolume* vol ) {
	if (quantContours[point] == NULL )
	{
		createContours(point,vol);
	}

	marAxisContours* axisCont = quantContours[point];

	for (int i = 0; i < axisCont->getSize(); i++)
	{
		double p1[2],p2[2];
		marContour *marcontour=getContour(point,vol,i);
		marContourVO* aContourVO = axisCont->getContour(i);

		marcontour->getMaximumLine(p1,p2);
		vtkPoints *_vtkPoints = vtkPoints::New();

		int sizeIma = getParameters( )->getSizeIma( );
		double dimIma = getParameters( )->getDimIma( );
		double factor=( float ) sizeIma / dimIma;
		p1[0]=p1[0]*factor+sizeIma/2;
		p1[1]=p1[1]*factor+sizeIma/2;
		p2[0]=p2[0]*factor+sizeIma/2;
		p2[1]=p2[1]*factor+sizeIma/2;
		_vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
		_vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
		
		if ( aContourVO->get2DDiameterMax()!=NULL ) { aContourVO->get2DDiameterMax()->Delete(); }
			aContourVO->set2DDiameterMax(_vtkPoints);
	}
}

// ----------------------------------------------------------------------------
void marAxisCT::createSignals (int point, kVolume* vol) {
	if (quantContours[point] == NULL )
	{
		create2Dcontours(point,vol); //AQUI
	//	filterSignal(vol);
	}
}

// ----------------------------------------------------------------------------
int	marAxisCT::getNumberOfContours(int point, kVolume* vol) {
	if (quantContours[point] == NULL )
	{
		create2Dcontours(point,vol);//AQUI
	//	filterSignal(vol);
		
	}

	return quantContours[point]->getSize();
}

// ----------------------------------------------------------------------------
int	marAxisCT::getContourType(int point, int index, kVolume* vol) {
	if (quantContours[point] == NULL )
	{
		create2Dcontours(point,vol);
	}
	
	return quantContours[point]->getContourType(index);

}

// ----------------------------------------------------------------------------
void marAxisCT::generatePoints(int point, kVolume* vol, std::vector<marIsocontour*> *list) {

	vtkImageData* imagedata;
	int dirs = 72;
	double aPoint[ 3 ];
	double *nPoint = new double[ 3 ];
	double x, y, prevIntensity, threshold;
	bool found, estStart, estDivided, isNeighbor;
	int actualZone, startZone, coordBase, maxPoint, lumenCount, calcCount;
	std::vector <marPoint  *> gradients;
	std::vector <double> average;
	std::vector <marIsocontour *> contoursLumen;
	std::vector <marIsocontour *> contoursCalc;

	contoursLumen.push_back(NULL);
	contoursCalc.push_back(NULL);
	
	
	
	//1. Check starting zone
	startZone = marAxisContours::LUMEN;
	actualZone = startZone;
	lumenCount = 0;
	calcCount = 0;

	imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));

	int limInfX, limInfY, limInfZ;
	int limSupX, limSupY, limSupZ;
	imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);

	vtkImageCast* cast = vtkImageCast::New();
	cast->SetInput(imagedata);
	cast->SetOutputScalarTypeToDouble();

	vtkImageGaussianSmooth* filter = vtkImageGaussianSmooth::New();
	filter->SetDimensionality(2);
	filter->SetInput(cast->GetOutput());
	filter->SetStandardDeviations((getParameters())->getStandardDeviation(),(getParameters())->getStandardDeviation());
	filter->SetRadiusFactors((getParameters())->getRadius(),(getParameters())->getRadius());
	filter->SetStandardDeviations(3,3);
	filter->SetRadiusFactors(3,3);

	//imagedata = filter->GetOutput();
	
	
	nPoint = getPoints(point);

	aPoint[ 0 ]	=  (limSupX - limInfX) / 2; //nPoint[0]; 
	aPoint[ 1 ]	=  (limSupY - limInfY) / 2; //nPoint[1]; 
	aPoint[ 2 ]	=  (limSupZ - limInfZ) / 2; // nPoint[2]; 
	estDivided = false;
	estStart = false;
	//3. Generate rays over all directions
	for (int dir = 0; dir < dirs; dir++)
	{
		
		x = aPoint[ 0 ];
		y = aPoint[ 1 ];
		coordBase = 1;
		prevIntensity = interpolate(x,y,imagedata);
		generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);

		//2.1 Evaluate gradient in every position of the ray 
		while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
		{
			
			//2.1.1 Interpolate
			double resp = interpolate(x,y,imagedata);
			//2.1.2 Calcultate gradient
			marPoint* p = new marPoint(x,y);
			p->setGradient(prevIntensity - resp);
			p->setDirection(dir);
			p->setIntensity(resp);
			
			prevIntensity = resp;

			gradients.push_back(p);

			//2.1.3 Increase ray for next evaluation
			generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);

		}// end while gradient

		int index = 0;
		found = false;
		x = aPoint[ 0 ];
		y = aPoint [ 1 ];
		threshold = ((marParameters *) getParameters())->getContourThresh();
		int mysize = gradients.size();
		
		//2.2 Compare every gradient value with threshold
		while (index < gradients.size() && !found)
		{
			//2.2.1 Evaluate Lumen
			if (actualZone == marAxisContours::LUMEN)
			{
				average.push_back(gradients[index]->getIntensity());
				
				//2.2.1.1 If value > threshold there is a contour
				if (fabs(gradients[index]->getGradient()) > threshold 
					&& fabs(gradients[index]->getGradient()) < 4000) {
						//|| (calibration && (abs(gradients[index]->getIntensity()) > abs(avgValue - 2.5*stdValue)
						//|| abs(gradients[index]->getIntensity()) > abs(avgValue + 2.5*stdValue)))) {

					

					if(gradients[index]->getGradient() > 0)
						maxPoint = getMaximumGrad(gradients,index,gradients.size(),threshold);
					else
						maxPoint = getMinimumGrad(gradients,index,gradients.size(),threshold);

					index = maxPoint;
					x = gradients[maxPoint]->getX();
					y = gradients[maxPoint]->getY();

					//2.2.1.1.1 Assign lumen structure
					isNeighbor = detectNeighbor(contoursLumen[lumenCount],dir,marAxisContours::LUMEN);
					if (!isNeighbor)
					{
						lumenCount++;
					}
					
					marIsocontour* cont = addPointToContour(contoursLumen[lumenCount],false,marAxisContours::LUMEN,gradients[maxPoint]);
					contoursLumen[lumenCount] = cont;

					//2.2.1.1.2 Check step lumen->calcification
					double actualInt = interpolate(x,y,imagedata);
					double auxX = x;
					double auxY = y;
					generateVector(dir,&coordBase,aPoint[0],aPoint[1],&auxX,&auxY);
				    double next = interpolate(auxX,auxY,imagedata);
					if (gradients[index]->getGradient() < 0 /*&& actualInt < next */&& (auxX > (limInfX+1) && auxX < (limSupX-1)) && (auxY > (limInfY+1) && auxY < (limSupY-1)))
					{
						//2.2.1.1.2.1 Assign calcification structure
						isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
						if (!isNeighbor)
						{
							calcCount++;
							contoursCalc.push_back(NULL);
						}

						cont = addPointToContour(contoursCalc[calcCount],true,marAxisContours::CALCIFICATION,gradients[maxPoint]);
						
						contoursCalc[calcCount] = cont;

						actualZone = marAxisContours::CALCIFICATION;
						//ASIGNACION varX???

						//2.2.1.1.2.3 Verify divided structure
						if (dir == 0)
						{
							estStart = true;
						}

						if ((dir == dirs - 1) && estStart)
						{
							estDivided = true;
						}

					}//2.2.1.1.3 Transition lumen->background
					else
					{
						//2.2.1.1.3.1 Asign border structure
						//ASIGNAR ESTRUCTURA
						
						//2.2.1.1.3.2 Generate stop condition
						found = true;
						actualZone = startZone;
						coordBase = 1;
						x = aPoint[ 0 ];
						y = aPoint[ 1 ];
						threshold = ((marParameters *) getParameters())->getContourThresh();
					}

				}
			}
			//2.2.2 Evaluate calcification (MODIFIED IF THERE IS HYPODENSE)
			else
			{
				coordBase = 1;
				int calcIndex;
				double actualInt = interpolate(x,y,imagedata);
				generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
				double next = interpolate(x,y,imagedata);

				//2.2.2.1 If actual intensity is smaller than previous there is a contour
				if (actualInt < next)
				{
					while (actualInt < interpolate(x,y,imagedata) && (x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
					{
						generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
						index++;
					}

				}

				generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);

				//2.2.2.2 Assign calcification structure
				isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
				if (!isNeighbor)
				{
					calcCount++;
					contoursCalc.push_back(NULL);
				}

				//Adjust index to avoiod crash
				if (index >= gradients.size())
				{
					calcIndex = gradients.size() - 1;
				}
				else
				{
					calcIndex = index;
				}

				marIsocontour* cont = addPointToContour(contoursCalc[calcCount],false,marAxisContours::CALCIFICATION,gradients[calcIndex]);
				contoursCalc[calcCount] = cont;
				
				//2.2.2.3 Asign border structure
				//ASIGNACION BORDE

				found = true;
				actualZone = startZone;
				x = aPoint[ 0 ];
				y = aPoint[ 1 ];
				threshold = ((marParameters *) getParameters())->getContourThresh();

			}//end if zones

			index++;

			if (index == gradients.size() && !found)
			{
				threshold = threshold - 1;
				index = 0;
			}

		}//end while evaluation

		//TODO - SACAR A UTILS
		for (int ind = 0; ind < gradients.size(); ind++)
		{
			marPoint* p = gradients[ind];
			gradients[ind] = NULL;
			delete p;
			
		}
		gradients.clear();

	}//end for directions

	FILE *archivo;
	char nomArchivo[30]; //, temp[3]

	strcpy(nomArchivo,"./debug");
//	itoa(point,temp,10);
//	strcat(nomArchivo,temp);
	strcat(nomArchivo,".txt");
	archivo = fopen(nomArchivo,"w");
	
  


	
//	lumenContour[point] = contoursLumen[marAxisContours::LUMEN];
	
	double  tempXi, tempYi;
	int polySides = contoursLumen[marAxisContours::LUMEN]->getSize();
	
    
	marIsocontour* alc = new marIsocontour();
	for (int m=0; m<polySides; m++) {
		contoursLumen[marAxisContours::LUMEN]->getPoint(m,&tempXi,&tempYi);          
		alc = addPointToContour(lumenContour[point],false,marAxisContours::LUMEN,new marPoint(tempXi,tempYi));
		lumenContour[point] = alc;
		fprintf(archivo,"X:%f,Y:%f\n",tempXi,tempYi);
	}
	

	
	//4. Verify divided unified structure
	if (estDivided && calcCount > 0)
	{
		unifyDividedStructure(&contoursCalc);
		calcCount--;
	}   
	
	//5. Obtain statistics
	double lumenInt = getStatistics(contoursLumen[marAxisContours::LUMEN],aPoint[ 0 ],aPoint[ 1 ], imagedata);
	contoursLumen[marAxisContours::LUMEN]->setMaxIntensity(lumenInt);

	if (contoursCalc[0] != NULL)
	{
		for (int i = 0; i < contoursCalc.size(); i++)
		{
			lumenCount++;
			double maxInt = getCalcStatistics(contoursCalc[i],aPoint[ 0 ],aPoint[1],imagedata,i-1);
			contoursCalc[i]->setMaxIntensity(maxInt);
			contoursLumen.push_back(NULL);
			contoursLumen[lumenCount] = contoursCalc[i];
		}
	}
	
	
	*list = contoursLumen;

	contoursCalc.clear();
    contoursLumen.clear();
	marUtils* util = new marUtils();
	double avg = util->obtainAverage(average);
	//FILE *archivo;
	//archivo = fopen("./debug.txt","w");
	fprintf(archivo,"tamaņo: %d umbra: %d",  (*list).size(),((marParameters *) getParameters())->getContourThresh());
	for (int k = 0; k < contoursLumen.size(); k++)
	{
		double cgx, cgy;

		(*list)[k]->getCG(&cgx,&cgy);
		fprintf(archivo,"cx: %f cy:%f Seņal: %f Type: %d\n", cgx, cgy, (*list)[k]->getMaxIntensity(), (*list)[k]->getType());
	}

	fprintf(archivo,"PROMEDIO %f - STD: %f\n",avg,util->obtainStandardDeviation(average,avg));
	fprintf(archivo,"%d\n",lumenContour.size());
	fclose(archivo);
	imagedata->Delete();
	
}

// ----------------------------------------------------------------------------
int marAxisCT::getMaximumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
	
	double max = list[iniPos]->getGradient();
	int coordMax = iniPos;
	int cont = iniPos + 1;

	while (cont < limit && list[cont]->getGradient()>=threshold){
		if (list[cont]->getGradient() > max){
			max = list[cont]->getGradient();
			coordMax = cont;
		}
		cont++;
	}

	return coordMax;
} 

// ----------------------------------------------------------------------------
int marAxisCT::getMinimumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
	
	double min = list[iniPos]->getGradient();
	int coordMin = iniPos;
	int cont = iniPos + 1;

	while (cont < limit && fabs(list[cont]->getGradient())>=threshold){
		if (list[cont]->getGradient() < min){
			min = list[cont]->getGradient();
			coordMin = cont;
		}
		cont++;
	}

	return coordMin;
} 

// ----------------------------------------------------------------------------
double marAxisCT::interpolate(double x, double y, vtkImageData* imagedata) {

	double a,b,c,d;
	double intensity;
	int i,j;

	i = (int)floor(x);
	j = (int)floor(y);

	
	c = imagedata->GetScalarComponentAsDouble(i,j,0,0) + (double) imagedata->GetScalarComponentAsDouble(i+1,j+1,0,0)
		- imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - imagedata->GetScalarComponentAsDouble(i,j+1,0,0);

	b = imagedata->GetScalarComponentAsDouble(i,j+1,0,0) - c*i - imagedata->GetScalarComponentAsDouble(i,j,0,0);

	a = imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - c*j - imagedata->GetScalarComponentAsDouble(i,j,0,0);

	d = imagedata->GetScalarComponentAsDouble(i,j,0,0) - a*i - b*j - c*i*j;

	intensity = a*x + b*y + c*x*y + d;

	return intensity;
}

// ----------------------------------------------------------------------------
void marAxisCT::generateVector(int dir,int *coordBase,double originX,double originY,double *x,double *y) {

	int dirs = 72;
	int angle = (dir * (360 / dirs));

	*x = 1 * cos(angle *(3.1415/180)) + (*x); 
	*y = 1 * sin(angle *(3.1415/180)) + (*y); 

	(*coordBase)++;	
}

// ----------------------------------------------------------------------------
bool marAxisCT::detectNeighbor(marIsocontour* contour,int dir, int type) {
	
	bool isNeighbor = false;

	if (contour == NULL){
		isNeighbor = true; 
	}else{
		for (int i = 0; i < contour->getSize(); i++) {
			
			if (((contour->getDir(i) == (dir - 1)) || (contour->getDir(i) == dir)) 
				    && (contour->getType() == type)){ 
				isNeighbor = true; 
			}	
		}
	}

	return isNeighbor;
}

// ----------------------------------------------------------------------------
marIsocontour*	marAxisCT::addPointToContour(marIsocontour* cont ,bool inside,int type,marPoint* point) {

	

	if (cont == NULL)
	{
		cont = new marIsocontour();
		cont->setType(type);
	}
	
	cont->insertPoint(point->getX(),point->getY());
	cont->setDir(cont->getSize() - 1,point->getDirection());
	cont->setInside(cont->getSize() - 1,inside);
	
	return cont;
	
}

// ----------------------------------------------------------------------------
void marAxisCT::unifyDividedStructure(std::vector<marIsocontour*> *contList) {

	marIsocontour* structOne = (*contList)[contList->size()-1];
	marIsocontour* structTwo = (*contList)[0];
	

	for (int i = 0; i < structOne->getSize(); i++)
	{
		double x;
		double y;
		structOne->getPoint(i,&x,&y);
		structTwo->insertPoint(x,y);
		structTwo->setDir(structTwo->getSize() - 1,structOne->getDir(i));
		structTwo->setInside(structTwo->getSize() - 1,structOne->getInside(i));

	}

	contList->pop_back();
}

// ----------------------------------------------------------------------------
double marAxisCT::getStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata) {

	double max, auxX, auxY;
	int basis;
	max = interpolate(x,y,imagedata);

	for (int i = 0; i < contour->getSize(); i++)
	{
		basis = 1;
		auxX = x;
		auxY = y;
		double pointX;
		double pointY;

		contour->getPoint(i,&pointX,&pointY);
		while (auxX != pointX || auxY != pointY)
		{
			double intensity = interpolate(auxX,auxY,imagedata);
			if (max < intensity)
			{
				max = intensity;
			}
			generateVector(contour->getDir(i),&basis,x,y,&auxX,&auxY);

		}
	}
	
	return max;
}


// ----------------------------------------------------------------------------
double marAxisCT::getCalcStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata, int dir) {

	
	double max, auxX,auxY;
	int basis;
	
	int i;
	FILE *archivo;
	char nomArchivo[30], temp[3];

	strcpy(nomArchivo,"./statistics");

// EED 06/06/2007 linux
//	itoa(dir,temp,10);
	sprintf(temp,"%d",dir);

	strcat(nomArchivo,temp);
	strcat(nomArchivo,".txt");


	archivo = fopen(nomArchivo,"w");
	i = 0; 
	max = 0;
	fprintf(archivo,"TAMAŅO %d\n",contour->getSize());

	while (i < contour->getSize()) 
	{
		basis = 1;
		double pointX, pointY;
		contour->getPoint(i+1,&pointX,&pointY);
		contour->getPoint(i,&auxX,&auxY);
		fprintf(archivo,"final: %f %f inicial %f %f \n",pointX,pointY,auxX,auxY);	
		
		while (auxX != pointX || auxY != pointY)
		{
			double intensity = interpolate(auxX,auxY,imagedata);
			if (max < intensity){
				max = intensity;
			}
			generateVector(contour->getDir(i),&basis,auxX,auxY,&auxX,&auxY);
		}
			
		i +=2;
	}
	fclose(archivo);
	return max;
}

// ----------------------------------------------------------------------------
int marAxisCT::round(double num){

	float decimal;
	int resp;
	
	decimal = num - floor(num) * 1.0000;

	if (decimal >= 0.50000)
		resp = (int)ceil(num);
	else
		resp = (int)floor(num);

	return resp;

}

// ----------------------------------------------------------------------------
void marAxisCT::createEmptyContours(){

	int nCnts = ( int ) getNumberOfSplinePoints( );
	

	for (int i = 0; i < nCnts; i++)
	{
		quantContours.push_back( NULL );
		lumenContour.push_back( NULL);
	}
}

// ----------------------------------------------------------------------------
void marAxisCT::updateLumenPercentage(int point, kVolume* vol)
{
	marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
	if (quantContours[point] == NULL || wallCont->get2DContour() == NULL)
	{
		create2Dcontours(point,vol); 
	}

	if (wallCont->getType() == marAxisContours::WALL
			&& wallCont->isReplaced() == false)
	{
		int a = ((marParameters *) getParameters())->getLumenPercentage();
		double oldValue = wallCont->getSignal();
		double lumenPercentage=((marParameters *) getParameters())->getLumenPercentage();
		double poww=pow(10.0,-2.0);
		double newValue = oldValue * lumenPercentage*poww;
		wallCont->getIsocontour()->SetValue(0,newValue);
		wallCont->getIsocontour()->Update();
		wallCont->get2DContour()->Update(); 
	}
	
}

// ----------------------------------------------------------------------------
void marAxisCT::updateCalcPercentage(int point, kVolume* vol)
{
	double oldValue;
	int posMax;

	if (quantContours[point] == NULL)
	{
		create2Dcontours(point,vol); 
	}

	posMax = getMaxIntensity(point);
	
	if (quantContours[point]->getSize() > 1)
	{
		oldValue = quantContours[point]->getContour(posMax)->getSignal();
	}
	
	for (int i = 1; i < quantContours[point]->getSize(); i++)
	{
		
		if (quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION 
				&& quantContours[point]->getContour(i)->isReplaced() == false)
		{

			double newValue = oldValue * ((((marParameters *) getParameters())->getCalcPercentage())*pow(10.0,-2.0));
			quantContours[point]->getContour(i)->getIsocontour()->SetValue(0,newValue);
			quantContours[point]->getContour(i)->getIsocontour()->Update();
			quantContours[point]->getContour(i)->get2DContour()->Update();
		}
	
	}
	
}

// ----------------------------------------------------------------------------
int	marAxisCT::getMaxIntensity(int point)
{
	int pos = 1; 
	double intMax;
	if (quantContours[point]->getSize() > 2)
	{
		intMax = quantContours[point]->getContour(0)->getSignal();
		for (int i = 0; i < quantContours[point]->getSize(); i++)
		{
			if (quantContours[point]->getContour(i)->getSignal() > intMax && quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION)
			{
				intMax = quantContours[point]->getContour(i)->getSignal() ;
				pos = i;
			}
		}
	}

	return pos;
	
}
// ----------------------------------------------------------------------------
void marAxisCT::eraseContours()
{

	
		int nCts = quantContours.size();
		vesselPoints.clear();

		int i = 0;

		while (i < nCts)
		{
			delete quantContours[i];
			quantContours[i] = NULL;
			lumenContour[i] = NULL;
			i++;
		}
		quantContours.clear();
		lumenContour.clear();

	
}

// ----------------------------------------------------------------------------
void marAxisCT::eraseContoursPartial(int start)
{

	int nCts = 0;
	int i = start;

	while (i >= nCts)
	{
		delete quantContours[i];
		quantContours[i] = NULL;
		lumenContour[i] = NULL;
		i--;
	}
	quantContours.clear();
	lumenContour.clear();

	//TODO Revisar esto que no me convence
}


// ----------------------------------------------------------------------------
bool marAxisCT::pointInPolygon(marIsocontour *c, double x, double y){
	int      i, j=0;
	bool  oddNODES=false;
	int polySides = c->getSize();
	double *polyX, *polyY, tempX, tempY;

	polyX = (double *) malloc(polySides*sizeof(double));
	polyY = (double *) malloc(polySides*sizeof(double));

	for (i=0; i<polySides; i++) {
		c->getPoint(i,&tempX,&tempY);          
		polyX[i] = tempX;
		polyY[i] = tempY;
	}

	for (i=0; i<polySides; i++) {
		j++; 
		if (j==polySides) 
			j=0;
		
		if (polyY[i]<y && polyY[j]>=y
			||  polyY[j]<y && polyY[i]>=y) {
			if (polyX[i]+(y-polyY[i])/(polyY[j]-polyY[i])*(polyX[j]-polyX[i])<x) {
				oddNODES=!oddNODES; 
			}
		}
	}

	free(polyX);
	free(polyY);
  return oddNODES; 

}

// ----------------------------------------------------------------------------
marIsocontour* marAxisCT::parsePolyDataToMarIsocontour(marContourVO* contourVO){

	vtkPolyData* polyDataCnt;
	vtkIdList* pointsId;
	double* pointPolyDataCnt;
	int numPointsPolyDataCnt, numTemp;
	marIsocontour *cont;

	if (contourVO->isReplaced())
	{
		polyDataCnt = contourVO->get2DContour();
	}
	else
	{
		polyDataCnt =    contourVO->getIsocontourStripped()->GetOutput();
	}
	
	numPointsPolyDataCnt = polyDataCnt->GetNumberOfPoints();
	pointsId = polyDataCnt->GetCell(0)->GetPointIds();
	numTemp = pointsId->GetNumberOfIds();
	
	cont = new marIsocontour();
	
	for (int i = 0; i < (numTemp -1); i++){
		pointPolyDataCnt = polyDataCnt->GetPoint(pointsId->GetId(i));
		cont->insertPoint(pointPolyDataCnt[0], pointPolyDataCnt[1]);
		
	}
		
	return cont;	
}



// ----------------------------------------------------------------------------
marIsocontour* marAxisCT::filterContour(marIsocontour* contExt, marIsocontour* contInt, double radio){
  
	double tmp_eed;
	double x1,y1,x2,y2,xInt,yInt, xAux, yAux;
	std::vector<marLine*> lineas;	
	std::vector<double> dists;
	int count, countAux, ultima;
	bool hayInterseccion;

	marIsocontour* newCont = new marIsocontour();

	for (int i = 0; i < contExt->getSize(); i++){
		contExt->getPoint(i,&x1,&y1);
		contInt->getPoint(i,&x2,&y2);
		lineas.push_back(new marLine(x1,y1,x2,y2));	
		tmp_eed = sqrt(pow((x2-x1),2.0) + pow((y2-y1),2.0)); 
		dists.push_back( tmp_eed );
	}

	count = 0;
	ultima = 0;
	while (count < lineas.size() - 1){
		countAux = 0;

		marLine* l1 = lineas[count];
	
		hayInterseccion = false;
		while (countAux < lineas.size()){
			marLine *l2 = lineas[countAux];

			l1->getIntersect(l2->getA(),l2->getB(),l2->getC(),&xInt,&yInt);

			
			if (xInt > 0 && yInt > 0){
				contExt->getPoint(countAux,&xAux,&yAux);
				if (dists[count] >= sqrt(pow((xInt-xAux),2.0) + pow((yInt-yAux),2.0))){
					contExt->getPoint(count,&x1,&y1);
					contInt->getPoint(count,&x2,&y2);

					//Distancia menor?
					if ((dists[count] > sqrt(pow((xInt-x1),2.0) + pow((yInt-y1),2.0))) && 
						(dists[count] > sqrt(pow((xInt-x2),2.0) + pow((yInt-y2),2.0))))
					{
						ultima = countAux;
						hayInterseccion = true;
					
					} 
				
				}
			}
			countAux++;
		}

	
		if (!hayInterseccion){
			contInt->getPoint(count,&x2,&y2);
			newCont->insertPoint(x2,y2);
		} 
			count++;
		
	}

	return newCont;
}

// ----------------------------------------------------------------------------
void marAxisCT::markUpLumen(int point, kVolume* vol)
{
	
	marContourVO* aVO = this->searchContour(marAxisContours::ELUMEN,point);

	if (aVO != NULL)
	{
		return;
	}

	int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
	vtkImageData* imagedata;
	vtkImageData* erodedata;
	std::vector <marIsocontour* > polygons;
	std::vector <marPoint* > points;
	double average = 0, deviation = 0;
	vtkImageData* image = vtkImageData::New();
	vesselPoints.clear();
	

	imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
	imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	image->SetDimensions(limSupX,limSupY,0);
	image->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	image->SetScalarTypeToUnsignedShort();
	image->AllocateScalars();
	image->Update();
	
	vtkImageThreshold *imageThreshold = vtkImageThreshold::New();
	imageThreshold->SetInput(imagedata);
	imageThreshold->ThresholdByUpper(maxSignal);
	imageThreshold->SetInValue(255);
	imageThreshold->SetOutValue(0.0);
	imageThreshold->Update();

	vtkImageContinuousErode3D *imageErode3D = vtkImageContinuousErode3D::New();
	imageErode3D->SetInput(imageThreshold->GetOutput());
	imageErode3D->SetKernelSize(4,4,1);
	imageErode3D->Update();
	erodedata = imageErode3D->GetOutput();

//	polygons.push_back(contExt);
	for (int i = 0; i < quantContours[point]->getSize(); i++)
	{
		polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
	}
	
	//7- Wipe image to detect lumen
	
        int x;
	for (x = limInfX + 30; x < (limSupX - 1); x++)
	{
		
		for (int y = limInfY; y < (limSupY - 1) ; y++)
		{
			bool inWall = false;
			bool outCalc = true;
			
			for (int numConts = 0; numConts < polygons.size(); numConts++)
			{
				if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
				{
					if (numConts == 0)
					{
						inWall = true;
					}
				
					else 
					{
						outCalc = false;
					}
					
				}
				
			}

			if (inWall && outCalc)
			{
				marPoint *p = new marPoint(x,y);
				p->setIntensity(imagedata->GetScalarComponentAsDouble(x,y,0,0));
				average += p->getIntensity();
				points.push_back(p);
			}

		/*	if (cont >= (limSupY - limSupX))
			{
				contin = false;
			}*/
		}
	}
	
	if (points.size() > 0)
	{
		average = average / points.size();
		for (int k = 0; k < points.size(); k++)
		{
			marPoint *p = points[k];
			double actualInt = p->getIntensity();
			double temp = average - actualInt;
			deviation += pow(temp,2.0);
		}
		deviation = sqrt(deviation / points.size());

	}


	polygons.push_back(lumenContour[point]);


	for (x = limInfX; x < (limSupX - 1); x++)
	{
		for (int y = limInfY; y < (limSupY - 1); y++)
		{
			bool inWall = false;
			bool outCalc = true;
			bool inLumen = false;
			unsigned short *refValue = (unsigned short *) image->GetScalarPointer(x,y,0);
			*refValue = 0;

			for (int numConts = 0; numConts < polygons.size(); numConts++)
			{
				bool adentro = pointInPolygon(polygons[numConts],x,y);
				if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
				{
					if (numConts == 0)
					{
						inWall = true;
					}
					else if (numConts == polygons.size() - 1)
					{
						inLumen = true;
						
					}
					else
					{
						outCalc = false;
					}
				}
			}

			double intensidad = imagedata->GetScalarComponentAsDouble(x,y,0,0);
			
			if (inWall && outCalc && inLumen)
			{

				if (imagedata->GetScalarComponentAsDouble(x,y,0,0) > (average - 1.5*deviation)
					 && imagedata->GetScalarComponentAsDouble(x,y,0,0) < (average + 2.5*deviation))
				{
					*refValue = 255;		
				}
			
			}

		}
	}

	

	extractLumen(image,lumenContour[point],point);
	
//	fclose(archivo);
}

// ----------------------------------------------------------------------------
double marAxisCT::obtainContourArea(marIsocontour* contour)
{

	double area = 0.0;
	double x1,y1,x2,y2;
	
	for (int i = 0; i < contour->getSize();i++)
	{
		contour->getPoint(i,&x1,&y1);
		if (i < (contour->getSize() - 1)) 
		{
			contour->getPoint(i+1,&x2,&y2);
		}
		else 
		{
			contour->getPoint(0,&x2,&y2);
		}

		area += (x1 * y2);
		area -= (y1 * x2);
	}

	area = area / 2;

	if (area < 0)
		area = -area;

	return area;
}

// ----------------------------------------------------------------------------
void marAxisCT::adjustWall(int point, kVolume* vol)
{
	//int actualPercentage = ((marParameters *) getParameters())->getLumenPercentage();
	marIsocontour* cont;
	double actualArea, previousArea;
	std::vector <double > grad;
	marContourVO* wallCont = searchContour(marAxisContours::WALL,point);

	//Obtain first area
	((marParameters *) getParameters())->setLumenPercentage(50);
	updateLumenPercentage(point,vol);
	cont = parsePolyDataToMarIsocontour(wallCont);
	previousArea = obtainContourArea(cont);

	for (int eval = 49; eval >= 0; eval--)
	{
		((marParameters *) getParameters())->setLumenPercentage(eval);
		updateLumenPercentage(point,vol);
		cont = parsePolyDataToMarIsocontour(wallCont);
		actualArea = obtainContourArea(cont);
		grad.push_back(fabs(actualArea - previousArea));
	}

	int newPercentage = 50 - maxValue(grad);
	((marParameters *) getParameters())->setLumenPercentage(newPercentage);
	updateLumenPercentage(point,vol);

	

}

// ----------------------------------------------------------------------------
void marAxisCT::adjustCalcification(int point, kVolume* vol)
{
}

// ----------------------------------------------------------------------------
int	marAxisCT::maxValue(std::vector <double> list)
{
	int max = 0;
	double maxVal = list[0];

	for (int i = 0; i < list.size(); i++)
	{
		if (list[i] > maxVal)
		{
			maxVal = list[i];
			max = i;
		}
	}

	return max;
}


// ----------------------------------------------------------------------------
int marAxisCT::getDiameter(marContourVO* contourVO, double x, double y, int limInfX, int limInfY, int limSupX, int limSupY)
{
	bool stop = false; 
	double originX = x;
	double originY = y;

	int coordBase = 1, diameterOne = 0, diameterTwo = 0;
	marIsocontour* cont = parsePolyDataToMarIsocontour(contourVO);
	generateVector(1,&coordBase,originX,originY,&x,&y);

	while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
	{
		generateVector(0,&coordBase,originX,originY,&x,&y);
		if (pointInPolygon(cont,x,y))
		{
			diameterOne++;
		}
		else
		{
			stop = true; 
		}
	}

	stop = false; 
	x = originX;
	y = originY; 
	while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
	{
		generateVector(35,&coordBase,originX,originY,&x,&y);
		if (pointInPolygon(cont,x,y))
		{
			diameterOne++;
		}
		else
		{
			stop = true; 
		}
	}

	stop = false; 
	x = originX;
	y = originY; 
	while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
	{
		generateVector(17,&coordBase,originX,originY,&x,&y);
		if (pointInPolygon(cont,x,y))
		{
			diameterTwo++;
		}
		else
		{
			stop = true; 
		}
	}

	stop = false; 
	x = originX;
	y = originY; 
	while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
	{
		generateVector(107,&coordBase,originX,originY,&x,&y);
		if (pointInPolygon(cont,x,y))
		{
			diameterTwo++;
		}
		else
		{
			stop = true; 
		}
	}

	delete cont;

	return (diameterOne + diameterTwo) / 2;
}

// ----------------------------------------------------------------------------
void marAxisCT::histogram(int point, kVolume* vol)
{

	vtkImageData* imagedata;
	int limInfX, limInfY, limInfZ;
	int limSupX, limSupY, limSupZ;
	int pos;
	
	
	FILE *archivo;
	char nomArchivo[30],temp[3];
	marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
	
		
//	for (int point = 0; point <  quantContours.size(); point++)
//	{
		imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
		imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
		((marParameters *) getParameters())->setLumenPercentage(86);
		updateLumenPercentage(point,vol);
		marIsocontour* contExt = parsePolyDataToMarIsocontour(wallCont);

		std::vector<double> intensity;
		std::vector<int> frecuencies;
	    
		for (int i = 0; i < limSupX; i++)
		{
			for (int j = 0; j < limSupY; j++)
			{
				if(pointInPolygon(contExt,i,j))
				{
					
					pos = searchData(intensity,imagedata->GetScalarComponentAsDouble(i,j,0,0));
					if (pos == -1)
					{
						intensity.push_back( imagedata->GetScalarComponentAsDouble(i,j,0,0));
						frecuencies.push_back(1);
					
						
					}
					else
					{
						frecuencies[pos]++;
					}
				}
				
			}
		}
			
		strcpy(nomArchivo,"./histogram");

// EED 06/06/2007 linux
//		itoa(point,temp,10);
		sprintf(temp,"%d",point);

		strcat(nomArchivo,temp);
		strcat(nomArchivo,".csv");
		
		

		archivo = fopen(nomArchivo,"w");
		for (int k = 0; k < frecuencies.size(); k++)
		{

			fprintf(archivo,"%f;%d\n",intensity[k],frecuencies[k]);

		}
		fclose(archivo);
		
	
		

		//imagedata->Delete();
		delete contExt;
//	}

}

// ----------------------------------------------------------------------------
int	marAxisCT::searchData(std::vector<double> vec, double value)
{
	int resp = -1;
	bool found = false;

	for (int i = 0; i< vec.size() && !found; i++)
	{
		if (vec[i] == value)
		{
			resp = i;
			found = true;
		}
	}

	return resp;
}

// ----------------------------------------------------------------------------
void marAxisCT::setCalibration(bool calib)
{
	calibration = calib;
}

// ----------------------------------------------------------------------------
bool marAxisCT::getCalibration()
{
	return calibration;
}

// ----------------------------------------------------------------------------
void marAxisCT::adjustContour(int point, kVolume* vol)
{

	double percentage = 0.05;
	double prevDifference = 0;
	if (quantContours[point + 1] == NULL)
	{
		create2Dcontours(point + 1 ,vol);
	}

	//1. Obtain actual signal 
	double signal = maxSignal;
	marContourVO* wallCont = searchContour(marAxisContours::WALL,point);

	//2. Obtain previous contour
	marIsocontour* prev = parsePolyDataToMarIsocontour(searchContour(marAxisContours::WALL, point + 1));
	
	//3. Adjust actual contour
	double newValue = signal; 
	wallCont->getIsocontour()->SetValue(0,newValue);
	wallCont->getIsocontour()->Update();
	wallCont->get2DContour()->Update();

	//4. Obtain points
	marIsocontour* actual = parsePolyDataToMarIsocontour(wallCont);

	//5. Compare areas
	double prevA = obtainContourArea(prev);
	double actuA = obtainContourArea(actual);

	prevDifference = fabs(prevA - actuA);


	if (prevA - actuA > 0 && prevA - actuA < percentage*prevA) //CASO 1
	{
		bool stop = false;
		while (!stop && newValue > 0)
		{
			newValue--;
			
			wallCont->getIsocontour()->SetValue(0,newValue);
			wallCont->getIsocontour()->Update();
			wallCont->get2DContour()->Update();

			actual = parsePolyDataToMarIsocontour(wallCont);
			actuA = obtainContourArea(actual);
 
			if (fabs(prevA - actuA) > percentage*prevA)
			{
				newValue++;
				wallCont->getIsocontour()->SetValue(0,newValue);
				wallCont->getIsocontour()->Update();
				wallCont->get2DContour()->Update();
				stop = true;
			}
		}
	
	}
	else if (prevA - actuA < 0 && prevA - actuA > -percentage*prevA) //CASO 2
	{

		bool stop = false;
		while (!stop && newValue > 0)
		{
			newValue--;
			
			wallCont->getIsocontour()->SetValue(0,newValue);
			wallCont->getIsocontour()->Update();
			wallCont->get2DContour()->Update();

			actual = parsePolyDataToMarIsocontour(wallCont);
			actuA = obtainContourArea(actual);
 
			if (prevA - actuA < -percentage*prevA)
			{
				newValue++;
				wallCont->getIsocontour()->SetValue(0,newValue);
				wallCont->getIsocontour()->Update();
				wallCont->get2DContour()->Update();
				stop = true;
			}
		}
	}
	else if (prevA - actuA < 0 && prevA - actuA < -percentage*prevA) //CASO 3
	{

		bool stop = false;
		
		while (!stop)
		{
			newValue++;
			
			wallCont->getIsocontour()->SetValue(0,newValue);
			wallCont->getIsocontour()->Update();
			wallCont->get2DContour()->Update();

			actual = parsePolyDataToMarIsocontour(wallCont);
			
			actuA = obtainContourArea(actual);
			if (prevA - actuA > -percentage*prevA)
			{
				if (prevDifference < fabs(prevA - actuA))
				{
					newValue--;
					wallCont->getIsocontour()->SetValue(0,newValue);
					wallCont->getIsocontour()->Update();
					wallCont->get2DContour()->Update();
				}
				
				stop = true;
				
			}
			prevDifference = fabs(prevA - actuA);
		}
	}
	else if (prevA - actuA > 0 && prevA - actuA > percentage*prevA) //CASO 4
	{
		bool stop = false;
		while (!stop && newValue > 0)
		{
			newValue--;
			
			wallCont->getIsocontour()->SetValue(0,newValue);
			wallCont->getIsocontour()->Update();
			wallCont->get2DContour()->Update();

			actual = parsePolyDataToMarIsocontour(wallCont);
			actuA = obtainContourArea(actual);
 
			if (prevA - actuA < percentage*prevA)
			{
				if (prevDifference < fabs(prevA - actuA))
				{
					newValue++;
					wallCont->getIsocontour()->SetValue(0,newValue);
					wallCont->getIsocontour()->Update();
					wallCont->get2DContour()->Update();
				}
				
				stop = true;
				
			}
			prevDifference = fabs(prevA - actuA);
		}
	}

	
	maxSignal = newValue;

	this->markUpLumen(point,vol);
}

// ----------------------------------------------------------------------------
int marAxisCT::getStartIndex()
{
	return startIndex;
}

// ----------------------------------------------------------------------------
void marAxisCT::setStartIndex(int start)
{
	startIndex = start;

	int per = ((marParameters *) getParameters())->getLumenPercentage();
	double oldValue = quantContours[start]->getContour(0)->getSignal()*per*pow(10.0,-2.0);

	maxSignal = oldValue; // / (per*pow(10.0,-2.0));
}

// ----------------------------------------------------------------------------
int	marAxisCT::getPointSize()
{
	return vesselPoints.size();
}

// ----------------------------------------------------------------------------
marPoint* marAxisCT::getPoint(int i)
{
	return vesselPoints[i];
}

// ----------------------------------------------------------------------------
void marAxisCT::extractLumen(vtkImageData *lumenImage, marIsocontour *lumenContour, int point)
{

	
	vtkImageCast* cast = vtkImageCast::New();
	cast->SetInput(lumenImage);
	cast->SetOutputScalarTypeToDouble();
	
	vtkImageGaussianSmooth* filter = vtkImageGaussianSmooth::New();
	filter->SetDimensionality(2);
	filter->SetInput(cast->GetOutput());
	filter->SetStandardDeviations((getParameters())->getStandardDeviation(),(getParameters())->getStandardDeviation());
	filter->SetRadiusFactors((getParameters())->getRadius(),(getParameters())->getRadius());
	filter->SetStandardDeviations(3,3);
	filter->SetRadiusFactors(3,3);
	


	marAxisContours* axisCont = quantContours[point];

	int pos = -1;
	
	for (int i = 0; i < axisCont->getSize()  && pos == -1; i++)
	{
		if (axisCont->getContourType(i) == marAxisContours::ELUMEN)
		{
			pos = i;
		}
	}

	bool stop = false;
	if (pos != -1)		//Verify that contour hasn't been replaced
	{
		stop = true;
		if (axisCont->isReplaced(pos))
		{
			stop = true;
		}
	
	}

	if (!stop)
	{
		marContourVO* contourVO = new marContourVO();
		marIsocontour* iso = lumenContour;
		
		contourVO->setIsocontour(vtkContourFilter::New());
		contourVO->getIsocontour()->SetInput( filter->GetOutput() );
		contourVO->getIsocontour()->SetNumberOfContours( 1 );
		contourVO->getIsocontour()->SetValue( 0, 128 );	
		contourVO->getIsocontour()->Update( );
			
		contourVO->setIsocontourCpd(vtkCleanPolyData::New());
		contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
		contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
		contourVO->getIsocontourCpd()->Update( );

		double cgx, cgy;

		iso->getCG(&cgx,&cgy);
			
		contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
		contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
			
		contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
		contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
		contourVO->getIsocontourDcf()->Update( );
			
		contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
		contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
		contourVO->getIsocontourCpd2()->Update();

		contourVO->setIsocontourStripped(vtkStripper::New( ));
		contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
		contourVO->getIsocontourStripped()->Update();
		
		contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
		contourVO->get2DContour()->Update();
		contourVO->setSignal(255);
		contourVO->setType(marAxisContours::ELUMEN); 
		
		
		
		marIsocontour* actual = parsePolyDataToMarIsocontour(contourVO);
		marContourVO* aVO = searchContour(marAxisContours::WALL,point);
		double wallArea = 0;
		if (aVO != NULL)
		{
			wallArea = obtainContourArea(parsePolyDataToMarIsocontour(aVO));
		}
	
		double prevA = obtainContourArea(actual);
		
		int newValue = 128;
		int finalValue = 128;
		double actuA = 1;
		while (newValue > 0 && actuA > 0)
		{
			
			newValue--;

			try
			{
				contourVO->getIsocontour()->SetValue(0,newValue);
				contourVO->getIsocontour()->Update();
				contourVO->get2DContour()->Update();
				actual = parsePolyDataToMarIsocontour(contourVO);
				actuA = obtainContourArea(actual);
				if ((actuA /1000) < 1 && (prevA / 1000) > 1)
				{
					finalValue = newValue;
					prevA = actuA;
				}
				else if (actuA > prevA && (prevA / 1000) < 1 && actuA < 0.1*wallArea)
				{
					finalValue = newValue;
					prevA = actuA;
				}
			}
			catch (std::exception e)
			{
				actuA = 0;
			}
			
		}

		contourVO->getIsocontour()->SetValue(0,finalValue);
		contourVO->getIsocontour()->Update();
		contourVO->get2DContour()->Update();



		/*iso->getType()*/

		if (pos != -1)
		{
			axisCont->replaceContour(contourVO,pos);
		}
		else
		{
			axisCont->addContour(contourVO);
		}
		//	axisCont->addContour(contourVO);


		quantContours[point] = axisCont;
	}

	
}

// ----------------------------------------------------------------------------
void marAxisCT::generateFile(int point,kVolume* vol)
{
	FILE *archivo;
	char nomArchivo[30], temp[3];

	strcpy(nomArchivo,"./point");

// EED 06/06/2007 linux
//	itoa(point,temp,10);
	sprintf(temp,"%d",point);

	strcat(nomArchivo,temp);
	strcat(nomArchivo,".txt");
	archivo = fopen(nomArchivo,"w");

	fprintf(archivo, "%d", quantContours[point]->getSize());

	int i;
	for (i = 0; i < quantContours[point]->getSize(); i++)
	{
		marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
		fprintf(archivo, " %d",iso->getSize());
		
	}
	fprintf(archivo,"\n");

	for ( i = 0; i < quantContours[point]->getSize(); i++)
	{
		marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
		fprintf(archivo,"Tipo: %d\n",quantContours[point]->getContourType(i));
		for (int j = 0; j < iso->getSize(); j++)
		{
			double tempX, tempY;
			iso->getPoint(j,&tempX,&tempY); 
			fprintf(archivo,"%f;%f;%d\n",tempX,tempY,point);
		}
	}

// EED 30 Janvier 2007
	vtkImageData *vtkimagedata = this->getSliceImage(point,vol);

	std::string directory = "./";
	std::string filename  = "xx";
	float rescalaSlope           =  1;
	float rescalaIntercept       =  0;
	int dim[3];
	vtkimagedata->GetDimensions(dim);
	int voi[6];
	voi[0]=0;
	voi[1]=dim[0];
	voi[2]=0;
	voi[3]=dim[1];
	voi[4]=0;
	voi[5]=dim[2];
	marRAW2Files marraw2;
	marraw2.saveVolume(directory,filename,vtkimagedata,voi,rescalaSlope,rescalaIntercept);



	fclose(archivo);

}

// ----------------------------------------------------------------------------
void marAxisCT::replaceContour2D(int point,int size,double *vx,double *vy, int type)
{

	vtkPoints *_pts = vtkPoints::New();
	_pts->SetNumberOfPoints(size);
	int j;

	for (j=0 ; j<size ; j++){
		_pts->SetPoint(j,	vx[j]	, vy[j]	, 0 );
	}
//	_pts->SetPoint(0,	vx[0]	, vy[0]	, 0 );

	vtkCellArray *lines = vtkCellArray::New();
	lines->InsertNextCell( size );
	for ( j=0 ; j<size+1 ; j++ ){
		lines->InsertCellPoint(j % size );
	}

	vtkPolyData *_pd = vtkPolyData::New();
	_pd->SetPoints( _pts );
	_pd->SetLines( lines );
	lines->Delete();  //do not delete lines ??
	_pts->Delete();  

	marContourVO* vo = new marContourVO();///this->searchContour(type,point); 
	vo->setReplaced(true);
	vo->set2DContour(_pd);
	vo->setType(type);
	quantContours[point]->addContour(vo);   //->replaceContour(vo,i);
/*	_2Dcontours[i]=_pd;
	createContour(i,NULL);*/
}

// ----------------------------------------------------------------------------
void marAxisCT::cleanContours(int type, int point)
{
	marAxisContours* newContours = new marAxisContours();
	if (quantContours[point] != NULL)
	{
		for (int i = 0; i < quantContours[point]->getSize(); i++)
		{
			if (quantContours[point]->getContourType(i) != type)
			{
				newContours->addContour(quantContours[point]->getContour(i));
			}
			else
			{
			/*	marContourVO* vo = quantContours[point]->getContour(i);
				vo->set2DContour(NULL);
				newContours->replaceContour(vo,i);*/
			}
		}
	}

	quantContours[point] = newContours;
}

// ----------------------------------------------------------------------------
marContourVO*   marAxisCT::searchContour(int type, int point)
{
	marContourVO* cont = NULL;

	for (int i = 0; i < quantContours[point]->getSize(); i++)
	{
		if (quantContours[point]->getContourType(i) == type)
		{
			cont = quantContours[point]->getContour(i);
		}
	}

	return cont;
}

// ----------------------------------------------------------------------------
double marAxisCT::performXOR(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
{
	if (manual.size() == 0)
	{
		return 0;
	}

	int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
	std::vector <marIsocontour* > polygons;
	vtkImageData* imagedata;

	imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
	imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	vtkImageData* imageA = vtkImageData::New();
	imageA->SetDimensions(limSupX,limSupY,0);
	imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	imageA->SetScalarTypeToUnsignedShort();
	imageA->AllocateScalars();
	imageA->Update();

	vtkImageData* imageM = vtkImageData::New();
	imageM->SetDimensions(limSupX,limSupY,0);
	imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	imageM->SetScalarTypeToUnsignedShort();
	imageM->AllocateScalars();
	imageM->Update();

	
	for (int i = 0; i < quantContours[point]->getSize(); i++)
	{
		if (quantContours[point]->getContourType(i) == type)
		{
			polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
		}
		
	}

	if (polygons.size() == 0)
	{
		return 0;
	}
	int x,y;
	for (x = limInfX; x < (limSupX - 1); x++)
	{
		for (y = limInfY; y < (limSupY - 1); y++)
		{
			unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
			unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
			*refValueA = 0;
			*refValueM = 0;

			int numConts;
			for (numConts = 0; numConts < polygons.size(); numConts++)
			{
				if (pointInPolygon(polygons[numConts],x,y))
				{
						*refValueA = 255;
				}
			}

			for (numConts = 0; numConts < manual.size(); numConts++)
			{
				if (pointInPolygon(manual[numConts],x,y))
				{
						*refValueM = 255;
				}
			}
		}
	}

/*	vtkImageLogic* xor = vtkImageLogic::New();
	xor->SetInput1(imageM);
	xor->SetInput2(imageA);
	xor->SetOutputTrueValue(255);
	xor->SetOperationToXor();

	vtkImageCast* cast = vtkImageCast::New();
	cast->SetInput(xor->GetOutput());
	cast->SetOutputScalarTypeToDouble();
*/
	int coincidencias = 0;
	for (x = limInfX; x < (limSupX - 1); x++)
	{
		for (y = limInfY; y < (limSupY - 1); y++)
		{
			int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
			int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);

			int xor_ = compA^compM;
			if (xor_ == 255)
			{
				coincidencias++;
			}
		}
	}

	double resp = (double)coincidencias ;

	return resp;

}

// ----------------------------------------------------------------------------
double  marAxisCT::performAND(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
{

	if (manual.size() == 0)
	{
		return 0;
	}

	int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
	std::vector <marIsocontour* > polygons;
	vtkImageData* imagedata;

	imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
	imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	vtkImageData* imageA = vtkImageData::New();
	imageA->SetDimensions(limSupX,limSupY,0);
	imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	imageA->SetScalarTypeToUnsignedShort();
	imageA->AllocateScalars();
	imageA->Update();

	vtkImageData* imageM = vtkImageData::New();
	imageM->SetDimensions(limSupX,limSupY,0);
	imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	imageM->SetScalarTypeToUnsignedShort();
	imageM->AllocateScalars();
	imageM->Update();

	for (int i = 0; i < quantContours[point]->getSize(); i++)
	{
		if (quantContours[point]->getContourType(i) == type)
		{
			polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
		}
		
	}

	if (polygons.size() == 0)
	{
		return 0;
	}

	int x,y;
	for (x = limInfX; x < (limSupX - 1); x++)
	{
		for (y = limInfY; y < (limSupY - 1); y++)
		{
			unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
			unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
			*refValueA = 0;
			*refValueM = 0;

			int numConts;
			for (numConts = 0; numConts < polygons.size(); numConts++)
			{
				if (pointInPolygon(polygons[numConts],x,y))
				{
						*refValueA = 255;
				}
			}

			for (numConts = 0; numConts < manual.size(); numConts++)
			{
				if (pointInPolygon(manual[numConts],x,y))
				{
						*refValueM = 255;
				}
			}
		}
	}

/*	vtkImageLogic* and = vtkImageLogic::New();
	and->SetInput1(imageM);
	and->SetInput2(imageA);
	and->SetOutputTrueValue(255);
	and->SetOperationToAnd();
*/
	int coincidencias = 0;
	for (x = limInfX; x < (limSupX - 1); x++)
	{
		for (y = limInfY; y < (limSupY - 1); y++)
		{
			int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
			int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);

			int and_ = compA & compM;
			if (and_ == 255)
			{
				coincidencias++;
			}
		}
	}

	double resp = (double)coincidencias ;


	return resp;
}

// ----------------------------------------------------------------------------
marIsocontour*	marAxisCT::loadMarIsocontour(int size, double *vx, double *vy)
{
	marIsocontour* cont = new marIsocontour();
		
	for ( int j=0 ; j<size ; j++ )
	{
		cont->insertPoint(vx[j], vy[j]);
	}

	return cont;
}

// ----------------------------------------------------------------------------
double	marAxisCT::performUnion(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
{
		if (manual.size() == 0)
	{
		return 0;
	}

	int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
	std::vector <marIsocontour* > polygons;
	vtkImageData* imagedata;

	imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
	imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	vtkImageData* imageA = vtkImageData::New();
	imageA->SetDimensions(limSupX,limSupY,0);
	imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	imageA->SetScalarTypeToUnsignedShort();
	imageA->AllocateScalars();
	imageA->Update();

	vtkImageData* imageM = vtkImageData::New();
	imageM->SetDimensions(limSupX,limSupY,0);
	imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
	imageM->SetScalarTypeToUnsignedShort();
	imageM->AllocateScalars();
	imageM->Update();

	int i;
	for (i = 0; i < quantContours[point]->getSize(); i++)
	{
		if (quantContours[point]->getContourType(i) == type)
		{
			polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
		}
		
	}

	if (polygons.size() == 0)
	{
		return 0;
	}

	int x,y;
	for ( x = limInfX; x < (limSupX - 1); x++)
	{
		for ( y = limInfY; y < (limSupY - 1); y++)
		{
			unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
			unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
			*refValueA = 0;
			*refValueM = 0;
			int  numConts; 
			for ( numConts = 0; numConts < polygons.size(); numConts++)
			{
				if (pointInPolygon(polygons[numConts],x,y))
				{
						*refValueA = 255;
				}
			}

			for (numConts = 0; numConts < manual.size(); numConts++)
			{
				if (pointInPolygon(manual[numConts],x,y))
				{
						*refValueM = 255;
				}
			}
		}
	}

/*	vtkImageLogic* and = vtkImageLogic::New();
	and->SetInput1(imageM);
	and->SetInput2(imageA);
	and->SetOutputTrueValue(255);
	and->SetOperationToAnd();
*/
	int coincidencias = 0;
	for (x = limInfX; x < (limSupX - 1); x++)
	{
		for ( y = limInfY; y < (limSupY - 1); y++)
		{
			int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
			int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);

			int or_ = compA | compM;
			if (or_ == 255)
			{
				coincidencias++;
			}
		}
	}

	double resp = (double)coincidencias;


	return resp;
}