segmentation/clitkTestFilter.cxx

   1
   2 /*=========================================================================
   3   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   4
   5   Authors belong to:
   6   - University of LYON              http://www.universite-lyon.fr/
   7   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
   8   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   9
  10   This software is distributed WITHOUT ANY WARRANTY; without even
  11   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12   PURPOSE.  See the copyright notices for more information.
  13
  14   It is distributed under dual licence
  15
  16   - BSD        See included LICENSE.txt file
  17   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  18   ======================================================================-====*/
  19
  20 // clitk
  21 #include "clitkTestFilter_ggo.h"
  22 #include "clitkImageCommon.h"
  23 #include "clitkBooleanOperatorLabelImageFilter.h"
  24 #include "clitkAutoCropFilter.h"
  25 //#include "clitkRelativePositionConstraintLabelImageFilter.h"
  26 #include "clitkResampleImageWithOptionsFilter.h"
  27 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
  28 #include "clitkExtractLungFilter.h"
  29 #include "clitkExtractPatientFilter.h"
  30 #include "clitkExtractMediastinumFilter.h"
  31 //#include "clitkTestStation7.h"
  32 #include "clitkSegmentationUtils.h"
  33 #include "clitkMorphoMathFilter.h"
  34
  35 // ITK ENST
  36 #include "RelativePositionPropImageFilter.h"
  37
  38 // VTK
  39 #include <vtkSmartPointer.h>
  40 #include <vtkPolyData.h>
  41 #include <vtkCellArray.h>
  42 #include <vtkPolyDataWriter.h>
  43 #include <vtkImageData.h>
  44 #include <vtkMetaImageWriter.h>
  45 #include <vtkPolyDataToImageStencil.h>
  46 #include <vtkLinearExtrusionFilter.h>
  47 #include <vtkImageStencil.h>
  48
  49 #include <algorithm>
  50
  51
  52 template<class PointType>
  53 class comparePointsX {
  54 public:
  55   bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
  56 };
  57
  58 template<int Dim>
  59 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
  60   std::vector<itk::Point<double, Dim-1> > previous;
  61   HypercubeCorners<Dim-1>(previous);
  62   out.resize(previous.size()*2);
  63   for(uint i=0; i<out.size(); i++) {
  64     itk::Point<double, Dim> p;
  65     if (i<previous.size()) p[0] = 0;
  66     else p[0] = 1;
  67     for(int j=0; j<Dim; j++) p[j+1] = previous[i%previous.size()][j]; //NON i p
  68     out[i] = p;
  69   }
  70 }
  71
  72 template<>
  73 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
  74   out.resize(2);
  75   out[0][0] = 0;
  76   out[1][0] = 1;
  77 }
  78
  79
  80 //--------------------------------------------------------------------
  81 int main(int argc, char * argv[]) {
  82
  83   // Init command line
  84   GGO(clitkTestFilter, args_info);
  85
  86   // Image types
  87   //typedef unsigned char PixelType;
  88   typedef short PixelType;
  89   static const int Dim=3;
  90   typedef itk::Image<PixelType, Dim> InputImageType;
  91
  92   // Read inputs
  93   InputImageType::Pointer input1;
  94   InputImageType::Pointer input2;
  95   InputImageType::Pointer input3;
  96   if (args_info.input1_given) input1 = clitk::readImage<InputImageType>(args_info.input1_arg, true);
  97   if (args_info.input2_given) input2 = clitk::readImage<InputImageType>(args_info.input2_arg, true);
  98   if (args_info.input3_given) input3 = clitk::readImage<InputImageType>(args_info.input3_arg, true);
  99
 100   // Declare output
 101   InputImageType::Pointer output;
 102
 103   //--------------------------------------------------------------------
 104   // Filter test BooleanOperatorLabelImageFilter
 105   if (0) {
 106     typedef clitk::BooleanOperatorLabelImageFilter<InputImageType> FilterType;
 107     FilterType::Pointer filter = FilterType::New();
 108     filter->SetInput1(input1);
 109     filter->SetInput2(input2);
 110     output = clitk::NewImageLike<InputImageType>(input1);
 111     filter->GraftOutput(output);  /// TO VERIFY !!!
 112     filter->Update();
 113     filter->SetInput2(input3);
 114     filter->Update();
 115     output = filter->GetOutput();
 116     clitk::writeImage<InputImageType>(output, args_info.output_arg);
 117   }
 118
 119   //--------------------------------------------------------------------
 120   // Filter test AutoCropLabelImageFilter
 121   if (0) {
 122     typedef clitk::AutoCropFilter<InputImageType> FilterType;
 123     FilterType::Pointer filter = FilterType::New();
 124     filter->SetInput(input1);
 125     filter->Update();
 126     output = filter->GetOutput();
 127     clitk::writeImage<InputImageType>(output, args_info.output_arg);
 128   }
 129
 130   //--------------------------------------------------------------------
 131   // Filter test RelativePositionPropImageFilter
 132   if (0) {
 133     typedef itk::Image<float, Dim> OutputImageType;
 134     OutputImageType::Pointer outputf;
 135     typedef itk::RelativePositionPropImageFilter<InputImageType, OutputImageType> FilterType;
 136     FilterType::Pointer filter = FilterType::New();
 137     filter->SetInput(input1);
 138
 139     filter->SetAlpha1(clitk::deg2rad(args_info.angle1_arg)); // xy plane
 140     filter->SetAlpha2(clitk::deg2rad(args_info.angle2_arg));
 141     filter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
 142     filter->SetFast(true);
 143     filter->SetRadius(2);
 144     filter->SetVerboseProgress(true);
 145
 146     /*         A1   A2
 147                Left      0    0
 148                Right   180    0
 149                Ant      90    0
 150                Post    -90    0
 151                Inf       0   90
 152                Sup       0  -90
 153     */
 154
 155     filter->Update();
 156     clitk::writeImage<OutputImageType>(filter->GetOutput(), args_info.output_arg);
 157   }
 158
 159   //--------------------------------------------------------------------
 160   // Filter test
 161   if (0) {
 162     typedef itk::Image<short, Dim> OutputImageType;
 163     typedef clitk::ResampleImageWithOptionsFilter<InputImageType, OutputImageType> FilterType;
 164     FilterType::Pointer filter = FilterType::New();
 165     filter->SetInput(input1);
 166     filter->SetOutputIsoSpacing(1);
 167     //filter->SetNumberOfThreads(4); // auto
 168     filter->SetGaussianFilteringEnabled(false);
 169     filter->Update();
 170     clitk::writeImage<OutputImageType>(filter->GetOutput(), args_info.output_arg);
 171   }
 172
 173   //--------------------------------------------------------------------
 174   // Filter AddRelativePositionConstraintToLabelImageFilter test
 175   if (0) {
 176     /*
 177     typedef clitk::AddRelativePositionConstraintToLabelImageFilter<InputImageType> FilterType;
 178     FilterType::Pointer filter = FilterType::New();
 179     filter->SetInput(input1);
 180     filter->SetInputObject(input2);
 181     filter->SetOrientationType(FilterType::LeftTo);
 182     filter->SetIntermediateSpacing(5);
 183     filter->SetFuzzyThreshold(0.5);
 184     filter->VerboseStepOn();
 185     filter->WriteStepOff();
 186     filter->Update();
 187
 188     filter->SetInput(filter->GetOutput());
 189     filter->SetInputObject(input2);
 190     filter->SetOrientationType(FilterType::RightTo);
 191     filter->SetIntermediateSpacing(5);
 192     filter->SetFuzzyThreshold(0.5);
 193     filter->Update();
 194
 195     clitk::writeImage<InputImageType>(filter->GetOutput(), args_info.output_arg);
 196     */
 197   }
 198
 199   //--------------------------------------------------------------------
 200   // Filter test ExtractPatientFilter
 201   if (0) {
 202     /*
 203     typedef itk::Image<char, Dim> OutputImageType;
 204     typedef clitk::ExtractPatientFilter<InputImageType, OutputImageType> FilterType;
 205     FilterType::Pointer filter = FilterType::New();
 206     filter->SetInput(input1);
 207     filter->VerboseStepOn();
 208     filter->WriteStepOn();
 209     // options (default)
 210     filter->SetUpperThreshold(-300);
 211     filter->FinalOpenCloseOff(); // default=on (not rezally needed ?)
 212     filter->Update();
 213     OutputImageType::Pointer output = filter->GetOutput();
 214     clitk::writeImage<OutputImageType>(output, args_info.output_arg);
 215     */
 216   }
 217
 218   //--------------------------------------------------------------------
 219   // Filter test ExtractLungsFilter
 220   if (0) {
 221     /*
 222     typedef itk::Image<PixelType, Dim> OutputImageType; // to change into char
 223     typedef clitk::ExtractLungFilter<InputImageType, OutputImageType> FilterType;
 224     FilterType::Pointer filter = FilterType::New();
 225     // DD(filter->GetNumberOfSteps());
 226     filter->SetInput(input1); // CT
 227     filter->SetInputPatientMask(input2, 0); // Patient mask and BG value
 228     filter->VerboseStepOn();
 229     filter->WriteStepOn();
 230     // options (default)
 231     //filter->SetBackgroundValue(0);
 232     filter->SetUpperThreshold(-300);
 233     // filter->SetMinimumComponentSize(100);
 234
 235     filter->Update();
 236     OutputImageType::Pointer output = filter->GetOutput();
 237     clitk::writeImage<OutputImageType>(output, args_info.output_arg);
 238     */
 239   }
 240
 241   //--------------------------------------------------------------------
 242   // Filter test ExtractMediastinumFilter
 243   if (0) {
 244     /*
 245     typedef clitk::ExtractMediastinumFilter<InputImageType> FilterType;
 246     FilterType::Pointer filter = FilterType::New();
 247     filter->SetInputPatientLabelImage(input1);
 248     filter->SetInputLungLabelImage(input2, 0, 1, 2); // BG, FG Left Lung, FG Right Lung
 249     filter->SetInputBonesLabelImage(input3);
 250     filter->VerboseStepOn();
 251     filter->WriteStepOn();
 252     filter->Update();
 253     output = filter->GetOutput();
 254     clitk::writeImage<InputImageType>(output, args_info.output_arg);
 255     */
 256   }
 257
 258   //--------------------------------------------------------------------
 259   // Test for auto register sub-task in a segmentation process
 260   if (0) {
 261     // ExtractLymphStation_7 * s7 = new ExtractLymphStation_7;
 262     //    s7->SetArgsInfo<args_info_clitkTestFilter>(args_info);
 263     // GetParent->SetArgsInfo<>
 264     //s7->StartSegmentation();
 265   }
 266
 267   //--------------------------------------------------------------------
 268   // Test for biinary image from a contour set
 269   if (0) {
 270     DD("here");
 271
 272     // Type of a slice
 273     typedef itk::Image<InputImageType::PixelType, InputImageType::ImageDimension-1> SliceType;
 274
 275     // Build the list of slices
 276     std::vector<SliceType::Pointer> slices;
 277     clitk::ExtractSlices<InputImageType>(input1, 2, slices);
 278     DD(slices.size());
 279
 280     // HOW TO DO SEVERAL BY SLICE !!! not a map ?
 281
 282     // Compute centroids 3D centroids by slices
 283     int BG = 0;
 284     std::map<int, std::vector<InputImageType::PointType> > centroids3D;
 285     for(uint i=0; i<slices.size(); i++) {
 286       // Labelize
 287       slices[i] = clitk::Labelize<SliceType>(slices[i], BG, true, 1);
 288       // ComputeCentroids
 289       std::vector<SliceType::PointType> temp;
 290       clitk::ComputeCentroids<SliceType>(slices[i], BG, temp);
 291       for(uint j=1; j<temp.size(); j++) {
 292         InputImageType::PointType a;
 293         clitk::PointsUtils<InputImageType>::Convert2DTo3D(temp[j], input1, i, a);
 294         centroids3D[i].push_back(a);
 295       }
 296     }
 297
 298     // REPRENDRE POUR TOUT FAIRE BY SLICE (pas de i);
 299
 300     // Take a given slice i=29
 301     int index=29;
 302     SliceType::Pointer slice = slices[index];
 303     std::vector<InputImageType::PointType> points = centroids3D[index];
 304
 305     // Sort points from R to L
 306     std::sort(points.begin(), points.end(),
 307               comparePointsX<InputImageType::PointType>());
 308
 309     // Slice corners (quel sens ?)
 310
 311     // Compute min and max coordinates
 312     const uint dim=3;
 313     typedef InputImageType::PointType PointType;
 314     typedef InputImageType::IndexType IndexType;
 315     PointType min_c, max_c;
 316     IndexType min_i, max_i;
 317     min_i = input1->GetLargestPossibleRegion().GetIndex();
 318     for(uint i=0; i<dim; i++) max_i[i] = input1->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
 319     input1->TransformIndexToPhysicalPoint(min_i, min_c);
 320     input1->TransformIndexToPhysicalPoint(max_i, max_c);
 321
 322     // Compute the corners coordinates
 323     std::vector<itk::Point<double, dim> > l;
 324     HypercubeCorners<dim>(l);
 325     for(uint i=0; i<l.size(); i++) {
 326       for(uint j=0; j<dim; j++) {
 327         if (l[i][j] == 0) l[i][j] = min_c[j];
 328         if (l[i][j] == 1) l[i][j] = max_c[j];
 329       }
 330     }
 331     DDV(l, 8);
 332
 333
 334     // Add first/last point, horizontally to the image boundaries
 335
 336     double sz = points[0][2]; // slice Z
 337     double margin = 2; // needed (if not forget to remove the first line)
 338     // Corners 1
 339     InputImageType::PointType p = min_c;
 340     p[0] -= margin; // margins
 341     p[1] -= margin; // margins
 342     p[2] = sz;
 343     points.insert(points.begin(), p);
 344     // vertical p
 345     p = points[0];
 346     p[0] = min_c[0] - margin; //margin
 347     points.insert(points.begin(), p);
 348     // last H point
 349     p = points.back();
 350     p[0] = max_c[0];
 351     points.push_back(p);
 352     // last corners
 353     p[0] = max_c[0];
 354     p[1] = min_c[1]-margin;
 355     p[2] = sz;
 356     points.push_back(p);
 357     // Same first point
 358     p = points[0];
 359     points.push_back(p);
 360
 361     DDV(points, points.size());
 362
 363     // create a contour, polydata.
 364     vtkSmartPointer<vtkPolyData> mesh = vtkSmartPointer<vtkPolyData>::New();
 365     mesh->Allocate(); //for cell structures
 366     mesh->SetPoints(vtkPoints::New());
 367     vtkIdType ids[2];
 368     int point_number = points.size();
 369     for (unsigned int i=0; i<points.size(); i++) {
 370       mesh->GetPoints()->InsertNextPoint(points[i][0],points[i][1],points[i][2]);
 371       ids[0]=i;
 372       ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
 373       mesh->GetLines()->InsertNextCell(2,ids);
 374     }
 375
 376     vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
 377     w->SetInput(mesh);
 378     w->SetFileName("bidon.vtk");
 379     w->Write();
 380
 381     // binarize
 382     DD("binarize");
 383     double *bounds=mesh->GetBounds();
 384     DDV(bounds, 6);
 385
 386     DD("create image");
 387     vtkSmartPointer<vtkImageData> binary_image=vtkSmartPointer<vtkImageData>::New();
 388     binary_image->SetScalarTypeToUnsignedChar();
 389     ///Use the smallest mask in which the mesh fits
 390     // Add two voxels on each side to make sure the mesh fits
 391     //    double * samp_origin=
 392     InputImageType::PointType samp_origin = input1->GetOrigin();
 393     //    double *
 394     InputImageType::SpacingType spacing=input1->GetSpacing();
 395     double * spacing2 = new double[3];
 396     spacing2[0] = spacing[0];
 397     spacing2[1] = spacing[1];
 398     spacing2[2] = spacing[2];
 399     DD(spacing2[0]);
 400     binary_image->SetSpacing(spacing2);
 401     /// Put the origin on a voxel to avoid small skips
 402     DD(floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]);
 403     DD(bounds[0]);
 404     DD(samp_origin[0]);
 405     DD(spacing[0]);
 406     binary_image->SetOrigin(//samp_origin[0], samp_origin[1], samp_origin[2]);
 407                             floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]+samp_origin[0],
 408                             floor((bounds[2]-samp_origin[1])/spacing[1]-2)*spacing[1]+samp_origin[1],
 409                             floor((bounds[4]-samp_origin[2])/spacing[2]-2)*spacing[2]+samp_origin[2]);
 410     double * origin=binary_image->GetOrigin();
 411     binary_image->SetExtent(0,ceil((bounds[1]-origin[0])/spacing[0]+4), // Joel used +4 here (?)
 412                             0,ceil((bounds[3]-origin[1])/spacing[1]+4),
 413                             0,ceil((bounds[5]-origin[2])/spacing[2]+4));
 414     binary_image->AllocateScalars();
 415     memset(binary_image->GetScalarPointer(),0,
 416            binary_image->GetDimensions()[0]*binary_image->GetDimensions()[1]*binary_image->GetDimensions()[2]*sizeof(unsigned char));
 417
 418     vtkSmartPointer<vtkPolyDataToImageStencil> sts=vtkSmartPointer<vtkPolyDataToImageStencil>::New();
 419     //The following line is extremely important
 420     //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933
 421     sts->SetTolerance(0);
 422     sts->SetInformationInput(binary_image);
 423
 424     bool extrude=true;
 425     if (extrude) {
 426       vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
 427       extrude->SetInput(mesh);
 428
 429       /// NO ????We extrude in the -slice_spacing direction to respect the FOCAL convention ???
 430
 431       extrude->SetVector(0, 0, input1->GetSpacing()[2]);//slice_spacing); // 2* ? yes use a single one
 432       sts->SetInput(extrude->GetOutput());
 433     } else
 434       sts->SetInput(mesh);
 435
 436     DD("stencil");
 437     vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
 438     stencil->SetStencil(sts->GetOutput());
 439     stencil->SetInput(binary_image);
 440     stencil->Update();
 441
 442     DD("write");
 443     vtkSmartPointer<vtkMetaImageWriter> ww = vtkSmartPointer<vtkMetaImageWriter>::New();
 444     ww->SetInput(stencil->GetOutput());
 445     ww->SetFileName("binary.mhd");
 446     ww->Write();
 447   }
 448   //--------------------------------------------------------------------
 449
 450   //--------------------------------------------------------------------
 451   // Test for vessels ReconstructionByDilatation
 452   if (1) {
 453     // Read input CT (already croped)
 454     // treshold 3D
 455     // erode n times (or in 2D ?)
 456     // slices extract
 457     // SBS
 458     //    - CCL (keep mask)
 459     //    - for all CCL, ReconstructionByDilatation in initial mask
 460     // joint for output
 461
 462     // input1
 463     DD("binarize")
 464     int BG = 0;
 465     int FG = 1;
 466     typedef itk::Image<unsigned char, InputImageType::ImageDimension> MaskImageType;
 467     typedef itk::BinaryThresholdImageFilter<InputImageType, MaskImageType> BinarizeFilterType;
 468     BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
 469     binarizeFilter->SetInput(input1);
 470     binarizeFilter->SetLowerThreshold(150);
 471     binarizeFilter->SetInsideValue(FG);
 472     binarizeFilter->SetOutsideValue(BG);
 473     binarizeFilter->Update();
 474     MaskImageType::Pointer mask = binarizeFilter->GetOutput();
 475     clitk::writeImage<MaskImageType>(mask, "m.mhd");
 476
 477     // Extract slices
 478     typedef itk::Image<MaskImageType::PixelType, MaskImageType::ImageDimension-1> SliceType;
 479     std::vector<SliceType> slices_mask;
 480     clitk::ExtractSlices<MaskImageType>(mask, 2, slices_mask);
 481     DD(slices_mask.size());
 482
 483     std::vector<SliceType> debug_eroded;
 484     std::vector<SliceType> debug_labeled;
 485
 486     // Loop
 487     for(uint i=0; i<slices_mask.size(); i++) {
 488       DD(i);
 489       // erode
 490       DD("erosion");
 491       clitk::MorphoMathFilter<SliceType>::Pointer f = clitk::MorphoMathFilter<SliceType>::New();
 492       f->SetInput(slices_mask[i]);
 493       f->SetBackgroundValue(BG);
 494       f->SetForegroundValue(FG);
 495       f->SetRadius(1);
 496       f->SetOperationType(0); // Erode
 497       f->Update();
 498       SliceType::Pointer eroded = f->GetOutput();
 499       debug_eroded.push_back(eroded);
 500
 501       // CCL
 502       DD("CCL");
 503       SliceType::Pointer labeled = Labelize<MaskSliceType>(slices_mask[i], 0, false, 10);
 504       debug_labeled.push_back(labeled);
 505
 506       // ReconstructionByDilatation
 507
 508
 509     } // end loop
 510
 511     MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, mask, 2);
 512     clitk::writeImage<MaskImageType>(eroded, "eroded.mhd");
 513
 514     MaskImageType::Pointer labeled = clitk::JoinSlices<MaskImageType>(debug_labeled, mask, 2);
 515     clitk::writeImage<MaskImageType>(labeled, "labeled.mhd");
 516
 517   }
 518   //--------------------------------------------------------------------
 519
 520   // This is the end my friend
 521   return EXIT_SUCCESS;
 522 }// end main
 523 //--------------------------------------------------------------------