1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
20 #include "clitkCommon.h"
21 #include "clitkBooleanOperatorLabelImageFilter.h"
22 #include "clitkAutoCropFilter.h"
23 #include "clitkResampleImageWithOptionsFilter.h"
24 #include "clitkBooleanOperatorLabelImageFilter.h"
28 #include "itkStatisticsLabelMapFilter.h"
29 #include "itkLabelImageToStatisticsLabelMapFilter.h"
30 #include "itkRegionOfInterestImageFilter.h"
31 #include "itkBinaryThresholdImageFilter.h"
32 #include "itkBinaryErodeImageFilter.h"
33 #include "itkBinaryBallStructuringElement.h"
36 #include "RelativePositionPropImageFilter.h"
38 //--------------------------------------------------------------------
39 template <class TImageType>
40 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
41 AddRelativePositionConstraintToLabelImageFilter():
43 itk::ImageToImageFilter<TImageType, TImageType>()
45 this->SetNumberOfRequiredInputs(2);
46 SetFuzzyThreshold(0.6);
47 SetBackgroundValue(0);
48 SetObjectBackgroundValue(0);
49 SetOrientationType(LeftTo);
50 ResampleBeforeRelativePositionFilterOn();
51 SetIntermediateSpacing(10);
53 // Step 1 : resample (option=sampling)
55 // Step 3 : relative position (option = angle)
57 // Step 5 : Erode for boundary
58 // Step 6 : resample to initial spacing
59 // Step 7 : pad if not the same size : it can occur when downsample and upsample
60 // Step 6: combine input+thresholded relpos
64 ResampleBeforeRelativePositionFilterOn();
65 SetIntermediateSpacing(10);
66 SetOrientationType(LeftTo);
67 // SegmentationStep * step = new SegmentationStep;
68 // step->SetName("Initial resampling and relative position map");
69 // SetStep(step, &Self::ResampleAndRelativePositionMap);
71 // Step 3 : threshold + postprocess
72 // SetFuzzyThreshold(0.4);
73 // step = new SegmentationStep;
74 // step->SetName("Map threshold and post-processing");
75 // SetStep(step, &Self::MapThresholdAndPostProcessing);
77 //--------------------------------------------------------------------
80 //--------------------------------------------------------------------
81 template <class TImageType>
83 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
84 SetInput(const ImageType * image) {
85 // Process object is not const-correct so the const casting is required.
86 this->SetNthInput(0, const_cast<ImageType *>(image));
88 //--------------------------------------------------------------------
91 //--------------------------------------------------------------------
92 template <class TImageType>
94 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
95 SetInputObject(const ImageType * image) {
96 // Process object is not const-correct so the const casting is required.
97 this->SetNthInput(1, const_cast<ImageType *>(image));
99 //--------------------------------------------------------------------
102 //--------------------------------------------------------------------
103 template <class TImageType>
105 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
106 GenerateOutputInformation() {
107 ImagePointer input = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
108 ImagePointer outputImage = this->GetOutput(0);
109 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
111 //--------------------------------------------------------------------
114 //--------------------------------------------------------------------
115 template <class TImageType>
117 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
118 GenerateInputRequestedRegion() {
120 itk::ImageToImageFilter<TImageType, TImageType>::GenerateInputRequestedRegion();
121 // Get input pointers and set requested region to common region
122 ImagePointer input1 = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
123 ImagePointer input2 = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
124 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
125 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
127 //--------------------------------------------------------------------
130 //--------------------------------------------------------------------
131 template <class TImageType>
133 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
134 SetAngle1(double a) {
135 SetOrientationType(Angle);
138 //--------------------------------------------------------------------
141 //--------------------------------------------------------------------
142 template <class TImageType>
144 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
145 SetAngle2(double a) {
146 SetOrientationType(Angle);
149 //--------------------------------------------------------------------
152 //--------------------------------------------------------------------
153 template <class TImageType>
155 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
156 SetOrientationType(OrientationTypeEnumeration orientation) {
157 m_OrientationType = orientation;
158 switch (m_OrientationType) {
159 case LeftTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(0); break;
160 case RightTo: m_Angle1 = clitk::deg2rad(180); m_Angle2 = clitk::deg2rad(0); break;
161 case AntTo: m_Angle1 = clitk::deg2rad(90); m_Angle2 = clitk::deg2rad(0); break;
162 case PostTo: m_Angle1 = clitk::deg2rad(-90); m_Angle2 = clitk::deg2rad(0); break;
163 case InfTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(90); break;
164 case SupTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(-90); break;
176 //--------------------------------------------------------------------
179 //--------------------------------------------------------------------
180 template <class TImageType>
182 clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType>::
185 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
186 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
188 //--------------------------------------------------------------------
189 //--------------------------------------------------------------------
190 StartNewStep("Resample And Relative Position Map");
192 static const unsigned int dim = ImageType::ImageDimension;
194 if (m_ResampleBeforeRelativePositionFilter) {
195 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
196 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
197 resampleFilter->SetInput(object);
198 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
199 resampleFilter->SetGaussianFilteringEnabled(false);
200 // resampleFilter->SetVerboseOptions(true);
201 resampleFilter->Update();
202 working_image = resampleFilter->GetOutput();
205 working_image = object;
207 // writeImage<ImageType>(working_image, "res.mhd");
209 // Step 2: object pad to input image -> we want to compute the
210 // relative position for each point belonging to the input image
211 // domain, so we have to extend (pad) the object image to fit the
213 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
214 typename ImageType::Pointer output = ImageType::New();
216 for(unsigned int i=0; i<dim; i++) {
217 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
220 region.SetSize(size);
221 // output->SetLargestPossibleRegion(region);
222 output->SetRegions(region);
223 output->SetSpacing(working_image->GetSpacing());
224 output->SetOrigin(input->GetOrigin());
226 output->FillBuffer(m_BackgroundValue);
227 typename PadFilterType::Pointer padFilter = PadFilterType::New();
228 typename PadFilterType::InputImageIndexType index;
229 for(unsigned int i=0; i<dim; i++) {
230 index[i] = lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/(double)m_IntermediateSpacing);
232 padFilter->SetSourceImage(working_image);
233 padFilter->SetDestinationImage(output);
234 padFilter->SetDestinationIndex(index);
235 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
237 working_image = padFilter->GetOutput();
240 DD("[debug] RelPos : same size and spacing : no padding");
243 // Keep object image (with resampline and pad)
244 object_resampled = working_image;
245 // writeImage<ImageType>(working_image, "pad.mhd");
247 // Step 3: compute rel pos in object
248 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
249 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
250 relPosFilter->SetInput(working_image);
251 relPosFilter->SetAlpha1(m_Angle1); // xy plane
252 relPosFilter->SetAlpha2(m_Angle2);
253 relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
254 relPosFilter->SetFast(true);
255 relPosFilter->SetRadius(1); // seems sufficient in this cas
256 // relPosFilter->SetVerboseProgress(true);
257 relPosFilter->Update();
258 relPos = relPosFilter->GetOutput();
259 this->template StopCurrentStep<FloatImageType>(relPos);
261 //--------------------------------------------------------------------
262 //--------------------------------------------------------------------
263 StartNewStep("Map Threshold And Post Processing");
266 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
267 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
268 thresholdFilter->SetInput(relPos);
269 thresholdFilter->SetOutsideValue(m_BackgroundValue);
270 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
271 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
272 thresholdFilter->Update();
273 working_image = thresholdFilter->GetOutput();
275 // Step 2 : erosion with initial mask to exclude pixels that were
276 // inside the resampled version and outside the original mask
277 typedef itk::BinaryBallStructuringElement<unsigned int, 3> StructuringElementType;
278 StructuringElementType kernel;
280 kernel.CreateStructuringElement();
281 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
282 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
283 erodeFilter->SetInput(working_image);
284 erodeFilter->SetKernel(kernel);
285 erodeFilter->SetBackgroundValue(m_BackgroundValue);
286 erodeFilter->SetErodeValue(m_BackgroundValue+1);
287 erodeFilter->Update();
288 working_image = erodeFilter->GetOutput();
290 // Step 5: resample to initial spacing
291 if (m_ResampleBeforeRelativePositionFilter) {
292 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
293 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
294 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
295 resampleFilter->SetInput(working_image);
296 resampleFilter->SetOutputSpacing(input->GetSpacing());
297 resampleFilter->SetGaussianFilteringEnabled(false);
298 // resampleFilter->SetVerboseOptions(true);
299 resampleFilter->Update();
300 working_image = resampleFilter->GetOutput();
303 // writeImage<ImageType>(working_image, "resinitial.mhd");
305 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
306 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
307 typename ImageType::Pointer temp = ImageType::New();
308 temp->CopyInformation(input);
309 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
311 temp->FillBuffer(m_BackgroundValue);
312 typename PadFilterType::Pointer padFilter2 = PadFilterType::New(); // if yes : redo relpos
313 padFilter2->SetSourceImage(working_image);
314 padFilter2->SetDestinationImage(temp);
315 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
316 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
317 padFilter2->Update();
318 working_image = padFilter2->GetOutput();
321 DD("[debug] Rel Pos : no padding after");
323 // writeImage<ImageType>(working_image, "pad2.mhd");
325 // Step 6: combine input+thresholded relpos
326 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
327 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
328 combineFilter->SetBackgroundValue(m_BackgroundValue);
329 combineFilter->SetBackgroundValue1(m_BackgroundValue);
330 combineFilter->SetBackgroundValue2(m_BackgroundValue);
331 combineFilter->SetForegroundValue(m_BackgroundValue+1);
332 combineFilter->SetInput1(input);
333 combineFilter->SetInput2(working_image);
334 combineFilter->SetOperationType(BoolFilterType::And);
335 combineFilter->InPlaceOn();
336 combineFilter->Update();
337 working_image = combineFilter->GetOutput();
339 combineFilter = BoolFilterType::New();
340 combineFilter->SetInput1(working_image);
341 combineFilter->SetInput2(object);
342 combineFilter->SetOperationType(BoolFilterType::AndNot);
343 combineFilter->InPlaceOn();
344 combineFilter->Update();
346 working_image = combineFilter->GetOutput();
347 // writeImage<ImageType>(working_image, "combine.mhd");
350 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
351 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
352 cropFilter->SetInput(working_image);
353 cropFilter->ReleaseDataFlagOff();
354 cropFilter->Update();
355 working_image = cropFilter->GetOutput();
356 this->template StopCurrentStep<ImageType>(working_image);
358 //--------------------------------------------------------------------
359 //--------------------------------------------------------------------
361 // Final Step -> set output
362 this->SetNthOutput(0, working_image);
365 //--------------------------------------------------------------------