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

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