1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
9 # This software is governed by the CeCILL-B license under French law and
10 # abiding by the rules of distribution of free software. You can use,
11 # modify and/ or redistribute the software under the terms of the CeCILL-B
12 # license as circulated by CEA, CNRS and INRIA at the following URL
13 # http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
14 # or in the file LICENSE.txt.
16 # As a counterpart to the access to the source code and rights to copy,
17 # modify and redistribute granted by the license, users are provided only
18 # with a limited warranty and the software's author, the holder of the
19 # economic rights, and the successive licensors have only limited
22 # The fact that you are presently reading this means that you have had
23 # knowledge of the CeCILL-B license and that you accept its terms.
24 # ------------------------------------------------------------------------ */
26 // marAxisCT.cpp: implementation of the marAxisCT class.
28 //////////////////////////////////////////////////////////////////////
31 #include <vtkImageCast.h>
32 #include <vtkImageGaussianSmooth.h>
33 #include <vtkCleanPolyData.h>
34 #include <vtkPolyDataConnectivityFilter.h>
35 #include <vtkStripper.h>
37 #include <vtkIdList.h>
38 #include <vtkMatrix4x4.h>
39 #include <vtkTransform.h>
40 #include <vtkPlaneSource.h>
41 #include <vtkContourFilter.h>
44 #include <vtkImageContinuousErode3D.h>
45 #include <vtkImageThreshold.h>
46 #include <vtkImageSobel2D.h>
47 #include <vtkCellArray.h>
48 #include <vtkImageLogic.h>
51 #include "marAxisCT.h"
52 #include "marGdcmDicom.h"
57 marAxisCT::marAxisCT( ) : marAxis( )
61 // ----------------------------------------------------------------------------
62 marContour* marAxisCT::getContour( int point , kVolume* vol, int index ) {
63 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->getContour() == NULL)
65 createContours(point,vol);
68 return quantContours[point]->getContour(index)->getContour();
71 // ----------------------------------------------------------------------------
72 vtkPoints* marAxisCT::get3Dcontour( int point , kVolume* vol, int index ) {
73 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get3DContour() == NULL)
75 create3Dcontours(point,vol);
78 marAxisContours* cont = quantContours[point];
80 return cont->getContour(index)->get3DContour();
83 // ----------------------------------------------------------------------------
84 vtkPolyData* marAxisCT::get2Dcontour( int point , kVolume* vol, int index ) {
85 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DContour() == NULL)
87 create2Dcontours(point,vol);
90 marAxisContours* cont = quantContours[point];
94 return cont->getContour(index)->get2DContour();
97 // ----------------------------------------------------------------------------
98 vtkPoints* marAxisCT::get2DDiameterMin( int point , kVolume* vol, int index ) {
99 // if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DDiameterMin() == NULL)
101 create2DDiametersMin(point,vol);
104 marAxisContours* cont = quantContours[point];
106 return cont->getContour(index)->get2DDiameterMin();
109 // ----------------------------------------------------------------------------
110 vtkPoints* marAxisCT::get2DDiameterMax( int point , kVolume* vol, int index ) {
111 // if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DDiameterMax() == NULL)
113 create2DDiametersMax(point,vol);
116 marAxisContours* cont = quantContours[point];
118 return cont->getContour(index)->get2DDiameterMax();
121 // ----------------------------------------------------------------------------
122 int marAxisCT::getSignal(int point, int index, kVolume* vol) {
123 if (quantContours[point] == NULL )
125 createSignals(point,vol);
128 return (int) ( quantContours[point]->getContour(index)->getSignal() );
131 // ----------------------------------------------------------------------------
132 void marAxisCT::createContours( int point, kVolume* vol ) {
135 int numberOfPoints,numberOfCells;
139 vtkPolyData* vtkPolyData_2DContour;
142 int sizeIma = getParameters( )->getSizeIma( );
143 double dimIma = getParameters( )->getDimIma( );
145 if (quantContours.size() < point - 1 )
147 create2Dcontours(point,vol);
150 marAxisContours* axisCont = quantContours[point];
152 for (int i = 0; i < axisCont->getSize(); i++)
154 marContourVO* aContourVO = axisCont->getContour(i);
155 aContourVO->setContour(new marContour( point, getParameters( ) ));
156 vtkPolyData_2DContour = aContourVO->get2DContour();
157 numberOfCells = vtkPolyData_2DContour->GetNumberOfCells( );
158 cell = vtkPolyData_2DContour->GetCell( 0 );
159 pids = cell->GetPointIds( );
160 numberOfPoints = pids->GetNumberOfIds( );
161 for( j = 0; j < numberOfPoints; j++ )
163 id = pids->GetId( j );
164 vtkPolyData_2DContour->GetPoint( id, p);
167 x=x * dimIma / ( float ) sizeIma;
168 y=y * dimIma / ( float ) sizeIma;
169 aContourVO->getContour()->addContourPoint( x , y );
172 aContourVO->getContour()->do_spline();
173 aContourVO->getContour()->calculateVariables();
178 // ----------------------------------------------------------------------------
179 void marAxisCT::create3Dcontours( int point , kVolume* vol ) {
181 if (quantContours[point] == NULL )
183 createContours(point,vol);
186 marAxisContours* axisCont = quantContours[point];
188 for (int i = 0; i < axisCont->getSize(); i++)
190 vtkPoints *_vtkPoints;
196 double c[3],p1[3],p2[3];
199 _vtkPoints = vtkPoints::New();
200 o = getPoints(point); //PENDIENTE _points[i];
202 mat = vtkMatrix4x4::New();
203 trans = vtkTransform::New();
206 nn[0]=n[0]; nn[1]=n[1]; nn[2]=n[2];
207 nC=sqrt( nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2] );
208 nn[0]=nn[0]/nC; nn[1]=nn[1]/nC; nn[2]=nn[2]/nC;
210 vtkPlaneSource* pSource = vtkPlaneSource::New( );
211 pSource->SetOrigin( 0, 0 , 0 );
212 pSource->SetPoint1( 1, 0 , 0 );
213 pSource->SetPoint2( 0, 0 , 1.0 );
214 pSource->SetCenter( 0,0,0 );
215 pSource->SetNormal( nn[ 0 ], nn[ 1 ], nn[ 2 ] );
218 pSource->GetOrigin(c);
219 pSource->GetPoint1(p1);
220 pSource->GetPoint2(p2);
222 p1[0]=p1[0]-c[0]; p1[1]=p1[1]-c[1]; p1[2]=p1[2]-c[2];
223 p2[0]=p2[0]-c[0]; p2[1]=p2[1]-c[1]; p2[2]=p2[2]-c[2];
224 nA=sqrt( p1[0]*p1[0] + p1[1]*p1[1] + p1[2]*p1[2] );
225 nB=sqrt( p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2] );
226 p1[0]=p1[0]/nA; p1[1]=p1[1]/nA; p1[2]=p1[2]/nA;
227 p2[0]=p2[0]/nB; p2[1]=p2[1]/nB; p2[2]=p2[2]/nB;
229 mat->SetElement (0,0, nn[0]);
230 mat->SetElement (1,0, nn[1]);
231 mat->SetElement (2,0, nn[2]);
232 mat->SetElement (3,0, 0);
233 mat->SetElement (0,1, p2[0]);
234 mat->SetElement (1,1, p2[1]);
235 mat->SetElement (2,1, p2[2]);
236 mat->SetElement (3,1, 0);
237 mat->SetElement (0,2, p1[0]);
238 mat->SetElement (1,2, p1[1]);
239 mat->SetElement (2,2, p1[2]);
240 mat->SetElement (3,2, 0);
241 mat->SetElement (0,3, 0);
242 mat->SetElement (1,3, 0);
243 mat->SetElement (2,3, 0);
244 mat->SetElement (3,3, 1);
246 double deter=mat->Determinant(mat);
248 trans->SetMatrix(mat);
251 trans->RotateWXYZ ( ang,0,1,0);
254 vtkMatrix4x4 *m=vtkMatrix4x4::New();
255 trans->GetMatrix ( m );
257 int j,numberOfPoints;
259 marContour* contour2D = getContour(point,vol,i);
260 marContourVO* aContourVO = axisCont->getContour(i);
261 numberOfPoints = contour2D->getNumberOfPoints( );
265 for( j = 0; j <= numberOfPoints; j++ ) {
266 contour2D->getPoint( fpA,j % numberOfPoints);
271 ppp[0] = m->GetElement(0,0)*pA[0] + m->GetElement(0,1)*pA[1] + m->GetElement(0,2)*pA[2] + m->GetElement(0,3);
272 ppp[1] = m->GetElement(1,0)*pA[0] + m->GetElement(1,1)*pA[1] + m->GetElement(1,2)*pA[2] + m->GetElement(1,3);
273 ppp[2] = m->GetElement(2,0)*pA[0] + m->GetElement(2,1)*pA[1] + m->GetElement(2,2)*pA[2] + m->GetElement(2,3);
275 ppp[0]=ppp[0]+o[0]-c[0]; ppp[1]=ppp[1]+o[1]-c[1]; ppp[2]=ppp[2]+o[2]-c[2];
276 _vtkPoints->InsertNextPoint( (float)ppp[0],(float)ppp[1],(float)ppp[2] );
286 if (aContourVO->get3DContour()!=NULL){ aContourVO->get3DContour()->Delete(); }
287 aContourVO->set3DContour(_vtkPoints);
296 // ----------------------------------------------------------------------------
297 void marAxisCT::create2Dcontours( int point , kVolume* vol ) {
298 std::vector <marIsocontour *> contours;
302 generatePoints(point,vol,&contours);
304 vtkImageData* imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
305 marAxisContours* axisCont = new marAxisContours();
307 for (int i = 0; i < contours.size(); i++)
309 marContourVO* contourVO = new marContourVO();
310 marIsocontour* iso = contours[i];
312 contourVO->setIsocontour(vtkContourFilter::New());
313 contourVO->getIsocontour()->SetInput( imagedata );
314 contourVO->getIsocontour()->SetNumberOfContours( 1 );
315 contourVO->getIsocontour()->SetValue( 0, iso->getMaxIntensity() );
316 contourVO->getIsocontour()->Update( );
318 contourVO->setIsocontourCpd(vtkCleanPolyData::New());
319 contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
320 contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
321 contourVO->getIsocontourCpd()->Update( );
325 iso->getCG(&cgx,&cgy);
327 contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
328 contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
330 contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
331 contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
332 contourVO->getIsocontourDcf()->Update( );
334 contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
335 contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
336 contourVO->getIsocontourCpd2()->Update();
338 contourVO->setIsocontourStripped(vtkStripper::New( ));
339 contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
340 contourVO->getIsocontourStripped()->Update();
342 contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
343 contourVO->get2DContour()->Update();
344 contourVO->setSignal(iso->getMaxIntensity());
346 if (iso->getType() == marAxisContours::LUMEN)
348 contourVO->setType(marAxisContours::WALL);
352 contourVO->setType(iso->getType());
357 axisCont->addContour(contourVO);
366 quantContours[point] = axisCont;
370 updateLumenPercentage(point,vol);
372 updateCalcPercentage(point,vol);
373 // adjustWall(point,vol);
374 if (calibration && point < startIndex)
376 this->adjustContour(point,vol);
380 adjustWall(point,vol);
386 // marIsocontour* cont = parsePolyDataToMarIsocontour(quantContours[point]->getContour(0));
389 char nomArchivo[30], temp[3];
391 strcpy(nomArchivo,"./points");
393 strcat(nomArchivo,temp);
394 strcat(nomArchivo,".csv");
396 archivo = fopen(nomArchivo,"w");
398 for (int h = 0; h < cont->getSize(); h++)
401 cont->getPoint(h,&x,&y);
402 fprintf(archivo,"%f;%f\n",x,y);
408 // ----------------------------------------------------------------------------
409 void marAxisCT::create2DDiametersMin(int point , kVolume* vol ) {
411 if (quantContours[point] == NULL )
413 createContours(point,vol);
416 marAxisContours* axisCont = quantContours[point];
418 for (int i = 0; i < axisCont->getSize(); i++)
421 marContour *marcontour=getContour(point,vol,i);
422 marContourVO* aContourVO = axisCont->getContour(i);
424 marcontour->getMinimumLine(p1,p2);
425 vtkPoints *_vtkPoints = vtkPoints::New();
427 int sizeIma = getParameters( )->getSizeIma( );
428 double dimIma = getParameters( )->getDimIma( );
429 double factor=( float ) sizeIma / dimIma;
431 p1[0]=p1[0]*factor+sizeIma/2;
432 p1[1]=p1[1]*factor+sizeIma/2;
433 p2[0]=p2[0]*factor+sizeIma/2;
434 p2[1]=p2[1]*factor+sizeIma/2;
435 _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
436 _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
438 if ( aContourVO->get2DDiameterMin()!=NULL ) { aContourVO->get2DDiameterMin()->Delete(); }
439 aContourVO->set2DDiameterMin(_vtkPoints);
445 // ----------------------------------------------------------------------------
446 void marAxisCT::create2DDiametersMax(int point , kVolume* vol ) {
447 if (quantContours[point] == NULL )
449 createContours(point,vol);
452 marAxisContours* axisCont = quantContours[point];
454 for (int i = 0; i < axisCont->getSize(); i++)
457 marContour *marcontour=getContour(point,vol,i);
458 marContourVO* aContourVO = axisCont->getContour(i);
460 marcontour->getMaximumLine(p1,p2);
461 vtkPoints *_vtkPoints = vtkPoints::New();
463 int sizeIma = getParameters( )->getSizeIma( );
464 double dimIma = getParameters( )->getDimIma( );
465 double factor=( float ) sizeIma / dimIma;
466 p1[0]=p1[0]*factor+sizeIma/2;
467 p1[1]=p1[1]*factor+sizeIma/2;
468 p2[0]=p2[0]*factor+sizeIma/2;
469 p2[1]=p2[1]*factor+sizeIma/2;
470 _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
471 _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
473 if ( aContourVO->get2DDiameterMax()!=NULL ) { aContourVO->get2DDiameterMax()->Delete(); }
474 aContourVO->set2DDiameterMax(_vtkPoints);
478 // ----------------------------------------------------------------------------
479 void marAxisCT::createSignals (int point, kVolume* vol) {
480 if (quantContours[point] == NULL )
482 create2Dcontours(point,vol); //AQUI
483 // filterSignal(vol);
487 // ----------------------------------------------------------------------------
488 int marAxisCT::getNumberOfContours(int point, kVolume* vol) {
489 if (quantContours[point] == NULL )
491 create2Dcontours(point,vol);//AQUI
492 // filterSignal(vol);
496 return quantContours[point]->getSize();
499 // ----------------------------------------------------------------------------
500 int marAxisCT::getContourType(int point, int index, kVolume* vol) {
501 if (quantContours[point] == NULL )
503 create2Dcontours(point,vol);
506 return quantContours[point]->getContourType(index);
510 // ----------------------------------------------------------------------------
511 void marAxisCT::generatePoints(int point, kVolume* vol, std::vector<marIsocontour*> *list) {
513 vtkImageData* imagedata;
516 double *nPoint = new double[ 3 ];
517 double x, y, prevIntensity, threshold;
518 bool found, estStart, estDivided, isNeighbor;
519 int actualZone, startZone, coordBase, maxPoint, lumenCount, calcCount;
520 std::vector <marPoint *> gradients;
521 std::vector <double> average;
522 std::vector <marIsocontour *> contoursLumen;
523 std::vector <marIsocontour *> contoursCalc;
525 contoursLumen.push_back(NULL);
526 contoursCalc.push_back(NULL);
530 //1. Check starting zone
531 startZone = marAxisContours::LUMEN;
532 actualZone = startZone;
536 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
538 int limInfX, limInfY, limInfZ;
539 int limSupX, limSupY, limSupZ;
540 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
542 vtkImageCast* cast = vtkImageCast::New();
543 cast->SetInput(imagedata);
544 cast->SetOutputScalarTypeToDouble();
546 vtkImageGaussianSmooth* filter = vtkImageGaussianSmooth::New();
547 filter->SetDimensionality(2);
548 filter->SetInput(cast->GetOutput());
549 filter->SetStandardDeviations((getParameters())->getStandardDeviation(),(getParameters())->getStandardDeviation());
550 filter->SetRadiusFactors((getParameters())->getRadius(),(getParameters())->getRadius());
551 filter->SetStandardDeviations(3,3);
552 filter->SetRadiusFactors(3,3);
554 //imagedata = filter->GetOutput();
557 nPoint = getPoints(point);
559 aPoint[ 0 ] = (limSupX - limInfX) / 2; //nPoint[0];
560 aPoint[ 1 ] = (limSupY - limInfY) / 2; //nPoint[1];
561 aPoint[ 2 ] = (limSupZ - limInfZ) / 2; // nPoint[2];
564 //3. Generate rays over all directions
565 for (int dir = 0; dir < dirs; dir++)
571 prevIntensity = interpolate(x,y,imagedata);
572 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
574 //2.1 Evaluate gradient in every position of the ray
575 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
579 double resp = interpolate(x,y,imagedata);
580 //2.1.2 Calcultate gradient
581 marPoint* p = new marPoint(x,y);
582 p->setGradient(prevIntensity - resp);
583 p->setDirection(dir);
584 p->setIntensity(resp);
586 prevIntensity = resp;
588 gradients.push_back(p);
590 //2.1.3 Increase ray for next evaluation
591 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
593 }// end while gradient
599 threshold = ((marParameters *) getParameters())->getContourThresh();
600 int mysize = gradients.size();
602 //2.2 Compare every gradient value with threshold
603 while (index < gradients.size() && !found)
605 //2.2.1 Evaluate Lumen
606 if (actualZone == marAxisContours::LUMEN)
608 average.push_back(gradients[index]->getIntensity());
610 //2.2.1.1 If value > threshold there is a contour
611 if (fabs(gradients[index]->getGradient()) > threshold
612 && fabs(gradients[index]->getGradient()) < 4000) {
613 //|| (calibration && (abs(gradients[index]->getIntensity()) > abs(avgValue - 2.5*stdValue)
614 //|| abs(gradients[index]->getIntensity()) > abs(avgValue + 2.5*stdValue)))) {
618 if(gradients[index]->getGradient() > 0)
619 maxPoint = getMaximumGrad(gradients,index,gradients.size(),threshold);
621 maxPoint = getMinimumGrad(gradients,index,gradients.size(),threshold);
624 x = gradients[maxPoint]->getX();
625 y = gradients[maxPoint]->getY();
627 //2.2.1.1.1 Assign lumen structure
628 isNeighbor = detectNeighbor(contoursLumen[lumenCount],dir,marAxisContours::LUMEN);
634 marIsocontour* cont = addPointToContour(contoursLumen[lumenCount],false,marAxisContours::LUMEN,gradients[maxPoint]);
635 contoursLumen[lumenCount] = cont;
637 //2.2.1.1.2 Check step lumen->calcification
638 double actualInt = interpolate(x,y,imagedata);
641 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&auxX,&auxY);
642 double next = interpolate(auxX,auxY,imagedata);
643 if (gradients[index]->getGradient() < 0 /*&& actualInt < next */&& (auxX > (limInfX+1) && auxX < (limSupX-1)) && (auxY > (limInfY+1) && auxY < (limSupY-1)))
645 //2.2.1.1.2.1 Assign calcification structure
646 isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
650 contoursCalc.push_back(NULL);
653 cont = addPointToContour(contoursCalc[calcCount],true,marAxisContours::CALCIFICATION,gradients[maxPoint]);
655 contoursCalc[calcCount] = cont;
657 actualZone = marAxisContours::CALCIFICATION;
660 //2.2.1.1.2.3 Verify divided structure
666 if ((dir == dirs - 1) && estStart)
671 }//2.2.1.1.3 Transition lumen->background
674 //2.2.1.1.3.1 Asign border structure
677 //2.2.1.1.3.2 Generate stop condition
679 actualZone = startZone;
683 threshold = ((marParameters *) getParameters())->getContourThresh();
688 //2.2.2 Evaluate calcification (MODIFIED IF THERE IS HYPODENSE)
693 double actualInt = interpolate(x,y,imagedata);
694 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
695 double next = interpolate(x,y,imagedata);
697 //2.2.2.1 If actual intensity is smaller than previous there is a contour
698 if (actualInt < next)
700 while (actualInt < interpolate(x,y,imagedata) && (x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
702 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
708 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
710 //2.2.2.2 Assign calcification structure
711 isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
715 contoursCalc.push_back(NULL);
718 //Adjust index to avoiod crash
719 if (index >= gradients.size())
721 calcIndex = gradients.size() - 1;
728 marIsocontour* cont = addPointToContour(contoursCalc[calcCount],false,marAxisContours::CALCIFICATION,gradients[calcIndex]);
729 contoursCalc[calcCount] = cont;
731 //2.2.2.3 Asign border structure
735 actualZone = startZone;
738 threshold = ((marParameters *) getParameters())->getContourThresh();
744 if (index == gradients.size() && !found)
746 threshold = threshold - 1;
750 }//end while evaluation
752 //TODO - SACAR A UTILS
753 for (int ind = 0; ind < gradients.size(); ind++)
755 marPoint* p = gradients[ind];
756 gradients[ind] = NULL;
762 }//end for directions
765 char nomArchivo[30]; //, temp[3]
767 strcpy(nomArchivo,"./debug");
768 // itoa(point,temp,10);
769 // strcat(nomArchivo,temp);
770 strcat(nomArchivo,".txt");
771 archivo = fopen(nomArchivo,"w");
777 // lumenContour[point] = contoursLumen[marAxisContours::LUMEN];
779 double tempXi, tempYi;
780 int polySides = contoursLumen[marAxisContours::LUMEN]->getSize();
783 marIsocontour* alc = new marIsocontour();
784 for (int m=0; m<polySides; m++) {
785 contoursLumen[marAxisContours::LUMEN]->getPoint(m,&tempXi,&tempYi);
786 alc = addPointToContour(lumenContour[point],false,marAxisContours::LUMEN,new marPoint(tempXi,tempYi));
787 lumenContour[point] = alc;
788 fprintf(archivo,"X:%f,Y:%f\n",tempXi,tempYi);
793 //4. Verify divided unified structure
794 if (estDivided && calcCount > 0)
796 unifyDividedStructure(&contoursCalc);
800 //5. Obtain statistics
801 double lumenInt = getStatistics(contoursLumen[marAxisContours::LUMEN],aPoint[ 0 ],aPoint[ 1 ], imagedata);
802 contoursLumen[marAxisContours::LUMEN]->setMaxIntensity(lumenInt);
804 if (contoursCalc[0] != NULL)
806 for (int i = 0; i < contoursCalc.size(); i++)
809 double maxInt = getCalcStatistics(contoursCalc[i],aPoint[ 0 ],aPoint[1],imagedata,i-1);
810 contoursCalc[i]->setMaxIntensity(maxInt);
811 contoursLumen.push_back(NULL);
812 contoursLumen[lumenCount] = contoursCalc[i];
817 *list = contoursLumen;
819 contoursCalc.clear();
820 contoursLumen.clear();
821 marUtils* util = new marUtils();
822 double avg = util->obtainAverage(average);
824 //archivo = fopen("./debug.txt","w");
825 fprintf(archivo,"tama�o: %d umbra: %d", (*list).size(),((marParameters *) getParameters())->getContourThresh());
826 for (int k = 0; k < contoursLumen.size(); k++)
830 (*list)[k]->getCG(&cgx,&cgy);
831 fprintf(archivo,"cx: %f cy:%f Se�al: %f Type: %d\n", cgx, cgy, (*list)[k]->getMaxIntensity(), (*list)[k]->getType());
834 fprintf(archivo,"PROMEDIO %f - STD: %f\n",avg,util->obtainStandardDeviation(average,avg));
835 fprintf(archivo,"%d\n",lumenContour.size());
841 // ----------------------------------------------------------------------------
842 int marAxisCT::getMaximumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
844 double max = list[iniPos]->getGradient();
845 int coordMax = iniPos;
846 int cont = iniPos + 1;
848 while (cont < limit && list[cont]->getGradient()>=threshold){
849 if (list[cont]->getGradient() > max){
850 max = list[cont]->getGradient();
859 // ----------------------------------------------------------------------------
860 int marAxisCT::getMinimumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
862 double min = list[iniPos]->getGradient();
863 int coordMin = iniPos;
864 int cont = iniPos + 1;
866 while (cont < limit && fabs(list[cont]->getGradient())>=threshold){
867 if (list[cont]->getGradient() < min){
868 min = list[cont]->getGradient();
877 // ----------------------------------------------------------------------------
878 double marAxisCT::interpolate(double x, double y, vtkImageData* imagedata) {
888 c = imagedata->GetScalarComponentAsDouble(i,j,0,0) + (double) imagedata->GetScalarComponentAsDouble(i+1,j+1,0,0)
889 - imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - imagedata->GetScalarComponentAsDouble(i,j+1,0,0);
891 b = imagedata->GetScalarComponentAsDouble(i,j+1,0,0) - c*i - imagedata->GetScalarComponentAsDouble(i,j,0,0);
893 a = imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - c*j - imagedata->GetScalarComponentAsDouble(i,j,0,0);
895 d = imagedata->GetScalarComponentAsDouble(i,j,0,0) - a*i - b*j - c*i*j;
897 intensity = a*x + b*y + c*x*y + d;
902 // ----------------------------------------------------------------------------
903 void marAxisCT::generateVector(int dir,int *coordBase,double originX,double originY,double *x,double *y) {
906 int angle = (dir * (360 / dirs));
908 *x = 1 * cos(angle *(3.1415/180)) + (*x);
909 *y = 1 * sin(angle *(3.1415/180)) + (*y);
914 // ----------------------------------------------------------------------------
915 bool marAxisCT::detectNeighbor(marIsocontour* contour,int dir, int type) {
917 bool isNeighbor = false;
919 if (contour == NULL){
922 for (int i = 0; i < contour->getSize(); i++) {
924 if (((contour->getDir(i) == (dir - 1)) || (contour->getDir(i) == dir))
925 && (contour->getType() == type)){
934 // ----------------------------------------------------------------------------
935 marIsocontour* marAxisCT::addPointToContour(marIsocontour* cont ,bool inside,int type,marPoint* point) {
941 cont = new marIsocontour();
945 cont->insertPoint(point->getX(),point->getY());
946 cont->setDir(cont->getSize() - 1,point->getDirection());
947 cont->setInside(cont->getSize() - 1,inside);
953 // ----------------------------------------------------------------------------
954 void marAxisCT::unifyDividedStructure(std::vector<marIsocontour*> *contList) {
956 marIsocontour* structOne = (*contList)[contList->size()-1];
957 marIsocontour* structTwo = (*contList)[0];
960 for (int i = 0; i < structOne->getSize(); i++)
964 structOne->getPoint(i,&x,&y);
965 structTwo->insertPoint(x,y);
966 structTwo->setDir(structTwo->getSize() - 1,structOne->getDir(i));
967 structTwo->setInside(structTwo->getSize() - 1,structOne->getInside(i));
971 contList->pop_back();
974 // ----------------------------------------------------------------------------
975 double marAxisCT::getStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata) {
977 double max, auxX, auxY;
979 max = interpolate(x,y,imagedata);
981 for (int i = 0; i < contour->getSize(); i++)
989 contour->getPoint(i,&pointX,&pointY);
990 while (auxX != pointX || auxY != pointY)
992 double intensity = interpolate(auxX,auxY,imagedata);
997 generateVector(contour->getDir(i),&basis,x,y,&auxX,&auxY);
1006 // ----------------------------------------------------------------------------
1007 double marAxisCT::getCalcStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata, int dir) {
1010 double max, auxX,auxY;
1015 char nomArchivo[30], temp[3];
1017 strcpy(nomArchivo,"./statistics");
1019 // EED 06/06/2007 linux
1020 // itoa(dir,temp,10);
1021 sprintf(temp,"%d",dir);
1023 strcat(nomArchivo,temp);
1024 strcat(nomArchivo,".txt");
1027 archivo = fopen(nomArchivo,"w");
1030 fprintf(archivo,"TAMA�O %d\n",contour->getSize());
1032 while (i < contour->getSize())
1035 double pointX, pointY;
1036 contour->getPoint(i+1,&pointX,&pointY);
1037 contour->getPoint(i,&auxX,&auxY);
1038 fprintf(archivo,"final: %f %f inicial %f %f \n",pointX,pointY,auxX,auxY);
1040 while (auxX != pointX || auxY != pointY)
1042 double intensity = interpolate(auxX,auxY,imagedata);
1043 if (max < intensity){
1046 generateVector(contour->getDir(i),&basis,auxX,auxY,&auxX,&auxY);
1055 // ----------------------------------------------------------------------------
1056 int marAxisCT::round(double num){
1061 decimal = num - floor(num) * 1.0000;
1063 if (decimal >= 0.50000)
1064 resp = (int)ceil(num);
1066 resp = (int)floor(num);
1072 // ----------------------------------------------------------------------------
1073 void marAxisCT::createEmptyContours(){
1075 int nCnts = ( int ) getNumberOfSplinePoints( );
1078 for (int i = 0; i < nCnts; i++)
1080 quantContours.push_back( NULL );
1081 lumenContour.push_back( NULL);
1085 // ----------------------------------------------------------------------------
1086 void marAxisCT::updateLumenPercentage(int point, kVolume* vol)
1088 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1089 if (quantContours[point] == NULL || wallCont->get2DContour() == NULL)
1091 create2Dcontours(point,vol);
1094 if (wallCont->getType() == marAxisContours::WALL
1095 && wallCont->isReplaced() == false)
1097 int a = ((marParameters *) getParameters())->getLumenPercentage();
1098 double oldValue = wallCont->getSignal();
1099 double lumenPercentage=((marParameters *) getParameters())->getLumenPercentage();
1100 double poww=pow(10.0,-2.0);
1101 double newValue = oldValue * lumenPercentage*poww;
1102 wallCont->getIsocontour()->SetValue(0,newValue);
1103 wallCont->getIsocontour()->Update();
1104 wallCont->get2DContour()->Update();
1109 // ----------------------------------------------------------------------------
1110 void marAxisCT::updateCalcPercentage(int point, kVolume* vol)
1115 if (quantContours[point] == NULL)
1117 create2Dcontours(point,vol);
1120 posMax = getMaxIntensity(point);
1122 if (quantContours[point]->getSize() > 1)
1124 oldValue = quantContours[point]->getContour(posMax)->getSignal();
1127 for (int i = 1; i < quantContours[point]->getSize(); i++)
1130 if (quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION
1131 && quantContours[point]->getContour(i)->isReplaced() == false)
1134 double newValue = oldValue * ((((marParameters *) getParameters())->getCalcPercentage())*pow(10.0,-2.0));
1135 quantContours[point]->getContour(i)->getIsocontour()->SetValue(0,newValue);
1136 quantContours[point]->getContour(i)->getIsocontour()->Update();
1137 quantContours[point]->getContour(i)->get2DContour()->Update();
1144 // ----------------------------------------------------------------------------
1145 int marAxisCT::getMaxIntensity(int point)
1149 if (quantContours[point]->getSize() > 2)
1151 intMax = quantContours[point]->getContour(0)->getSignal();
1152 for (int i = 0; i < quantContours[point]->getSize(); i++)
1154 if (quantContours[point]->getContour(i)->getSignal() > intMax && quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION)
1156 intMax = quantContours[point]->getContour(i)->getSignal() ;
1165 // ----------------------------------------------------------------------------
1166 void marAxisCT::eraseContours()
1170 int nCts = quantContours.size();
1171 vesselPoints.clear();
1177 delete quantContours[i];
1178 quantContours[i] = NULL;
1179 lumenContour[i] = NULL;
1182 quantContours.clear();
1183 lumenContour.clear();
1188 // ----------------------------------------------------------------------------
1189 void marAxisCT::eraseContoursPartial(int start)
1197 delete quantContours[i];
1198 quantContours[i] = NULL;
1199 lumenContour[i] = NULL;
1202 quantContours.clear();
1203 lumenContour.clear();
1205 //TODO Revisar esto que no me convence
1209 // ----------------------------------------------------------------------------
1210 bool marAxisCT::pointInPolygon(marIsocontour *c, double x, double y){
1212 bool oddNODES=false;
1213 int polySides = c->getSize();
1214 double *polyX, *polyY, tempX, tempY;
1216 polyX = (double *) malloc(polySides*sizeof(double));
1217 polyY = (double *) malloc(polySides*sizeof(double));
1219 for (i=0; i<polySides; i++) {
1220 c->getPoint(i,&tempX,&tempY);
1225 for (i=0; i<polySides; i++) {
1230 if (polyY[i]<y && polyY[j]>=y
1231 || polyY[j]<y && polyY[i]>=y) {
1232 if (polyX[i]+(y-polyY[i])/(polyY[j]-polyY[i])*(polyX[j]-polyX[i])<x) {
1244 // ----------------------------------------------------------------------------
1245 marIsocontour* marAxisCT::parsePolyDataToMarIsocontour(marContourVO* contourVO){
1247 vtkPolyData* polyDataCnt;
1248 vtkIdList* pointsId;
1249 double* pointPolyDataCnt;
1250 int numPointsPolyDataCnt, numTemp;
1251 marIsocontour *cont;
1253 if (contourVO->isReplaced())
1255 polyDataCnt = contourVO->get2DContour();
1259 polyDataCnt = contourVO->getIsocontourStripped()->GetOutput();
1262 numPointsPolyDataCnt = polyDataCnt->GetNumberOfPoints();
1263 pointsId = polyDataCnt->GetCell(0)->GetPointIds();
1264 numTemp = pointsId->GetNumberOfIds();
1266 cont = new marIsocontour();
1268 for (int i = 0; i < (numTemp -1); i++){
1269 pointPolyDataCnt = polyDataCnt->GetPoint(pointsId->GetId(i));
1270 cont->insertPoint(pointPolyDataCnt[0], pointPolyDataCnt[1]);
1279 // ----------------------------------------------------------------------------
1280 marIsocontour* marAxisCT::filterContour(marIsocontour* contExt, marIsocontour* contInt, double radio){
1283 double x1,y1,x2,y2,xInt,yInt, xAux, yAux;
1284 std::vector<marLine*> lineas;
1285 std::vector<double> dists;
1286 int count, countAux, ultima;
1287 bool hayInterseccion;
1289 marIsocontour* newCont = new marIsocontour();
1291 for (int i = 0; i < contExt->getSize(); i++){
1292 contExt->getPoint(i,&x1,&y1);
1293 contInt->getPoint(i,&x2,&y2);
1294 lineas.push_back(new marLine(x1,y1,x2,y2));
1295 tmp_eed = sqrt(pow((x2-x1),2.0) + pow((y2-y1),2.0));
1296 dists.push_back( tmp_eed );
1301 while (count < lineas.size() - 1){
1304 marLine* l1 = lineas[count];
1306 hayInterseccion = false;
1307 while (countAux < lineas.size()){
1308 marLine *l2 = lineas[countAux];
1310 l1->getIntersect(l2->getA(),l2->getB(),l2->getC(),&xInt,&yInt);
1313 if (xInt > 0 && yInt > 0){
1314 contExt->getPoint(countAux,&xAux,&yAux);
1315 if (dists[count] >= sqrt(pow((xInt-xAux),2.0) + pow((yInt-yAux),2.0))){
1316 contExt->getPoint(count,&x1,&y1);
1317 contInt->getPoint(count,&x2,&y2);
1320 if ((dists[count] > sqrt(pow((xInt-x1),2.0) + pow((yInt-y1),2.0))) &&
1321 (dists[count] > sqrt(pow((xInt-x2),2.0) + pow((yInt-y2),2.0))))
1324 hayInterseccion = true;
1334 if (!hayInterseccion){
1335 contInt->getPoint(count,&x2,&y2);
1336 newCont->insertPoint(x2,y2);
1345 // ----------------------------------------------------------------------------
1346 void marAxisCT::markUpLumen(int point, kVolume* vol)
1349 marContourVO* aVO = this->searchContour(marAxisContours::ELUMEN,point);
1356 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
1357 vtkImageData* imagedata;
1358 vtkImageData* erodedata;
1359 std::vector <marIsocontour* > polygons;
1360 std::vector <marPoint* > points;
1361 double average = 0, deviation = 0;
1362 vtkImageData* image = vtkImageData::New();
1363 vesselPoints.clear();
1366 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
1367 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
1368 image->SetDimensions(limSupX,limSupY,0);
1369 image->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
1370 image->SetScalarTypeToUnsignedShort();
1371 image->AllocateScalars();
1374 vtkImageThreshold *imageThreshold = vtkImageThreshold::New();
1375 imageThreshold->SetInput(imagedata);
1376 imageThreshold->ThresholdByUpper(maxSignal);
1377 imageThreshold->SetInValue(255);
1378 imageThreshold->SetOutValue(0.0);
1379 imageThreshold->Update();
1381 vtkImageContinuousErode3D *imageErode3D = vtkImageContinuousErode3D::New();
1382 imageErode3D->SetInput(imageThreshold->GetOutput());
1383 imageErode3D->SetKernelSize(4,4,1);
1384 imageErode3D->Update();
1385 erodedata = imageErode3D->GetOutput();
1387 // polygons.push_back(contExt);
1388 for (int i = 0; i < quantContours[point]->getSize(); i++)
1390 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
1393 //7- Wipe image to detect lumen
1396 for (x = limInfX + 30; x < (limSupX - 1); x++)
1399 for (int y = limInfY; y < (limSupY - 1) ; y++)
1401 bool inWall = false;
1402 bool outCalc = true;
1404 for (int numConts = 0; numConts < polygons.size(); numConts++)
1406 if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
1422 if (inWall && outCalc)
1424 marPoint *p = new marPoint(x,y);
1425 p->setIntensity(imagedata->GetScalarComponentAsDouble(x,y,0,0));
1426 average += p->getIntensity();
1427 points.push_back(p);
1430 /* if (cont >= (limSupY - limSupX))
1437 if (points.size() > 0)
1439 average = average / points.size();
1440 for (int k = 0; k < points.size(); k++)
1442 marPoint *p = points[k];
1443 double actualInt = p->getIntensity();
1444 double temp = average - actualInt;
1445 deviation += pow(temp,2.0);
1447 deviation = sqrt(deviation / points.size());
1452 polygons.push_back(lumenContour[point]);
1455 for (x = limInfX; x < (limSupX - 1); x++)
1457 for (int y = limInfY; y < (limSupY - 1); y++)
1459 bool inWall = false;
1460 bool outCalc = true;
1461 bool inLumen = false;
1462 unsigned short *refValue = (unsigned short *) image->GetScalarPointer(x,y,0);
1465 for (int numConts = 0; numConts < polygons.size(); numConts++)
1467 bool adentro = pointInPolygon(polygons[numConts],x,y);
1468 if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
1474 else if (numConts == polygons.size() - 1)
1486 double intensidad = imagedata->GetScalarComponentAsDouble(x,y,0,0);
1488 if (inWall && outCalc && inLumen)
1491 if (imagedata->GetScalarComponentAsDouble(x,y,0,0) > (average - 1.5*deviation)
1492 && imagedata->GetScalarComponentAsDouble(x,y,0,0) < (average + 2.5*deviation))
1504 extractLumen(image,lumenContour[point],point);
1509 // ----------------------------------------------------------------------------
1510 double marAxisCT::obtainContourArea(marIsocontour* contour)
1516 for (int i = 0; i < contour->getSize();i++)
1518 contour->getPoint(i,&x1,&y1);
1519 if (i < (contour->getSize() - 1))
1521 contour->getPoint(i+1,&x2,&y2);
1525 contour->getPoint(0,&x2,&y2);
1540 // ----------------------------------------------------------------------------
1541 void marAxisCT::adjustWall(int point, kVolume* vol)
1543 //int actualPercentage = ((marParameters *) getParameters())->getLumenPercentage();
1544 marIsocontour* cont;
1545 double actualArea, previousArea;
1546 std::vector <double > grad;
1547 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1550 ((marParameters *) getParameters())->setLumenPercentage(50);
1551 updateLumenPercentage(point,vol);
1552 cont = parsePolyDataToMarIsocontour(wallCont);
1553 previousArea = obtainContourArea(cont);
1555 for (int eval = 49; eval >= 0; eval--)
1557 ((marParameters *) getParameters())->setLumenPercentage(eval);
1558 updateLumenPercentage(point,vol);
1559 cont = parsePolyDataToMarIsocontour(wallCont);
1560 actualArea = obtainContourArea(cont);
1561 grad.push_back(fabs(actualArea - previousArea));
1564 int newPercentage = 50 - maxValue(grad);
1565 ((marParameters *) getParameters())->setLumenPercentage(newPercentage);
1566 updateLumenPercentage(point,vol);
1572 // ----------------------------------------------------------------------------
1573 void marAxisCT::adjustCalcification(int point, kVolume* vol)
1577 // ----------------------------------------------------------------------------
1578 int marAxisCT::maxValue(std::vector <double> list)
1581 double maxVal = list[0];
1583 for (int i = 0; i < list.size(); i++)
1585 if (list[i] > maxVal)
1596 // ----------------------------------------------------------------------------
1597 int marAxisCT::getDiameter(marContourVO* contourVO, double x, double y, int limInfX, int limInfY, int limSupX, int limSupY)
1603 int coordBase = 1, diameterOne = 0, diameterTwo = 0;
1604 marIsocontour* cont = parsePolyDataToMarIsocontour(contourVO);
1605 generateVector(1,&coordBase,originX,originY,&x,&y);
1607 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1609 generateVector(0,&coordBase,originX,originY,&x,&y);
1610 if (pointInPolygon(cont,x,y))
1623 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1625 generateVector(35,&coordBase,originX,originY,&x,&y);
1626 if (pointInPolygon(cont,x,y))
1639 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1641 generateVector(17,&coordBase,originX,originY,&x,&y);
1642 if (pointInPolygon(cont,x,y))
1655 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1657 generateVector(107,&coordBase,originX,originY,&x,&y);
1658 if (pointInPolygon(cont,x,y))
1670 return (diameterOne + diameterTwo) / 2;
1673 // ----------------------------------------------------------------------------
1674 void marAxisCT::histogram(int point, kVolume* vol)
1677 vtkImageData* imagedata;
1678 int limInfX, limInfY, limInfZ;
1679 int limSupX, limSupY, limSupZ;
1684 char nomArchivo[30],temp[3];
1685 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1688 // for (int point = 0; point < quantContours.size(); point++)
1690 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
1691 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
1692 ((marParameters *) getParameters())->setLumenPercentage(86);
1693 updateLumenPercentage(point,vol);
1694 marIsocontour* contExt = parsePolyDataToMarIsocontour(wallCont);
1696 std::vector<double> intensity;
1697 std::vector<int> frecuencies;
1699 for (int i = 0; i < limSupX; i++)
1701 for (int j = 0; j < limSupY; j++)
1703 if(pointInPolygon(contExt,i,j))
1706 pos = searchData(intensity,imagedata->GetScalarComponentAsDouble(i,j,0,0));
1709 intensity.push_back( imagedata->GetScalarComponentAsDouble(i,j,0,0));
1710 frecuencies.push_back(1);
1723 strcpy(nomArchivo,"./histogram");
1725 // EED 06/06/2007 linux
1726 // itoa(point,temp,10);
1727 sprintf(temp,"%d",point);
1729 strcat(nomArchivo,temp);
1730 strcat(nomArchivo,".csv");
1734 archivo = fopen(nomArchivo,"w");
1735 for (int k = 0; k < frecuencies.size(); k++)
1738 fprintf(archivo,"%f;%d\n",intensity[k],frecuencies[k]);
1746 //imagedata->Delete();
1752 // ----------------------------------------------------------------------------
1753 int marAxisCT::searchData(std::vector<double> vec, double value)
1758 for (int i = 0; i< vec.size() && !found; i++)
1760 if (vec[i] == value)
1770 // ----------------------------------------------------------------------------
1771 void marAxisCT::setCalibration(bool calib)
1773 calibration = calib;
1776 // ----------------------------------------------------------------------------
1777 bool marAxisCT::getCalibration()
1782 // ----------------------------------------------------------------------------
1783 void marAxisCT::adjustContour(int point, kVolume* vol)
1786 double percentage = 0.05;
1787 double prevDifference = 0;
1788 if (quantContours[point + 1] == NULL)
1790 create2Dcontours(point + 1 ,vol);
1793 //1. Obtain actual signal
1794 double signal = maxSignal;
1795 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1797 //2. Obtain previous contour
1798 marIsocontour* prev = parsePolyDataToMarIsocontour(searchContour(marAxisContours::WALL, point + 1));
1800 //3. Adjust actual contour
1801 double newValue = signal;
1802 wallCont->getIsocontour()->SetValue(0,newValue);
1803 wallCont->getIsocontour()->Update();
1804 wallCont->get2DContour()->Update();
1807 marIsocontour* actual = parsePolyDataToMarIsocontour(wallCont);
1810 double prevA = obtainContourArea(prev);
1811 double actuA = obtainContourArea(actual);
1813 prevDifference = fabs(prevA - actuA);
1816 if (prevA - actuA > 0 && prevA - actuA < percentage*prevA) //CASO 1
1819 while (!stop && newValue > 0)
1823 wallCont->getIsocontour()->SetValue(0,newValue);
1824 wallCont->getIsocontour()->Update();
1825 wallCont->get2DContour()->Update();
1827 actual = parsePolyDataToMarIsocontour(wallCont);
1828 actuA = obtainContourArea(actual);
1830 if (fabs(prevA - actuA) > percentage*prevA)
1833 wallCont->getIsocontour()->SetValue(0,newValue);
1834 wallCont->getIsocontour()->Update();
1835 wallCont->get2DContour()->Update();
1841 else if (prevA - actuA < 0 && prevA - actuA > -percentage*prevA) //CASO 2
1845 while (!stop && newValue > 0)
1849 wallCont->getIsocontour()->SetValue(0,newValue);
1850 wallCont->getIsocontour()->Update();
1851 wallCont->get2DContour()->Update();
1853 actual = parsePolyDataToMarIsocontour(wallCont);
1854 actuA = obtainContourArea(actual);
1856 if (prevA - actuA < -percentage*prevA)
1859 wallCont->getIsocontour()->SetValue(0,newValue);
1860 wallCont->getIsocontour()->Update();
1861 wallCont->get2DContour()->Update();
1866 else if (prevA - actuA < 0 && prevA - actuA < -percentage*prevA) //CASO 3
1875 wallCont->getIsocontour()->SetValue(0,newValue);
1876 wallCont->getIsocontour()->Update();
1877 wallCont->get2DContour()->Update();
1879 actual = parsePolyDataToMarIsocontour(wallCont);
1881 actuA = obtainContourArea(actual);
1882 if (prevA - actuA > -percentage*prevA)
1884 if (prevDifference < fabs(prevA - actuA))
1887 wallCont->getIsocontour()->SetValue(0,newValue);
1888 wallCont->getIsocontour()->Update();
1889 wallCont->get2DContour()->Update();
1895 prevDifference = fabs(prevA - actuA);
1898 else if (prevA - actuA > 0 && prevA - actuA > percentage*prevA) //CASO 4
1901 while (!stop && newValue > 0)
1905 wallCont->getIsocontour()->SetValue(0,newValue);
1906 wallCont->getIsocontour()->Update();
1907 wallCont->get2DContour()->Update();
1909 actual = parsePolyDataToMarIsocontour(wallCont);
1910 actuA = obtainContourArea(actual);
1912 if (prevA - actuA < percentage*prevA)
1914 if (prevDifference < fabs(prevA - actuA))
1917 wallCont->getIsocontour()->SetValue(0,newValue);
1918 wallCont->getIsocontour()->Update();
1919 wallCont->get2DContour()->Update();
1925 prevDifference = fabs(prevA - actuA);
1930 maxSignal = newValue;
1932 this->markUpLumen(point,vol);
1935 // ----------------------------------------------------------------------------
1936 int marAxisCT::getStartIndex()
1941 // ----------------------------------------------------------------------------
1942 void marAxisCT::setStartIndex(int start)
1946 int per = ((marParameters *) getParameters())->getLumenPercentage();
1947 double oldValue = quantContours[start]->getContour(0)->getSignal()*per*pow(10.0,-2.0);
1949 maxSignal = oldValue; // / (per*pow(10.0,-2.0));
1952 // ----------------------------------------------------------------------------
1953 int marAxisCT::getPointSize()
1955 return vesselPoints.size();
1958 // ----------------------------------------------------------------------------
1959 marPoint* marAxisCT::getPoint(int i)
1961 return vesselPoints[i];
1964 // ----------------------------------------------------------------------------
1965 void marAxisCT::extractLumen(vtkImageData *lumenImage, marIsocontour *lumenContour, int point)
1969 vtkImageCast* cast = vtkImageCast::New();
1970 cast->SetInput(lumenImage);
1971 cast->SetOutputScalarTypeToDouble();
1973 vtkImageGaussianSmooth* filter = vtkImageGaussianSmooth::New();
1974 filter->SetDimensionality(2);
1975 filter->SetInput(cast->GetOutput());
1976 filter->SetStandardDeviations((getParameters())->getStandardDeviation(),(getParameters())->getStandardDeviation());
1977 filter->SetRadiusFactors((getParameters())->getRadius(),(getParameters())->getRadius());
1978 filter->SetStandardDeviations(3,3);
1979 filter->SetRadiusFactors(3,3);
1983 marAxisContours* axisCont = quantContours[point];
1987 for (int i = 0; i < axisCont->getSize() && pos == -1; i++)
1989 if (axisCont->getContourType(i) == marAxisContours::ELUMEN)
1996 if (pos != -1) //Verify that contour hasn't been replaced
1999 if (axisCont->isReplaced(pos))
2008 marContourVO* contourVO = new marContourVO();
2009 marIsocontour* iso = lumenContour;
2011 contourVO->setIsocontour(vtkContourFilter::New());
2012 contourVO->getIsocontour()->SetInput( filter->GetOutput() );
2013 contourVO->getIsocontour()->SetNumberOfContours( 1 );
2014 contourVO->getIsocontour()->SetValue( 0, 128 );
2015 contourVO->getIsocontour()->Update( );
2017 contourVO->setIsocontourCpd(vtkCleanPolyData::New());
2018 contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
2019 contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
2020 contourVO->getIsocontourCpd()->Update( );
2024 iso->getCG(&cgx,&cgy);
2026 contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
2027 contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
2029 contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
2030 contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
2031 contourVO->getIsocontourDcf()->Update( );
2033 contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
2034 contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
2035 contourVO->getIsocontourCpd2()->Update();
2037 contourVO->setIsocontourStripped(vtkStripper::New( ));
2038 contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
2039 contourVO->getIsocontourStripped()->Update();
2041 contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
2042 contourVO->get2DContour()->Update();
2043 contourVO->setSignal(255);
2044 contourVO->setType(marAxisContours::ELUMEN);
2048 marIsocontour* actual = parsePolyDataToMarIsocontour(contourVO);
2049 marContourVO* aVO = searchContour(marAxisContours::WALL,point);
2050 double wallArea = 0;
2053 wallArea = obtainContourArea(parsePolyDataToMarIsocontour(aVO));
2056 double prevA = obtainContourArea(actual);
2059 int finalValue = 128;
2061 while (newValue > 0 && actuA > 0)
2068 contourVO->getIsocontour()->SetValue(0,newValue);
2069 contourVO->getIsocontour()->Update();
2070 contourVO->get2DContour()->Update();
2071 actual = parsePolyDataToMarIsocontour(contourVO);
2072 actuA = obtainContourArea(actual);
2073 if ((actuA /1000) < 1 && (prevA / 1000) > 1)
2075 finalValue = newValue;
2078 else if (actuA > prevA && (prevA / 1000) < 1 && actuA < 0.1*wallArea)
2080 finalValue = newValue;
2084 catch (std::exception e)
2091 contourVO->getIsocontour()->SetValue(0,finalValue);
2092 contourVO->getIsocontour()->Update();
2093 contourVO->get2DContour()->Update();
2101 axisCont->replaceContour(contourVO,pos);
2105 axisCont->addContour(contourVO);
2107 // axisCont->addContour(contourVO);
2110 quantContours[point] = axisCont;
2116 // ----------------------------------------------------------------------------
2117 void marAxisCT::generateFile(int point,kVolume* vol)
2120 char nomArchivo[30], temp[3];
2122 strcpy(nomArchivo,"./point");
2124 // EED 06/06/2007 linux
2125 // itoa(point,temp,10);
2126 sprintf(temp,"%d",point);
2128 strcat(nomArchivo,temp);
2129 strcat(nomArchivo,".txt");
2130 archivo = fopen(nomArchivo,"w");
2132 fprintf(archivo, "%d", quantContours[point]->getSize());
2135 for (i = 0; i < quantContours[point]->getSize(); i++)
2137 marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
2138 fprintf(archivo, " %d",iso->getSize());
2141 fprintf(archivo,"\n");
2143 for ( i = 0; i < quantContours[point]->getSize(); i++)
2145 marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
2146 fprintf(archivo,"Tipo: %d\n",quantContours[point]->getContourType(i));
2147 for (int j = 0; j < iso->getSize(); j++)
2149 double tempX, tempY;
2150 iso->getPoint(j,&tempX,&tempY);
2151 fprintf(archivo,"%f;%f;%d\n",tempX,tempY,point);
2155 // EED 30 Janvier 2007
2156 vtkImageData *vtkimagedata = this->getSliceImage(point,vol);
2158 std::string directory = "./";
2159 std::string filename = "xx";
2160 float rescalaSlope = 1;
2161 float rescalaIntercept = 0;
2163 vtkimagedata->GetDimensions(dim);
2171 marRAW2Files marraw2;
2172 marraw2.saveVolume(directory,filename,vtkimagedata,voi,rescalaSlope,rescalaIntercept);
2180 // ----------------------------------------------------------------------------
2181 void marAxisCT::replaceContour2D(int point,int size,double *vx,double *vy, int type)
2184 vtkPoints *_pts = vtkPoints::New();
2185 _pts->SetNumberOfPoints(size);
2188 for (j=0 ; j<size ; j++){
2189 _pts->SetPoint(j, vx[j] , vy[j] , 0 );
2191 // _pts->SetPoint(0, vx[0] , vy[0] , 0 );
2193 vtkCellArray *lines = vtkCellArray::New();
2194 lines->InsertNextCell( size );
2195 for ( j=0 ; j<size+1 ; j++ ){
2196 lines->InsertCellPoint(j % size );
2199 vtkPolyData *_pd = vtkPolyData::New();
2200 _pd->SetPoints( _pts );
2201 _pd->SetLines( lines );
2202 lines->Delete(); //do not delete lines ??
2205 marContourVO* vo = new marContourVO();///this->searchContour(type,point);
2206 vo->setReplaced(true);
2207 vo->set2DContour(_pd);
2209 quantContours[point]->addContour(vo); //->replaceContour(vo,i);
2210 /* _2Dcontours[i]=_pd;
2211 createContour(i,NULL);*/
2214 // ----------------------------------------------------------------------------
2215 void marAxisCT::cleanContours(int type, int point)
2217 marAxisContours* newContours = new marAxisContours();
2218 if (quantContours[point] != NULL)
2220 for (int i = 0; i < quantContours[point]->getSize(); i++)
2222 if (quantContours[point]->getContourType(i) != type)
2224 newContours->addContour(quantContours[point]->getContour(i));
2228 /* marContourVO* vo = quantContours[point]->getContour(i);
2229 vo->set2DContour(NULL);
2230 newContours->replaceContour(vo,i);*/
2235 quantContours[point] = newContours;
2238 // ----------------------------------------------------------------------------
2239 marContourVO* marAxisCT::searchContour(int type, int point)
2241 marContourVO* cont = NULL;
2243 for (int i = 0; i < quantContours[point]->getSize(); i++)
2245 if (quantContours[point]->getContourType(i) == type)
2247 cont = quantContours[point]->getContour(i);
2254 // ----------------------------------------------------------------------------
2255 double marAxisCT::performXOR(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2257 if (manual.size() == 0)
2262 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2263 std::vector <marIsocontour* > polygons;
2264 vtkImageData* imagedata;
2266 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
2267 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2268 vtkImageData* imageA = vtkImageData::New();
2269 imageA->SetDimensions(limSupX,limSupY,0);
2270 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2271 imageA->SetScalarTypeToUnsignedShort();
2272 imageA->AllocateScalars();
2275 vtkImageData* imageM = vtkImageData::New();
2276 imageM->SetDimensions(limSupX,limSupY,0);
2277 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2278 imageM->SetScalarTypeToUnsignedShort();
2279 imageM->AllocateScalars();
2283 for (int i = 0; i < quantContours[point]->getSize(); i++)
2285 if (quantContours[point]->getContourType(i) == type)
2287 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2292 if (polygons.size() == 0)
2297 for (x = limInfX; x < (limSupX - 1); x++)
2299 for (y = limInfY; y < (limSupY - 1); y++)
2301 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2302 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2307 for (numConts = 0; numConts < polygons.size(); numConts++)
2309 if (pointInPolygon(polygons[numConts],x,y))
2315 for (numConts = 0; numConts < manual.size(); numConts++)
2317 if (pointInPolygon(manual[numConts],x,y))
2325 /* vtkImageLogic* xor = vtkImageLogic::New();
2326 xor->SetInput1(imageM);
2327 xor->SetInput2(imageA);
2328 xor->SetOutputTrueValue(255);
2329 xor->SetOperationToXor();
2331 vtkImageCast* cast = vtkImageCast::New();
2332 cast->SetInput(xor->GetOutput());
2333 cast->SetOutputScalarTypeToDouble();
2335 int coincidencias = 0;
2336 for (x = limInfX; x < (limSupX - 1); x++)
2338 for (y = limInfY; y < (limSupY - 1); y++)
2340 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2341 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2343 int xor_ = compA^compM;
2351 double resp = (double)coincidencias ;
2357 // ----------------------------------------------------------------------------
2358 double marAxisCT::performAND(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2361 if (manual.size() == 0)
2366 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2367 std::vector <marIsocontour* > polygons;
2368 vtkImageData* imagedata;
2370 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
2371 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2372 vtkImageData* imageA = vtkImageData::New();
2373 imageA->SetDimensions(limSupX,limSupY,0);
2374 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2375 imageA->SetScalarTypeToUnsignedShort();
2376 imageA->AllocateScalars();
2379 vtkImageData* imageM = vtkImageData::New();
2380 imageM->SetDimensions(limSupX,limSupY,0);
2381 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2382 imageM->SetScalarTypeToUnsignedShort();
2383 imageM->AllocateScalars();
2386 for (int i = 0; i < quantContours[point]->getSize(); i++)
2388 if (quantContours[point]->getContourType(i) == type)
2390 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2395 if (polygons.size() == 0)
2401 for (x = limInfX; x < (limSupX - 1); x++)
2403 for (y = limInfY; y < (limSupY - 1); y++)
2405 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2406 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2411 for (numConts = 0; numConts < polygons.size(); numConts++)
2413 if (pointInPolygon(polygons[numConts],x,y))
2419 for (numConts = 0; numConts < manual.size(); numConts++)
2421 if (pointInPolygon(manual[numConts],x,y))
2429 /* vtkImageLogic* and = vtkImageLogic::New();
2430 and->SetInput1(imageM);
2431 and->SetInput2(imageA);
2432 and->SetOutputTrueValue(255);
2433 and->SetOperationToAnd();
2435 int coincidencias = 0;
2436 for (x = limInfX; x < (limSupX - 1); x++)
2438 for (y = limInfY; y < (limSupY - 1); y++)
2440 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2441 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2443 int and_ = compA & compM;
2451 double resp = (double)coincidencias ;
2457 // ----------------------------------------------------------------------------
2458 marIsocontour* marAxisCT::loadMarIsocontour(int size, double *vx, double *vy)
2460 marIsocontour* cont = new marIsocontour();
2462 for ( int j=0 ; j<size ; j++ )
2464 cont->insertPoint(vx[j], vy[j]);
2470 // ----------------------------------------------------------------------------
2471 double marAxisCT::performUnion(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2473 if (manual.size() == 0)
2478 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2479 std::vector <marIsocontour* > polygons;
2480 vtkImageData* imagedata;
2482 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
2483 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2484 vtkImageData* imageA = vtkImageData::New();
2485 imageA->SetDimensions(limSupX,limSupY,0);
2486 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2487 imageA->SetScalarTypeToUnsignedShort();
2488 imageA->AllocateScalars();
2491 vtkImageData* imageM = vtkImageData::New();
2492 imageM->SetDimensions(limSupX,limSupY,0);
2493 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2494 imageM->SetScalarTypeToUnsignedShort();
2495 imageM->AllocateScalars();
2499 for (i = 0; i < quantContours[point]->getSize(); i++)
2501 if (quantContours[point]->getContourType(i) == type)
2503 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2508 if (polygons.size() == 0)
2514 for ( x = limInfX; x < (limSupX - 1); x++)
2516 for ( y = limInfY; y < (limSupY - 1); y++)
2518 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2519 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2523 for ( numConts = 0; numConts < polygons.size(); numConts++)
2525 if (pointInPolygon(polygons[numConts],x,y))
2531 for (numConts = 0; numConts < manual.size(); numConts++)
2533 if (pointInPolygon(manual[numConts],x,y))
2541 /* vtkImageLogic* and = vtkImageLogic::New();
2542 and->SetInput1(imageM);
2543 and->SetInput2(imageA);
2544 and->SetOutputTrueValue(255);
2545 and->SetOperationToAnd();
2547 int coincidencias = 0;
2548 for (x = limInfX; x < (limSupX - 1); x++)
2550 for ( y = limInfY; y < (limSupY - 1); y++)
2552 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2553 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2555 int or_ = compA | compM;
2563 double resp = (double)coincidencias;