segmentation/clitkFillMaskGenericFilter.txx

   1 #ifndef clitkFillMaskGenericFilter_txx
   2 #define clitkFillMaskGenericFilter_txx
   3
   4 /* =================================================
   5  * @file   clitkFillMaskGenericFilter.txx
   6  * @author
   7  * @date
   8  *
   9  * @brief
  10  *
  11  ===================================================*/
  12
  13
  14 namespace clitk
  15 {
  16
  17   //-------------------------------------------------------------------
  18   // Update with the pixeltype
  19   //-------------------------------------------------------------------
  20   template <class  PixelType>
  21   void
  22   FillMaskGenericFilter::UpdateWithPixelType()
  23   {
  24     // Dim & Pix
  25     const unsigned int Dimension=3;
  26     typedef int InternalPixelType;
  27
  28     // ImageTypes
  29     typedef itk::Image<PixelType, Dimension> InputImageType;
  30     typedef itk::Image<InternalPixelType, Dimension> InternalImageType;
  31     typedef itk::Image<PixelType, Dimension> OutputImageType;
  32
  33     // Read the input
  34     typedef itk::ImageFileReader<InputImageType> InputReaderType;
  35     typename InputReaderType::Pointer reader = InputReaderType::New();
  36     reader->SetFileName( m_InputFileName);
  37     reader->Update();
  38     typename InputImageType::Pointer input= reader->GetOutput();
  39
  40     // Read the directions over which to fill holes
  41     std::vector<unsigned int> direction;
  42     if (m_ArgsInfo.dir_given)
  43       for ( unsigned int i=0;i<m_ArgsInfo.dir_given;i++)
  44         direction.push_back(m_ArgsInfo.dir_arg[i]);
  45     else
  46       for ( unsigned int i=0;i<Dimension;i++)
  47         direction.push_back(Dimension-i-1);
  48
  49     //----------------------------------------
  50     // Cast to internal type
  51     //----------------------------------------
  52     typedef itk::CastImageFilter<InputImageType,InternalImageType> InputCastImageFilterType;
  53     typename InputCastImageFilterType::Pointer inputCaster= InputCastImageFilterType::New();
  54     inputCaster->SetInput(input);
  55     inputCaster->Update();
  56
  57     //----------------------------------------
  58     // Loop over directions
  59     //----------------------------------------
  60     typename InternalImageType::Pointer output=inputCaster->GetOutput();
  61     for (unsigned int i=0; i<direction.size();i++)
  62       {
  63
  64         //----------------------------------------
  65         // Fill the holes of a mask in 2D
  66         //----------------------------------------
  67         if(m_Verbose) std::cout<<"Fill holes in the mask slice by slice in direction "<<direction[i]<<"..."<<std::endl;
  68
  69         // We define the region to be extracted.
  70         typedef  itk::Image<InternalPixelType, Dimension-1> ImageSliceType;
  71         typedef  itk::Image<InternalPixelType, Dimension-1> MaskSliceType;
  72         typename InternalImageType::RegionType region3D= input->GetLargestPossibleRegion();
  73         typename InternalImageType::RegionType::SizeType size3D= region3D.GetSize();
  74         typename InternalImageType::RegionType::SizeType size2D=size3D;
  75         size2D[direction[i]]=0;
  76         typename InternalImageType::IndexType start2D;
  77         start2D.Fill(0);
  78         typename InternalImageType::RegionType desiredRegion;
  79         desiredRegion.SetSize( size2D );
  80         desiredRegion.SetIndex( start2D );
  81
  82         // Extract and Join
  83         typedef itk::ExtractImageFilter<InternalImageType, ImageSliceType> ExtractImageFilterType;
  84         typedef itk::JoinSeriesImageFilter<ImageSliceType, InternalImageType> JoinSeriesFilterType;
  85         typename JoinSeriesFilterType::Pointer joinFilter=JoinSeriesFilterType::New();
  86         joinFilter->SetSpacing(input->GetSpacing()[direction[i]]);
  87
  88         //----------------------------------------
  89         // Run over the sliceIndexs
  90         // ----------------------------------------
  91         for(unsigned int sliceIndex=0; sliceIndex <size3D[direction[i]]; sliceIndex++)
  92           {
  93             //----------------------------------------
  94             // Extract mask sliceIndex
  95             //----------------------------------------
  96             typename ExtractImageFilterType::Pointer extractFilter=ExtractImageFilterType::New();
  97             extractFilter->SetInput(output);
  98             start2D[direction[i]]=sliceIndex;
  99             desiredRegion.SetIndex( start2D );
 100             extractFilter->SetExtractionRegion( desiredRegion );
 101             extractFilter->Update( );
 102             typename ImageSliceType::Pointer slice= extractFilter->GetOutput();
 103
 104             // Binarize the image (Before: OBJECT!=0, rest=0, After: object=1, rest=0 )
 105             typedef itk::BinaryThresholdImageFilter<ImageSliceType,ImageSliceType> BinarizeFilterType;
 106             typename BinarizeFilterType::Pointer binarizeFilter=BinarizeFilterType::New();
 107             binarizeFilter->SetInput(slice);
 108             binarizeFilter->SetUpperThreshold(0);
 109             binarizeFilter->SetOutsideValue(0);
 110             binarizeFilter->SetInsideValue(1);
 111             // writeImage<ImageSliceType>(binarizeFilter->GetOutput(),"/home/jef/tmp/input.mhd");
 112
 113             // Perform connected labelling on the slice (body+air=0 )
 114             typedef itk::ConnectedComponentImageFilter<ImageSliceType, ImageSliceType> ConnectFilterType;
 115             typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 116             connectFilter->SetInput(binarizeFilter->GetOutput());
 117             connectFilter->SetBackgroundValue(0);
 118             connectFilter->SetFullyConnected(false);
 119             //connectFilter->Update();
 120             //writeImage<ImageSliceType>(connectFilter->GetOutput(),"/home/jef/tmp/connect.mhd");
 121
 122             // Sort the labels
 123             typedef itk::RelabelComponentImageFilter<ImageSliceType, ImageSliceType> RelabelFilterType;
 124             typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 125             relabelFilter->SetInput(connectFilter->GetOutput());
 126             //relabelFilter->Update();
 127             //writeImage<ImageSliceType>(relabelFilter->GetOutput(),"/home/jef/tmp/label.mhd");
 128
 129             // Keep the first
 130             typedef itk::ThresholdImageFilter<ImageSliceType> ThresholdFilterType;
 131             typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
 132             thresholdFilter->SetInput(relabelFilter->GetOutput());
 133             thresholdFilter->SetUpper(1);
 134             thresholdFilter->SetOutsideValue(0);
 135             // thresholdFilter->Update();
 136             // writeImage<ImageSliceType>(thresholdFilter->GetOutput(),"/home/jef/tmp/bin.mhd");
 137
 138             // Invert the labels (lung 1, rest 0)
 139             typename BinarizeFilterType::Pointer switchFilter=BinarizeFilterType::New();
 140             switchFilter->SetInput(thresholdFilter->GetOutput());
 141             switchFilter->SetUpperThreshold(0);
 142             switchFilter->SetOutsideValue(0);
 143             switchFilter->SetInsideValue(1);
 144             switchFilter->Update();
 145             //writeImage<ImageSliceType>(switchFilter->GetOutput(),"/home/jef/tmp/inv_bin.mhd");
 146
 147             //Join
 148             joinFilter->SetInput( sliceIndex, switchFilter->GetOutput());
 149           }
 150
 151         // Join to a 3D image
 152         if (m_Verbose) std::cout<<"Joining the slices..."<<std::endl;
 153         joinFilter->Update();
 154
 155         // Permute the axes to reset to orientation
 156         typedef itk::PermuteAxesImageFilter<InternalImageType> PermuteFilterType;
 157         typename PermuteFilterType::Pointer permuteFilter=PermuteFilterType::New();
 158         permuteFilter->SetInput(joinFilter->GetOutput());
 159         typename PermuteFilterType::PermuteOrderArrayType order;
 160         order[direction[i]]=2;
 161         if( direction[i]==2)
 162           {
 163             order[0]=0;
 164             order[1]=1;
 165           }
 166         else if ( direction[i]==1)
 167           {
 168             order[0]=0;
 169             order[2]=1;
 170           }
 171         else if (direction[i]==0)
 172           {
 173             order[1]=0;
 174             order[2]=1;
 175           }
 176         permuteFilter->SetOrder(order);
 177         permuteFilter->Update();
 178         output =permuteFilter->GetOutput();
 179
 180         // Set the image direction to the input one
 181         output->SetDirection(input->GetDirection());
 182         output->SetOrigin(input->GetOrigin());
 183       }
 184
 185
 186     // Cast
 187     typedef itk::CastImageFilter<InternalImageType,OutputImageType> OutputCastImageFilterType;
 188     typename OutputCastImageFilterType::Pointer outputCaster =OutputCastImageFilterType::New();
 189     outputCaster->SetInput(output);
 190
 191     // Output
 192     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 193     typename WriterType::Pointer writer = WriterType::New();
 194     writer->SetFileName(m_ArgsInfo.output_arg);
 195     writer->SetInput(outputCaster->GetOutput());
 196     writer->Update();
 197   }
 198
 199 }//end clitk
 200
 201 #endif //#define clitkFillMaskGenericFilter_txx