1 // marAxisCT.cpp: implementation of the marAxisCT class.
3 //////////////////////////////////////////////////////////////////////
6 #include <vtkImageCast.h>
7 #include <vtkImageGaussianSmooth.h>
8 #include <vtkCleanPolyData.h>
9 #include <vtkPolyDataConnectivityFilter.h>
10 #include <vtkStripper.h>
12 #include <vtkIdList.h>
13 #include <vtkMatrix4x4.h>
14 #include <vtkTransform.h>
15 #include <vtkPlaneSource.h>
16 #include <vtkContourFilter.h>
19 #include <vtkImageContinuousErode3D.h>
20 #include <vtkImageThreshold.h>
21 #include <vtkImageSobel2D.h>
22 #include <vtkCellArray.h>
23 #include <vtkImageLogic.h>
26 #include "marAxisCT.h"
27 #include "marGdcmDicom.h"
32 marAxisCT::marAxisCT( ) : marAxis( )
36 // ----------------------------------------------------------------------------
37 marContour* marAxisCT::getContour( int point , kVolume* vol, int index ) {
38 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->getContour() == NULL)
40 createContours(point,vol);
43 return quantContours[point]->getContour(index)->getContour();
46 // ----------------------------------------------------------------------------
47 vtkPoints* marAxisCT::get3Dcontour( int point , kVolume* vol, int index ) {
48 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get3DContour() == NULL)
50 create3Dcontours(point,vol);
53 marAxisContours* cont = quantContours[point];
55 return cont->getContour(index)->get3DContour();
58 // ----------------------------------------------------------------------------
59 vtkPolyData* marAxisCT::get2Dcontour( int point , kVolume* vol, int index ) {
60 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DContour() == NULL)
62 create2Dcontours(point,vol);
65 marAxisContours* cont = quantContours[point];
69 return cont->getContour(index)->get2DContour();
72 // ----------------------------------------------------------------------------
73 vtkPoints* marAxisCT::get2DDiameterMin( int point , kVolume* vol, int index ) {
74 // if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DDiameterMin() == NULL)
76 create2DDiametersMin(point,vol);
79 marAxisContours* cont = quantContours[point];
81 return cont->getContour(index)->get2DDiameterMin();
84 // ----------------------------------------------------------------------------
85 vtkPoints* marAxisCT::get2DDiameterMax( int point , kVolume* vol, int index ) {
86 // if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DDiameterMax() == NULL)
88 create2DDiametersMax(point,vol);
91 marAxisContours* cont = quantContours[point];
93 return cont->getContour(index)->get2DDiameterMax();
96 // ----------------------------------------------------------------------------
97 int marAxisCT::getSignal(int point, int index, kVolume* vol) {
98 if (quantContours[point] == NULL )
100 createSignals(point,vol);
103 return (int) ( quantContours[point]->getContour(index)->getSignal() );
106 // ----------------------------------------------------------------------------
107 void marAxisCT::createContours( int point, kVolume* vol ) {
110 int numberOfPoints,numberOfCells;
114 vtkPolyData* vtkPolyData_2DContour;
117 int sizeIma = getParameters( )->getSizeIma( );
118 double dimIma = getParameters( )->getDimIma( );
120 if (quantContours.size() < point - 1 )
122 create2Dcontours(point,vol);
125 marAxisContours* axisCont = quantContours[point];
127 for (int i = 0; i < axisCont->getSize(); i++)
129 marContourVO* aContourVO = axisCont->getContour(i);
130 aContourVO->setContour(new marContour( point, getParameters( ) ));
131 vtkPolyData_2DContour = aContourVO->get2DContour();
132 numberOfCells = vtkPolyData_2DContour->GetNumberOfCells( );
133 cell = vtkPolyData_2DContour->GetCell( 0 );
134 pids = cell->GetPointIds( );
135 numberOfPoints = pids->GetNumberOfIds( );
136 for( j = 0; j < numberOfPoints; j++ )
138 id = pids->GetId( j );
139 vtkPolyData_2DContour->GetPoint( id, p);
142 x=x * dimIma / ( float ) sizeIma;
143 y=y * dimIma / ( float ) sizeIma;
144 aContourVO->getContour()->addContourPoint( x , y );
147 aContourVO->getContour()->do_spline();
148 aContourVO->getContour()->calculateVariables();
153 // ----------------------------------------------------------------------------
154 void marAxisCT::create3Dcontours( int point , kVolume* vol ) {
156 if (quantContours[point] == NULL )
158 createContours(point,vol);
161 marAxisContours* axisCont = quantContours[point];
163 for (int i = 0; i < axisCont->getSize(); i++)
165 vtkPoints *_vtkPoints;
171 double c[3],p1[3],p2[3];
174 _vtkPoints = vtkPoints::New();
175 o = getPoints(point); //PENDIENTE _points[i];
177 mat = vtkMatrix4x4::New();
178 trans = vtkTransform::New();
181 nn[0]=n[0]; nn[1]=n[1]; nn[2]=n[2];
182 nC=sqrt( nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2] );
183 nn[0]=nn[0]/nC; nn[1]=nn[1]/nC; nn[2]=nn[2]/nC;
185 vtkPlaneSource* pSource = vtkPlaneSource::New( );
186 pSource->SetOrigin( 0, 0 , 0 );
187 pSource->SetPoint1( 1, 0 , 0 );
188 pSource->SetPoint2( 0, 0 , 1.0 );
189 pSource->SetCenter( 0,0,0 );
190 pSource->SetNormal( nn[ 0 ], nn[ 1 ], nn[ 2 ] );
193 pSource->GetOrigin(c);
194 pSource->GetPoint1(p1);
195 pSource->GetPoint2(p2);
197 p1[0]=p1[0]-c[0]; p1[1]=p1[1]-c[1]; p1[2]=p1[2]-c[2];
198 p2[0]=p2[0]-c[0]; p2[1]=p2[1]-c[1]; p2[2]=p2[2]-c[2];
199 nA=sqrt( p1[0]*p1[0] + p1[1]*p1[1] + p1[2]*p1[2] );
200 nB=sqrt( p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2] );
201 p1[0]=p1[0]/nA; p1[1]=p1[1]/nA; p1[2]=p1[2]/nA;
202 p2[0]=p2[0]/nB; p2[1]=p2[1]/nB; p2[2]=p2[2]/nB;
204 mat->SetElement (0,0, nn[0]);
205 mat->SetElement (1,0, nn[1]);
206 mat->SetElement (2,0, nn[2]);
207 mat->SetElement (3,0, 0);
208 mat->SetElement (0,1, p2[0]);
209 mat->SetElement (1,1, p2[1]);
210 mat->SetElement (2,1, p2[2]);
211 mat->SetElement (3,1, 0);
212 mat->SetElement (0,2, p1[0]);
213 mat->SetElement (1,2, p1[1]);
214 mat->SetElement (2,2, p1[2]);
215 mat->SetElement (3,2, 0);
216 mat->SetElement (0,3, 0);
217 mat->SetElement (1,3, 0);
218 mat->SetElement (2,3, 0);
219 mat->SetElement (3,3, 1);
221 double deter=mat->Determinant(mat);
223 trans->SetMatrix(mat);
226 trans->RotateWXYZ ( ang,0,1,0);
229 vtkMatrix4x4 *m=vtkMatrix4x4::New();
230 trans->GetMatrix ( m );
232 int j,numberOfPoints;
234 marContour* contour2D = getContour(point,vol,i);
235 marContourVO* aContourVO = axisCont->getContour(i);
236 numberOfPoints = contour2D->getNumberOfPoints( );
240 for( j = 0; j <= numberOfPoints; j++ ) {
241 contour2D->getPoint( fpA,j % numberOfPoints);
246 ppp[0] = m->GetElement(0,0)*pA[0] + m->GetElement(0,1)*pA[1] + m->GetElement(0,2)*pA[2] + m->GetElement(0,3);
247 ppp[1] = m->GetElement(1,0)*pA[0] + m->GetElement(1,1)*pA[1] + m->GetElement(1,2)*pA[2] + m->GetElement(1,3);
248 ppp[2] = m->GetElement(2,0)*pA[0] + m->GetElement(2,1)*pA[1] + m->GetElement(2,2)*pA[2] + m->GetElement(2,3);
250 ppp[0]=ppp[0]+o[0]-c[0]; ppp[1]=ppp[1]+o[1]-c[1]; ppp[2]=ppp[2]+o[2]-c[2];
251 _vtkPoints->InsertNextPoint( (float)ppp[0],(float)ppp[1],(float)ppp[2] );
261 if (aContourVO->get3DContour()!=NULL){ aContourVO->get3DContour()->Delete(); }
262 aContourVO->set3DContour(_vtkPoints);
271 // ----------------------------------------------------------------------------
272 void marAxisCT::create2Dcontours( int point , kVolume* vol ) {
273 std::vector <marIsocontour *> contours;
277 generatePoints(point,vol,&contours);
279 vtkImageData* imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
280 marAxisContours* axisCont = new marAxisContours();
282 for (int i = 0; i < contours.size(); i++)
284 marContourVO* contourVO = new marContourVO();
285 marIsocontour* iso = contours[i];
287 contourVO->setIsocontour(vtkContourFilter::New());
288 contourVO->getIsocontour()->SetInput( imagedata );
289 contourVO->getIsocontour()->SetNumberOfContours( 1 );
290 contourVO->getIsocontour()->SetValue( 0, iso->getMaxIntensity() );
291 contourVO->getIsocontour()->Update( );
293 contourVO->setIsocontourCpd(vtkCleanPolyData::New());
294 contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
295 contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
296 contourVO->getIsocontourCpd()->Update( );
300 iso->getCG(&cgx,&cgy);
302 contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
303 contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
305 contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
306 contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
307 contourVO->getIsocontourDcf()->Update( );
309 contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
310 contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
311 contourVO->getIsocontourCpd2()->Update();
313 contourVO->setIsocontourStripped(vtkStripper::New( ));
314 contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
315 contourVO->getIsocontourStripped()->Update();
317 contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
318 contourVO->get2DContour()->Update();
319 contourVO->setSignal(iso->getMaxIntensity());
321 if (iso->getType() == marAxisContours::LUMEN)
323 contourVO->setType(marAxisContours::WALL);
327 contourVO->setType(iso->getType());
332 axisCont->addContour(contourVO);
341 quantContours[point] = axisCont;
345 updateLumenPercentage(point,vol);
347 updateCalcPercentage(point,vol);
348 // adjustWall(point,vol);
349 if (calibration && point < startIndex)
351 this->adjustContour(point,vol);
355 adjustWall(point,vol);
361 // marIsocontour* cont = parsePolyDataToMarIsocontour(quantContours[point]->getContour(0));
364 char nomArchivo[30], temp[3];
366 strcpy(nomArchivo,"./points");
368 strcat(nomArchivo,temp);
369 strcat(nomArchivo,".csv");
371 archivo = fopen(nomArchivo,"w");
373 for (int h = 0; h < cont->getSize(); h++)
376 cont->getPoint(h,&x,&y);
377 fprintf(archivo,"%f;%f\n",x,y);
383 // ----------------------------------------------------------------------------
384 void marAxisCT::create2DDiametersMin(int point , kVolume* vol ) {
386 if (quantContours[point] == NULL )
388 createContours(point,vol);
391 marAxisContours* axisCont = quantContours[point];
393 for (int i = 0; i < axisCont->getSize(); i++)
396 marContour *marcontour=getContour(point,vol,i);
397 marContourVO* aContourVO = axisCont->getContour(i);
399 marcontour->getMinimumLine(p1,p2);
400 vtkPoints *_vtkPoints = vtkPoints::New();
402 int sizeIma = getParameters( )->getSizeIma( );
403 double dimIma = getParameters( )->getDimIma( );
404 double factor=( float ) sizeIma / dimIma;
406 p1[0]=p1[0]*factor+sizeIma/2;
407 p1[1]=p1[1]*factor+sizeIma/2;
408 p2[0]=p2[0]*factor+sizeIma/2;
409 p2[1]=p2[1]*factor+sizeIma/2;
410 _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
411 _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
413 if ( aContourVO->get2DDiameterMin()!=NULL ) { aContourVO->get2DDiameterMin()->Delete(); }
414 aContourVO->set2DDiameterMin(_vtkPoints);
420 // ----------------------------------------------------------------------------
421 void marAxisCT::create2DDiametersMax(int point , kVolume* vol ) {
422 if (quantContours[point] == NULL )
424 createContours(point,vol);
427 marAxisContours* axisCont = quantContours[point];
429 for (int i = 0; i < axisCont->getSize(); i++)
432 marContour *marcontour=getContour(point,vol,i);
433 marContourVO* aContourVO = axisCont->getContour(i);
435 marcontour->getMaximumLine(p1,p2);
436 vtkPoints *_vtkPoints = vtkPoints::New();
438 int sizeIma = getParameters( )->getSizeIma( );
439 double dimIma = getParameters( )->getDimIma( );
440 double factor=( float ) sizeIma / dimIma;
441 p1[0]=p1[0]*factor+sizeIma/2;
442 p1[1]=p1[1]*factor+sizeIma/2;
443 p2[0]=p2[0]*factor+sizeIma/2;
444 p2[1]=p2[1]*factor+sizeIma/2;
445 _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
446 _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
448 if ( aContourVO->get2DDiameterMax()!=NULL ) { aContourVO->get2DDiameterMax()->Delete(); }
449 aContourVO->set2DDiameterMax(_vtkPoints);
453 // ----------------------------------------------------------------------------
454 void marAxisCT::createSignals (int point, kVolume* vol) {
455 if (quantContours[point] == NULL )
457 create2Dcontours(point,vol); //AQUI
458 // filterSignal(vol);
462 // ----------------------------------------------------------------------------
463 int marAxisCT::getNumberOfContours(int point, kVolume* vol) {
464 if (quantContours[point] == NULL )
466 create2Dcontours(point,vol);//AQUI
467 // filterSignal(vol);
471 return quantContours[point]->getSize();
474 // ----------------------------------------------------------------------------
475 int marAxisCT::getContourType(int point, int index, kVolume* vol) {
476 if (quantContours[point] == NULL )
478 create2Dcontours(point,vol);
481 return quantContours[point]->getContourType(index);
485 // ----------------------------------------------------------------------------
486 void marAxisCT::generatePoints(int point, kVolume* vol, std::vector<marIsocontour*> *list) {
488 vtkImageData* imagedata;
491 double *nPoint = new double[ 3 ];
492 double x, y, prevIntensity, threshold;
493 bool found, estStart, estDivided, isNeighbor;
494 int actualZone, startZone, coordBase, maxPoint, lumenCount, calcCount;
495 std::vector <marPoint *> gradients;
496 std::vector <double> average;
497 std::vector <marIsocontour *> contoursLumen;
498 std::vector <marIsocontour *> contoursCalc;
500 contoursLumen.push_back(NULL);
501 contoursCalc.push_back(NULL);
505 //1. Check starting zone
506 startZone = marAxisContours::LUMEN;
507 actualZone = startZone;
511 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
513 int limInfX, limInfY, limInfZ;
514 int limSupX, limSupY, limSupZ;
515 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
517 vtkImageCast* cast = vtkImageCast::New();
518 cast->SetInput(imagedata);
519 cast->SetOutputScalarTypeToDouble();
521 vtkImageGaussianSmooth* filter = vtkImageGaussianSmooth::New();
522 filter->SetDimensionality(2);
523 filter->SetInput(cast->GetOutput());
524 filter->SetStandardDeviations((getParameters())->getStandardDeviation(),(getParameters())->getStandardDeviation());
525 filter->SetRadiusFactors((getParameters())->getRadius(),(getParameters())->getRadius());
526 filter->SetStandardDeviations(3,3);
527 filter->SetRadiusFactors(3,3);
529 //imagedata = filter->GetOutput();
532 nPoint = getPoints(point);
534 aPoint[ 0 ] = (limSupX - limInfX) / 2; //nPoint[0];
535 aPoint[ 1 ] = (limSupY - limInfY) / 2; //nPoint[1];
536 aPoint[ 2 ] = (limSupZ - limInfZ) / 2; // nPoint[2];
539 //3. Generate rays over all directions
540 for (int dir = 0; dir < dirs; dir++)
546 prevIntensity = interpolate(x,y,imagedata);
547 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
549 //2.1 Evaluate gradient in every position of the ray
550 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
554 double resp = interpolate(x,y,imagedata);
555 //2.1.2 Calcultate gradient
556 marPoint* p = new marPoint(x,y);
557 p->setGradient(prevIntensity - resp);
558 p->setDirection(dir);
559 p->setIntensity(resp);
561 prevIntensity = resp;
563 gradients.push_back(p);
565 //2.1.3 Increase ray for next evaluation
566 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
568 }// end while gradient
574 threshold = ((marParameters *) getParameters())->getContourThresh();
575 int mysize = gradients.size();
577 //2.2 Compare every gradient value with threshold
578 while (index < gradients.size() && !found)
580 //2.2.1 Evaluate Lumen
581 if (actualZone == marAxisContours::LUMEN)
583 average.push_back(gradients[index]->getIntensity());
585 //2.2.1.1 If value > threshold there is a contour
586 if (fabs(gradients[index]->getGradient()) > threshold
587 && fabs(gradients[index]->getGradient()) < 4000) {
588 //|| (calibration && (abs(gradients[index]->getIntensity()) > abs(avgValue - 2.5*stdValue)
589 //|| abs(gradients[index]->getIntensity()) > abs(avgValue + 2.5*stdValue)))) {
593 if(gradients[index]->getGradient() > 0)
594 maxPoint = getMaximumGrad(gradients,index,gradients.size(),threshold);
596 maxPoint = getMinimumGrad(gradients,index,gradients.size(),threshold);
599 x = gradients[maxPoint]->getX();
600 y = gradients[maxPoint]->getY();
602 //2.2.1.1.1 Assign lumen structure
603 isNeighbor = detectNeighbor(contoursLumen[lumenCount],dir,marAxisContours::LUMEN);
609 marIsocontour* cont = addPointToContour(contoursLumen[lumenCount],false,marAxisContours::LUMEN,gradients[maxPoint]);
610 contoursLumen[lumenCount] = cont;
612 //2.2.1.1.2 Check step lumen->calcification
613 double actualInt = interpolate(x,y,imagedata);
616 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&auxX,&auxY);
617 double next = interpolate(auxX,auxY,imagedata);
618 if (gradients[index]->getGradient() < 0 /*&& actualInt < next */&& (auxX > (limInfX+1) && auxX < (limSupX-1)) && (auxY > (limInfY+1) && auxY < (limSupY-1)))
620 //2.2.1.1.2.1 Assign calcification structure
621 isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
625 contoursCalc.push_back(NULL);
628 cont = addPointToContour(contoursCalc[calcCount],true,marAxisContours::CALCIFICATION,gradients[maxPoint]);
630 contoursCalc[calcCount] = cont;
632 actualZone = marAxisContours::CALCIFICATION;
635 //2.2.1.1.2.3 Verify divided structure
641 if ((dir == dirs - 1) && estStart)
646 }//2.2.1.1.3 Transition lumen->background
649 //2.2.1.1.3.1 Asign border structure
652 //2.2.1.1.3.2 Generate stop condition
654 actualZone = startZone;
658 threshold = ((marParameters *) getParameters())->getContourThresh();
663 //2.2.2 Evaluate calcification (MODIFIED IF THERE IS HYPODENSE)
668 double actualInt = interpolate(x,y,imagedata);
669 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
670 double next = interpolate(x,y,imagedata);
672 //2.2.2.1 If actual intensity is smaller than previous there is a contour
673 if (actualInt < next)
675 while (actualInt < interpolate(x,y,imagedata) && (x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
677 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
683 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
685 //2.2.2.2 Assign calcification structure
686 isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
690 contoursCalc.push_back(NULL);
693 //Adjust index to avoiod crash
694 if (index >= gradients.size())
696 calcIndex = gradients.size() - 1;
703 marIsocontour* cont = addPointToContour(contoursCalc[calcCount],false,marAxisContours::CALCIFICATION,gradients[calcIndex]);
704 contoursCalc[calcCount] = cont;
706 //2.2.2.3 Asign border structure
710 actualZone = startZone;
713 threshold = ((marParameters *) getParameters())->getContourThresh();
719 if (index == gradients.size() && !found)
721 threshold = threshold - 1;
725 }//end while evaluation
727 //TODO - SACAR A UTILS
728 for (int ind = 0; ind < gradients.size(); ind++)
730 marPoint* p = gradients[ind];
731 gradients[ind] = NULL;
737 }//end for directions
740 char nomArchivo[30]; //, temp[3]
742 strcpy(nomArchivo,"./debug");
743 // itoa(point,temp,10);
744 // strcat(nomArchivo,temp);
745 strcat(nomArchivo,".txt");
746 archivo = fopen(nomArchivo,"w");
752 // lumenContour[point] = contoursLumen[marAxisContours::LUMEN];
754 double tempXi, tempYi;
755 int polySides = contoursLumen[marAxisContours::LUMEN]->getSize();
758 marIsocontour* alc = new marIsocontour();
759 for (int m=0; m<polySides; m++) {
760 contoursLumen[marAxisContours::LUMEN]->getPoint(m,&tempXi,&tempYi);
761 alc = addPointToContour(lumenContour[point],false,marAxisContours::LUMEN,new marPoint(tempXi,tempYi));
762 lumenContour[point] = alc;
763 fprintf(archivo,"X:%f,Y:%f\n",tempXi,tempYi);
768 //4. Verify divided unified structure
769 if (estDivided && calcCount > 0)
771 unifyDividedStructure(&contoursCalc);
775 //5. Obtain statistics
776 double lumenInt = getStatistics(contoursLumen[marAxisContours::LUMEN],aPoint[ 0 ],aPoint[ 1 ], imagedata);
777 contoursLumen[marAxisContours::LUMEN]->setMaxIntensity(lumenInt);
779 if (contoursCalc[0] != NULL)
781 for (int i = 0; i < contoursCalc.size(); i++)
784 double maxInt = getCalcStatistics(contoursCalc[i],aPoint[ 0 ],aPoint[1],imagedata,i-1);
785 contoursCalc[i]->setMaxIntensity(maxInt);
786 contoursLumen.push_back(NULL);
787 contoursLumen[lumenCount] = contoursCalc[i];
792 *list = contoursLumen;
794 contoursCalc.clear();
795 contoursLumen.clear();
796 marUtils* util = new marUtils();
797 double avg = util->obtainAverage(average);
799 //archivo = fopen("./debug.txt","w");
800 fprintf(archivo,"tamaño: %d umbra: %d", (*list).size(),((marParameters *) getParameters())->getContourThresh());
801 for (int k = 0; k < contoursLumen.size(); k++)
805 (*list)[k]->getCG(&cgx,&cgy);
806 fprintf(archivo,"cx: %f cy:%f Señal: %f Type: %d\n", cgx, cgy, (*list)[k]->getMaxIntensity(), (*list)[k]->getType());
809 fprintf(archivo,"PROMEDIO %f - STD: %f\n",avg,util->obtainStandardDeviation(average,avg));
810 fprintf(archivo,"%d\n",lumenContour.size());
816 // ----------------------------------------------------------------------------
817 int marAxisCT::getMaximumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
819 double max = list[iniPos]->getGradient();
820 int coordMax = iniPos;
821 int cont = iniPos + 1;
823 while (cont < limit && list[cont]->getGradient()>=threshold){
824 if (list[cont]->getGradient() > max){
825 max = list[cont]->getGradient();
834 // ----------------------------------------------------------------------------
835 int marAxisCT::getMinimumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
837 double min = list[iniPos]->getGradient();
838 int coordMin = iniPos;
839 int cont = iniPos + 1;
841 while (cont < limit && fabs(list[cont]->getGradient())>=threshold){
842 if (list[cont]->getGradient() < min){
843 min = list[cont]->getGradient();
852 // ----------------------------------------------------------------------------
853 double marAxisCT::interpolate(double x, double y, vtkImageData* imagedata) {
863 c = imagedata->GetScalarComponentAsDouble(i,j,0,0) + (double) imagedata->GetScalarComponentAsDouble(i+1,j+1,0,0)
864 - imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - imagedata->GetScalarComponentAsDouble(i,j+1,0,0);
866 b = imagedata->GetScalarComponentAsDouble(i,j+1,0,0) - c*i - imagedata->GetScalarComponentAsDouble(i,j,0,0);
868 a = imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - c*j - imagedata->GetScalarComponentAsDouble(i,j,0,0);
870 d = imagedata->GetScalarComponentAsDouble(i,j,0,0) - a*i - b*j - c*i*j;
872 intensity = a*x + b*y + c*x*y + d;
877 // ----------------------------------------------------------------------------
878 void marAxisCT::generateVector(int dir,int *coordBase,double originX,double originY,double *x,double *y) {
881 int angle = (dir * (360 / dirs));
883 *x = 1 * cos(angle *(3.1415/180)) + (*x);
884 *y = 1 * sin(angle *(3.1415/180)) + (*y);
889 // ----------------------------------------------------------------------------
890 bool marAxisCT::detectNeighbor(marIsocontour* contour,int dir, int type) {
892 bool isNeighbor = false;
894 if (contour == NULL){
897 for (int i = 0; i < contour->getSize(); i++) {
899 if (((contour->getDir(i) == (dir - 1)) || (contour->getDir(i) == dir))
900 && (contour->getType() == type)){
909 // ----------------------------------------------------------------------------
910 marIsocontour* marAxisCT::addPointToContour(marIsocontour* cont ,bool inside,int type,marPoint* point) {
916 cont = new marIsocontour();
920 cont->insertPoint(point->getX(),point->getY());
921 cont->setDir(cont->getSize() - 1,point->getDirection());
922 cont->setInside(cont->getSize() - 1,inside);
928 // ----------------------------------------------------------------------------
929 void marAxisCT::unifyDividedStructure(std::vector<marIsocontour*> *contList) {
931 marIsocontour* structOne = (*contList)[contList->size()-1];
932 marIsocontour* structTwo = (*contList)[0];
935 for (int i = 0; i < structOne->getSize(); i++)
939 structOne->getPoint(i,&x,&y);
940 structTwo->insertPoint(x,y);
941 structTwo->setDir(structTwo->getSize() - 1,structOne->getDir(i));
942 structTwo->setInside(structTwo->getSize() - 1,structOne->getInside(i));
946 contList->pop_back();
949 // ----------------------------------------------------------------------------
950 double marAxisCT::getStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata) {
952 double max, auxX, auxY;
954 max = interpolate(x,y,imagedata);
956 for (int i = 0; i < contour->getSize(); i++)
964 contour->getPoint(i,&pointX,&pointY);
965 while (auxX != pointX || auxY != pointY)
967 double intensity = interpolate(auxX,auxY,imagedata);
972 generateVector(contour->getDir(i),&basis,x,y,&auxX,&auxY);
981 // ----------------------------------------------------------------------------
982 double marAxisCT::getCalcStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata, int dir) {
985 double max, auxX,auxY;
990 char nomArchivo[30], temp[3];
992 strcpy(nomArchivo,"./statistics");
994 // EED 06/06/2007 linux
995 // itoa(dir,temp,10);
996 sprintf(temp,"%d",dir);
998 strcat(nomArchivo,temp);
999 strcat(nomArchivo,".txt");
1002 archivo = fopen(nomArchivo,"w");
1005 fprintf(archivo,"TAMAÑO %d\n",contour->getSize());
1007 while (i < contour->getSize())
1010 double pointX, pointY;
1011 contour->getPoint(i+1,&pointX,&pointY);
1012 contour->getPoint(i,&auxX,&auxY);
1013 fprintf(archivo,"final: %f %f inicial %f %f \n",pointX,pointY,auxX,auxY);
1015 while (auxX != pointX || auxY != pointY)
1017 double intensity = interpolate(auxX,auxY,imagedata);
1018 if (max < intensity){
1021 generateVector(contour->getDir(i),&basis,auxX,auxY,&auxX,&auxY);
1030 // ----------------------------------------------------------------------------
1031 int marAxisCT::round(double num){
1036 decimal = num - floor(num) * 1.0000;
1038 if (decimal >= 0.50000)
1039 resp = (int)ceil(num);
1041 resp = (int)floor(num);
1047 // ----------------------------------------------------------------------------
1048 void marAxisCT::createEmptyContours(){
1050 int nCnts = ( int ) getNumberOfSplinePoints( );
1053 for (int i = 0; i < nCnts; i++)
1055 quantContours.push_back( NULL );
1056 lumenContour.push_back( NULL);
1060 // ----------------------------------------------------------------------------
1061 void marAxisCT::updateLumenPercentage(int point, kVolume* vol)
1063 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1064 if (quantContours[point] == NULL || wallCont->get2DContour() == NULL)
1066 create2Dcontours(point,vol);
1069 if (wallCont->getType() == marAxisContours::WALL
1070 && wallCont->isReplaced() == false)
1072 int a = ((marParameters *) getParameters())->getLumenPercentage();
1073 double oldValue = wallCont->getSignal();
1074 double lumenPercentage=((marParameters *) getParameters())->getLumenPercentage();
1075 double poww=pow(10.0,-2.0);
1076 double newValue = oldValue * lumenPercentage*poww;
1077 wallCont->getIsocontour()->SetValue(0,newValue);
1078 wallCont->getIsocontour()->Update();
1079 wallCont->get2DContour()->Update();
1084 // ----------------------------------------------------------------------------
1085 void marAxisCT::updateCalcPercentage(int point, kVolume* vol)
1090 if (quantContours[point] == NULL)
1092 create2Dcontours(point,vol);
1095 posMax = getMaxIntensity(point);
1097 if (quantContours[point]->getSize() > 1)
1099 oldValue = quantContours[point]->getContour(posMax)->getSignal();
1102 for (int i = 1; i < quantContours[point]->getSize(); i++)
1105 if (quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION
1106 && quantContours[point]->getContour(i)->isReplaced() == false)
1109 double newValue = oldValue * ((((marParameters *) getParameters())->getCalcPercentage())*pow(10.0,-2.0));
1110 quantContours[point]->getContour(i)->getIsocontour()->SetValue(0,newValue);
1111 quantContours[point]->getContour(i)->getIsocontour()->Update();
1112 quantContours[point]->getContour(i)->get2DContour()->Update();
1119 // ----------------------------------------------------------------------------
1120 int marAxisCT::getMaxIntensity(int point)
1124 if (quantContours[point]->getSize() > 2)
1126 intMax = quantContours[point]->getContour(0)->getSignal();
1127 for (int i = 0; i < quantContours[point]->getSize(); i++)
1129 if (quantContours[point]->getContour(i)->getSignal() > intMax && quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION)
1131 intMax = quantContours[point]->getContour(i)->getSignal() ;
1140 // ----------------------------------------------------------------------------
1141 void marAxisCT::eraseContours()
1145 int nCts = quantContours.size();
1146 vesselPoints.clear();
1152 delete quantContours[i];
1153 quantContours[i] = NULL;
1154 lumenContour[i] = NULL;
1157 quantContours.clear();
1158 lumenContour.clear();
1163 // ----------------------------------------------------------------------------
1164 void marAxisCT::eraseContoursPartial(int start)
1172 delete quantContours[i];
1173 quantContours[i] = NULL;
1174 lumenContour[i] = NULL;
1177 quantContours.clear();
1178 lumenContour.clear();
1180 //TODO Revisar esto que no me convence
1184 // ----------------------------------------------------------------------------
1185 bool marAxisCT::pointInPolygon(marIsocontour *c, double x, double y){
1187 bool oddNODES=false;
1188 int polySides = c->getSize();
1189 double *polyX, *polyY, tempX, tempY;
1191 polyX = (double *) malloc(polySides*sizeof(double));
1192 polyY = (double *) malloc(polySides*sizeof(double));
1194 for (i=0; i<polySides; i++) {
1195 c->getPoint(i,&tempX,&tempY);
1200 for (i=0; i<polySides; i++) {
1205 if (polyY[i]<y && polyY[j]>=y
1206 || polyY[j]<y && polyY[i]>=y) {
1207 if (polyX[i]+(y-polyY[i])/(polyY[j]-polyY[i])*(polyX[j]-polyX[i])<x) {
1219 // ----------------------------------------------------------------------------
1220 marIsocontour* marAxisCT::parsePolyDataToMarIsocontour(marContourVO* contourVO){
1222 vtkPolyData* polyDataCnt;
1223 vtkIdList* pointsId;
1224 double* pointPolyDataCnt;
1225 int numPointsPolyDataCnt, numTemp;
1226 marIsocontour *cont;
1228 if (contourVO->isReplaced())
1230 polyDataCnt = contourVO->get2DContour();
1234 polyDataCnt = contourVO->getIsocontourStripped()->GetOutput();
1237 numPointsPolyDataCnt = polyDataCnt->GetNumberOfPoints();
1238 pointsId = polyDataCnt->GetCell(0)->GetPointIds();
1239 numTemp = pointsId->GetNumberOfIds();
1241 cont = new marIsocontour();
1243 for (int i = 0; i < (numTemp -1); i++){
1244 pointPolyDataCnt = polyDataCnt->GetPoint(pointsId->GetId(i));
1245 cont->insertPoint(pointPolyDataCnt[0], pointPolyDataCnt[1]);
1254 // ----------------------------------------------------------------------------
1255 marIsocontour* marAxisCT::filterContour(marIsocontour* contExt, marIsocontour* contInt, double radio){
1258 double x1,y1,x2,y2,xInt,yInt, xAux, yAux;
1259 std::vector<marLine*> lineas;
1260 std::vector<double> dists;
1261 int count, countAux, ultima;
1262 bool hayInterseccion;
1264 marIsocontour* newCont = new marIsocontour();
1266 for (int i = 0; i < contExt->getSize(); i++){
1267 contExt->getPoint(i,&x1,&y1);
1268 contInt->getPoint(i,&x2,&y2);
1269 lineas.push_back(new marLine(x1,y1,x2,y2));
1270 tmp_eed = sqrt(pow((x2-x1),2.0) + pow((y2-y1),2.0));
1271 dists.push_back( tmp_eed );
1276 while (count < lineas.size() - 1){
1279 marLine* l1 = lineas[count];
1281 hayInterseccion = false;
1282 while (countAux < lineas.size()){
1283 marLine *l2 = lineas[countAux];
1285 l1->getIntersect(l2->getA(),l2->getB(),l2->getC(),&xInt,&yInt);
1288 if (xInt > 0 && yInt > 0){
1289 contExt->getPoint(countAux,&xAux,&yAux);
1290 if (dists[count] >= sqrt(pow((xInt-xAux),2.0) + pow((yInt-yAux),2.0))){
1291 contExt->getPoint(count,&x1,&y1);
1292 contInt->getPoint(count,&x2,&y2);
1295 if ((dists[count] > sqrt(pow((xInt-x1),2.0) + pow((yInt-y1),2.0))) &&
1296 (dists[count] > sqrt(pow((xInt-x2),2.0) + pow((yInt-y2),2.0))))
1299 hayInterseccion = true;
1309 if (!hayInterseccion){
1310 contInt->getPoint(count,&x2,&y2);
1311 newCont->insertPoint(x2,y2);
1320 // ----------------------------------------------------------------------------
1321 void marAxisCT::markUpLumen(int point, kVolume* vol)
1324 marContourVO* aVO = this->searchContour(marAxisContours::ELUMEN,point);
1331 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
1332 vtkImageData* imagedata;
1333 vtkImageData* erodedata;
1334 std::vector <marIsocontour* > polygons;
1335 std::vector <marPoint* > points;
1336 double average = 0, deviation = 0;
1337 vtkImageData* image = vtkImageData::New();
1338 vesselPoints.clear();
1341 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
1342 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
1343 image->SetDimensions(limSupX,limSupY,0);
1344 image->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
1345 image->SetScalarTypeToUnsignedShort();
1346 image->AllocateScalars();
1349 vtkImageThreshold *imageThreshold = vtkImageThreshold::New();
1350 imageThreshold->SetInput(imagedata);
1351 imageThreshold->ThresholdByUpper(maxSignal);
1352 imageThreshold->SetInValue(255);
1353 imageThreshold->SetOutValue(0.0);
1354 imageThreshold->Update();
1356 vtkImageContinuousErode3D *imageErode3D = vtkImageContinuousErode3D::New();
1357 imageErode3D->SetInput(imageThreshold->GetOutput());
1358 imageErode3D->SetKernelSize(4,4,1);
1359 imageErode3D->Update();
1360 erodedata = imageErode3D->GetOutput();
1362 // polygons.push_back(contExt);
1363 for (int i = 0; i < quantContours[point]->getSize(); i++)
1365 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
1368 //7- Wipe image to detect lumen
1371 for (x = limInfX + 30; x < (limSupX - 1); x++)
1374 for (int y = limInfY; y < (limSupY - 1) ; y++)
1376 bool inWall = false;
1377 bool outCalc = true;
1379 for (int numConts = 0; numConts < polygons.size(); numConts++)
1381 if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
1397 if (inWall && outCalc)
1399 marPoint *p = new marPoint(x,y);
1400 p->setIntensity(imagedata->GetScalarComponentAsDouble(x,y,0,0));
1401 average += p->getIntensity();
1402 points.push_back(p);
1405 /* if (cont >= (limSupY - limSupX))
1412 if (points.size() > 0)
1414 average = average / points.size();
1415 for (int k = 0; k < points.size(); k++)
1417 marPoint *p = points[k];
1418 double actualInt = p->getIntensity();
1419 double temp = average - actualInt;
1420 deviation += pow(temp,2.0);
1422 deviation = sqrt(deviation / points.size());
1427 polygons.push_back(lumenContour[point]);
1430 for (x = limInfX; x < (limSupX - 1); x++)
1432 for (int y = limInfY; y < (limSupY - 1); y++)
1434 bool inWall = false;
1435 bool outCalc = true;
1436 bool inLumen = false;
1437 unsigned short *refValue = (unsigned short *) image->GetScalarPointer(x,y,0);
1440 for (int numConts = 0; numConts < polygons.size(); numConts++)
1442 bool adentro = pointInPolygon(polygons[numConts],x,y);
1443 if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
1449 else if (numConts == polygons.size() - 1)
1461 double intensidad = imagedata->GetScalarComponentAsDouble(x,y,0,0);
1463 if (inWall && outCalc && inLumen)
1466 if (imagedata->GetScalarComponentAsDouble(x,y,0,0) > (average - 1.5*deviation)
1467 && imagedata->GetScalarComponentAsDouble(x,y,0,0) < (average + 2.5*deviation))
1479 extractLumen(image,lumenContour[point],point);
1484 // ----------------------------------------------------------------------------
1485 double marAxisCT::obtainContourArea(marIsocontour* contour)
1491 for (int i = 0; i < contour->getSize();i++)
1493 contour->getPoint(i,&x1,&y1);
1494 if (i < (contour->getSize() - 1))
1496 contour->getPoint(i+1,&x2,&y2);
1500 contour->getPoint(0,&x2,&y2);
1515 // ----------------------------------------------------------------------------
1516 void marAxisCT::adjustWall(int point, kVolume* vol)
1518 //int actualPercentage = ((marParameters *) getParameters())->getLumenPercentage();
1519 marIsocontour* cont;
1520 double actualArea, previousArea;
1521 std::vector <double > grad;
1522 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1525 ((marParameters *) getParameters())->setLumenPercentage(50);
1526 updateLumenPercentage(point,vol);
1527 cont = parsePolyDataToMarIsocontour(wallCont);
1528 previousArea = obtainContourArea(cont);
1530 for (int eval = 49; eval >= 0; eval--)
1532 ((marParameters *) getParameters())->setLumenPercentage(eval);
1533 updateLumenPercentage(point,vol);
1534 cont = parsePolyDataToMarIsocontour(wallCont);
1535 actualArea = obtainContourArea(cont);
1536 grad.push_back(fabs(actualArea - previousArea));
1539 int newPercentage = 50 - maxValue(grad);
1540 ((marParameters *) getParameters())->setLumenPercentage(newPercentage);
1541 updateLumenPercentage(point,vol);
1547 // ----------------------------------------------------------------------------
1548 void marAxisCT::adjustCalcification(int point, kVolume* vol)
1552 // ----------------------------------------------------------------------------
1553 int marAxisCT::maxValue(std::vector <double> list)
1556 double maxVal = list[0];
1558 for (int i = 0; i < list.size(); i++)
1560 if (list[i] > maxVal)
1571 // ----------------------------------------------------------------------------
1572 int marAxisCT::getDiameter(marContourVO* contourVO, double x, double y, int limInfX, int limInfY, int limSupX, int limSupY)
1578 int coordBase = 1, diameterOne = 0, diameterTwo = 0;
1579 marIsocontour* cont = parsePolyDataToMarIsocontour(contourVO);
1580 generateVector(1,&coordBase,originX,originY,&x,&y);
1582 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1584 generateVector(0,&coordBase,originX,originY,&x,&y);
1585 if (pointInPolygon(cont,x,y))
1598 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1600 generateVector(35,&coordBase,originX,originY,&x,&y);
1601 if (pointInPolygon(cont,x,y))
1614 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1616 generateVector(17,&coordBase,originX,originY,&x,&y);
1617 if (pointInPolygon(cont,x,y))
1630 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1632 generateVector(107,&coordBase,originX,originY,&x,&y);
1633 if (pointInPolygon(cont,x,y))
1645 return (diameterOne + diameterTwo) / 2;
1648 // ----------------------------------------------------------------------------
1649 void marAxisCT::histogram(int point, kVolume* vol)
1652 vtkImageData* imagedata;
1653 int limInfX, limInfY, limInfZ;
1654 int limSupX, limSupY, limSupZ;
1659 char nomArchivo[30],temp[3];
1660 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1663 // for (int point = 0; point < quantContours.size(); point++)
1665 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
1666 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
1667 ((marParameters *) getParameters())->setLumenPercentage(86);
1668 updateLumenPercentage(point,vol);
1669 marIsocontour* contExt = parsePolyDataToMarIsocontour(wallCont);
1671 std::vector<double> intensity;
1672 std::vector<int> frecuencies;
1674 for (int i = 0; i < limSupX; i++)
1676 for (int j = 0; j < limSupY; j++)
1678 if(pointInPolygon(contExt,i,j))
1681 pos = searchData(intensity,imagedata->GetScalarComponentAsDouble(i,j,0,0));
1684 intensity.push_back( imagedata->GetScalarComponentAsDouble(i,j,0,0));
1685 frecuencies.push_back(1);
1698 strcpy(nomArchivo,"./histogram");
1700 // EED 06/06/2007 linux
1701 // itoa(point,temp,10);
1702 sprintf(temp,"%d",point);
1704 strcat(nomArchivo,temp);
1705 strcat(nomArchivo,".csv");
1709 archivo = fopen(nomArchivo,"w");
1710 for (int k = 0; k < frecuencies.size(); k++)
1713 fprintf(archivo,"%f;%d\n",intensity[k],frecuencies[k]);
1721 //imagedata->Delete();
1727 // ----------------------------------------------------------------------------
1728 int marAxisCT::searchData(std::vector<double> vec, double value)
1733 for (int i = 0; i< vec.size() && !found; i++)
1735 if (vec[i] == value)
1745 // ----------------------------------------------------------------------------
1746 void marAxisCT::setCalibration(bool calib)
1748 calibration = calib;
1751 // ----------------------------------------------------------------------------
1752 bool marAxisCT::getCalibration()
1757 // ----------------------------------------------------------------------------
1758 void marAxisCT::adjustContour(int point, kVolume* vol)
1761 double percentage = 0.05;
1762 double prevDifference = 0;
1763 if (quantContours[point + 1] == NULL)
1765 create2Dcontours(point + 1 ,vol);
1768 //1. Obtain actual signal
1769 double signal = maxSignal;
1770 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1772 //2. Obtain previous contour
1773 marIsocontour* prev = parsePolyDataToMarIsocontour(searchContour(marAxisContours::WALL, point + 1));
1775 //3. Adjust actual contour
1776 double newValue = signal;
1777 wallCont->getIsocontour()->SetValue(0,newValue);
1778 wallCont->getIsocontour()->Update();
1779 wallCont->get2DContour()->Update();
1782 marIsocontour* actual = parsePolyDataToMarIsocontour(wallCont);
1785 double prevA = obtainContourArea(prev);
1786 double actuA = obtainContourArea(actual);
1788 prevDifference = fabs(prevA - actuA);
1791 if (prevA - actuA > 0 && prevA - actuA < percentage*prevA) //CASO 1
1794 while (!stop && newValue > 0)
1798 wallCont->getIsocontour()->SetValue(0,newValue);
1799 wallCont->getIsocontour()->Update();
1800 wallCont->get2DContour()->Update();
1802 actual = parsePolyDataToMarIsocontour(wallCont);
1803 actuA = obtainContourArea(actual);
1805 if (fabs(prevA - actuA) > percentage*prevA)
1808 wallCont->getIsocontour()->SetValue(0,newValue);
1809 wallCont->getIsocontour()->Update();
1810 wallCont->get2DContour()->Update();
1816 else if (prevA - actuA < 0 && prevA - actuA > -percentage*prevA) //CASO 2
1820 while (!stop && newValue > 0)
1824 wallCont->getIsocontour()->SetValue(0,newValue);
1825 wallCont->getIsocontour()->Update();
1826 wallCont->get2DContour()->Update();
1828 actual = parsePolyDataToMarIsocontour(wallCont);
1829 actuA = obtainContourArea(actual);
1831 if (prevA - actuA < -percentage*prevA)
1834 wallCont->getIsocontour()->SetValue(0,newValue);
1835 wallCont->getIsocontour()->Update();
1836 wallCont->get2DContour()->Update();
1841 else if (prevA - actuA < 0 && prevA - actuA < -percentage*prevA) //CASO 3
1850 wallCont->getIsocontour()->SetValue(0,newValue);
1851 wallCont->getIsocontour()->Update();
1852 wallCont->get2DContour()->Update();
1854 actual = parsePolyDataToMarIsocontour(wallCont);
1856 actuA = obtainContourArea(actual);
1857 if (prevA - actuA > -percentage*prevA)
1859 if (prevDifference < fabs(prevA - actuA))
1862 wallCont->getIsocontour()->SetValue(0,newValue);
1863 wallCont->getIsocontour()->Update();
1864 wallCont->get2DContour()->Update();
1870 prevDifference = fabs(prevA - actuA);
1873 else if (prevA - actuA > 0 && prevA - actuA > percentage*prevA) //CASO 4
1876 while (!stop && newValue > 0)
1880 wallCont->getIsocontour()->SetValue(0,newValue);
1881 wallCont->getIsocontour()->Update();
1882 wallCont->get2DContour()->Update();
1884 actual = parsePolyDataToMarIsocontour(wallCont);
1885 actuA = obtainContourArea(actual);
1887 if (prevA - actuA < percentage*prevA)
1889 if (prevDifference < fabs(prevA - actuA))
1892 wallCont->getIsocontour()->SetValue(0,newValue);
1893 wallCont->getIsocontour()->Update();
1894 wallCont->get2DContour()->Update();
1900 prevDifference = fabs(prevA - actuA);
1905 maxSignal = newValue;
1907 this->markUpLumen(point,vol);
1910 // ----------------------------------------------------------------------------
1911 int marAxisCT::getStartIndex()
1916 // ----------------------------------------------------------------------------
1917 void marAxisCT::setStartIndex(int start)
1921 int per = ((marParameters *) getParameters())->getLumenPercentage();
1922 double oldValue = quantContours[start]->getContour(0)->getSignal()*per*pow(10.0,-2.0);
1924 maxSignal = oldValue; // / (per*pow(10.0,-2.0));
1927 // ----------------------------------------------------------------------------
1928 int marAxisCT::getPointSize()
1930 return vesselPoints.size();
1933 // ----------------------------------------------------------------------------
1934 marPoint* marAxisCT::getPoint(int i)
1936 return vesselPoints[i];
1939 // ----------------------------------------------------------------------------
1940 void marAxisCT::extractLumen(vtkImageData *lumenImage, marIsocontour *lumenContour, int point)
1944 vtkImageCast* cast = vtkImageCast::New();
1945 cast->SetInput(lumenImage);
1946 cast->SetOutputScalarTypeToDouble();
1948 vtkImageGaussianSmooth* filter = vtkImageGaussianSmooth::New();
1949 filter->SetDimensionality(2);
1950 filter->SetInput(cast->GetOutput());
1951 filter->SetStandardDeviations((getParameters())->getStandardDeviation(),(getParameters())->getStandardDeviation());
1952 filter->SetRadiusFactors((getParameters())->getRadius(),(getParameters())->getRadius());
1953 filter->SetStandardDeviations(3,3);
1954 filter->SetRadiusFactors(3,3);
1958 marAxisContours* axisCont = quantContours[point];
1962 for (int i = 0; i < axisCont->getSize() && pos == -1; i++)
1964 if (axisCont->getContourType(i) == marAxisContours::ELUMEN)
1971 if (pos != -1) //Verify that contour hasn't been replaced
1974 if (axisCont->isReplaced(pos))
1983 marContourVO* contourVO = new marContourVO();
1984 marIsocontour* iso = lumenContour;
1986 contourVO->setIsocontour(vtkContourFilter::New());
1987 contourVO->getIsocontour()->SetInput( filter->GetOutput() );
1988 contourVO->getIsocontour()->SetNumberOfContours( 1 );
1989 contourVO->getIsocontour()->SetValue( 0, 128 );
1990 contourVO->getIsocontour()->Update( );
1992 contourVO->setIsocontourCpd(vtkCleanPolyData::New());
1993 contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
1994 contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
1995 contourVO->getIsocontourCpd()->Update( );
1999 iso->getCG(&cgx,&cgy);
2001 contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
2002 contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
2004 contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
2005 contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
2006 contourVO->getIsocontourDcf()->Update( );
2008 contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
2009 contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
2010 contourVO->getIsocontourCpd2()->Update();
2012 contourVO->setIsocontourStripped(vtkStripper::New( ));
2013 contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
2014 contourVO->getIsocontourStripped()->Update();
2016 contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
2017 contourVO->get2DContour()->Update();
2018 contourVO->setSignal(255);
2019 contourVO->setType(marAxisContours::ELUMEN);
2023 marIsocontour* actual = parsePolyDataToMarIsocontour(contourVO);
2024 marContourVO* aVO = searchContour(marAxisContours::WALL,point);
2025 double wallArea = 0;
2028 wallArea = obtainContourArea(parsePolyDataToMarIsocontour(aVO));
2031 double prevA = obtainContourArea(actual);
2034 int finalValue = 128;
2036 while (newValue > 0 && actuA > 0)
2043 contourVO->getIsocontour()->SetValue(0,newValue);
2044 contourVO->getIsocontour()->Update();
2045 contourVO->get2DContour()->Update();
2046 actual = parsePolyDataToMarIsocontour(contourVO);
2047 actuA = obtainContourArea(actual);
2048 if ((actuA /1000) < 1 && (prevA / 1000) > 1)
2050 finalValue = newValue;
2053 else if (actuA > prevA && (prevA / 1000) < 1 && actuA < 0.1*wallArea)
2055 finalValue = newValue;
2059 catch (std::exception e)
2066 contourVO->getIsocontour()->SetValue(0,finalValue);
2067 contourVO->getIsocontour()->Update();
2068 contourVO->get2DContour()->Update();
2076 axisCont->replaceContour(contourVO,pos);
2080 axisCont->addContour(contourVO);
2082 // axisCont->addContour(contourVO);
2085 quantContours[point] = axisCont;
2091 // ----------------------------------------------------------------------------
2092 void marAxisCT::generateFile(int point,kVolume* vol)
2095 char nomArchivo[30], temp[3];
2097 strcpy(nomArchivo,"./point");
2099 // EED 06/06/2007 linux
2100 // itoa(point,temp,10);
2101 sprintf(temp,"%d",point);
2103 strcat(nomArchivo,temp);
2104 strcat(nomArchivo,".txt");
2105 archivo = fopen(nomArchivo,"w");
2107 fprintf(archivo, "%d", quantContours[point]->getSize());
2110 for (i = 0; i < quantContours[point]->getSize(); i++)
2112 marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
2113 fprintf(archivo, " %d",iso->getSize());
2116 fprintf(archivo,"\n");
2118 for ( i = 0; i < quantContours[point]->getSize(); i++)
2120 marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
2121 fprintf(archivo,"Tipo: %d\n",quantContours[point]->getContourType(i));
2122 for (int j = 0; j < iso->getSize(); j++)
2124 double tempX, tempY;
2125 iso->getPoint(j,&tempX,&tempY);
2126 fprintf(archivo,"%f;%f;%d\n",tempX,tempY,point);
2130 // EED 30 Janvier 2007
2131 vtkImageData *vtkimagedata = this->getSliceImage(point,vol);
2133 std::string directory = "./";
2134 std::string filename = "xx";
2135 float rescalaSlope = 1;
2136 float rescalaIntercept = 0;
2138 vtkimagedata->GetDimensions(dim);
2146 marRAW2Files marraw2;
2147 marraw2.saveVolume(directory,filename,vtkimagedata,voi,rescalaSlope,rescalaIntercept);
2155 // ----------------------------------------------------------------------------
2156 void marAxisCT::replaceContour2D(int point,int size,double *vx,double *vy, int type)
2159 vtkPoints *_pts = vtkPoints::New();
2160 _pts->SetNumberOfPoints(size);
2163 for (j=0 ; j<size ; j++){
2164 _pts->SetPoint(j, vx[j] , vy[j] , 0 );
2166 // _pts->SetPoint(0, vx[0] , vy[0] , 0 );
2168 vtkCellArray *lines = vtkCellArray::New();
2169 lines->InsertNextCell( size );
2170 for ( j=0 ; j<size+1 ; j++ ){
2171 lines->InsertCellPoint(j % size );
2174 vtkPolyData *_pd = vtkPolyData::New();
2175 _pd->SetPoints( _pts );
2176 _pd->SetLines( lines );
2177 lines->Delete(); //do not delete lines ??
2180 marContourVO* vo = new marContourVO();///this->searchContour(type,point);
2181 vo->setReplaced(true);
2182 vo->set2DContour(_pd);
2184 quantContours[point]->addContour(vo); //->replaceContour(vo,i);
2185 /* _2Dcontours[i]=_pd;
2186 createContour(i,NULL);*/
2189 // ----------------------------------------------------------------------------
2190 void marAxisCT::cleanContours(int type, int point)
2192 marAxisContours* newContours = new marAxisContours();
2193 if (quantContours[point] != NULL)
2195 for (int i = 0; i < quantContours[point]->getSize(); i++)
2197 if (quantContours[point]->getContourType(i) != type)
2199 newContours->addContour(quantContours[point]->getContour(i));
2203 /* marContourVO* vo = quantContours[point]->getContour(i);
2204 vo->set2DContour(NULL);
2205 newContours->replaceContour(vo,i);*/
2210 quantContours[point] = newContours;
2213 // ----------------------------------------------------------------------------
2214 marContourVO* marAxisCT::searchContour(int type, int point)
2216 marContourVO* cont = NULL;
2218 for (int i = 0; i < quantContours[point]->getSize(); i++)
2220 if (quantContours[point]->getContourType(i) == type)
2222 cont = quantContours[point]->getContour(i);
2229 // ----------------------------------------------------------------------------
2230 double marAxisCT::performXOR(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2232 if (manual.size() == 0)
2237 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2238 std::vector <marIsocontour* > polygons;
2239 vtkImageData* imagedata;
2241 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
2242 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2243 vtkImageData* imageA = vtkImageData::New();
2244 imageA->SetDimensions(limSupX,limSupY,0);
2245 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2246 imageA->SetScalarTypeToUnsignedShort();
2247 imageA->AllocateScalars();
2250 vtkImageData* imageM = vtkImageData::New();
2251 imageM->SetDimensions(limSupX,limSupY,0);
2252 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2253 imageM->SetScalarTypeToUnsignedShort();
2254 imageM->AllocateScalars();
2258 for (int i = 0; i < quantContours[point]->getSize(); i++)
2260 if (quantContours[point]->getContourType(i) == type)
2262 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2267 if (polygons.size() == 0)
2272 for (x = limInfX; x < (limSupX - 1); x++)
2274 for (y = limInfY; y < (limSupY - 1); y++)
2276 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2277 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2282 for (numConts = 0; numConts < polygons.size(); numConts++)
2284 if (pointInPolygon(polygons[numConts],x,y))
2290 for (numConts = 0; numConts < manual.size(); numConts++)
2292 if (pointInPolygon(manual[numConts],x,y))
2300 /* vtkImageLogic* xor = vtkImageLogic::New();
2301 xor->SetInput1(imageM);
2302 xor->SetInput2(imageA);
2303 xor->SetOutputTrueValue(255);
2304 xor->SetOperationToXor();
2306 vtkImageCast* cast = vtkImageCast::New();
2307 cast->SetInput(xor->GetOutput());
2308 cast->SetOutputScalarTypeToDouble();
2310 int coincidencias = 0;
2311 for (x = limInfX; x < (limSupX - 1); x++)
2313 for (y = limInfY; y < (limSupY - 1); y++)
2315 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2316 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2318 int xor_ = compA^compM;
2326 double resp = (double)coincidencias ;
2332 // ----------------------------------------------------------------------------
2333 double marAxisCT::performAND(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2336 if (manual.size() == 0)
2341 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2342 std::vector <marIsocontour* > polygons;
2343 vtkImageData* imagedata;
2345 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
2346 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2347 vtkImageData* imageA = vtkImageData::New();
2348 imageA->SetDimensions(limSupX,limSupY,0);
2349 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2350 imageA->SetScalarTypeToUnsignedShort();
2351 imageA->AllocateScalars();
2354 vtkImageData* imageM = vtkImageData::New();
2355 imageM->SetDimensions(limSupX,limSupY,0);
2356 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2357 imageM->SetScalarTypeToUnsignedShort();
2358 imageM->AllocateScalars();
2361 for (int i = 0; i < quantContours[point]->getSize(); i++)
2363 if (quantContours[point]->getContourType(i) == type)
2365 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2370 if (polygons.size() == 0)
2376 for (x = limInfX; x < (limSupX - 1); x++)
2378 for (y = limInfY; y < (limSupY - 1); y++)
2380 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2381 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2386 for (numConts = 0; numConts < polygons.size(); numConts++)
2388 if (pointInPolygon(polygons[numConts],x,y))
2394 for (numConts = 0; numConts < manual.size(); numConts++)
2396 if (pointInPolygon(manual[numConts],x,y))
2404 /* vtkImageLogic* and = vtkImageLogic::New();
2405 and->SetInput1(imageM);
2406 and->SetInput2(imageA);
2407 and->SetOutputTrueValue(255);
2408 and->SetOperationToAnd();
2410 int coincidencias = 0;
2411 for (x = limInfX; x < (limSupX - 1); x++)
2413 for (y = limInfY; y < (limSupY - 1); y++)
2415 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2416 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2418 int and_ = compA & compM;
2426 double resp = (double)coincidencias ;
2432 // ----------------------------------------------------------------------------
2433 marIsocontour* marAxisCT::loadMarIsocontour(int size, double *vx, double *vy)
2435 marIsocontour* cont = new marIsocontour();
2437 for ( int j=0 ; j<size ; j++ )
2439 cont->insertPoint(vx[j], vy[j]);
2445 // ----------------------------------------------------------------------------
2446 double marAxisCT::performUnion(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2448 if (manual.size() == 0)
2453 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2454 std::vector <marIsocontour* > polygons;
2455 vtkImageData* imagedata;
2457 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
2458 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2459 vtkImageData* imageA = vtkImageData::New();
2460 imageA->SetDimensions(limSupX,limSupY,0);
2461 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2462 imageA->SetScalarTypeToUnsignedShort();
2463 imageA->AllocateScalars();
2466 vtkImageData* imageM = vtkImageData::New();
2467 imageM->SetDimensions(limSupX,limSupY,0);
2468 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
2469 imageM->SetScalarTypeToUnsignedShort();
2470 imageM->AllocateScalars();
2474 for (i = 0; i < quantContours[point]->getSize(); i++)
2476 if (quantContours[point]->getContourType(i) == type)
2478 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2483 if (polygons.size() == 0)
2489 for ( x = limInfX; x < (limSupX - 1); x++)
2491 for ( y = limInfY; y < (limSupY - 1); y++)
2493 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2494 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2498 for ( numConts = 0; numConts < polygons.size(); numConts++)
2500 if (pointInPolygon(polygons[numConts],x,y))
2506 for (numConts = 0; numConts < manual.size(); numConts++)
2508 if (pointInPolygon(manual[numConts],x,y))
2516 /* vtkImageLogic* and = vtkImageLogic::New();
2517 and->SetInput1(imageM);
2518 and->SetInput2(imageA);
2519 and->SetOutputTrueValue(255);
2520 and->SetOperationToAnd();
2522 int coincidencias = 0;
2523 for (x = limInfX; x < (limSupX - 1); x++)
2525 for ( y = limInfY; y < (limSupY - 1); y++)
2527 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2528 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2530 int or_ = compA | compM;
2538 double resp = (double)coincidencias;