lib/maracasVisuLib/src/kernel/include/marAxisCT.cpp

   1 // marAxisCT.cpp: implementation of the marAxisCT class.
   2 //
   3 //////////////////////////////////////////////////////////////////////
   4
   5
   6 #include <vtkImageCast.h>
   7 #include <vtkImageGaussianSmooth.h>
   8 #include <vtkCleanPolyData.h>
   9 #include <vtkPolyDataConnectivityFilter.h>
  10 #include <vtkStripper.h>
  11 #include <vtkCell.h>
  12 #include <vtkIdList.h>
  13 #include <vtkMatrix4x4.h>
  14 #include <vtkTransform.h>
  15 #include <vtkPlaneSource.h>
  16 #include <vtkContourFilter.h>
  17 #include <marLine.h>
  18 #include <marUtils.h>
  19 #include <vtkImageContinuousErode3D.h>
  20 #include <vtkImageThreshold.h>
  21 #include <vtkImageSobel2D.h>
  22 #include <vtkCellArray.h>
  23 #include <vtkImageLogic.h>
  24
  25
  26 #include "marAxisCT.h"
  27 #include "marGdcmDicom.h"
  28
  29
  30
  31
  32 marAxisCT::marAxisCT(  ) : marAxis(  )
  33 {
  34 }
  35
  36 // ----------------------------------------------------------------------------
  37 marContour* marAxisCT::getContour( int point , kVolume* vol, int index ) {
  38         if (quantContours[point] == NULL || quantContours[point]->getContour(index)->getContour() == NULL)
  39         {
  40                 createContours(point,vol);
  41         }
  42
  43         return quantContours[point]->getContour(index)->getContour();
  44 }
  45
  46 // ----------------------------------------------------------------------------
  47 vtkPoints* marAxisCT::get3Dcontour( int point , kVolume* vol, int index ) {
  48         if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get3DContour() == NULL)
  49         {
  50                 create3Dcontours(point,vol);
  51         }
  52
  53         marAxisContours* cont = quantContours[point];
  54
  55         return cont->getContour(index)->get3DContour();
  56 }
  57
  58 // ----------------------------------------------------------------------------
  59 vtkPolyData* marAxisCT::get2Dcontour( int point , kVolume* vol, int index ) {
  60         if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DContour() == NULL)
  61         {
  62                 create2Dcontours(point,vol);
  63         }
  64
  65         marAxisContours* cont = quantContours[point];
  66
  67
  68
  69         return cont->getContour(index)->get2DContour();
  70 }
  71
  72 // ----------------------------------------------------------------------------
  73 vtkPoints* marAxisCT::get2DDiameterMin( int point , kVolume* vol, int index ) {
  74 //      if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DDiameterMin() == NULL)
  75 //      {
  76                 create2DDiametersMin(point,vol);
  77 //      }
  78
  79         marAxisContours* cont = quantContours[point];
  80
  81         return cont->getContour(index)->get2DDiameterMin();
  82 }
  83
  84 // ----------------------------------------------------------------------------
  85 vtkPoints* marAxisCT::get2DDiameterMax( int point , kVolume* vol, int index ) {
  86 //      if (quantContours[point] == NULL || quantContours[point]->getContour(index)->get2DDiameterMax() == NULL)
  87 //      {
  88                 create2DDiametersMax(point,vol);
  89 //      }
  90
  91         marAxisContours* cont = quantContours[point];
  92
  93         return cont->getContour(index)->get2DDiameterMax();
  94 }
  95
  96 // ----------------------------------------------------------------------------
  97 int marAxisCT::getSignal(int point, int index, kVolume* vol) {
  98         if (quantContours[point] == NULL )
  99         {
 100                 createSignals(point,vol);
 101         }
 102
 103         return (int) ( quantContours[point]->getContour(index)->getSignal() );
 104 }
 105
 106 // ----------------------------------------------------------------------------
 107 void marAxisCT::createContours( int point, kVolume* vol ) {
 108
 109         int j,id;
 110         int numberOfPoints,numberOfCells;
 111         vtkCell* cell;
 112         vtkIdList* pids;
 113         double p[ 3 ];
 114         vtkPolyData* vtkPolyData_2DContour;
 115         double x,y;
 116
 117         int sizeIma                             =       getParameters( )->getSizeIma( );
 118         double dimIma                   =       getParameters( )->getDimIma( );
 119
 120         if (quantContours.size() < point - 1 )
 121         {
 122                 create2Dcontours(point,vol);
 123         }
 124
 125         marAxisContours* axisCont = quantContours[point];
 126
 127         for (int i = 0; i < axisCont->getSize(); i++)
 128         {
 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++ )
 137                 {
 138                         id = pids->GetId( j );
 139                         vtkPolyData_2DContour->GetPoint( id, p);
 140                         x=p[0]-64.0;
 141                         y=p[1]-64.0;
 142                         x=x * dimIma / ( float ) sizeIma;
 143                         y=y * dimIma / ( float ) sizeIma;
 144                         aContourVO->getContour()->addContourPoint( x , y );
 145                 }
 146
 147                 aContourVO->getContour()->do_spline();
 148                 aContourVO->getContour()->calculateVariables();
 149         }
 150
 151 }
 152
 153 // ----------------------------------------------------------------------------
 154 void marAxisCT::create3Dcontours(       int point       , kVolume* vol ) {
 155
 156         if (quantContours[point] == NULL )
 157         {
 158                 createContours(point,vol);
 159         }
 160
 161         marAxisContours* axisCont = quantContours[point];
 162
 163         for (int i = 0; i < axisCont->getSize(); i++)
 164         {
 165                 vtkPoints *_vtkPoints;
 166                 double *o;
 167                 double *n;
 168                 vtkMatrix4x4* mat;
 169                 vtkTransform* trans;
 170                 double nn[3];
 171                 double c[3],p1[3],p2[3];
 172                 double nA,nB,nC;
 173
 174                 _vtkPoints                                              = vtkPoints::New();
 175                 o                                                               = getPoints(point); //PENDIENTE  _points[i];
 176                 n                                                               = getNormal( i );
 177                 mat                                                             = vtkMatrix4x4::New();
 178                 trans                                                   = vtkTransform::New();
 179
 180
 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;
 184
 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 ] );
 191                 pSource->Update( );
 192                 pSource->Update( );
 193                 pSource->GetOrigin(c);
 194                 pSource->GetPoint1(p1);
 195                 pSource->GetPoint2(p2);
 196                 pSource->Delete();
 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;
 203
 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);
 220
 221                 double deter=mat->Determinant(mat);
 222                 trans->Identity();
 223                 trans->SetMatrix(mat);
 224                 float ang;
 225                 ang=-90;
 226                 trans->RotateWXYZ  ( ang,0,1,0);
 227                 trans->Update();
 228
 229                 vtkMatrix4x4 *m=vtkMatrix4x4::New();
 230                 trans->GetMatrix  (  m  );
 231
 232                 int j,numberOfPoints;
 233
 234                 marContour* contour2D   = getContour(point,vol,i);
 235                 marContourVO* aContourVO = axisCont->getContour(i);
 236                 numberOfPoints                  = contour2D->getNumberOfPoints( );
 237                 double fpA[3];
 238                 double pA[3],ppp[3];
 239
 240                 for( j = 0; j <= numberOfPoints; j++ ) {
 241                         contour2D->getPoint( fpA,j % numberOfPoints);
 242                         pA[0]=fpA[0];
 243                         pA[1]=fpA[1];
 244                         pA[2]=0;
 245
 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);
 249
 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] );
 252
 253
 254                 }
 255
 256
 257                 m->Delete();
 258                 mat->Delete();
 259                 trans->Delete();
 260
 261                 if (aContourVO->get3DContour()!=NULL){ aContourVO->get3DContour()->Delete(); }
 262                         aContourVO->set3DContour(_vtkPoints);
 263
 264         }
 265
 266
 267
 268
 269 }
 270
 271 // ----------------------------------------------------------------------------
 272 void marAxisCT::create2Dcontours(       int point       , kVolume* vol ) {
 273         std::vector <marIsocontour *> contours;
 274
 275
 276
 277         generatePoints(point,vol,&contours);
 278
 279         vtkImageData* imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
 280         marAxisContours* axisCont = new marAxisContours();
 281
 282         for (int i = 0; i < contours.size(); i++)
 283         {
 284                 marContourVO* contourVO = new marContourVO();
 285                 marIsocontour* iso = contours[i];
 286
 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( );
 292
 293                 contourVO->setIsocontourCpd(vtkCleanPolyData::New());
 294                 contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
 295                 contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
 296                 contourVO->getIsocontourCpd()->Update( );
 297
 298                 double cgx, cgy;
 299
 300                 iso->getCG(&cgx,&cgy);
 301
 302                 contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
 303                 contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
 304
 305                 contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
 306                 contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
 307                 contourVO->getIsocontourDcf()->Update( );
 308
 309                 contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
 310                 contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
 311                 contourVO->getIsocontourCpd2()->Update();
 312
 313                 contourVO->setIsocontourStripped(vtkStripper::New( ));
 314                 contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
 315                 contourVO->getIsocontourStripped()->Update();
 316
 317                 contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
 318                 contourVO->get2DContour()->Update();
 319                 contourVO->setSignal(iso->getMaxIntensity());
 320
 321                 if (iso->getType() == marAxisContours::LUMEN)
 322                 {
 323                         contourVO->setType(marAxisContours::WALL);
 324                 }
 325                 else
 326                 {
 327                         contourVO->setType(iso->getType());
 328                 }
 329
 330
 331
 332                 axisCont->addContour(contourVO);
 333
 334
 335                 delete iso;
 336                 contours[i] = NULL;
 337         }
 338
 339
 340
 341         quantContours[point] = axisCont;
 342
 343
 344
 345         updateLumenPercentage(point,vol);
 346
 347         updateCalcPercentage(point,vol);
 348 //      adjustWall(point,vol);
 349         if (calibration && point < startIndex)
 350         {
 351                 this->adjustContour(point,vol);
 352         }
 353         else
 354         {
 355                 adjustWall(point,vol);
 356         }
 357         contours.clear();
 358         imagedata->Delete();
 359
 360
 361 //      marIsocontour* cont = parsePolyDataToMarIsocontour(quantContours[point]->getContour(0));
 362
 363         /*FILE *archivo;
 364         char nomArchivo[30], temp[3];
 365
 366         strcpy(nomArchivo,"./points");
 367         itoa(point,temp,10);
 368         strcat(nomArchivo,temp);
 369         strcat(nomArchivo,".csv");
 370
 371         archivo = fopen(nomArchivo,"w");
 372
 373         for (int h = 0; h < cont->getSize(); h++)
 374         {
 375                 double x, y;
 376                 cont->getPoint(h,&x,&y);
 377                 fprintf(archivo,"%f;%f\n",x,y);
 378         }
 379
 380         fclose(archivo);*/
 381 }
 382
 383 // ----------------------------------------------------------------------------
 384 void marAxisCT::create2DDiametersMin(int point  , kVolume* vol ) {
 385
 386 if (quantContours[point] == NULL )
 387         {
 388                 createContours(point,vol);
 389         }
 390
 391         marAxisContours* axisCont = quantContours[point];
 392
 393         for (int i = 0; i < axisCont->getSize(); i++)
 394         {
 395                 double p1[2],p2[2];
 396                 marContour *marcontour=getContour(point,vol,i);
 397                 marContourVO* aContourVO = axisCont->getContour(i);
 398
 399                 marcontour->getMinimumLine(p1,p2);
 400                 vtkPoints *_vtkPoints = vtkPoints::New();
 401
 402                 int sizeIma = getParameters( )->getSizeIma( );
 403                 double dimIma = getParameters( )->getDimIma( );
 404                 double factor=( float ) sizeIma / dimIma;
 405
 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 );
 412
 413                 if ( aContourVO->get2DDiameterMin()!=NULL ) { aContourVO->get2DDiameterMin()->Delete(); }
 414                         aContourVO->set2DDiameterMin(_vtkPoints);
 415
 416         }
 417
 418 }
 419
 420 // ----------------------------------------------------------------------------
 421 void marAxisCT::create2DDiametersMax(int point  , kVolume* vol ) {
 422         if (quantContours[point] == NULL )
 423         {
 424                 createContours(point,vol);
 425         }
 426
 427         marAxisContours* axisCont = quantContours[point];
 428
 429         for (int i = 0; i < axisCont->getSize(); i++)
 430         {
 431                 double p1[2],p2[2];
 432                 marContour *marcontour=getContour(point,vol,i);
 433                 marContourVO* aContourVO = axisCont->getContour(i);
 434
 435                 marcontour->getMaximumLine(p1,p2);
 436                 vtkPoints *_vtkPoints = vtkPoints::New();
 437
 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 );
 447
 448                 if ( aContourVO->get2DDiameterMax()!=NULL ) { aContourVO->get2DDiameterMax()->Delete(); }
 449                         aContourVO->set2DDiameterMax(_vtkPoints);
 450         }
 451 }
 452
 453 // ----------------------------------------------------------------------------
 454 void marAxisCT::createSignals (int point, kVolume* vol) {
 455         if (quantContours[point] == NULL )
 456         {
 457                 create2Dcontours(point,vol); //AQUI
 458         //      filterSignal(vol);
 459         }
 460 }
 461
 462 // ----------------------------------------------------------------------------
 463 int     marAxisCT::getNumberOfContours(int point, kVolume* vol) {
 464         if (quantContours[point] == NULL )
 465         {
 466                 create2Dcontours(point,vol);//AQUI
 467         //      filterSignal(vol);
 468
 469         }
 470
 471         return quantContours[point]->getSize();
 472 }
 473
 474 // ----------------------------------------------------------------------------
 475 int     marAxisCT::getContourType(int point, int index, kVolume* vol) {
 476         if (quantContours[point] == NULL )
 477         {
 478                 create2Dcontours(point,vol);
 479         }
 480
 481         return quantContours[point]->getContourType(index);
 482
 483 }
 484
 485 // ----------------------------------------------------------------------------
 486 void marAxisCT::generatePoints(int point, kVolume* vol, std::vector<marIsocontour*> *list) {
 487
 488         vtkImageData* imagedata;
 489         int dirs = 72;
 490         double aPoint[ 3 ];
 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;
 499
 500         contoursLumen.push_back(NULL);
 501         contoursCalc.push_back(NULL);
 502
 503
 504
 505         //1. Check starting zone
 506         startZone = marAxisContours::LUMEN;
 507         actualZone = startZone;
 508         lumenCount = 0;
 509         calcCount = 0;
 510
 511         imagedata = (vtkImageData*) ((getSlice( point , vol ))->castVtk( ));
 512
 513         int limInfX, limInfY, limInfZ;
 514         int limSupX, limSupY, limSupZ;
 515         imagedata->GetExtent(limInfX,limSupX,limInfY,limSupY,limInfZ,limSupZ);
 516
 517         vtkImageCast* cast = vtkImageCast::New();
 518         cast->SetInput(imagedata);
 519         cast->SetOutputScalarTypeToDouble();
 520
 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);
 528
 529         //imagedata = filter->GetOutput();
 530
 531
 532         nPoint = getPoints(point);
 533
 534         aPoint[ 0 ]     =  (limSupX - limInfX) / 2; //nPoint[0];
 535         aPoint[ 1 ]     =  (limSupY - limInfY) / 2; //nPoint[1];
 536         aPoint[ 2 ]     =  (limSupZ - limInfZ) / 2; // nPoint[2];
 537         estDivided = false;
 538         estStart = false;
 539         //3. Generate rays over all directions
 540         for (int dir = 0; dir < dirs; dir++)
 541         {
 542
 543                 x = aPoint[ 0 ];
 544                 y = aPoint[ 1 ];
 545                 coordBase = 1;
 546                 prevIntensity = interpolate(x,y,imagedata);
 547                 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
 548
 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)))
 551                 {
 552
 553                         //2.1.1 Interpolate
 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);
 560
 561                         prevIntensity = resp;
 562
 563                         gradients.push_back(p);
 564
 565                         //2.1.3 Increase ray for next evaluation
 566                         generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
 567
 568                 }// end while gradient
 569
 570                 int index = 0;
 571                 found = false;
 572                 x = aPoint[ 0 ];
 573                 y = aPoint [ 1 ];
 574                 threshold = ((marParameters *) getParameters())->getContourThresh();
 575                 int mysize = gradients.size();
 576
 577                 //2.2 Compare every gradient value with threshold
 578                 while (index < gradients.size() && !found)
 579                 {
 580                         //2.2.1 Evaluate Lumen
 581                         if (actualZone == marAxisContours::LUMEN)
 582                         {
 583                                 average.push_back(gradients[index]->getIntensity());
 584
 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)))) {
 590
 591
 592
 593                                         if(gradients[index]->getGradient() > 0)
 594                                                 maxPoint = getMaximumGrad(gradients,index,gradients.size(),threshold);
 595                                         else
 596                                                 maxPoint = getMinimumGrad(gradients,index,gradients.size(),threshold);
 597
 598                                         index = maxPoint;
 599                                         x = gradients[maxPoint]->getX();
 600                                         y = gradients[maxPoint]->getY();
 601
 602                                         //2.2.1.1.1 Assign lumen structure
 603                                         isNeighbor = detectNeighbor(contoursLumen[lumenCount],dir,marAxisContours::LUMEN);
 604                                         if (!isNeighbor)
 605                                         {
 606                                                 lumenCount++;
 607                                         }
 608
 609                                         marIsocontour* cont = addPointToContour(contoursLumen[lumenCount],false,marAxisContours::LUMEN,gradients[maxPoint]);
 610                                         contoursLumen[lumenCount] = cont;
 611
 612                                         //2.2.1.1.2 Check step lumen->calcification
 613                                         double actualInt = interpolate(x,y,imagedata);
 614                                         double auxX = x;
 615                                         double auxY = y;
 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)))
 619                                         {
 620                                                 //2.2.1.1.2.1 Assign calcification structure
 621                                                 isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
 622                                                 if (!isNeighbor)
 623                                                 {
 624                                                         calcCount++;
 625                                                         contoursCalc.push_back(NULL);
 626                                                 }
 627
 628                                                 cont = addPointToContour(contoursCalc[calcCount],true,marAxisContours::CALCIFICATION,gradients[maxPoint]);
 629
 630                                                 contoursCalc[calcCount] = cont;
 631
 632                                                 actualZone = marAxisContours::CALCIFICATION;
 633                                                 //ASIGNACION varX???
 634
 635                                                 //2.2.1.1.2.3 Verify divided structure
 636                                                 if (dir == 0)
 637                                                 {
 638                                                         estStart = true;
 639                                                 }
 640
 641                                                 if ((dir == dirs - 1) && estStart)
 642                                                 {
 643                                                         estDivided = true;
 644                                                 }
 645
 646                                         }//2.2.1.1.3 Transition lumen->background
 647                                         else
 648                                         {
 649                                                 //2.2.1.1.3.1 Asign border structure
 650                                                 //ASIGNAR ESTRUCTURA
 651
 652                                                 //2.2.1.1.3.2 Generate stop condition
 653                                                 found = true;
 654                                                 actualZone = startZone;
 655                                                 coordBase = 1;
 656                                                 x = aPoint[ 0 ];
 657                                                 y = aPoint[ 1 ];
 658                                                 threshold = ((marParameters *) getParameters())->getContourThresh();
 659                                         }
 660
 661                                 }
 662                         }
 663                         //2.2.2 Evaluate calcification (MODIFIED IF THERE IS HYPODENSE)
 664                         else
 665                         {
 666                                 coordBase = 1;
 667                                 int calcIndex;
 668                                 double actualInt = interpolate(x,y,imagedata);
 669                                 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
 670                                 double next = interpolate(x,y,imagedata);
 671
 672                                 //2.2.2.1 If actual intensity is smaller than previous there is a contour
 673                                 if (actualInt < next)
 674                                 {
 675                                         while (actualInt < interpolate(x,y,imagedata) && (x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)))
 676                                         {
 677                                                 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
 678                                                 index++;
 679                                         }
 680
 681                                 }
 682
 683                                 generateVector(dir,&coordBase,aPoint[0],aPoint[1],&x,&y);
 684
 685                                 //2.2.2.2 Assign calcification structure
 686                                 isNeighbor = detectNeighbor(contoursCalc[calcCount],dir,marAxisContours::CALCIFICATION);
 687                                 if (!isNeighbor)
 688                                 {
 689                                         calcCount++;
 690                                         contoursCalc.push_back(NULL);
 691                                 }
 692
 693                                 //Adjust index to avoiod crash
 694                                 if (index >= gradients.size())
 695                                 {
 696                                         calcIndex = gradients.size() - 1;
 697                                 }
 698                                 else
 699                                 {
 700                                         calcIndex = index;
 701                                 }
 702
 703                                 marIsocontour* cont = addPointToContour(contoursCalc[calcCount],false,marAxisContours::CALCIFICATION,gradients[calcIndex]);
 704                                 contoursCalc[calcCount] = cont;
 705
 706                                 //2.2.2.3 Asign border structure
 707                                 //ASIGNACION BORDE
 708
 709                                 found = true;
 710                                 actualZone = startZone;
 711                                 x = aPoint[ 0 ];
 712                                 y = aPoint[ 1 ];
 713                                 threshold = ((marParameters *) getParameters())->getContourThresh();
 714
 715                         }//end if zones
 716
 717                         index++;
 718
 719                         if (index == gradients.size() && !found)
 720                         {
 721                                 threshold = threshold - 1;
 722                                 index = 0;
 723                         }
 724
 725                 }//end while evaluation
 726
 727                 //TODO - SACAR A UTILS
 728                 for (int ind = 0; ind < gradients.size(); ind++)
 729                 {
 730                         marPoint* p = gradients[ind];
 731                         gradients[ind] = NULL;
 732                         delete p;
 733
 734                 }
 735                 gradients.clear();
 736
 737         }//end for directions
 738
 739         FILE *archivo;
 740         char nomArchivo[30]; //, temp[3]
 741
 742         strcpy(nomArchivo,"./debug");
 743 //      itoa(point,temp,10);
 744 //      strcat(nomArchivo,temp);
 745         strcat(nomArchivo,".txt");
 746         archivo = fopen(nomArchivo,"w");
 747
 748
 749
 750
 751
 752 //      lumenContour[point] = contoursLumen[marAxisContours::LUMEN];
 753
 754         double  tempXi, tempYi;
 755         int polySides = contoursLumen[marAxisContours::LUMEN]->getSize();
 756
 757
 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);
 764         }
 765
 766
 767
 768         //4. Verify divided unified structure
 769         if (estDivided && calcCount > 0)
 770         {
 771                 unifyDividedStructure(&contoursCalc);
 772                 calcCount--;
 773         }
 774
 775         //5. Obtain statistics
 776         double lumenInt = getStatistics(contoursLumen[marAxisContours::LUMEN],aPoint[ 0 ],aPoint[ 1 ], imagedata);
 777         contoursLumen[marAxisContours::LUMEN]->setMaxIntensity(lumenInt);
 778
 779         if (contoursCalc[0] != NULL)
 780         {
 781                 for (int i = 0; i < contoursCalc.size(); i++)
 782                 {
 783                         lumenCount++;
 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];
 788                 }
 789         }
 790
 791
 792         *list = contoursLumen;
 793
 794         contoursCalc.clear();
 795     contoursLumen.clear();
 796         marUtils* util = new marUtils();
 797         double avg = util->obtainAverage(average);
 798         //FILE *archivo;
 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++)
 802         {
 803                 double cgx, cgy;
 804
 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());
 807         }
 808
 809         fprintf(archivo,"PROMEDIO %f - STD: %f\n",avg,util->obtainStandardDeviation(average,avg));
 810         fprintf(archivo,"%d\n",lumenContour.size());
 811         fclose(archivo);
 812         imagedata->Delete();
 813
 814 }
 815
 816 // ----------------------------------------------------------------------------
 817 int marAxisCT::getMaximumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
 818
 819         double max = list[iniPos]->getGradient();
 820         int coordMax = iniPos;
 821         int cont = iniPos + 1;
 822
 823         while (cont < limit && list[cont]->getGradient()>=threshold){
 824                 if (list[cont]->getGradient() > max){
 825                         max = list[cont]->getGradient();
 826                         coordMax = cont;
 827                 }
 828                 cont++;
 829         }
 830
 831         return coordMax;
 832 }
 833
 834 // ----------------------------------------------------------------------------
 835 int marAxisCT::getMinimumGrad(std::vector <marPoint *> list, int iniPos, int limit, double threshold) {
 836
 837         double min = list[iniPos]->getGradient();
 838         int coordMin = iniPos;
 839         int cont = iniPos + 1;
 840
 841         while (cont < limit && fabs(list[cont]->getGradient())>=threshold){
 842                 if (list[cont]->getGradient() < min){
 843                         min = list[cont]->getGradient();
 844                         coordMin = cont;
 845                 }
 846                 cont++;
 847         }
 848
 849         return coordMin;
 850 }
 851
 852 // ----------------------------------------------------------------------------
 853 double marAxisCT::interpolate(double x, double y, vtkImageData* imagedata) {
 854
 855         double a,b,c,d;
 856         double intensity;
 857         int i,j;
 858
 859         i = (int)floor(x);
 860         j = (int)floor(y);
 861
 862
 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);
 865
 866         b = imagedata->GetScalarComponentAsDouble(i,j+1,0,0) - c*i - imagedata->GetScalarComponentAsDouble(i,j,0,0);
 867
 868         a = imagedata->GetScalarComponentAsDouble(i+1,j,0,0) - c*j - imagedata->GetScalarComponentAsDouble(i,j,0,0);
 869
 870         d = imagedata->GetScalarComponentAsDouble(i,j,0,0) - a*i - b*j - c*i*j;
 871
 872         intensity = a*x + b*y + c*x*y + d;
 873
 874         return intensity;
 875 }
 876
 877 // ----------------------------------------------------------------------------
 878 void marAxisCT::generateVector(int dir,int *coordBase,double originX,double originY,double *x,double *y) {
 879
 880         int dirs = 72;
 881         int angle = (dir * (360 / dirs));
 882
 883         *x = 1 * cos(angle *(3.1415/180)) + (*x);
 884         *y = 1 * sin(angle *(3.1415/180)) + (*y);
 885
 886         (*coordBase)++;
 887 }
 888
 889 // ----------------------------------------------------------------------------
 890 bool marAxisCT::detectNeighbor(marIsocontour* contour,int dir, int type) {
 891
 892         bool isNeighbor = false;
 893
 894         if (contour == NULL){
 895                 isNeighbor = true;
 896         }else{
 897                 for (int i = 0; i < contour->getSize(); i++) {
 898
 899                         if (((contour->getDir(i) == (dir - 1)) || (contour->getDir(i) == dir))
 900                                     && (contour->getType() == type)){
 901                                 isNeighbor = true;
 902                         }
 903                 }
 904         }
 905
 906         return isNeighbor;
 907 }
 908
 909 // ----------------------------------------------------------------------------
 910 marIsocontour*  marAxisCT::addPointToContour(marIsocontour* cont ,bool inside,int type,marPoint* point) {
 911
 912
 913
 914         if (cont == NULL)
 915         {
 916                 cont = new marIsocontour();
 917                 cont->setType(type);
 918         }
 919
 920         cont->insertPoint(point->getX(),point->getY());
 921         cont->setDir(cont->getSize() - 1,point->getDirection());
 922         cont->setInside(cont->getSize() - 1,inside);
 923
 924         return cont;
 925
 926 }
 927
 928 // ----------------------------------------------------------------------------
 929 void marAxisCT::unifyDividedStructure(std::vector<marIsocontour*> *contList) {
 930
 931         marIsocontour* structOne = (*contList)[contList->size()-1];
 932         marIsocontour* structTwo = (*contList)[0];
 933
 934
 935         for (int i = 0; i < structOne->getSize(); i++)
 936         {
 937                 double x;
 938                 double y;
 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));
 943
 944         }
 945
 946         contList->pop_back();
 947 }
 948
 949 // ----------------------------------------------------------------------------
 950 double marAxisCT::getStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata) {
 951
 952         double max, auxX, auxY;
 953         int basis;
 954         max = interpolate(x,y,imagedata);
 955
 956         for (int i = 0; i < contour->getSize(); i++)
 957         {
 958                 basis = 1;
 959                 auxX = x;
 960                 auxY = y;
 961                 double pointX;
 962                 double pointY;
 963
 964                 contour->getPoint(i,&pointX,&pointY);
 965                 while (auxX != pointX || auxY != pointY)
 966                 {
 967                         double intensity = interpolate(auxX,auxY,imagedata);
 968                         if (max < intensity)
 969                         {
 970                                 max = intensity;
 971                         }
 972                         generateVector(contour->getDir(i),&basis,x,y,&auxX,&auxY);
 973
 974                 }
 975         }
 976
 977         return max;
 978 }
 979
 980
 981 // ----------------------------------------------------------------------------
 982 double marAxisCT::getCalcStatistics(marIsocontour* contour, double x, double y, vtkImageData* imagedata, int dir) {
 983
 984
 985         double max, auxX,auxY;
 986         int basis;
 987
 988         int i;
 989         FILE *archivo;
 990         char nomArchivo[30], temp[3];
 991
 992         strcpy(nomArchivo,"./statistics");
 993
 994 // EED 06/06/2007 linux
 995 //      itoa(dir,temp,10);
 996         sprintf(temp,"%d",dir);
 997
 998         strcat(nomArchivo,temp);
 999         strcat(nomArchivo,".txt");
1000
1001
1002         archivo = fopen(nomArchivo,"w");
1003         i = 0;
1004         max = 0;
1005         fprintf(archivo,"TAMAÑO %d\n",contour->getSize());
1006
1007         while (i < contour->getSize())
1008         {
1009                 basis = 1;
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);
1014
1015                 while (auxX != pointX || auxY != pointY)
1016                 {
1017                         double intensity = interpolate(auxX,auxY,imagedata);
1018                         if (max < intensity){
1019                                 max = intensity;
1020                         }
1021                         generateVector(contour->getDir(i),&basis,auxX,auxY,&auxX,&auxY);
1022                 }
1023
1024                 i +=2;
1025         }
1026         fclose(archivo);
1027         return max;
1028 }
1029
1030 // ----------------------------------------------------------------------------
1031 int marAxisCT::round(double num){
1032
1033         float decimal;
1034         int resp;
1035
1036         decimal = num - floor(num) * 1.0000;
1037
1038         if (decimal >= 0.50000)
1039                 resp = (int)ceil(num);
1040         else
1041                 resp = (int)floor(num);
1042
1043         return resp;
1044
1045 }
1046
1047 // ----------------------------------------------------------------------------
1048 void marAxisCT::createEmptyContours(){
1049
1050         int nCnts = ( int ) getNumberOfSplinePoints( );
1051
1052
1053         for (int i = 0; i < nCnts; i++)
1054         {
1055                 quantContours.push_back( NULL );
1056                 lumenContour.push_back( NULL);
1057         }
1058 }
1059
1060 // ----------------------------------------------------------------------------
1061 void marAxisCT::updateLumenPercentage(int point, kVolume* vol)
1062 {
1063         marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1064         if (quantContours[point] == NULL || wallCont->get2DContour() == NULL)
1065         {
1066                 create2Dcontours(point,vol);
1067         }
1068
1069         if (wallCont->getType() == marAxisContours::WALL
1070                         && wallCont->isReplaced() == false)
1071         {
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();
1080         }
1081
1082 }
1083
1084 // ----------------------------------------------------------------------------
1085 void marAxisCT::updateCalcPercentage(int point, kVolume* vol)
1086 {
1087         double oldValue;
1088         int posMax;
1089
1090         if (quantContours[point] == NULL)
1091         {
1092                 create2Dcontours(point,vol);
1093         }
1094
1095         posMax = getMaxIntensity(point);
1096
1097         if (quantContours[point]->getSize() > 1)
1098         {
1099                 oldValue = quantContours[point]->getContour(posMax)->getSignal();
1100         }
1101
1102         for (int i = 1; i < quantContours[point]->getSize(); i++)
1103         {
1104
1105                 if (quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION
1106                                 && quantContours[point]->getContour(i)->isReplaced() == false)
1107                 {
1108
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();
1113                 }
1114
1115         }
1116
1117 }
1118
1119 // ----------------------------------------------------------------------------
1120 int     marAxisCT::getMaxIntensity(int point)
1121 {
1122         int pos = 1;
1123         double intMax;
1124         if (quantContours[point]->getSize() > 2)
1125         {
1126                 intMax = quantContours[point]->getContour(0)->getSignal();
1127                 for (int i = 0; i < quantContours[point]->getSize(); i++)
1128                 {
1129                         if (quantContours[point]->getContour(i)->getSignal() > intMax && quantContours[point]->getContour(i)->getType() == marAxisContours::CALCIFICATION)
1130                         {
1131                                 intMax = quantContours[point]->getContour(i)->getSignal() ;
1132                                 pos = i;
1133                         }
1134                 }
1135         }
1136
1137         return pos;
1138
1139 }
1140 // ----------------------------------------------------------------------------
1141 void marAxisCT::eraseContours()
1142 {
1143
1144
1145                 int nCts = quantContours.size();
1146                 vesselPoints.clear();
1147
1148                 int i = 0;
1149
1150                 while (i < nCts)
1151                 {
1152                         delete quantContours[i];
1153                         quantContours[i] = NULL;
1154                         lumenContour[i] = NULL;
1155                         i++;
1156                 }
1157                 quantContours.clear();
1158                 lumenContour.clear();
1159
1160
1161 }
1162
1163 // ----------------------------------------------------------------------------
1164 void marAxisCT::eraseContoursPartial(int start)
1165 {
1166
1167         int nCts = 0;
1168         int i = start;
1169
1170         while (i >= nCts)
1171         {
1172                 delete quantContours[i];
1173                 quantContours[i] = NULL;
1174                 lumenContour[i] = NULL;
1175                 i--;
1176         }
1177         quantContours.clear();
1178         lumenContour.clear();
1179
1180         //TODO Revisar esto que no me convence
1181 }
1182
1183
1184 // ----------------------------------------------------------------------------
1185 bool marAxisCT::pointInPolygon(marIsocontour *c, double x, double y){
1186         int      i, j=0;
1187         bool  oddNODES=false;
1188         int polySides = c->getSize();
1189         double *polyX, *polyY, tempX, tempY;
1190
1191         polyX = (double *) malloc(polySides*sizeof(double));
1192         polyY = (double *) malloc(polySides*sizeof(double));
1193
1194         for (i=0; i<polySides; i++) {
1195                 c->getPoint(i,&tempX,&tempY);
1196                 polyX[i] = tempX;
1197                 polyY[i] = tempY;
1198         }
1199
1200         for (i=0; i<polySides; i++) {
1201                 j++;
1202                 if (j==polySides)
1203                         j=0;
1204
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) {
1208                                 oddNODES=!oddNODES;
1209                         }
1210                 }
1211         }
1212
1213         free(polyX);
1214         free(polyY);
1215   return oddNODES;
1216
1217 }
1218
1219 // ----------------------------------------------------------------------------
1220 marIsocontour* marAxisCT::parsePolyDataToMarIsocontour(marContourVO* contourVO){
1221
1222         vtkPolyData* polyDataCnt;
1223         vtkIdList* pointsId;
1224         double* pointPolyDataCnt;
1225         int numPointsPolyDataCnt, numTemp;
1226         marIsocontour *cont;
1227
1228         if (contourVO->isReplaced())
1229         {
1230                 polyDataCnt = contourVO->get2DContour();
1231         }
1232         else
1233         {
1234                 polyDataCnt =    contourVO->getIsocontourStripped()->GetOutput();
1235         }
1236
1237         numPointsPolyDataCnt = polyDataCnt->GetNumberOfPoints();
1238         pointsId = polyDataCnt->GetCell(0)->GetPointIds();
1239         numTemp = pointsId->GetNumberOfIds();
1240
1241         cont = new marIsocontour();
1242
1243         for (int i = 0; i < (numTemp -1); i++){
1244                 pointPolyDataCnt = polyDataCnt->GetPoint(pointsId->GetId(i));
1245                 cont->insertPoint(pointPolyDataCnt[0], pointPolyDataCnt[1]);
1246
1247         }
1248
1249         return cont;
1250 }
1251
1252
1253
1254 // ----------------------------------------------------------------------------
1255 marIsocontour* marAxisCT::filterContour(marIsocontour* contExt, marIsocontour* contInt, double radio){
1256
1257         double tmp_eed;
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;
1263
1264         marIsocontour* newCont = new marIsocontour();
1265
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 );
1272         }
1273
1274         count = 0;
1275         ultima = 0;
1276         while (count < lineas.size() - 1){
1277                 countAux = 0;
1278
1279                 marLine* l1 = lineas[count];
1280
1281                 hayInterseccion = false;
1282                 while (countAux < lineas.size()){
1283                         marLine *l2 = lineas[countAux];
1284
1285                         l1->getIntersect(l2->getA(),l2->getB(),l2->getC(),&xInt,&yInt);
1286
1287
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);
1293
1294                                         //Distancia menor?
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))))
1297                                         {
1298                                                 ultima = countAux;
1299                                                 hayInterseccion = true;
1300
1301                                         }
1302
1303                                 }
1304                         }
1305                         countAux++;
1306                 }
1307
1308
1309                 if (!hayInterseccion){
1310                         contInt->getPoint(count,&x2,&y2);
1311                         newCont->insertPoint(x2,y2);
1312                 }
1313                         count++;
1314
1315         }
1316
1317         return newCont;
1318 }
1319
1320 // ----------------------------------------------------------------------------
1321 void marAxisCT::markUpLumen(int point, kVolume* vol)
1322 {
1323
1324         marContourVO* aVO = this->searchContour(marAxisContours::ELUMEN,point);
1325
1326         if (aVO != NULL)
1327         {
1328                 return;
1329         }
1330
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();
1339
1340
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();
1347         image->Update();
1348
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();
1355
1356         vtkImageContinuousErode3D *imageErode3D = vtkImageContinuousErode3D::New();
1357         imageErode3D->SetInput(imageThreshold->GetOutput());
1358         imageErode3D->SetKernelSize(4,4,1);
1359         imageErode3D->Update();
1360         erodedata = imageErode3D->GetOutput();
1361
1362 //      polygons.push_back(contExt);
1363         for (int i = 0; i < quantContours[point]->getSize(); i++)
1364         {
1365                 polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
1366         }
1367
1368         //7- Wipe image to detect lumen
1369
1370         int x;
1371         for (x = limInfX + 30; x < (limSupX - 1); x++)
1372         {
1373
1374                 for (int y = limInfY; y < (limSupY - 1) ; y++)
1375                 {
1376                         bool inWall = false;
1377                         bool outCalc = true;
1378
1379                         for (int numConts = 0; numConts < polygons.size(); numConts++)
1380                         {
1381                                 if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
1382                                 {
1383                                         if (numConts == 0)
1384                                         {
1385                                                 inWall = true;
1386                                         }
1387
1388                                         else
1389                                         {
1390                                                 outCalc = false;
1391                                         }
1392
1393                                 }
1394
1395                         }
1396
1397                         if (inWall && outCalc)
1398                         {
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);
1403                         }
1404
1405                 /*      if (cont >= (limSupY - limSupX))
1406                         {
1407                                 contin = false;
1408                         }*/
1409                 }
1410         }
1411
1412         if (points.size() > 0)
1413         {
1414                 average = average / points.size();
1415                 for (int k = 0; k < points.size(); k++)
1416                 {
1417                         marPoint *p = points[k];
1418                         double actualInt = p->getIntensity();
1419                         double temp = average - actualInt;
1420                         deviation += pow(temp,2.0);
1421                 }
1422                 deviation = sqrt(deviation / points.size());
1423
1424         }
1425
1426
1427         polygons.push_back(lumenContour[point]);
1428
1429
1430         for (x = limInfX; x < (limSupX - 1); x++)
1431         {
1432                 for (int y = limInfY; y < (limSupY - 1); y++)
1433                 {
1434                         bool inWall = false;
1435                         bool outCalc = true;
1436                         bool inLumen = false;
1437                         unsigned short *refValue = (unsigned short *) image->GetScalarPointer(x,y,0);
1438                         *refValue = 0;
1439
1440                         for (int numConts = 0; numConts < polygons.size(); numConts++)
1441                         {
1442                                 bool adentro = pointInPolygon(polygons[numConts],x,y);
1443                                 if (pointInPolygon(polygons[numConts],x,y) && erodedata->GetScalarComponentAsDouble(x,y,0,0) > 0)
1444                                 {
1445                                         if (numConts == 0)
1446                                         {
1447                                                 inWall = true;
1448                                         }
1449                                         else if (numConts == polygons.size() - 1)
1450                                         {
1451                                                 inLumen = true;
1452
1453                                         }
1454                                         else
1455                                         {
1456                                                 outCalc = false;
1457                                         }
1458                                 }
1459                         }
1460
1461                         double intensidad = imagedata->GetScalarComponentAsDouble(x,y,0,0);
1462
1463                         if (inWall && outCalc && inLumen)
1464                         {
1465
1466                                 if (imagedata->GetScalarComponentAsDouble(x,y,0,0) > (average - 1.5*deviation)
1467                                          && imagedata->GetScalarComponentAsDouble(x,y,0,0) < (average + 2.5*deviation))
1468                                 {
1469                                         *refValue = 255;
1470                                 }
1471
1472                         }
1473
1474                 }
1475         }
1476
1477
1478
1479         extractLumen(image,lumenContour[point],point);
1480
1481 //      fclose(archivo);
1482 }
1483
1484 // ----------------------------------------------------------------------------
1485 double marAxisCT::obtainContourArea(marIsocontour* contour)
1486 {
1487
1488         double area = 0.0;
1489         double x1,y1,x2,y2;
1490
1491         for (int i = 0; i < contour->getSize();i++)
1492         {
1493                 contour->getPoint(i,&x1,&y1);
1494                 if (i < (contour->getSize() - 1))
1495                 {
1496                         contour->getPoint(i+1,&x2,&y2);
1497                 }
1498                 else
1499                 {
1500                         contour->getPoint(0,&x2,&y2);
1501                 }
1502
1503                 area += (x1 * y2);
1504                 area -= (y1 * x2);
1505         }
1506
1507         area = area / 2;
1508
1509         if (area < 0)
1510                 area = -area;
1511
1512         return area;
1513 }
1514
1515 // ----------------------------------------------------------------------------
1516 void marAxisCT::adjustWall(int point, kVolume* vol)
1517 {
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);
1523
1524         //Obtain first area
1525         ((marParameters *) getParameters())->setLumenPercentage(50);
1526         updateLumenPercentage(point,vol);
1527         cont = parsePolyDataToMarIsocontour(wallCont);
1528         previousArea = obtainContourArea(cont);
1529
1530         for (int eval = 49; eval >= 0; eval--)
1531         {
1532                 ((marParameters *) getParameters())->setLumenPercentage(eval);
1533                 updateLumenPercentage(point,vol);
1534                 cont = parsePolyDataToMarIsocontour(wallCont);
1535                 actualArea = obtainContourArea(cont);
1536                 grad.push_back(fabs(actualArea - previousArea));
1537         }
1538
1539         int newPercentage = 50 - maxValue(grad);
1540         ((marParameters *) getParameters())->setLumenPercentage(newPercentage);
1541         updateLumenPercentage(point,vol);
1542
1543
1544
1545 }
1546
1547 // ----------------------------------------------------------------------------
1548 void marAxisCT::adjustCalcification(int point, kVolume* vol)
1549 {
1550 }
1551
1552 // ----------------------------------------------------------------------------
1553 int     marAxisCT::maxValue(std::vector <double> list)
1554 {
1555         int max = 0;
1556         double maxVal = list[0];
1557
1558         for (int i = 0; i < list.size(); i++)
1559         {
1560                 if (list[i] > maxVal)
1561                 {
1562                         maxVal = list[i];
1563                         max = i;
1564                 }
1565         }
1566
1567         return max;
1568 }
1569
1570
1571 // ----------------------------------------------------------------------------
1572 int marAxisCT::getDiameter(marContourVO* contourVO, double x, double y, int limInfX, int limInfY, int limSupX, int limSupY)
1573 {
1574         bool stop = false;
1575         double originX = x;
1576         double originY = y;
1577
1578         int coordBase = 1, diameterOne = 0, diameterTwo = 0;
1579         marIsocontour* cont = parsePolyDataToMarIsocontour(contourVO);
1580         generateVector(1,&coordBase,originX,originY,&x,&y);
1581
1582         while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1583         {
1584                 generateVector(0,&coordBase,originX,originY,&x,&y);
1585                 if (pointInPolygon(cont,x,y))
1586                 {
1587                         diameterOne++;
1588                 }
1589                 else
1590                 {
1591                         stop = true;
1592                 }
1593         }
1594
1595         stop = false;
1596         x = originX;
1597         y = originY;
1598         while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1599         {
1600                 generateVector(35,&coordBase,originX,originY,&x,&y);
1601                 if (pointInPolygon(cont,x,y))
1602                 {
1603                         diameterOne++;
1604                 }
1605                 else
1606                 {
1607                         stop = true;
1608                 }
1609         }
1610
1611         stop = false;
1612         x = originX;
1613         y = originY;
1614         while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1615         {
1616                 generateVector(17,&coordBase,originX,originY,&x,&y);
1617                 if (pointInPolygon(cont,x,y))
1618                 {
1619                         diameterTwo++;
1620                 }
1621                 else
1622                 {
1623                         stop = true;
1624                 }
1625         }
1626
1627         stop = false;
1628         x = originX;
1629         y = originY;
1630         while ((x > (limInfX+1) && x < (limSupX-1)) && (y > (limInfY+1) && y < (limSupY-1)) && !stop)
1631         {
1632                 generateVector(107,&coordBase,originX,originY,&x,&y);
1633                 if (pointInPolygon(cont,x,y))
1634                 {
1635                         diameterTwo++;
1636                 }
1637                 else
1638                 {
1639                         stop = true;
1640                 }
1641         }
1642
1643         delete cont;
1644
1645         return (diameterOne + diameterTwo) / 2;
1646 }
1647
1648 // ----------------------------------------------------------------------------
1649 void marAxisCT::histogram(int point, kVolume* vol)
1650 {
1651
1652         vtkImageData* imagedata;
1653         int limInfX, limInfY, limInfZ;
1654         int limSupX, limSupY, limSupZ;
1655         int pos;
1656
1657
1658         FILE *archivo;
1659         char nomArchivo[30],temp[3];
1660         marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1661
1662
1663 //      for (int point = 0; point <  quantContours.size(); point++)
1664 //      {
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);
1670
1671                 std::vector<double> intensity;
1672                 std::vector<int> frecuencies;
1673
1674                 for (int i = 0; i < limSupX; i++)
1675                 {
1676                         for (int j = 0; j < limSupY; j++)
1677                         {
1678                                 if(pointInPolygon(contExt,i,j))
1679                                 {
1680
1681                                         pos = searchData(intensity,imagedata->GetScalarComponentAsDouble(i,j,0,0));
1682                                         if (pos == -1)
1683                                         {
1684                                                 intensity.push_back( imagedata->GetScalarComponentAsDouble(i,j,0,0));
1685                                                 frecuencies.push_back(1);
1686
1687
1688                                         }
1689                                         else
1690                                         {
1691                                                 frecuencies[pos]++;
1692                                         }
1693                                 }
1694
1695                         }
1696                 }
1697
1698                 strcpy(nomArchivo,"./histogram");
1699
1700 // EED 06/06/2007 linux
1701 //              itoa(point,temp,10);
1702                 sprintf(temp,"%d",point);
1703
1704                 strcat(nomArchivo,temp);
1705                 strcat(nomArchivo,".csv");
1706
1707
1708
1709                 archivo = fopen(nomArchivo,"w");
1710                 for (int k = 0; k < frecuencies.size(); k++)
1711                 {
1712
1713                         fprintf(archivo,"%f;%d\n",intensity[k],frecuencies[k]);
1714
1715                 }
1716                 fclose(archivo);
1717
1718
1719
1720
1721                 //imagedata->Delete();
1722                 delete contExt;
1723 //      }
1724
1725 }
1726
1727 // ----------------------------------------------------------------------------
1728 int     marAxisCT::searchData(std::vector<double> vec, double value)
1729 {
1730         int resp = -1;
1731         bool found = false;
1732
1733         for (int i = 0; i< vec.size() && !found; i++)
1734         {
1735                 if (vec[i] == value)
1736                 {
1737                         resp = i;
1738                         found = true;
1739                 }
1740         }
1741
1742         return resp;
1743 }
1744
1745 // ----------------------------------------------------------------------------
1746 void marAxisCT::setCalibration(bool calib)
1747 {
1748         calibration = calib;
1749 }
1750
1751 // ----------------------------------------------------------------------------
1752 bool marAxisCT::getCalibration()
1753 {
1754         return calibration;
1755 }
1756
1757 // ----------------------------------------------------------------------------
1758 void marAxisCT::adjustContour(int point, kVolume* vol)
1759 {
1760
1761         double percentage = 0.05;
1762         double prevDifference = 0;
1763         if (quantContours[point + 1] == NULL)
1764         {
1765                 create2Dcontours(point + 1 ,vol);
1766         }
1767
1768         //1. Obtain actual signal
1769         double signal = maxSignal;
1770         marContourVO* wallCont = searchContour(marAxisContours::WALL,point);
1771
1772         //2. Obtain previous contour
1773         marIsocontour* prev = parsePolyDataToMarIsocontour(searchContour(marAxisContours::WALL, point + 1));
1774
1775         //3. Adjust actual contour
1776         double newValue = signal;
1777         wallCont->getIsocontour()->SetValue(0,newValue);
1778         wallCont->getIsocontour()->Update();
1779         wallCont->get2DContour()->Update();
1780
1781         //4. Obtain points
1782         marIsocontour* actual = parsePolyDataToMarIsocontour(wallCont);
1783
1784         //5. Compare areas
1785         double prevA = obtainContourArea(prev);
1786         double actuA = obtainContourArea(actual);
1787
1788         prevDifference = fabs(prevA - actuA);
1789
1790
1791         if (prevA - actuA > 0 && prevA - actuA < percentage*prevA) //CASO 1
1792         {
1793                 bool stop = false;
1794                 while (!stop && newValue > 0)
1795                 {
1796                         newValue--;
1797
1798                         wallCont->getIsocontour()->SetValue(0,newValue);
1799                         wallCont->getIsocontour()->Update();
1800                         wallCont->get2DContour()->Update();
1801
1802                         actual = parsePolyDataToMarIsocontour(wallCont);
1803                         actuA = obtainContourArea(actual);
1804
1805                         if (fabs(prevA - actuA) > percentage*prevA)
1806                         {
1807                                 newValue++;
1808                                 wallCont->getIsocontour()->SetValue(0,newValue);
1809                                 wallCont->getIsocontour()->Update();
1810                                 wallCont->get2DContour()->Update();
1811                                 stop = true;
1812                         }
1813                 }
1814
1815         }
1816         else if (prevA - actuA < 0 && prevA - actuA > -percentage*prevA) //CASO 2
1817         {
1818
1819                 bool stop = false;
1820                 while (!stop && newValue > 0)
1821                 {
1822                         newValue--;
1823
1824                         wallCont->getIsocontour()->SetValue(0,newValue);
1825                         wallCont->getIsocontour()->Update();
1826                         wallCont->get2DContour()->Update();
1827
1828                         actual = parsePolyDataToMarIsocontour(wallCont);
1829                         actuA = obtainContourArea(actual);
1830
1831                         if (prevA - actuA < -percentage*prevA)
1832                         {
1833                                 newValue++;
1834                                 wallCont->getIsocontour()->SetValue(0,newValue);
1835                                 wallCont->getIsocontour()->Update();
1836                                 wallCont->get2DContour()->Update();
1837                                 stop = true;
1838                         }
1839                 }
1840         }
1841         else if (prevA - actuA < 0 && prevA - actuA < -percentage*prevA) //CASO 3
1842         {
1843
1844                 bool stop = false;
1845
1846                 while (!stop)
1847                 {
1848                         newValue++;
1849
1850                         wallCont->getIsocontour()->SetValue(0,newValue);
1851                         wallCont->getIsocontour()->Update();
1852                         wallCont->get2DContour()->Update();
1853
1854                         actual = parsePolyDataToMarIsocontour(wallCont);
1855
1856                         actuA = obtainContourArea(actual);
1857                         if (prevA - actuA > -percentage*prevA)
1858                         {
1859                                 if (prevDifference < fabs(prevA - actuA))
1860                                 {
1861                                         newValue--;
1862                                         wallCont->getIsocontour()->SetValue(0,newValue);
1863                                         wallCont->getIsocontour()->Update();
1864                                         wallCont->get2DContour()->Update();
1865                                 }
1866
1867                                 stop = true;
1868
1869                         }
1870                         prevDifference = fabs(prevA - actuA);
1871                 }
1872         }
1873         else if (prevA - actuA > 0 && prevA - actuA > percentage*prevA) //CASO 4
1874         {
1875                 bool stop = false;
1876                 while (!stop && newValue > 0)
1877                 {
1878                         newValue--;
1879
1880                         wallCont->getIsocontour()->SetValue(0,newValue);
1881                         wallCont->getIsocontour()->Update();
1882                         wallCont->get2DContour()->Update();
1883
1884                         actual = parsePolyDataToMarIsocontour(wallCont);
1885                         actuA = obtainContourArea(actual);
1886
1887                         if (prevA - actuA < percentage*prevA)
1888                         {
1889                                 if (prevDifference < fabs(prevA - actuA))
1890                                 {
1891                                         newValue++;
1892                                         wallCont->getIsocontour()->SetValue(0,newValue);
1893                                         wallCont->getIsocontour()->Update();
1894                                         wallCont->get2DContour()->Update();
1895                                 }
1896
1897                                 stop = true;
1898
1899                         }
1900                         prevDifference = fabs(prevA - actuA);
1901                 }
1902         }
1903
1904
1905         maxSignal = newValue;
1906
1907         this->markUpLumen(point,vol);
1908 }
1909
1910 // ----------------------------------------------------------------------------
1911 int marAxisCT::getStartIndex()
1912 {
1913         return startIndex;
1914 }
1915
1916 // ----------------------------------------------------------------------------
1917 void marAxisCT::setStartIndex(int start)
1918 {
1919         startIndex = start;
1920
1921         int per = ((marParameters *) getParameters())->getLumenPercentage();
1922         double oldValue = quantContours[start]->getContour(0)->getSignal()*per*pow(10.0,-2.0);
1923
1924         maxSignal = oldValue; // / (per*pow(10.0,-2.0));
1925 }
1926
1927 // ----------------------------------------------------------------------------
1928 int     marAxisCT::getPointSize()
1929 {
1930         return vesselPoints.size();
1931 }
1932
1933 // ----------------------------------------------------------------------------
1934 marPoint* marAxisCT::getPoint(int i)
1935 {
1936         return vesselPoints[i];
1937 }
1938
1939 // ----------------------------------------------------------------------------
1940 void marAxisCT::extractLumen(vtkImageData *lumenImage, marIsocontour *lumenContour, int point)
1941 {
1942
1943
1944         vtkImageCast* cast = vtkImageCast::New();
1945         cast->SetInput(lumenImage);
1946         cast->SetOutputScalarTypeToDouble();
1947
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);
1955
1956
1957
1958         marAxisContours* axisCont = quantContours[point];
1959
1960         int pos = -1;
1961
1962         for (int i = 0; i < axisCont->getSize()  && pos == -1; i++)
1963         {
1964                 if (axisCont->getContourType(i) == marAxisContours::ELUMEN)
1965                 {
1966                         pos = i;
1967                 }
1968         }
1969
1970         bool stop = false;
1971         if (pos != -1)          //Verify that contour hasn't been replaced
1972         {
1973                 stop = true;
1974                 if (axisCont->isReplaced(pos))
1975                 {
1976                         stop = true;
1977                 }
1978
1979         }
1980
1981         if (!stop)
1982         {
1983                 marContourVO* contourVO = new marContourVO();
1984                 marIsocontour* iso = lumenContour;
1985
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( );
1991
1992                 contourVO->setIsocontourCpd(vtkCleanPolyData::New());
1993                 contourVO->getIsocontourCpd()->SetInput( contourVO->getIsocontour()->GetOutput( ) );
1994                 contourVO->getIsocontourCpd()->ConvertLinesToPointsOff( );
1995                 contourVO->getIsocontourCpd()->Update( );
1996
1997                 double cgx, cgy;
1998
1999                 iso->getCG(&cgx,&cgy);
2000
2001                 contourVO->setIsocontourDcf(vtkPolyDataConnectivityFilter::New( ));
2002                 contourVO->getIsocontourDcf()->SetInput( contourVO->getIsocontourCpd()->GetOutput( ) );
2003
2004                 contourVO->getIsocontourDcf()->SetClosestPoint( cgx, cgy, 0 );
2005                 contourVO->getIsocontourDcf()->SetExtractionModeToClosestPointRegion( );
2006                 contourVO->getIsocontourDcf()->Update( );
2007
2008                 contourVO->setIsocontourCpd2(vtkCleanPolyData::New());
2009                 contourVO->getIsocontourCpd2()->SetInput( contourVO->getIsocontourDcf()->GetOutput( ) );
2010                 contourVO->getIsocontourCpd2()->Update();
2011
2012                 contourVO->setIsocontourStripped(vtkStripper::New( ));
2013                 contourVO->getIsocontourStripped()->SetInput( contourVO->getIsocontourCpd2()->GetOutput() );
2014                 contourVO->getIsocontourStripped()->Update();
2015
2016                 contourVO->set2DContour(contourVO->getIsocontourStripped()->GetOutput());
2017                 contourVO->get2DContour()->Update();
2018                 contourVO->setSignal(255);
2019                 contourVO->setType(marAxisContours::ELUMEN);
2020
2021
2022
2023                 marIsocontour* actual = parsePolyDataToMarIsocontour(contourVO);
2024                 marContourVO* aVO = searchContour(marAxisContours::WALL,point);
2025                 double wallArea = 0;
2026                 if (aVO != NULL)
2027                 {
2028                         wallArea = obtainContourArea(parsePolyDataToMarIsocontour(aVO));
2029                 }
2030
2031                 double prevA = obtainContourArea(actual);
2032
2033                 int newValue = 128;
2034                 int finalValue = 128;
2035                 double actuA = 1;
2036                 while (newValue > 0 && actuA > 0)
2037                 {
2038
2039                         newValue--;
2040
2041                         try
2042                         {
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)
2049                                 {
2050                                         finalValue = newValue;
2051                                         prevA = actuA;
2052                                 }
2053                                 else if (actuA > prevA && (prevA / 1000) < 1 && actuA < 0.1*wallArea)
2054                                 {
2055                                         finalValue = newValue;
2056                                         prevA = actuA;
2057                                 }
2058                         }
2059                         catch (std::exception e)
2060                         {
2061                                 actuA = 0;
2062                         }
2063
2064                 }
2065
2066                 contourVO->getIsocontour()->SetValue(0,finalValue);
2067                 contourVO->getIsocontour()->Update();
2068                 contourVO->get2DContour()->Update();
2069
2070
2071
2072                 /*iso->getType()*/
2073
2074                 if (pos != -1)
2075                 {
2076                         axisCont->replaceContour(contourVO,pos);
2077                 }
2078                 else
2079                 {
2080                         axisCont->addContour(contourVO);
2081                 }
2082                 //      axisCont->addContour(contourVO);
2083
2084
2085                 quantContours[point] = axisCont;
2086         }
2087
2088
2089 }
2090
2091 // ----------------------------------------------------------------------------
2092 void marAxisCT::generateFile(int point,kVolume* vol)
2093 {
2094         FILE *archivo;
2095         char nomArchivo[30], temp[3];
2096
2097         strcpy(nomArchivo,"./point");
2098
2099 // EED 06/06/2007 linux
2100 //      itoa(point,temp,10);
2101         sprintf(temp,"%d",point);
2102
2103         strcat(nomArchivo,temp);
2104         strcat(nomArchivo,".txt");
2105         archivo = fopen(nomArchivo,"w");
2106
2107         fprintf(archivo, "%d", quantContours[point]->getSize());
2108
2109         int i;
2110         for (i = 0; i < quantContours[point]->getSize(); i++)
2111         {
2112                 marIsocontour* iso =parsePolyDataToMarIsocontour(quantContours[point]->getContour(i));
2113                 fprintf(archivo, " %d",iso->getSize());
2114
2115         }
2116         fprintf(archivo,"\n");
2117
2118         for ( i = 0; i < quantContours[point]->getSize(); i++)
2119         {
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++)
2123                 {
2124                         double tempX, tempY;
2125                         iso->getPoint(j,&tempX,&tempY);
2126                         fprintf(archivo,"%f;%f;%d\n",tempX,tempY,point);
2127                 }
2128         }
2129
2130 // EED 30 Janvier 2007
2131         vtkImageData *vtkimagedata = this->getSliceImage(point,vol);
2132
2133         std::string directory = "./";
2134         std::string filename  = "xx";
2135         float rescalaSlope           =  1;
2136         float rescalaIntercept       =  0;
2137         int dim[3];
2138         vtkimagedata->GetDimensions(dim);
2139         int voi[6];
2140         voi[0]=0;
2141         voi[1]=dim[0];
2142         voi[2]=0;
2143         voi[3]=dim[1];
2144         voi[4]=0;
2145         voi[5]=dim[2];
2146         marRAW2Files marraw2;
2147         marraw2.saveVolume(directory,filename,vtkimagedata,voi,rescalaSlope,rescalaIntercept);
2148
2149
2150
2151         fclose(archivo);
2152
2153 }
2154
2155 // ----------------------------------------------------------------------------
2156 void marAxisCT::replaceContour2D(int point,int size,double *vx,double *vy, int type)
2157 {
2158
2159         vtkPoints *_pts = vtkPoints::New();
2160         _pts->SetNumberOfPoints(size);
2161         int j;
2162
2163         for (j=0 ; j<size ; j++){
2164                 _pts->SetPoint(j,       vx[j]   , vy[j] , 0 );
2165         }
2166 //      _pts->SetPoint(0,       vx[0]   , vy[0] , 0 );
2167
2168         vtkCellArray *lines = vtkCellArray::New();
2169         lines->InsertNextCell( size );
2170         for ( j=0 ; j<size+1 ; j++ ){
2171                 lines->InsertCellPoint(j % size );
2172         }
2173
2174         vtkPolyData *_pd = vtkPolyData::New();
2175         _pd->SetPoints( _pts );
2176         _pd->SetLines( lines );
2177         lines->Delete();  //do not delete lines ??
2178         _pts->Delete();
2179
2180         marContourVO* vo = new marContourVO();///this->searchContour(type,point);
2181         vo->setReplaced(true);
2182         vo->set2DContour(_pd);
2183         vo->setType(type);
2184         quantContours[point]->addContour(vo);   //->replaceContour(vo,i);
2185 /*      _2Dcontours[i]=_pd;
2186         createContour(i,NULL);*/
2187 }
2188
2189 // ----------------------------------------------------------------------------
2190 void marAxisCT::cleanContours(int type, int point)
2191 {
2192         marAxisContours* newContours = new marAxisContours();
2193         if (quantContours[point] != NULL)
2194         {
2195                 for (int i = 0; i < quantContours[point]->getSize(); i++)
2196                 {
2197                         if (quantContours[point]->getContourType(i) != type)
2198                         {
2199                                 newContours->addContour(quantContours[point]->getContour(i));
2200                         }
2201                         else
2202                         {
2203                         /*      marContourVO* vo = quantContours[point]->getContour(i);
2204                                 vo->set2DContour(NULL);
2205                                 newContours->replaceContour(vo,i);*/
2206                         }
2207                 }
2208         }
2209
2210         quantContours[point] = newContours;
2211 }
2212
2213 // ----------------------------------------------------------------------------
2214 marContourVO*   marAxisCT::searchContour(int type, int point)
2215 {
2216         marContourVO* cont = NULL;
2217
2218         for (int i = 0; i < quantContours[point]->getSize(); i++)
2219         {
2220                 if (quantContours[point]->getContourType(i) == type)
2221                 {
2222                         cont = quantContours[point]->getContour(i);
2223                 }
2224         }
2225
2226         return cont;
2227 }
2228
2229 // ----------------------------------------------------------------------------
2230 double marAxisCT::performXOR(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2231 {
2232         if (manual.size() == 0)
2233         {
2234                 return 0;
2235         }
2236
2237         int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2238         std::vector <marIsocontour* > polygons;
2239         vtkImageData* imagedata;
2240
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();
2248         imageA->Update();
2249
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();
2255         imageM->Update();
2256
2257
2258         for (int i = 0; i < quantContours[point]->getSize(); i++)
2259         {
2260                 if (quantContours[point]->getContourType(i) == type)
2261                 {
2262                         polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2263                 }
2264
2265         }
2266
2267         if (polygons.size() == 0)
2268         {
2269                 return 0;
2270         }
2271         int x,y;
2272         for (x = limInfX; x < (limSupX - 1); x++)
2273         {
2274                 for (y = limInfY; y < (limSupY - 1); y++)
2275                 {
2276                         unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2277                         unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2278                         *refValueA = 0;
2279                         *refValueM = 0;
2280
2281                         int numConts;
2282                         for (numConts = 0; numConts < polygons.size(); numConts++)
2283                         {
2284                                 if (pointInPolygon(polygons[numConts],x,y))
2285                                 {
2286                                                 *refValueA = 255;
2287                                 }
2288                         }
2289
2290                         for (numConts = 0; numConts < manual.size(); numConts++)
2291                         {
2292                                 if (pointInPolygon(manual[numConts],x,y))
2293                                 {
2294                                                 *refValueM = 255;
2295                                 }
2296                         }
2297                 }
2298         }
2299
2300 /*      vtkImageLogic* xor = vtkImageLogic::New();
2301         xor->SetInput1(imageM);
2302         xor->SetInput2(imageA);
2303         xor->SetOutputTrueValue(255);
2304         xor->SetOperationToXor();
2305
2306         vtkImageCast* cast = vtkImageCast::New();
2307         cast->SetInput(xor->GetOutput());
2308         cast->SetOutputScalarTypeToDouble();
2309 */
2310         int coincidencias = 0;
2311         for (x = limInfX; x < (limSupX - 1); x++)
2312         {
2313                 for (y = limInfY; y < (limSupY - 1); y++)
2314                 {
2315                         int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2316                         int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2317
2318                         int xor_ = compA^compM;
2319                         if (xor_ == 255)
2320                         {
2321                                 coincidencias++;
2322                         }
2323                 }
2324         }
2325
2326         double resp = (double)coincidencias ;
2327
2328         return resp;
2329
2330 }
2331
2332 // ----------------------------------------------------------------------------
2333 double  marAxisCT::performAND(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2334 {
2335
2336         if (manual.size() == 0)
2337         {
2338                 return 0;
2339         }
2340
2341         int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2342         std::vector <marIsocontour* > polygons;
2343         vtkImageData* imagedata;
2344
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();
2352         imageA->Update();
2353
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();
2359         imageM->Update();
2360
2361         for (int i = 0; i < quantContours[point]->getSize(); i++)
2362         {
2363                 if (quantContours[point]->getContourType(i) == type)
2364                 {
2365                         polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2366                 }
2367
2368         }
2369
2370         if (polygons.size() == 0)
2371         {
2372                 return 0;
2373         }
2374
2375         int x,y;
2376         for (x = limInfX; x < (limSupX - 1); x++)
2377         {
2378                 for (y = limInfY; y < (limSupY - 1); y++)
2379                 {
2380                         unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2381                         unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2382                         *refValueA = 0;
2383                         *refValueM = 0;
2384
2385                         int numConts;
2386                         for (numConts = 0; numConts < polygons.size(); numConts++)
2387                         {
2388                                 if (pointInPolygon(polygons[numConts],x,y))
2389                                 {
2390                                                 *refValueA = 255;
2391                                 }
2392                         }
2393
2394                         for (numConts = 0; numConts < manual.size(); numConts++)
2395                         {
2396                                 if (pointInPolygon(manual[numConts],x,y))
2397                                 {
2398                                                 *refValueM = 255;
2399                                 }
2400                         }
2401                 }
2402         }
2403
2404 /*      vtkImageLogic* and = vtkImageLogic::New();
2405         and->SetInput1(imageM);
2406         and->SetInput2(imageA);
2407         and->SetOutputTrueValue(255);
2408         and->SetOperationToAnd();
2409 */
2410         int coincidencias = 0;
2411         for (x = limInfX; x < (limSupX - 1); x++)
2412         {
2413                 for (y = limInfY; y < (limSupY - 1); y++)
2414                 {
2415                         int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2416                         int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2417
2418                         int and_ = compA & compM;
2419                         if (and_ == 255)
2420                         {
2421                                 coincidencias++;
2422                         }
2423                 }
2424         }
2425
2426         double resp = (double)coincidencias ;
2427
2428
2429         return resp;
2430 }
2431
2432 // ----------------------------------------------------------------------------
2433 marIsocontour*  marAxisCT::loadMarIsocontour(int size, double *vx, double *vy)
2434 {
2435         marIsocontour* cont = new marIsocontour();
2436
2437         for ( int j=0 ; j<size ; j++ )
2438         {
2439                 cont->insertPoint(vx[j], vy[j]);
2440         }
2441
2442         return cont;
2443 }
2444
2445 // ----------------------------------------------------------------------------
2446 double  marAxisCT::performUnion(int type, int point, std::vector<marIsocontour *> manual, kVolume* vol)
2447 {
2448                 if (manual.size() == 0)
2449         {
2450                 return 0;
2451         }
2452
2453         int limInfX,limInfY,limInfZ,limSupX,limSupY,limSupZ;
2454         std::vector <marIsocontour* > polygons;
2455         vtkImageData* imagedata;
2456
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();
2464         imageA->Update();
2465
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();
2471         imageM->Update();
2472
2473         int i;
2474         for (i = 0; i < quantContours[point]->getSize(); i++)
2475         {
2476                 if (quantContours[point]->getContourType(i) == type)
2477                 {
2478                         polygons.push_back(parsePolyDataToMarIsocontour(quantContours[point]->getContour(i)));
2479                 }
2480
2481         }
2482
2483         if (polygons.size() == 0)
2484         {
2485                 return 0;
2486         }
2487
2488         int x,y;
2489         for ( x = limInfX; x < (limSupX - 1); x++)
2490         {
2491                 for ( y = limInfY; y < (limSupY - 1); y++)
2492                 {
2493                         unsigned short *refValueA = (unsigned short *) imageA->GetScalarPointer(x,y,0);
2494                         unsigned short *refValueM = (unsigned short *) imageM->GetScalarPointer(x,y,0);
2495                         *refValueA = 0;
2496                         *refValueM = 0;
2497                         int  numConts;
2498                         for ( numConts = 0; numConts < polygons.size(); numConts++)
2499                         {
2500                                 if (pointInPolygon(polygons[numConts],x,y))
2501                                 {
2502                                                 *refValueA = 255;
2503                                 }
2504                         }
2505
2506                         for (numConts = 0; numConts < manual.size(); numConts++)
2507                         {
2508                                 if (pointInPolygon(manual[numConts],x,y))
2509                                 {
2510                                                 *refValueM = 255;
2511                                 }
2512                         }
2513                 }
2514         }
2515
2516 /*      vtkImageLogic* and = vtkImageLogic::New();
2517         and->SetInput1(imageM);
2518         and->SetInput2(imageA);
2519         and->SetOutputTrueValue(255);
2520         and->SetOperationToAnd();
2521 */
2522         int coincidencias = 0;
2523         for (x = limInfX; x < (limSupX - 1); x++)
2524         {
2525                 for ( y = limInfY; y < (limSupY - 1); y++)
2526                 {
2527                         int compA = (int) imageA->GetScalarComponentAsDouble(x,y,0,0);
2528                         int compM = (int) imageM->GetScalarComponentAsDouble(x,y,0,0);
2529
2530                         int or_ = compA | compM;
2531                         if (or_ == 255)
2532                         {
2533                                 coincidencias++;
2534                         }
2535                 }
2536         }
2537
2538         double resp = (double)coincidencias;
2539
2540
2541         return resp;
2542 }
2543