00001
00002
00004
00005
00006 #include <vtkImageCast.h>
00007 #include <vtkImageGaussianSmooth.h>
00008 #include <vtkCleanPolyData.h>
00009 #include <vtkPolyDataConnectivityFilter.h>
00010 #include <vtkStripper.h>
00011 #include <vtkCell.h>
00012 #include <vtkIdList.h>
00013 #include <vtkMatrix4x4.h>
00014 #include <vtkTransform.h>
00015 #include <vtkPlaneSource.h>
00016 #include <vtkContourFilter.h>
00017 #include <marLine.h>
00018 #include <marUtils.h>
00019 #include <vtkImageContinuousErode3D.h>
00020 #include <vtkImageThreshold.h>
00021 #include <vtkImageSobel2D.h>
00022 #include <vtkCellArray.h>
00023 #include <vtkImageLogic.h>
00024
00025
00026 #include "marAxisCT.h"
00027 #include "marGdcmDicom.h"
00028
00029
00030
00031
00032 marAxisCT::marAxisCT( ) : marAxis( )
00033 {
00034 }
00035
00036
00037 marContour* marAxisCT::getContour( int point , kVolume* vol, int index ) {
00038 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->getContour() == NULL)
00039 {
00040 createContours(point,vol);
00041 }
00042
00043 return quantContours[point]->getContour(index)->getContour();
00044 }
00045
00046
00047 vtkPoints* marAxisCT::get3Dcontour( int point , kVolume* vol, int index ) {
00048 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get3DContour() == NULL)
00049 {
00050 create3Dcontours(point,vol);
00051 }
00052
00053 marAxisContours* cont = quantContours[point];
00054
00055 return cont->getContour(index)->get3DContour();
00056 }
00057
00058
00059 vtkPolyData* marAxisCT::get2Dcontour( int point , kVolume* vol, int index ) {
00060 if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DContour() == NULL)
00061 {
00062 create2Dcontours(point,vol);
00063 }
00064
00065 marAxisContours* cont = quantContours[point];
00066
00067
00068
00069 return cont->getContour(index)->get2DContour();
00070 }
00071
00072
00073 vtkPoints* marAxisCT::get2DDiameterMin( int point , kVolume* vol, int index ) {
00074
00075
00076 create2DDiametersMin(point,vol);
00077
00078
00079 marAxisContours* cont = quantContours[point];
00080
00081 return cont->getContour(index)->get2DDiameterMin();
00082 }
00083
00084
00085 vtkPoints* marAxisCT::get2DDiameterMax( int point , kVolume* vol, int index ) {
00086
00087
00088 create2DDiametersMax(point,vol);
00089
00090
00091 marAxisContours* cont = quantContours[point];
00092
00093 return cont->getContour(index)->get2DDiameterMax();
00094 }
00095
00096
00097 int marAxisCT::getSignal(int point, int index, kVolume* vol) {
00098 if (quantContours[point] == NULL )
00099 {
00100 createSignals(point,vol);
00101 }
00102
00103 return (int) ( quantContours[point]->getContour(index)->getSignal() );
00104 }
00105
00106
00107 void marAxisCT::createContours( int point, kVolume* vol ) {
00108
00109 int j,id;
00110 int numberOfPoints,numberOfCells;
00111 vtkCell* cell;
00112 vtkIdList* pids;
00113 double p[ 3 ];
00114 vtkPolyData* vtkPolyData_2DContour;
00115 double x,y;
00116
00117 int sizeIma = getParameters( )->getSizeIma( );
00118 double dimIma = getParameters( )->getDimIma( );
00119
00120 if (quantContours.size() < point - 1 )
00121 {
00122 create2Dcontours(point,vol);
00123 }
00124
00125 marAxisContours* axisCont = quantContours[point];
00126
00127 for (int i = 0; i < axisCont->getSize(); i++)
00128 {
00129 marContourVO* aContourVO = axisCont->getContour(i);
00130 aContourVO->setContour(new marContour( point, getParameters( ) ));
00131 vtkPolyData_2DContour = aContourVO->get2DContour();
00132 numberOfCells = vtkPolyData_2DContour->GetNumberOfCells( );
00133 cell = vtkPolyData_2DContour->GetCell( 0 );
00134 pids = cell->GetPointIds( );
00135 numberOfPoints = pids->GetNumberOfIds( );
00136 for( j = 0; j < numberOfPoints; j++ )
00137 {
00138 id = pids->GetId( j );
00139 vtkPolyData_2DContour->GetPoint( id, p);
00140 x=p[0]-64.0;
00141 y=p[1]-64.0;
00142 x=x * dimIma / ( float ) sizeIma;
00143 y=y * dimIma / ( float ) sizeIma;
00144 aContourVO->getContour()->addContourPoint( x , y );
00145 }
00146
00147 aContourVO->getContour()->do_spline();
00148 aContourVO->getContour()->calculateVariables();
00149 }
00150
00151 }
00152
00153
00154 void marAxisCT::create3Dcontours( int point , kVolume* vol ) {
00155
00156 if (quantContours[point] == NULL )
00157 {
00158 createContours(point,vol);
00159 }
00160
00161 marAxisContours* axisCont = quantContours[point];
00162
00163 for (int i = 0; i < axisCont->getSize(); i++)
00164 {
00165 vtkPoints *_vtkPoints;
00166 double *o;
00167 double *n;
00168 vtkMatrix4x4* mat;
00169 vtkTransform* trans;
00170 double nn[3];
00171 double c[3],p1[3],p2[3];
00172 double nA,nB,nC;
00173
00174 _vtkPoints = vtkPoints::New();
00175 o = getPoints(point);
00176 n = getNormal( i );
00177 mat = vtkMatrix4x4::New();
00178 trans = vtkTransform::New();
00179
00180
00181 nn[0]=n[0]; nn[1]=n[1]; nn[2]=n[2];
00182 nC=sqrt( nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2] );
00183 nn[0]=nn[0]/nC; nn[1]=nn[1]/nC; nn[2]=nn[2]/nC;
00184
00185 vtkPlaneSource* pSource = vtkPlaneSource::New( );
00186 pSource->SetOrigin( 0, 0 , 0 );
00187 pSource->SetPoint1( 1, 0 , 0 );
00188 pSource->SetPoint2( 0, 0 , 1.0 );
00189 pSource->SetCenter( 0,0,0 );
00190 pSource->SetNormal( nn[ 0 ], nn[ 1 ], nn[ 2 ] );
00191 pSource->Update( );
00192 pSource->Update( );
00193 pSource->GetOrigin(c);
00194 pSource->GetPoint1(p1);
00195 pSource->GetPoint2(p2);
00196 pSource->Delete();
00197 p1[0]=p1[0]-c[0]; p1[1]=p1[1]-c[1]; p1[2]=p1[2]-c[2];
00198 p2[0]=p2[0]-c[0]; p2[1]=p2[1]-c[1]; p2[2]=p2[2]-c[2];
00199 nA=sqrt( p1[0]*p1[0] + p1[1]*p1[1] + p1[2]*p1[2] );
00200 nB=sqrt( p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2] );
00201 p1[0]=p1[0]/nA; p1[1]=p1[1]/nA; p1[2]=p1[2]/nA;
00202 p2[0]=p2[0]/nB; p2[1]=p2[1]/nB; p2[2]=p2[2]/nB;
00203
00204 mat->SetElement (0,0, nn[0]);
00205 mat->SetElement (1,0, nn[1]);
00206 mat->SetElement (2,0, nn[2]);
00207 mat->SetElement (3,0, 0);
00208 mat->SetElement (0,1, p2[0]);
00209 mat->SetElement (1,1, p2[1]);
00210 mat->SetElement (2,1, p2[2]);
00211 mat->SetElement (3,1, 0);
00212 mat->SetElement (0,2, p1[0]);
00213 mat->SetElement (1,2, p1[1]);
00214 mat->SetElement (2,2, p1[2]);
00215 mat->SetElement (3,2, 0);
00216 mat->SetElement (0,3, 0);
00217 mat->SetElement (1,3, 0);
00218 mat->SetElement (2,3, 0);
00219 mat->SetElement (3,3, 1);
00220
00221 double deter=mat->Determinant(mat);
00222 trans->Identity();
00223 trans->SetMatrix(mat);
00224 float ang;
00225 ang=-90;
00226 trans->RotateWXYZ ( ang,0,1,0);
00227 trans->Update();
00228
00229 vtkMatrix4x4 *m=vtkMatrix4x4::New();
00230 trans->GetMatrix ( m );
00231
00232 int j,numberOfPoints;
00233
00234 marContour* contour2D = getContour(point,vol,i);
00235 marContourVO* aContourVO = axisCont->getContour(i);
00236 numberOfPoints = contour2D->getNumberOfPoints( );
00237 double fpA[3];
00238 double pA[3],ppp[3];
00239
00240 for( j = 0; j <= numberOfPoints; j++ ) {
00241 contour2D->getPoint( fpA,j % numberOfPoints);
00242 pA[0]=fpA[0];
00243 pA[1]=fpA[1];
00244 pA[2]=0;
00245
00246 ppp[0] = m->GetElement(0,0)*pA[0] + m->GetElement(0,1)*pA[1] + m->GetElement(0,2)*pA[2] + m->GetElement(0,3);
00247 ppp[1] = m->GetElement(1,0)*pA[0] + m->GetElement(1,1)*pA[1] + m->GetElement(1,2)*pA[2] + m->GetElement(1,3);
00248 ppp[2] = m->GetElement(2,0)*pA[0] + m->GetElement(2,1)*pA[1] + m->GetElement(2,2)*pA[2] + m->GetElement(2,3);
00249
00250 ppp[0]=ppp[0]+o[0]-c[0]; ppp[1]=ppp[1]+o[1]-c[1]; ppp[2]=ppp[2]+o[2]-c[2];
00251 _vtkPoints->InsertNextPoint( (float)ppp[0],(float)ppp[1],(float)ppp[2] );
00252
00253
00254 }
00255
00256
00257 m->Delete();
00258 mat->Delete();
00259 trans->Delete();
00260
00261 if (aContourVO->get3DContour()!=NULL){ aContourVO->get3DContour()->Delete(); }
00262 aContourVO->set3DContour(_vtkPoints);
00263
00264 }
00265
00266
00267
00268
00269 }
00270
00271
00272 void marAxisCT::create2Dcontours( int point , kVolume* vol ) {
00273 std::vector <marIsocontour *> contours;
00274
00275
00276
00277 generatePoints(point,vol,&contours);
00278
00279 vtkImageData* imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
00280 marAxisContours* axisCont = new marAxisContours();
00281
00282 for (int i = 0; i < contours.size(); i++)
00283 {
00284 marContourVO* contourVO = new marContourVO();
00285 marIsocontour* iso = contours[i];
00286
00287 contourVO->setIsocontour(vtkContourFilter::New());
00288 contourVO->getIsocontour()->SetInput( imagedata );
00289 contourVO->getIsocontour()->SetNumberOfContours( 1 );
00290 contourVO->getIsocontour()->SetValue( 0, iso->getMaxIntensity() );
00291 contourVO->getIsocontour()->Update( );
00292
00293 contourVO->setIsocontourCpd(vtkCleanPolyData::New());
00294 contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
00295 contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
00296 contourVO->getIsocontourCpd()->Update( );
00297
00298 double cgx, cgy;
00299
00300 iso->getCG(&cgx,&cgy);
00301
00302 contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
00303 contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
00304
00305 contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
00306 contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
00307 contourVO->getIsocontourDcf()->Update( );
00308
00309 contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
00310 contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
00311 contourVO->getIsocontourCpd2()->Update();
00312
00313 contourVO->setIsocontourStripped(vtkStripper::New( ));
00314 contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
00315 contourVO->getIsocontourStripped()->Update();
00316
00317 contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
00318 contourVO->get2DContour()->Update();
00319 contourVO->setSignal(iso->getMaxIntensity());
00320
00321 if (iso->getType() == marAxisContours::LUMEN)
00322 {
00323 contourVO->setType(marAxisContours::WALL);
00324 }
00325 else
00326 {
00327 contourVO->setType(iso->getType());
00328 }
00329
00330
00331
00332 axisCont->addContour(contourVO);
00333
00334
00335 delete iso;
00336 contours[i] = NULL;
00337 }
00338
00339
00340
00341 quantContours[point] = axisCont;
00342
00343
00344
00345 updateLumenPercentage(point,vol);
00346
00347 updateCalcPercentage(point,vol);
00348
00349 if (calibration && point < startIndex)
00350 {
00351 this->adjustContour(point,vol);
00352 }
00353 else
00354 {
00355 adjustWall(point,vol);
00356 }
00357 contours.clear();
00358 imagedata->Delete();
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381 }
00382
00383
00384 void marAxisCT::create2DDiametersMin(int point , kVolume* vol ) {
00385
00386 if (quantContours[point] == NULL )
00387 {
00388 createContours(point,vol);
00389 }
00390
00391 marAxisContours* axisCont = quantContours[point];
00392
00393 for (int i = 0; i < axisCont->getSize(); i++)
00394 {
00395 double p1[2],p2[2];
00396 marContour *marcontour=getContour(point,vol,i);
00397 marContourVO* aContourVO = axisCont->getContour(i);
00398
00399 marcontour->getMinimumLine(p1,p2);
00400 vtkPoints *_vtkPoints = vtkPoints::New();
00401
00402 int sizeIma = getParameters( )->getSizeIma( );
00403 double dimIma = getParameters( )->getDimIma( );
00404 double factor=( float ) sizeIma / dimIma;
00405
00406 p1[0]=p1[0]*factor+sizeIma/2;
00407 p1[1]=p1[1]*factor+sizeIma/2;
00408 p2[0]=p2[0]*factor+sizeIma/2;
00409 p2[1]=p2[1]*factor+sizeIma/2;
00410 _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
00411 _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
00412
00413 if ( aContourVO->get2DDiameterMin()!=NULL ) { aContourVO->get2DDiameterMin()->Delete(); }
00414 aContourVO->set2DDiameterMin(_vtkPoints);
00415
00416 }
00417
00418 }
00419
00420
00421 void marAxisCT::create2DDiametersMax(int point , kVolume* vol ) {
00422 if (quantContours[point] == NULL )
00423 {
00424 createContours(point,vol);
00425 }
00426
00427 marAxisContours* axisCont = quantContours[point];
00428
00429 for (int i = 0; i < axisCont->getSize(); i++)
00430 {
00431 double p1[2],p2[2];
00432 marContour *marcontour=getContour(point,vol,i);
00433 marContourVO* aContourVO = axisCont->getContour(i);
00434
00435 marcontour->getMaximumLine(p1,p2);
00436 vtkPoints *_vtkPoints = vtkPoints::New();
00437
00438 int sizeIma = getParameters( )->getSizeIma( );
00439 double dimIma = getParameters( )->getDimIma( );
00440 double factor=( float ) sizeIma / dimIma;
00441 p1[0]=p1[0]*factor+sizeIma/2;
00442 p1[1]=p1[1]*factor+sizeIma/2;
00443 p2[0]=p2[0]*factor+sizeIma/2;
00444 p2[1]=p2[1]*factor+sizeIma/2;
00445 _vtkPoints->InsertNextPoint( p1[0] , p1[1] , 1 );
00446 _vtkPoints->InsertNextPoint( p2[0] , p2[1] , 1 );
00447
00448 if ( aContourVO->get2DDiameterMax()!=NULL ) { aContourVO->get2DDiameterMax()->Delete(); }
00449 aContourVO->set2DDiameterMax(_vtkPoints);
00450 }
00451 }
00452
00453
00454 void marAxisCT::createSignals (int point, kVolume* vol) {
00455 if (quantContours[point] == NULL )
00456 {
00457 create2Dcontours(point,vol);
00458
00459 }
00460 }
00461
00462
00463 int marAxisCT::getNumberOfContours(int point, kVolume* vol) {
00464 if (quantContours[point] == NULL )
00465 {
00466 create2Dcontours(point,vol);
00467
00468
00469 }
00470
00471 return quantContours[point]->getSize();
00472 }
00473
00474
00475 int marAxisCT::getContourType(int point, int index, kVolume* vol) {
00476 if (quantContours[point] == NULL )
00477 {
00478 create2Dcontours(point,vol);
00479 }
00480
00481 return quantContours[point]->getContourType(index);
00482
00483 }
00484
00485
00486 void marAxisCT::generatePoints(int point, kVolume* vol, std::vector<marIsocontour*> *list) {
00487
00488 vtkImageData* imagedata;
00489 int dirs = 72;
00490 double aPoint[ 3 ];
00491 double *nPoint = new double[ 3 ];
00492 double x, y, prevIntensity, threshold;
00493 bool found, estStart, estDivided, isNeighbor;
00494 int actualZone, startZone, coordBase, maxPoint, lumenCount, calcCount;
00495 std::vector <marPoint *> gradients;
00496 std::vector <double> average;
00497 std::vector <marIsocontour *> contoursLumen;
00498 std::vector <marIsocontour *> contoursCalc;
00499
00500 contoursLumen.push_back(NULL);
00501 contoursCalc.push_back(NULL);
00502
00503
00504
00505
00506 startZone = marAxisContours::LUMEN;
00507 actualZone = startZone;
00508 lumenCount = 0;
00509 calcCount = 0;
00510
00511 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
00512
00513 int limInfX, limInfY, limInfZ;
00514 int limSupX, limSupY, limSupZ;
00515 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
00516
00517 vtkImageCast* cast = vtkImageCast::New();
00518 cast->SetInput(imagedata);
00519 cast->SetOutputScalarTypeToDouble();
00520
00521 vtkImageGaussianSmooth* filter = vtkImageGaussianSmooth::New();
00522 filter->SetDimensionality(2);
00523 filter->SetInput(cast->GetOutput());
00524 filter->SetStandardDeviations((getParameters())->getStandardDeviation(),(getParameters())->getStandardDeviation());
00525 filter->SetRadiusFactors((getParameters())->getRadius(),(getParameters())->getRadius());
00526 filter->SetStandardDeviations(3,3);
00527 filter->SetRadiusFactors(3,3);
00528
00529
00530
00531
00532 nPoint = getPoints(point);
00533
00534 aPoint[ 0 ] = (limSupX - limInfX) / 2;
00535 aPoint[ 1 ] = (limSupY - limInfY) / 2;
00536 aPoint[ 2 ] = (limSupZ - limInfZ) / 2;
00537 estDivided = false;
00538 estStart = false;
00539
00540 for (int dir = 0; dir < dirs; dir++)
00541 {
00542
00543 x = aPoint[ 0 ];
00544 y = aPoint[ 1 ];
00545 coordBase = 1;
00546 prevIntensity = interpolate(x,y,imagedata);
00547 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
00548
00549
00550 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
00551 {
00552
00553
00554 double resp = interpolate(x,y,imagedata);
00555
00556 marPoint* p = new marPoint(x,y);
00557 p->setGradient(prevIntensity - resp);
00558 p->setDirection(dir);
00559 p->setIntensity(resp);
00560
00561 prevIntensity = resp;
00562
00563 gradients.push_back(p);
00564
00565
00566 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
00567
00568 }
00569
00570 int index = 0;
00571 found = false;
00572 x = aPoint[ 0 ];
00573 y = aPoint [ 1 ];
00574 threshold = ((marParameters *) getParameters())->getContourThresh();
00575 int mysize = gradients.size();
00576
00577
00578 while (index < gradients.size() && !found)
00579 {
00580
00581 if (actualZone == marAxisContours::LUMEN)
00582 {
00583 average.push_back(gradients[index]->getIntensity());
00584
00585
00586 if (fabs(gradients[index]->getGradient()) > threshold
00587 && fabs(gradients[index]->getGradient()) < 4000) {
00588
00589
00590
00591
00592
00593 if(gradients[index]->getGradient() > 0)
00594 maxPoint = getMaximumGrad(gradients,index,gradients.size(),threshold);
00595 else
00596 maxPoint = getMinimumGrad(gradients,index,gradients.size(),threshold);
00597
00598 index = maxPoint;
00599 x = gradients[maxPoint]->getX();
00600 y = gradients[maxPoint]->getY();
00601
00602
00603 isNeighbor = detectNeighbor(contoursLumen[lumenCount],dir,marAxisContours::LUMEN);
00604 if (!isNeighbor)
00605 {
00606 lumenCount++;
00607 }
00608
00609 marIsocontour* cont = addPointToContour(contoursLumen[lumenCount],false,marAxisContours::LUMEN,gradients[maxPoint]);
00610 contoursLumen[lumenCount] = cont;
00611
00612
00613 double actualInt = interpolate(x,y,imagedata);
00614 double auxX = x;
00615 double auxY = y;
00616 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&auxX,&auxY);
00617 double next = interpolate(auxX,auxY,imagedata);
00618 if (gradients[index]->getGradient() < 0 && (auxX > (limInfX+1) && auxX < (limSupX-1)) && (auxY > (limInfY+1) && auxY < (limSupY-1)))
00619 {
00620
00621 isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
00622 if (!isNeighbor)
00623 {
00624 calcCount++;
00625 contoursCalc.push_back(NULL);
00626 }
00627
00628 cont = addPointToContour(contoursCalc[calcCount],true,marAxisContours::CALCIFICATION,gradients[maxPoint]);
00629
00630 contoursCalc[calcCount] = cont;
00631
00632 actualZone = marAxisContours::CALCIFICATION;
00633
00634
00635
00636 if (dir == 0)
00637 {
00638 estStart = true;
00639 }
00640
00641 if ((dir == dirs - 1) && estStart)
00642 {
00643 estDivided = true;
00644 }
00645
00646 }
00647 else
00648 {
00649
00650
00651
00652
00653 found = true;
00654 actualZone = startZone;
00655 coordBase = 1;
00656 x = aPoint[ 0 ];
00657 y = aPoint[ 1 ];
00658 threshold = ((marParameters *) getParameters())->getContourThresh();
00659 }
00660
00661 }
00662 }
00663
00664 else
00665 {
00666 coordBase = 1;
00667 int calcIndex;
00668 double actualInt = interpolate(x,y,imagedata);
00669 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
00670 double next = interpolate(x,y,imagedata);
00671
00672
00673 if (actualInt < next)
00674 {
00675 while (actualInt < interpolate(x,y,imagedata) && (x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
00676 {
00677 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
00678 index++;
00679 }
00680
00681 }
00682
00683 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
00684
00685
00686 isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
00687 if (!isNeighbor)
00688 {
00689 calcCount++;
00690 contoursCalc.push_back(NULL);
00691 }
00692
00693
00694 if (index >= gradients.size())
00695 {
00696 calcIndex = gradients.size() - 1;
00697 }
00698 else
00699 {
00700 calcIndex = index;
00701 }
00702
00703 marIsocontour* cont = addPointToContour(contoursCalc[calcCount],false,marAxisContours::CALCIFICATION,gradients[calcIndex]);
00704 contoursCalc[calcCount] = cont;
00705
00706
00707
00708
00709 found = true;
00710 actualZone = startZone;
00711 x = aPoint[ 0 ];
00712 y = aPoint[ 1 ];
00713 threshold = ((marParameters *) getParameters())->getContourThresh();
00714
00715 }
00716
00717 index++;
00718
00719 if (index == gradients.size() && !found)
00720 {
00721 threshold = threshold - 1;
00722 index = 0;
00723 }
00724
00725 }
00726
00727
00728 for (int ind = 0; ind < gradients.size(); ind++)
00729 {
00730 marPoint* p = gradients[ind];
00731 gradients[ind] = NULL;
00732 delete p;
00733
00734 }
00735 gradients.clear();
00736
00737 }
00738
00739 FILE *archivo;
00740 char nomArchivo[30];
00741
00742 strcpy(nomArchivo,"./debug");
00743
00744
00745 strcat(nomArchivo,".txt");
00746 archivo = fopen(nomArchivo,"w");
00747
00748
00749
00750
00751
00752
00753
00754 double tempXi, tempYi;
00755 int polySides = contoursLumen[marAxisContours::LUMEN]->getSize();
00756
00757
00758 marIsocontour* alc = new marIsocontour();
00759 for (int m=0; m<polySides; m++) {
00760 contoursLumen[marAxisContours::LUMEN]->getPoint(m,&tempXi,&tempYi);
00761 alc = addPointToContour(lumenContour[point],false,marAxisContours::LUMEN,new marPoint(tempXi,tempYi));
00762 lumenContour[point] = alc;
00763 fprintf(archivo,"X:%f,Y:%f\n",tempXi,tempYi);
00764 }
00765
00766
00767
00768
00769 if (estDivided && calcCount > 0)
00770 {
00771 unifyDividedStructure(&contoursCalc);
00772 calcCount--;
00773 }
00774
00775
00776 double lumenInt = getStatistics(contoursLumen[marAxisContours::LUMEN],aPoint[ 0 ],aPoint[ 1 ], imagedata);
00777 contoursLumen[marAxisContours::LUMEN]->setMaxIntensity(lumenInt);
00778
00779 if (contoursCalc[0] != NULL)
00780 {
00781 for (int i = 0; i < contoursCalc.size(); i++)
00782 {
00783 lumenCount++;
00784 double maxInt = getCalcStatistics(contoursCalc[i],aPoint[ 0 ],aPoint[1],imagedata,i-1);
00785 contoursCalc[i]->setMaxIntensity(maxInt);
00786 contoursLumen.push_back(NULL);
00787 contoursLumen[lumenCount] = contoursCalc[i];
00788 }
00789 }
00790
00791
00792 *list = contoursLumen;
00793
00794 contoursCalc.clear();
00795 contoursLumen.clear();
00796 marUtils* util = new marUtils();
00797 double avg = util->obtainAverage(average);
00798
00799
00800 fprintf(archivo,"tamaño: %d umbra: %d", (*list).size(),((marParameters *) getParameters())->getContourThresh());
00801 for (int k = 0; k < contoursLumen.size(); k++)
00802 {
00803 double cgx, cgy;
00804
00805 (*list)[k]->getCG(&cgx,&cgy);
00806 fprintf(archivo,"cx: %f cy:%f Señal: %f Type: %d\n", cgx, cgy, (*list)[k]->getMaxIntensity(), (*list)[k]->getType());
00807 }
00808
00809 fprintf(archivo,"PROMEDIO %f - STD: %f\n",avg,util->obtainStandardDeviation(average,avg));
00810 fprintf(archivo,"%d\n",lumenContour.size());
00811 fclose(archivo);
00812 imagedata->Delete();
00813
00814 }
00815
00816
00817 int marAxisCT::getMaximumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
00818
00819 double max = list[iniPos]->getGradient();
00820 int coordMax = iniPos;
00821 int cont = iniPos + 1;
00822
00823 while (cont < limit && list[cont]->getGradient()>=threshold){
00824 if (list[cont]->getGradient() > max){
00825 max = list[cont]->getGradient();
00826 coordMax = cont;
00827 }
00828 cont++;
00829 }
00830
00831 return coordMax;
00832 }
00833
00834
00835 int marAxisCT::getMinimumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
00836
00837 double min = list[iniPos]->getGradient();
00838 int coordMin = iniPos;
00839 int cont = iniPos + 1;
00840
00841 while (cont < limit && fabs(list[cont]->getGradient())>=threshold){
00842 if (list[cont]->getGradient() < min){
00843 min = list[cont]->getGradient();
00844 coordMin = cont;
00845 }
00846 cont++;
00847 }
00848
00849 return coordMin;
00850 }
00851
00852
00853 double marAxisCT::interpolate(double x, double y, vtkImageData* imagedata) {
00854
00855 double a,b,c,d;
00856 double intensity;
00857 int i,j;
00858
00859 i = (int)floor(x);
00860 j = (int)floor(y);
00861
00862
00863 c = imagedata->GetScalarComponentAsDouble(i,j,0,0) + (double) imagedata->GetScalarComponentAsDouble(i+1,j+1,0,0)
00864 - imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - imagedata->GetScalarComponentAsDouble(i,j+1,0,0);
00865
00866 b = imagedata->GetScalarComponentAsDouble(i,j+1,0,0) - c*i - imagedata->GetScalarComponentAsDouble(i,j,0,0);
00867
00868 a = imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - c*j - imagedata->GetScalarComponentAsDouble(i,j,0,0);
00869
00870 d = imagedata->GetScalarComponentAsDouble(i,j,0,0) - a*i - b*j - c*i*j;
00871
00872 intensity = a*x + b*y + c*x*y + d;
00873
00874 return intensity;
00875 }
00876
00877
00878 void marAxisCT::generateVector(int dir,int *coordBase,double originX,double originY,double *x,double *y) {
00879
00880 int dirs = 72;
00881 int angle = (dir * (360 / dirs));
00882
00883 *x = 1 * cos(angle *(3.1415/180)) + (*x);
00884 *y = 1 * sin(angle *(3.1415/180)) + (*y);
00885
00886 (*coordBase)++;
00887 }
00888
00889
00890 bool marAxisCT::detectNeighbor(marIsocontour* contour,int dir, int type) {
00891
00892 bool isNeighbor = false;
00893
00894 if (contour == NULL){
00895 isNeighbor = true;
00896 }else{
00897 for (int i = 0; i < contour->getSize(); i++) {
00898
00899 if (((contour->getDir(i) == (dir - 1)) || (contour->getDir(i) == dir))
00900 && (contour->getType() == type)){
00901 isNeighbor = true;
00902 }
00903 }
00904 }
00905
00906 return isNeighbor;
00907 }
00908
00909
00910 marIsocontour* marAxisCT::addPointToContour(marIsocontour* cont ,bool inside,int type,marPoint* point) {
00911
00912
00913
00914 if (cont == NULL)
00915 {
00916 cont = new marIsocontour();
00917 cont->setType(type);
00918 }
00919
00920 cont->insertPoint(point->getX(),point->getY());
00921 cont->setDir(cont->getSize() - 1,point->getDirection());
00922 cont->setInside(cont->getSize() - 1,inside);
00923
00924 return cont;
00925
00926 }
00927
00928
00929 void marAxisCT::unifyDividedStructure(std::vector<marIsocontour*> *contList) {
00930
00931 marIsocontour* structOne = (*contList)[contList->size()-1];
00932 marIsocontour* structTwo = (*contList)[0];
00933
00934
00935 for (int i = 0; i < structOne->getSize(); i++)
00936 {
00937 double x;
00938 double y;
00939 structOne->getPoint(i,&x,&y);
00940 structTwo->insertPoint(x,y);
00941 structTwo->setDir(structTwo->getSize() - 1,structOne->getDir(i));
00942 structTwo->setInside(structTwo->getSize() - 1,structOne->getInside(i));
00943
00944 }
00945
00946 contList->pop_back();
00947 }
00948
00949
00950 double marAxisCT::getStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata) {
00951
00952 double max, auxX, auxY;
00953 int basis;
00954 max = interpolate(x,y,imagedata);
00955
00956 for (int i = 0; i < contour->getSize(); i++)
00957 {
00958 basis = 1;
00959 auxX = x;
00960 auxY = y;
00961 double pointX;
00962 double pointY;
00963
00964 contour->getPoint(i,&pointX,&pointY);
00965 while (auxX != pointX || auxY != pointY)
00966 {
00967 double intensity = interpolate(auxX,auxY,imagedata);
00968 if (max < intensity)
00969 {
00970 max = intensity;
00971 }
00972 generateVector(contour->getDir(i),&basis,x,y,&auxX,&auxY);
00973
00974 }
00975 }
00976
00977 return max;
00978 }
00979
00980
00981
00982 double marAxisCT::getCalcStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata, int dir) {
00983
00984
00985 double max, auxX,auxY;
00986 int basis;
00987
00988 int i;
00989 FILE *archivo;
00990 char nomArchivo[30], temp[3];
00991
00992 strcpy(nomArchivo,"./statistics");
00993
00994
00995
00996 sprintf(temp,"%d",dir);
00997
00998 strcat(nomArchivo,temp);
00999 strcat(nomArchivo,".txt");
01000
01001
01002 archivo = fopen(nomArchivo,"w");
01003 i = 0;
01004 max = 0;
01005 fprintf(archivo,"TAMAÑO %d\n",contour->getSize());
01006
01007 while (i < contour->getSize())
01008 {
01009 basis = 1;
01010 double pointX, pointY;
01011 contour->getPoint(i+1,&pointX,&pointY);
01012 contour->getPoint(i,&auxX,&auxY);
01013 fprintf(archivo,"final: %f %f inicial %f %f \n",pointX,pointY,auxX,auxY);
01014
01015 while (auxX != pointX || auxY != pointY)
01016 {
01017 double intensity = interpolate(auxX,auxY,imagedata);
01018 if (max < intensity){
01019 max = intensity;
01020 }
01021 generateVector(contour->getDir(i),&basis,auxX,auxY,&auxX,&auxY);
01022 }
01023
01024 i +=2;
01025 }
01026 fclose(archivo);
01027 return max;
01028 }
01029
01030
01031 int marAxisCT::round(double num){
01032
01033 float decimal;
01034 int resp;
01035
01036 decimal = num - floor(num) * 1.0000;
01037
01038 if (decimal >= 0.50000)
01039 resp = (int)ceil(num);
01040 else
01041 resp = (int)floor(num);
01042
01043 return resp;
01044
01045 }
01046
01047
01048 void marAxisCT::createEmptyContours(){
01049
01050 int nCnts = ( int ) getNumberOfSplinePoints( );
01051
01052
01053 for (int i = 0; i < nCnts; i++)
01054 {
01055 quantContours.push_back( NULL );
01056 lumenContour.push_back( NULL);
01057 }
01058 }
01059
01060
01061 void marAxisCT::updateLumenPercentage(int point, kVolume* vol)
01062 {
01063 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
01064 if (quantContours[point] == NULL || wallCont->get2DContour() == NULL)
01065 {
01066 create2Dcontours(point,vol);
01067 }
01068
01069 if (wallCont->getType() == marAxisContours::WALL
01070 && wallCont->isReplaced() == false)
01071 {
01072 int a = ((marParameters *) getParameters())->getLumenPercentage();
01073 double oldValue = wallCont->getSignal();
01074 double lumenPercentage=((marParameters *) getParameters())->getLumenPercentage();
01075 double poww=pow(10.0,-2.0);
01076 double newValue = oldValue * lumenPercentage*poww;
01077 wallCont->getIsocontour()->SetValue(0,newValue);
01078 wallCont->getIsocontour()->Update();
01079 wallCont->get2DContour()->Update();
01080 }
01081
01082 }
01083
01084
01085 void marAxisCT::updateCalcPercentage(int point, kVolume* vol)
01086 {
01087 double oldValue;
01088 int posMax;
01089
01090 if (quantContours[point] == NULL)
01091 {
01092 create2Dcontours(point,vol);
01093 }
01094
01095 posMax = getMaxIntensity(point);
01096
01097 if (quantContours[point]->getSize() > 1)
01098 {
01099 oldValue = quantContours[point]->getContour(posMax)->getSignal();
01100 }
01101
01102 for (int i = 1; i < quantContours[point]->getSize(); i++)
01103 {
01104
01105 if (quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION
01106 && quantContours[point]->getContour(i)->isReplaced() == false)
01107 {
01108
01109 double newValue = oldValue * ((((marParameters *) getParameters())->getCalcPercentage())*pow(10.0,-2.0));
01110 quantContours[point]->getContour(i)->getIsocontour()->SetValue(0,newValue);
01111 quantContours[point]->getContour(i)->getIsocontour()->Update();
01112 quantContours[point]->getContour(i)->get2DContour()->Update();
01113 }
01114
01115 }
01116
01117 }
01118
01119
01120 int marAxisCT::getMaxIntensity(int point)
01121 {
01122 int pos = 1;
01123 double intMax;
01124 if (quantContours[point]->getSize() > 2)
01125 {
01126 intMax = quantContours[point]->getContour(0)->getSignal();
01127 for (int i = 0; i < quantContours[point]->getSize(); i++)
01128 {
01129 if (quantContours[point]->getContour(i)->getSignal() > intMax && quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION)
01130 {
01131 intMax = quantContours[point]->getContour(i)->getSignal() ;
01132 pos = i;
01133 }
01134 }
01135 }
01136
01137 return pos;
01138
01139 }
01140
01141 void marAxisCT::eraseContours()
01142 {
01143
01144
01145 int nCts = quantContours.size();
01146 vesselPoints.clear();
01147
01148 int i = 0;
01149
01150 while (i < nCts)
01151 {
01152 delete quantContours[i];
01153 quantContours[i] = NULL;
01154 lumenContour[i] = NULL;
01155 i++;
01156 }
01157 quantContours.clear();
01158 lumenContour.clear();
01159
01160
01161 }
01162
01163
01164 void marAxisCT::eraseContoursPartial(int start)
01165 {
01166
01167 int nCts = 0;
01168 int i = start;
01169
01170 while (i >= nCts)
01171 {
01172 delete quantContours[i];
01173 quantContours[i] = NULL;
01174 lumenContour[i] = NULL;
01175 i--;
01176 }
01177 quantContours.clear();
01178 lumenContour.clear();
01179
01180
01181 }
01182
01183
01184
01185 bool marAxisCT::pointInPolygon(marIsocontour *c, double x, double y){
01186 int i, j=0;
01187 bool oddNODES=false;
01188 int polySides = c->getSize();
01189 double *polyX, *polyY, tempX, tempY;
01190
01191 polyX = (double *) malloc(polySides*sizeof(double));
01192 polyY = (double *) malloc(polySides*sizeof(double));
01193
01194 for (i=0; i<polySides; i++) {
01195 c->getPoint(i,&tempX,&tempY);
01196 polyX[i] = tempX;
01197 polyY[i] = tempY;
01198 }
01199
01200 for (i=0; i<polySides; i++) {
01201 j++;
01202 if (j==polySides)
01203 j=0;
01204
01205 if (polyY[i]<y && polyY[j]>=y
01206 || polyY[j]<y && polyY[i]>=y) {
01207 if (polyX[i]+(y-polyY[i])/(polyY[j]-polyY[i])*(polyX[j]-polyX[i])<x) {
01208 oddNODES=!oddNODES;
01209 }
01210 }
01211 }
01212
01213 free(polyX);
01214 free(polyY);
01215 return oddNODES;
01216
01217 }
01218
01219
01220 marIsocontour* marAxisCT::parsePolyDataToMarIsocontour(marContourVO* contourVO){
01221
01222 vtkPolyData* polyDataCnt;
01223 vtkIdList* pointsId;
01224 double* pointPolyDataCnt;
01225 int numPointsPolyDataCnt, numTemp;
01226 marIsocontour *cont;
01227
01228 if (contourVO->isReplaced())
01229 {
01230 polyDataCnt = contourVO->get2DContour();
01231 }
01232 else
01233 {
01234 polyDataCnt = contourVO->getIsocontourStripped()->GetOutput();
01235 }
01236
01237 numPointsPolyDataCnt = polyDataCnt->GetNumberOfPoints();
01238 pointsId = polyDataCnt->GetCell(0)->GetPointIds();
01239 numTemp = pointsId->GetNumberOfIds();
01240
01241 cont = new marIsocontour();
01242
01243 for (int i = 0; i < (numTemp -1); i++){
01244 pointPolyDataCnt = polyDataCnt->GetPoint(pointsId->GetId(i));
01245 cont->insertPoint(pointPolyDataCnt[0], pointPolyDataCnt[1]);
01246
01247 }
01248
01249 return cont;
01250 }
01251
01252
01253
01254
01255 marIsocontour* marAxisCT::filterContour(marIsocontour* contExt, marIsocontour* contInt, double radio){
01256
01257 double tmp_eed;
01258 double x1,y1,x2,y2,xInt,yInt, xAux, yAux;
01259 std::vector<marLine*> lineas;
01260 std::vector<double> dists;
01261 int count, countAux, ultima;
01262 bool hayInterseccion;
01263
01264 marIsocontour* newCont = new marIsocontour();
01265
01266 for (int i = 0; i < contExt->getSize(); i++){
01267 contExt->getPoint(i,&x1,&y1);
01268 contInt->getPoint(i,&x2,&y2);
01269 lineas.push_back(new marLine(x1,y1,x2,y2));
01270 tmp_eed = sqrt(pow((x2-x1),2.0) + pow((y2-y1),2.0));
01271 dists.push_back( tmp_eed );
01272 }
01273
01274 count = 0;
01275 ultima = 0;
01276 while (count < lineas.size() - 1){
01277 countAux = 0;
01278
01279 marLine* l1 = lineas[count];
01280
01281 hayInterseccion = false;
01282 while (countAux < lineas.size()){
01283 marLine *l2 = lineas[countAux];
01284
01285 l1->getIntersect(l2->getA(),l2->getB(),l2->getC(),&xInt,&yInt);
01286
01287
01288 if (xInt > 0 && yInt > 0){
01289 contExt->getPoint(countAux,&xAux,&yAux);
01290 if (dists[count] >= sqrt(pow((xInt-xAux),2.0) + pow((yInt-yAux),2.0))){
01291 contExt->getPoint(count,&x1,&y1);
01292 contInt->getPoint(count,&x2,&y2);
01293
01294
01295 if ((dists[count] > sqrt(pow((xInt-x1),2.0) + pow((yInt-y1),2.0))) &&
01296 (dists[count] > sqrt(pow((xInt-x2),2.0) + pow((yInt-y2),2.0))))
01297 {
01298 ultima = countAux;
01299 hayInterseccion = true;
01300
01301 }
01302
01303 }
01304 }
01305 countAux++;
01306 }
01307
01308
01309 if (!hayInterseccion){
01310 contInt->getPoint(count,&x2,&y2);
01311 newCont->insertPoint(x2,y2);
01312 }
01313 count++;
01314
01315 }
01316
01317 return newCont;
01318 }
01319
01320
01321 void marAxisCT::markUpLumen(int point, kVolume* vol)
01322 {
01323
01324 marContourVO* aVO = this->searchContour(marAxisContours::ELUMEN,point);
01325
01326 if (aVO != NULL)
01327 {
01328 return;
01329 }
01330
01331 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
01332 vtkImageData* imagedata;
01333 vtkImageData* erodedata;
01334 std::vector <marIsocontour* > polygons;
01335 std::vector <marPoint* > points;
01336 double average = 0, deviation = 0;
01337 vtkImageData* image = vtkImageData::New();
01338 vesselPoints.clear();
01339
01340
01341 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
01342 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
01343 image->SetDimensions(limSupX,limSupY,0);
01344 image->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
01345 image->SetScalarTypeToUnsignedShort();
01346 image->AllocateScalars();
01347 image->Update();
01348
01349 vtkImageThreshold *imageThreshold = vtkImageThreshold::New();
01350 imageThreshold->SetInput(imagedata);
01351 imageThreshold->ThresholdByUpper(maxSignal);
01352 imageThreshold->SetInValue(255);
01353 imageThreshold->SetOutValue(0.0);
01354 imageThreshold->Update();
01355
01356 vtkImageContinuousErode3D *imageErode3D = vtkImageContinuousErode3D::New();
01357 imageErode3D->SetInput(imageThreshold->GetOutput());
01358 imageErode3D->SetKernelSize(4,4,1);
01359 imageErode3D->Update();
01360 erodedata = imageErode3D->GetOutput();
01361
01362
01363 for (int i = 0; i < quantContours[point]->getSize(); i++)
01364 {
01365 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
01366 }
01367
01368
01369
01370 int x;
01371 for (x = limInfX + 30; x < (limSupX - 1); x++)
01372 {
01373
01374 for (int y = limInfY; y < (limSupY - 1) ; y++)
01375 {
01376 bool inWall = false;
01377 bool outCalc = true;
01378
01379 for (int numConts = 0; numConts < polygons.size(); numConts++)
01380 {
01381 if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
01382 {
01383 if (numConts == 0)
01384 {
01385 inWall = true;
01386 }
01387
01388 else
01389 {
01390 outCalc = false;
01391 }
01392
01393 }
01394
01395 }
01396
01397 if (inWall && outCalc)
01398 {
01399 marPoint *p = new marPoint(x,y);
01400 p->setIntensity(imagedata->GetScalarComponentAsDouble(x,y,0,0));
01401 average += p->getIntensity();
01402 points.push_back(p);
01403 }
01404
01405
01406
01407
01408
01409 }
01410 }
01411
01412 if (points.size() > 0)
01413 {
01414 average = average / points.size();
01415 for (int k = 0; k < points.size(); k++)
01416 {
01417 marPoint *p = points[k];
01418 double actualInt = p->getIntensity();
01419 double temp = average - actualInt;
01420 deviation += pow(temp,2.0);
01421 }
01422 deviation = sqrt(deviation / points.size());
01423
01424 }
01425
01426
01427 polygons.push_back(lumenContour[point]);
01428
01429
01430 for (x = limInfX; x < (limSupX - 1); x++)
01431 {
01432 for (int y = limInfY; y < (limSupY - 1); y++)
01433 {
01434 bool inWall = false;
01435 bool outCalc = true;
01436 bool inLumen = false;
01437 unsigned short *refValue = (unsigned short *) image->GetScalarPointer(x,y,0);
01438 *refValue = 0;
01439
01440 for (int numConts = 0; numConts < polygons.size(); numConts++)
01441 {
01442 bool adentro = pointInPolygon(polygons[numConts],x,y);
01443 if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
01444 {
01445 if (numConts == 0)
01446 {
01447 inWall = true;
01448 }
01449 else if (numConts == polygons.size() - 1)
01450 {
01451 inLumen = true;
01452
01453 }
01454 else
01455 {
01456 outCalc = false;
01457 }
01458 }
01459 }
01460
01461 double intensidad = imagedata->GetScalarComponentAsDouble(x,y,0,0);
01462
01463 if (inWall && outCalc && inLumen)
01464 {
01465
01466 if (imagedata->GetScalarComponentAsDouble(x,y,0,0) > (average - 1.5*deviation)
01467 && imagedata->GetScalarComponentAsDouble(x,y,0,0) < (average + 2.5*deviation))
01468 {
01469 *refValue = 255;
01470 }
01471
01472 }
01473
01474 }
01475 }
01476
01477
01478
01479 extractLumen(image,lumenContour[point],point);
01480
01481
01482 }
01483
01484
01485 double marAxisCT::obtainContourArea(marIsocontour* contour)
01486 {
01487
01488 double area = 0.0;
01489 double x1,y1,x2,y2;
01490
01491 for (int i = 0; i < contour->getSize();i++)
01492 {
01493 contour->getPoint(i,&x1,&y1);
01494 if (i < (contour->getSize() - 1))
01495 {
01496 contour->getPoint(i+1,&x2,&y2);
01497 }
01498 else
01499 {
01500 contour->getPoint(0,&x2,&y2);
01501 }
01502
01503 area += (x1 * y2);
01504 area -= (y1 * x2);
01505 }
01506
01507 area = area / 2;
01508
01509 if (area < 0)
01510 area = -area;
01511
01512 return area;
01513 }
01514
01515
01516 void marAxisCT::adjustWall(int point, kVolume* vol)
01517 {
01518
01519 marIsocontour* cont;
01520 double actualArea, previousArea;
01521 std::vector <double > grad;
01522 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
01523
01524
01525 ((marParameters *) getParameters())->setLumenPercentage(50);
01526 updateLumenPercentage(point,vol);
01527 cont = parsePolyDataToMarIsocontour(wallCont);
01528 previousArea = obtainContourArea(cont);
01529
01530 for (int eval = 49; eval >= 0; eval--)
01531 {
01532 ((marParameters *) getParameters())->setLumenPercentage(eval);
01533 updateLumenPercentage(point,vol);
01534 cont = parsePolyDataToMarIsocontour(wallCont);
01535 actualArea = obtainContourArea(cont);
01536 grad.push_back(fabs(actualArea - previousArea));
01537 }
01538
01539 int newPercentage = 50 - maxValue(grad);
01540 ((marParameters *) getParameters())->setLumenPercentage(newPercentage);
01541 updateLumenPercentage(point,vol);
01542
01543
01544
01545 }
01546
01547
01548 void marAxisCT::adjustCalcification(int point, kVolume* vol)
01549 {
01550 }
01551
01552
01553 int marAxisCT::maxValue(std::vector <double> list)
01554 {
01555 int max = 0;
01556 double maxVal = list[0];
01557
01558 for (int i = 0; i < list.size(); i++)
01559 {
01560 if (list[i] > maxVal)
01561 {
01562 maxVal = list[i];
01563 max = i;
01564 }
01565 }
01566
01567 return max;
01568 }
01569
01570
01571
01572 int marAxisCT::getDiameter(marContourVO* contourVO, double x, double y, int limInfX, int limInfY, int limSupX, int limSupY)
01573 {
01574 bool stop = false;
01575 double originX = x;
01576 double originY = y;
01577
01578 int coordBase = 1, diameterOne = 0, diameterTwo = 0;
01579 marIsocontour* cont = parsePolyDataToMarIsocontour(contourVO);
01580 generateVector(1,&coordBase,originX,originY,&x,&y);
01581
01582 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
01583 {
01584 generateVector(0,&coordBase,originX,originY,&x,&y);
01585 if (pointInPolygon(cont,x,y))
01586 {
01587 diameterOne++;
01588 }
01589 else
01590 {
01591 stop = true;
01592 }
01593 }
01594
01595 stop = false;
01596 x = originX;
01597 y = originY;
01598 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
01599 {
01600 generateVector(35,&coordBase,originX,originY,&x,&y);
01601 if (pointInPolygon(cont,x,y))
01602 {
01603 diameterOne++;
01604 }
01605 else
01606 {
01607 stop = true;
01608 }
01609 }
01610
01611 stop = false;
01612 x = originX;
01613 y = originY;
01614 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
01615 {
01616 generateVector(17,&coordBase,originX,originY,&x,&y);
01617 if (pointInPolygon(cont,x,y))
01618 {
01619 diameterTwo++;
01620 }
01621 else
01622 {
01623 stop = true;
01624 }
01625 }
01626
01627 stop = false;
01628 x = originX;
01629 y = originY;
01630 while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
01631 {
01632 generateVector(107,&coordBase,originX,originY,&x,&y);
01633 if (pointInPolygon(cont,x,y))
01634 {
01635 diameterTwo++;
01636 }
01637 else
01638 {
01639 stop = true;
01640 }
01641 }
01642
01643 delete cont;
01644
01645 return (diameterOne + diameterTwo) / 2;
01646 }
01647
01648
01649 void marAxisCT::histogram(int point, kVolume* vol)
01650 {
01651
01652 vtkImageData* imagedata;
01653 int limInfX, limInfY, limInfZ;
01654 int limSupX, limSupY, limSupZ;
01655 int pos;
01656
01657
01658 FILE *archivo;
01659 char nomArchivo[30],temp[3];
01660 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
01661
01662
01663
01664
01665 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
01666 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
01667 ((marParameters *) getParameters())->setLumenPercentage(86);
01668 updateLumenPercentage(point,vol);
01669 marIsocontour* contExt = parsePolyDataToMarIsocontour(wallCont);
01670
01671 std::vector<double> intensity;
01672 std::vector<int> frecuencies;
01673
01674 for (int i = 0; i < limSupX; i++)
01675 {
01676 for (int j = 0; j < limSupY; j++)
01677 {
01678 if(pointInPolygon(contExt,i,j))
01679 {
01680
01681 pos = searchData(intensity,imagedata->GetScalarComponentAsDouble(i,j,0,0));
01682 if (pos == -1)
01683 {
01684 intensity.push_back( imagedata->GetScalarComponentAsDouble(i,j,0,0));
01685 frecuencies.push_back(1);
01686
01687
01688 }
01689 else
01690 {
01691 frecuencies[pos]++;
01692 }
01693 }
01694
01695 }
01696 }
01697
01698 strcpy(nomArchivo,"./histogram");
01699
01700
01701
01702 sprintf(temp,"%d",point);
01703
01704 strcat(nomArchivo,temp);
01705 strcat(nomArchivo,".csv");
01706
01707
01708
01709 archivo = fopen(nomArchivo,"w");
01710 for (int k = 0; k < frecuencies.size(); k++)
01711 {
01712
01713 fprintf(archivo,"%f;%d\n",intensity[k],frecuencies[k]);
01714
01715 }
01716 fclose(archivo);
01717
01718
01719
01720
01721
01722 delete contExt;
01723
01724
01725 }
01726
01727
01728 int marAxisCT::searchData(std::vector<double> vec, double value)
01729 {
01730 int resp = -1;
01731 bool found = false;
01732
01733 for (int i = 0; i< vec.size() && !found; i++)
01734 {
01735 if (vec[i] == value)
01736 {
01737 resp = i;
01738 found = true;
01739 }
01740 }
01741
01742 return resp;
01743 }
01744
01745
01746 void marAxisCT::setCalibration(bool calib)
01747 {
01748 calibration = calib;
01749 }
01750
01751
01752 bool marAxisCT::getCalibration()
01753 {
01754 return calibration;
01755 }
01756
01757
01758 void marAxisCT::adjustContour(int point, kVolume* vol)
01759 {
01760
01761 double percentage = 0.05;
01762 double prevDifference = 0;
01763 if (quantContours[point + 1] == NULL)
01764 {
01765 create2Dcontours(point + 1 ,vol);
01766 }
01767
01768
01769 double signal = maxSignal;
01770 marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
01771
01772
01773 marIsocontour* prev = parsePolyDataToMarIsocontour(searchContour(marAxisContours::WALL, point + 1));
01774
01775
01776 double newValue = signal;
01777 wallCont->getIsocontour()->SetValue(0,newValue);
01778 wallCont->getIsocontour()->Update();
01779 wallCont->get2DContour()->Update();
01780
01781
01782 marIsocontour* actual = parsePolyDataToMarIsocontour(wallCont);
01783
01784
01785 double prevA = obtainContourArea(prev);
01786 double actuA = obtainContourArea(actual);
01787
01788 prevDifference = fabs(prevA - actuA);
01789
01790
01791 if (prevA - actuA > 0 && prevA - actuA < percentage*prevA)
01792 {
01793 bool stop = false;
01794 while (!stop && newValue > 0)
01795 {
01796 newValue--;
01797
01798 wallCont->getIsocontour()->SetValue(0,newValue);
01799 wallCont->getIsocontour()->Update();
01800 wallCont->get2DContour()->Update();
01801
01802 actual = parsePolyDataToMarIsocontour(wallCont);
01803 actuA = obtainContourArea(actual);
01804
01805 if (fabs(prevA - actuA) > percentage*prevA)
01806 {
01807 newValue++;
01808 wallCont->getIsocontour()->SetValue(0,newValue);
01809 wallCont->getIsocontour()->Update();
01810 wallCont->get2DContour()->Update();
01811 stop = true;
01812 }
01813 }
01814
01815 }
01816 else if (prevA - actuA < 0 && prevA - actuA > -percentage*prevA)
01817 {
01818
01819 bool stop = false;
01820 while (!stop && newValue > 0)
01821 {
01822 newValue--;
01823
01824 wallCont->getIsocontour()->SetValue(0,newValue);
01825 wallCont->getIsocontour()->Update();
01826 wallCont->get2DContour()->Update();
01827
01828 actual = parsePolyDataToMarIsocontour(wallCont);
01829 actuA = obtainContourArea(actual);
01830
01831 if (prevA - actuA < -percentage*prevA)
01832 {
01833 newValue++;
01834 wallCont->getIsocontour()->SetValue(0,newValue);
01835 wallCont->getIsocontour()->Update();
01836 wallCont->get2DContour()->Update();
01837 stop = true;
01838 }
01839 }
01840 }
01841 else if (prevA - actuA < 0 && prevA - actuA < -percentage*prevA)
01842 {
01843
01844 bool stop = false;
01845
01846 while (!stop)
01847 {
01848 newValue++;
01849
01850 wallCont->getIsocontour()->SetValue(0,newValue);
01851 wallCont->getIsocontour()->Update();
01852 wallCont->get2DContour()->Update();
01853
01854 actual = parsePolyDataToMarIsocontour(wallCont);
01855
01856 actuA = obtainContourArea(actual);
01857 if (prevA - actuA > -percentage*prevA)
01858 {
01859 if (prevDifference < fabs(prevA - actuA))
01860 {
01861 newValue--;
01862 wallCont->getIsocontour()->SetValue(0,newValue);
01863 wallCont->getIsocontour()->Update();
01864 wallCont->get2DContour()->Update();
01865 }
01866
01867 stop = true;
01868
01869 }
01870 prevDifference = fabs(prevA - actuA);
01871 }
01872 }
01873 else if (prevA - actuA > 0 && prevA - actuA > percentage*prevA)
01874 {
01875 bool stop = false;
01876 while (!stop && newValue > 0)
01877 {
01878 newValue--;
01879
01880 wallCont->getIsocontour()->SetValue(0,newValue);
01881 wallCont->getIsocontour()->Update();
01882 wallCont->get2DContour()->Update();
01883
01884 actual = parsePolyDataToMarIsocontour(wallCont);
01885 actuA = obtainContourArea(actual);
01886
01887 if (prevA - actuA < percentage*prevA)
01888 {
01889 if (prevDifference < fabs(prevA - actuA))
01890 {
01891 newValue++;
01892 wallCont->getIsocontour()->SetValue(0,newValue);
01893 wallCont->getIsocontour()->Update();
01894 wallCont->get2DContour()->Update();
01895 }
01896
01897 stop = true;
01898
01899 }
01900 prevDifference = fabs(prevA - actuA);
01901 }
01902 }
01903
01904
01905 maxSignal = newValue;
01906
01907 this->markUpLumen(point,vol);
01908 }
01909
01910
01911 int marAxisCT::getStartIndex()
01912 {
01913 return startIndex;
01914 }
01915
01916
01917 void marAxisCT::setStartIndex(int start)
01918 {
01919 startIndex = start;
01920
01921 int per = ((marParameters *) getParameters())->getLumenPercentage();
01922 double oldValue = quantContours[start]->getContour(0)->getSignal()*per*pow(10.0,-2.0);
01923
01924 maxSignal = oldValue;
01925 }
01926
01927
01928 int marAxisCT::getPointSize()
01929 {
01930 return vesselPoints.size();
01931 }
01932
01933
01934 marPoint* marAxisCT::getPoint(int i)
01935 {
01936 return vesselPoints[i];
01937 }
01938
01939
01940 void marAxisCT::extractLumen(vtkImageData *lumenImage, marIsocontour *lumenContour, int point)
01941 {
01942
01943
01944 vtkImageCast* cast = vtkImageCast::New();
01945 cast->SetInput(lumenImage);
01946 cast->SetOutputScalarTypeToDouble();
01947
01948 vtkImageGaussianSmooth* filter = vtkImageGaussianSmooth::New();
01949 filter->SetDimensionality(2);
01950 filter->SetInput(cast->GetOutput());
01951 filter->SetStandardDeviations((getParameters())->getStandardDeviation(),(getParameters())->getStandardDeviation());
01952 filter->SetRadiusFactors((getParameters())->getRadius(),(getParameters())->getRadius());
01953 filter->SetStandardDeviations(3,3);
01954 filter->SetRadiusFactors(3,3);
01955
01956
01957
01958 marAxisContours* axisCont = quantContours[point];
01959
01960 int pos = -1;
01961
01962 for (int i = 0; i < axisCont->getSize() && pos == -1; i++)
01963 {
01964 if (axisCont->getContourType(i) == marAxisContours::ELUMEN)
01965 {
01966 pos = i;
01967 }
01968 }
01969
01970 bool stop = false;
01971 if (pos != -1)
01972 {
01973 stop = true;
01974 if (axisCont->isReplaced(pos))
01975 {
01976 stop = true;
01977 }
01978
01979 }
01980
01981 if (!stop)
01982 {
01983 marContourVO* contourVO = new marContourVO();
01984 marIsocontour* iso = lumenContour;
01985
01986 contourVO->setIsocontour(vtkContourFilter::New());
01987 contourVO->getIsocontour()->SetInput( filter->GetOutput() );
01988 contourVO->getIsocontour()->SetNumberOfContours( 1 );
01989 contourVO->getIsocontour()->SetValue( 0, 128 );
01990 contourVO->getIsocontour()->Update( );
01991
01992 contourVO->setIsocontourCpd(vtkCleanPolyData::New());
01993 contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
01994 contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
01995 contourVO->getIsocontourCpd()->Update( );
01996
01997 double cgx, cgy;
01998
01999 iso->getCG(&cgx,&cgy);
02000
02001 contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
02002 contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
02003
02004 contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
02005 contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
02006 contourVO->getIsocontourDcf()->Update( );
02007
02008 contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
02009 contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
02010 contourVO->getIsocontourCpd2()->Update();
02011
02012 contourVO->setIsocontourStripped(vtkStripper::New( ));
02013 contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
02014 contourVO->getIsocontourStripped()->Update();
02015
02016 contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
02017 contourVO->get2DContour()->Update();
02018 contourVO->setSignal(255);
02019 contourVO->setType(marAxisContours::ELUMEN);
02020
02021
02022
02023 marIsocontour* actual = parsePolyDataToMarIsocontour(contourVO);
02024 marContourVO* aVO = searchContour(marAxisContours::WALL,point);
02025 double wallArea = 0;
02026 if (aVO != NULL)
02027 {
02028 wallArea = obtainContourArea(parsePolyDataToMarIsocontour(aVO));
02029 }
02030
02031 double prevA = obtainContourArea(actual);
02032
02033 int newValue = 128;
02034 int finalValue = 128;
02035 double actuA = 1;
02036 while (newValue > 0 && actuA > 0)
02037 {
02038
02039 newValue--;
02040
02041 try
02042 {
02043 contourVO->getIsocontour()->SetValue(0,newValue);
02044 contourVO->getIsocontour()->Update();
02045 contourVO->get2DContour()->Update();
02046 actual = parsePolyDataToMarIsocontour(contourVO);
02047 actuA = obtainContourArea(actual);
02048 if ((actuA /1000) < 1 && (prevA / 1000) > 1)
02049 {
02050 finalValue = newValue;
02051 prevA = actuA;
02052 }
02053 else if (actuA > prevA && (prevA / 1000) < 1 && actuA < 0.1*wallArea)
02054 {
02055 finalValue = newValue;
02056 prevA = actuA;
02057 }
02058 }
02059 catch (std::exception e)
02060 {
02061 actuA = 0;
02062 }
02063
02064 }
02065
02066 contourVO->getIsocontour()->SetValue(0,finalValue);
02067 contourVO->getIsocontour()->Update();
02068 contourVO->get2DContour()->Update();
02069
02070
02071
02072
02073
02074 if (pos != -1)
02075 {
02076 axisCont->replaceContour(contourVO,pos);
02077 }
02078 else
02079 {
02080 axisCont->addContour(contourVO);
02081 }
02082
02083
02084
02085 quantContours[point] = axisCont;
02086 }
02087
02088
02089 }
02090
02091
02092 void marAxisCT::generateFile(int point,kVolume* vol)
02093 {
02094 FILE *archivo;
02095 char nomArchivo[30], temp[3];
02096
02097 strcpy(nomArchivo,"./point");
02098
02099
02100
02101 sprintf(temp,"%d",point);
02102
02103 strcat(nomArchivo,temp);
02104 strcat(nomArchivo,".txt");
02105 archivo = fopen(nomArchivo,"w");
02106
02107 fprintf(archivo, "%d", quantContours[point]->getSize());
02108
02109 int i;
02110 for (i = 0; i < quantContours[point]->getSize(); i++)
02111 {
02112 marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
02113 fprintf(archivo, " %d",iso->getSize());
02114
02115 }
02116 fprintf(archivo,"\n");
02117
02118 for ( i = 0; i < quantContours[point]->getSize(); i++)
02119 {
02120 marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
02121 fprintf(archivo,"Tipo: %d\n",quantContours[point]->getContourType(i));
02122 for (int j = 0; j < iso->getSize(); j++)
02123 {
02124 double tempX, tempY;
02125 iso->getPoint(j,&tempX,&tempY);
02126 fprintf(archivo,"%f;%f;%d\n",tempX,tempY,point);
02127 }
02128 }
02129
02130
02131 vtkImageData *vtkimagedata = this->getSliceImage(point,vol);
02132
02133 std::string directory = "./";
02134 std::string filename = "xx";
02135 float rescalaSlope = 1;
02136 float rescalaIntercept = 0;
02137 int dim[3];
02138 vtkimagedata->GetDimensions(dim);
02139 int voi[6];
02140 voi[0]=0;
02141 voi[1]=dim[0];
02142 voi[2]=0;
02143 voi[3]=dim[1];
02144 voi[4]=0;
02145 voi[5]=dim[2];
02146 marRAW2Files marraw2;
02147 marraw2.saveVolume(directory,filename,vtkimagedata,voi,rescalaSlope,rescalaIntercept);
02148
02149
02150
02151 fclose(archivo);
02152
02153 }
02154
02155
02156 void marAxisCT::replaceContour2D(int point,int size,double *vx,double *vy, int type)
02157 {
02158
02159 vtkPoints *_pts = vtkPoints::New();
02160 _pts->SetNumberOfPoints(size);
02161 int j;
02162
02163 for (j=0 ; j<size ; j++){
02164 _pts->SetPoint(j, vx[j] , vy[j] , 0 );
02165 }
02166
02167
02168 vtkCellArray *lines = vtkCellArray::New();
02169 lines->InsertNextCell( size );
02170 for ( j=0 ; j<size+1 ; j++ ){
02171 lines->InsertCellPoint(j % size );
02172 }
02173
02174 vtkPolyData *_pd = vtkPolyData::New();
02175 _pd->SetPoints( _pts );
02176 _pd->SetLines( lines );
02177 lines->Delete();
02178 _pts->Delete();
02179
02180 marContourVO* vo = new marContourVO();
02181 vo->setReplaced(true);
02182 vo->set2DContour(_pd);
02183 vo->setType(type);
02184 quantContours[point]->addContour(vo);
02185
02186
02187 }
02188
02189
02190 void marAxisCT::cleanContours(int type, int point)
02191 {
02192 marAxisContours* newContours = new marAxisContours();
02193 if (quantContours[point] != NULL)
02194 {
02195 for (int i = 0; i < quantContours[point]->getSize(); i++)
02196 {
02197 if (quantContours[point]->getContourType(i) != type)
02198 {
02199 newContours->addContour(quantContours[point]->getContour(i));
02200 }
02201 else
02202 {
02203
02204
02205
02206 }
02207 }
02208 }
02209
02210 quantContours[point] = newContours;
02211 }
02212
02213
02214 marContourVO* marAxisCT::searchContour(int type, int point)
02215 {
02216 marContourVO* cont = NULL;
02217
02218 for (int i = 0; i < quantContours[point]->getSize(); i++)
02219 {
02220 if (quantContours[point]->getContourType(i) == type)
02221 {
02222 cont = quantContours[point]->getContour(i);
02223 }
02224 }
02225
02226 return cont;
02227 }
02228
02229
02230 double marAxisCT::performXOR(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
02231 {
02232 if (manual.size() == 0)
02233 {
02234 return 0;
02235 }
02236
02237 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
02238 std::vector <marIsocontour* > polygons;
02239 vtkImageData* imagedata;
02240
02241 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
02242 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02243 vtkImageData* imageA = vtkImageData::New();
02244 imageA->SetDimensions(limSupX,limSupY,0);
02245 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02246 imageA->SetScalarTypeToUnsignedShort();
02247 imageA->AllocateScalars();
02248 imageA->Update();
02249
02250 vtkImageData* imageM = vtkImageData::New();
02251 imageM->SetDimensions(limSupX,limSupY,0);
02252 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02253 imageM->SetScalarTypeToUnsignedShort();
02254 imageM->AllocateScalars();
02255 imageM->Update();
02256
02257
02258 for (int i = 0; i < quantContours[point]->getSize(); i++)
02259 {
02260 if (quantContours[point]->getContourType(i) == type)
02261 {
02262 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
02263 }
02264
02265 }
02266
02267 if (polygons.size() == 0)
02268 {
02269 return 0;
02270 }
02271 int x,y;
02272 for (x = limInfX; x < (limSupX - 1); x++)
02273 {
02274 for (y = limInfY; y < (limSupY - 1); y++)
02275 {
02276 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
02277 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
02278 *refValueA = 0;
02279 *refValueM = 0;
02280
02281 int numConts;
02282 for (numConts = 0; numConts < polygons.size(); numConts++)
02283 {
02284 if (pointInPolygon(polygons[numConts],x,y))
02285 {
02286 *refValueA = 255;
02287 }
02288 }
02289
02290 for (numConts = 0; numConts < manual.size(); numConts++)
02291 {
02292 if (pointInPolygon(manual[numConts],x,y))
02293 {
02294 *refValueM = 255;
02295 }
02296 }
02297 }
02298 }
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310 int coincidencias = 0;
02311 for (x = limInfX; x < (limSupX - 1); x++)
02312 {
02313 for (y = limInfY; y < (limSupY - 1); y++)
02314 {
02315 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
02316 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
02317
02318 int xor_ = compA^compM;
02319 if (xor_ == 255)
02320 {
02321 coincidencias++;
02322 }
02323 }
02324 }
02325
02326 double resp = (double)coincidencias ;
02327
02328 return resp;
02329
02330 }
02331
02332
02333 double marAxisCT::performAND(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
02334 {
02335
02336 if (manual.size() == 0)
02337 {
02338 return 0;
02339 }
02340
02341 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
02342 std::vector <marIsocontour* > polygons;
02343 vtkImageData* imagedata;
02344
02345 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
02346 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02347 vtkImageData* imageA = vtkImageData::New();
02348 imageA->SetDimensions(limSupX,limSupY,0);
02349 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02350 imageA->SetScalarTypeToUnsignedShort();
02351 imageA->AllocateScalars();
02352 imageA->Update();
02353
02354 vtkImageData* imageM = vtkImageData::New();
02355 imageM->SetDimensions(limSupX,limSupY,0);
02356 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02357 imageM->SetScalarTypeToUnsignedShort();
02358 imageM->AllocateScalars();
02359 imageM->Update();
02360
02361 for (int i = 0; i < quantContours[point]->getSize(); i++)
02362 {
02363 if (quantContours[point]->getContourType(i) == type)
02364 {
02365 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
02366 }
02367
02368 }
02369
02370 if (polygons.size() == 0)
02371 {
02372 return 0;
02373 }
02374
02375 int x,y;
02376 for (x = limInfX; x < (limSupX - 1); x++)
02377 {
02378 for (y = limInfY; y < (limSupY - 1); y++)
02379 {
02380 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
02381 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
02382 *refValueA = 0;
02383 *refValueM = 0;
02384
02385 int numConts;
02386 for (numConts = 0; numConts < polygons.size(); numConts++)
02387 {
02388 if (pointInPolygon(polygons[numConts],x,y))
02389 {
02390 *refValueA = 255;
02391 }
02392 }
02393
02394 for (numConts = 0; numConts < manual.size(); numConts++)
02395 {
02396 if (pointInPolygon(manual[numConts],x,y))
02397 {
02398 *refValueM = 255;
02399 }
02400 }
02401 }
02402 }
02403
02404
02405
02406
02407
02408
02409
02410 int coincidencias = 0;
02411 for (x = limInfX; x < (limSupX - 1); x++)
02412 {
02413 for (y = limInfY; y < (limSupY - 1); y++)
02414 {
02415 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
02416 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
02417
02418 int and_ = compA & compM;
02419 if (and_ == 255)
02420 {
02421 coincidencias++;
02422 }
02423 }
02424 }
02425
02426 double resp = (double)coincidencias ;
02427
02428
02429 return resp;
02430 }
02431
02432
02433 marIsocontour* marAxisCT::loadMarIsocontour(int size, double *vx, double *vy)
02434 {
02435 marIsocontour* cont = new marIsocontour();
02436
02437 for ( int j=0 ; j<size ; j++ )
02438 {
02439 cont->insertPoint(vx[j], vy[j]);
02440 }
02441
02442 return cont;
02443 }
02444
02445
02446 double marAxisCT::performUnion(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
02447 {
02448 if (manual.size() == 0)
02449 {
02450 return 0;
02451 }
02452
02453 int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
02454 std::vector <marIsocontour* > polygons;
02455 vtkImageData* imagedata;
02456
02457 imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
02458 imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02459 vtkImageData* imageA = vtkImageData::New();
02460 imageA->SetDimensions(limSupX,limSupY,0);
02461 imageA->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02462 imageA->SetScalarTypeToUnsignedShort();
02463 imageA->AllocateScalars();
02464 imageA->Update();
02465
02466 vtkImageData* imageM = vtkImageData::New();
02467 imageM->SetDimensions(limSupX,limSupY,0);
02468 imageM->SetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
02469 imageM->SetScalarTypeToUnsignedShort();
02470 imageM->AllocateScalars();
02471 imageM->Update();
02472
02473 int i;
02474 for (i = 0; i < quantContours[point]->getSize(); i++)
02475 {
02476 if (quantContours[point]->getContourType(i) == type)
02477 {
02478 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
02479 }
02480
02481 }
02482
02483 if (polygons.size() == 0)
02484 {
02485 return 0;
02486 }
02487
02488 int x,y;
02489 for ( x = limInfX; x < (limSupX - 1); x++)
02490 {
02491 for ( y = limInfY; y < (limSupY - 1); y++)
02492 {
02493 unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
02494 unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
02495 *refValueA = 0;
02496 *refValueM = 0;
02497 int numConts;
02498 for ( numConts = 0; numConts < polygons.size(); numConts++)
02499 {
02500 if (pointInPolygon(polygons[numConts],x,y))
02501 {
02502 *refValueA = 255;
02503 }
02504 }
02505
02506 for (numConts = 0; numConts < manual.size(); numConts++)
02507 {
02508 if (pointInPolygon(manual[numConts],x,y))
02509 {
02510 *refValueM = 255;
02511 }
02512 }
02513 }
02514 }
02515
02516
02517
02518
02519
02520
02521
02522 int coincidencias = 0;
02523 for (x = limInfX; x < (limSupX - 1); x++)
02524 {
02525 for ( y = limInfY; y < (limSupY - 1); y++)
02526 {
02527 int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
02528 int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
02529
02530 int or_ = compA | compM;
02531 if (or_ == 255)
02532 {
02533 coincidencias++;
02534 }
02535 }
02536 }
02537
02538 double resp = (double)coincidencias;
02539
02540
02541 return resp;
02542 }
02543
