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 ImageType>
40 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
41 AddRelativePositionConstraintToLabelImageFilter():
43 itk::ImageToImageFilter<ImageType, ImageType>()
45 this->SetNumberOfRequiredInputs(2);
46 SetFuzzyThreshold(0.6);
47 SetBackgroundValue(0);
48 SetObjectBackgroundValue(0);
49 SetOrientationType(LeftTo);
50 ResampleBeforeRelativePositionFilterOn();
51 SetIntermediateSpacing(10);
54 //--------------------------------------------------------------------
57 //--------------------------------------------------------------------
58 template <class ImageType>
60 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
61 SetInput(const ImageType * image)
63 // Process object is not const-correct so the const casting is required.
64 this->SetNthInput(0, const_cast<ImageType *>(image));
66 //--------------------------------------------------------------------
69 //--------------------------------------------------------------------
70 template <class ImageType>
72 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
73 SetInputObject(const ImageType * image)
75 // Process object is not const-correct so the const casting is required.
76 this->SetNthInput(1, const_cast<ImageType *>(image));
78 //--------------------------------------------------------------------
81 //--------------------------------------------------------------------
82 template <class ImageType>
84 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
85 GenerateOutputInformation()
87 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
88 ImagePointer outputImage = this->GetOutput(0);
89 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
91 //--------------------------------------------------------------------
94 //--------------------------------------------------------------------
95 template <class ImageType>
97 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
98 GenerateInputRequestedRegion()
101 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
102 // Get input pointers and set requested region to common region
103 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
104 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
105 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
106 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
108 //--------------------------------------------------------------------
111 //--------------------------------------------------------------------
112 template <class ImageType>
114 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
117 SetOrientationType(Angle);
120 //--------------------------------------------------------------------
123 //--------------------------------------------------------------------
124 template <class ImageType>
126 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
129 SetOrientationType(Angle);
132 //--------------------------------------------------------------------
135 //--------------------------------------------------------------------
136 template <class ImageType>
138 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
139 SetOrientationType(OrientationTypeEnumeration orientation)
141 m_OrientationType = orientation;
142 switch (m_OrientationType) {
143 case LeftTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(0); break;
144 case RightTo: m_Angle1 = clitk::deg2rad(180); m_Angle2 = clitk::deg2rad(0); break;
145 case AntTo: m_Angle1 = clitk::deg2rad(90); m_Angle2 = clitk::deg2rad(0); break;
146 case PostTo: m_Angle1 = clitk::deg2rad(-90); m_Angle2 = clitk::deg2rad(0); break;
147 case InfTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(90); break;
148 case SupTo: m_Angle1 = clitk::deg2rad(0); m_Angle2 = clitk::deg2rad(-90); break;
160 //--------------------------------------------------------------------
163 //--------------------------------------------------------------------
164 template <class ImageType>
166 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
170 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
171 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
173 //--------------------------------------------------------------------
174 //--------------------------------------------------------------------
175 static const unsigned int dim = ImageType::ImageDimension;
176 StartNewStep("Initial resample");
178 if (m_ResampleBeforeRelativePositionFilter) {
179 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
180 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
181 resampleFilter->SetInput(object);
182 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
183 resampleFilter->SetGaussianFilteringEnabled(false);
184 // resampleFilter->SetVerboseOptions(true);
185 resampleFilter->Update();
186 working_image = resampleFilter->GetOutput();
189 working_image = object;
191 StopCurrentStep<ImageType>(working_image);
193 // Step 2: object pad to input image -> we want to compute the
194 // relative position for each point belonging to the input image
195 // domain, so we have to extend (pad) the object image to fit the
197 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
198 StartNewStep("Pad object to image size");
199 typename ImageType::Pointer output = ImageType::New();
201 for(unsigned int i=0; i<dim; i++) {
202 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
205 // The index of the input is not necessarily zero, so we have to
206 // take it into account (not done)
208 IndexType index = input->GetLargestPossibleRegion().GetIndex();
209 region.SetSize(size);
210 for(unsigned int i=0; i<dim; i++) {
212 std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
216 // output->SetLargestPossibleRegion(region);
217 output->SetRegions(region);
218 output->SetSpacing(working_image->GetSpacing());
219 PointType origin = input->GetOrigin();
220 for(unsigned int i=0; i<dim; i++) {
221 origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
223 output->SetOrigin(origin);
224 // output->SetOrigin(input->GetOrigin());
227 output->FillBuffer(m_BackgroundValue);
228 typename PadFilterType::Pointer padFilter = PadFilterType::New();
229 // typename PadFilterType::InputImageIndexType index;
230 for(unsigned int i=0; i<dim; i++) {
231 index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
232 + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
234 padFilter->SetSourceImage(working_image);
235 padFilter->SetDestinationImage(output);
236 padFilter->SetDestinationIndex(index);
237 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
239 working_image = padFilter->GetOutput();
240 StopCurrentStep<ImageType>(working_image);
243 // DD("[debug] RelPos : same size and spacing : no padding");
245 // Keep object image (with resampline and pad)
246 object_resampled = working_image;
247 // StopCurrentStep<ImageType>(working_image);
249 // Step 3: compute rel pos in object
250 StartNewStep("Relative Position Map");
251 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
252 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
253 relPosFilter->SetInput(working_image);
254 relPosFilter->SetAlpha1(m_Angle1); // xy plane
255 relPosFilter->SetAlpha2(m_Angle2);
256 relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
257 relPosFilter->SetFast(true);
258 relPosFilter->SetRadius(1); // seems sufficient in this cas
259 // relPosFilter->SetVerboseProgress(true);
260 relPosFilter->Update();
261 relPos = relPosFilter->GetOutput();
262 StopCurrentStep<FloatImageType>(relPos);
264 //--------------------------------------------------------------------
265 //--------------------------------------------------------------------
266 StartNewStep("Map Threshold");
268 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
269 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
270 thresholdFilter->SetInput(relPos);
271 thresholdFilter->SetOutsideValue(m_BackgroundValue);
272 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
273 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
274 thresholdFilter->Update();
275 working_image = thresholdFilter->GetOutput();
276 StopCurrentStep<ImageType>(working_image);
278 //--------------------------------------------------------------------
279 //--------------------------------------------------------------------
280 StartNewStep("Post Processing: erosion with initial mask");
281 // Step 2 : erosion with initial mask to exclude pixels that were
282 // inside the resampled version and outside the original mask
283 typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
284 StructuringElementType kernel;
286 kernel.CreateStructuringElement();
287 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
288 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
289 erodeFilter->SetInput(working_image);
290 erodeFilter->SetKernel(kernel);
291 erodeFilter->SetBackgroundValue(m_BackgroundValue);
292 erodeFilter->SetErodeValue(m_BackgroundValue+1);
293 erodeFilter->Update();
294 working_image = erodeFilter->GetOutput();
295 StopCurrentStep<ImageType>(working_image);
297 //--------------------------------------------------------------------
298 //--------------------------------------------------------------------
299 // Step 5: resample to initial spacing
300 if (m_ResampleBeforeRelativePositionFilter) {
301 StartNewStep("Resample to get the same sampling than input");
302 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
303 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
304 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
305 resampleFilter->SetInput(working_image);
306 resampleFilter->SetOutputSpacing(input->GetSpacing());
307 resampleFilter->SetGaussianFilteringEnabled(false);
308 // resampleFilter->SetVerboseOptions(true);
309 resampleFilter->Update();
310 working_image = resampleFilter->GetOutput();
311 StopCurrentStep<ImageType>(working_image);
314 //--------------------------------------------------------------------
315 //--------------------------------------------------------------------
316 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
317 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
318 StartNewStep("Pad to get the same size than input");
319 typename ImageType::Pointer temp = ImageType::New();
320 temp->CopyInformation(input);
321 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
323 temp->FillBuffer(m_BackgroundValue);
324 typename PadFilterType::Pointer padFilter2 = PadFilterType::New();
325 padFilter2->SetSourceImage(working_image);
326 padFilter2->SetDestinationImage(temp);
327 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
328 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
329 padFilter2->Update();
330 working_image = padFilter2->GetOutput();
331 StopCurrentStep<ImageType>(working_image);
334 //DD("[debug] Rel Pos : no padding after");
337 //--------------------------------------------------------------------
338 //--------------------------------------------------------------------
339 // Step 6: combine input+thresholded relpos
340 StartNewStep("Combine with initial input (boolean And)");
341 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
342 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
343 combineFilter->SetBackgroundValue(m_BackgroundValue);
344 combineFilter->SetBackgroundValue1(m_BackgroundValue);
345 combineFilter->SetBackgroundValue2(m_BackgroundValue);
346 combineFilter->SetForegroundValue(m_BackgroundValue+1);
347 combineFilter->SetInput1(input);
348 combineFilter->SetInput2(working_image);
349 combineFilter->SetOperationType(BoolFilterType::And);
350 combineFilter->InPlaceOn();
351 combineFilter->Update();
352 working_image = combineFilter->GetOutput();
354 combineFilter = BoolFilterType::New();
355 combineFilter->SetInput1(working_image);
356 combineFilter->SetInput2(object);
357 combineFilter->SetOperationType(BoolFilterType::AndNot);
358 combineFilter->InPlaceOn();
359 combineFilter->Update();
361 working_image = combineFilter->GetOutput();
362 StopCurrentStep<ImageType>(working_image);
364 //--------------------------------------------------------------------
365 //--------------------------------------------------------------------
368 StartNewStep("Final AutoCrop");
369 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
370 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
371 cropFilter->SetInput(working_image);
372 cropFilter->ReleaseDataFlagOff();
373 cropFilter->Update();
374 working_image = cropFilter->GetOutput();
375 StopCurrentStep<ImageType>(working_image);
378 //--------------------------------------------------------------------
379 //--------------------------------------------------------------------
381 // Final Step -> set output
382 this->SetNthOutput(0, working_image);
385 //--------------------------------------------------------------------