+ template<class InputImageType>
+ void BLUTDIRGenericFilter::UpdateWithInputImageType()
+ {
+ if (m_Verbose) std::cout << "BLUTDIRGenericFilter::UpdateWithInputImageType()" << std::endl;
+
+ //=============================================================================
+ //Input
+ //=============================================================================
+ bool threadsGiven=m_ArgsInfo.threads_given;
+ int threads=m_ArgsInfo.threads_arg;
+ typedef typename InputImageType::PixelType PixelType;
+
+ typedef double TCoordRep;
+
+ typename InputImageType::Pointer fixedImage = this->template GetInput<InputImageType>(0);
+
+ typename InputImageType::Pointer inputFixedImage = this->template GetInput<InputImageType>(0);
+
+ // typedef input2
+ typename InputImageType::Pointer movingImage = this->template GetInput<InputImageType>(1);
+
+ typename InputImageType::Pointer inputMovingImage = this->template GetInput<InputImageType>(1);
+
+ typedef itk::Image< PixelType,InputImageType::ImageDimension > FixedImageType;
+ typedef itk::Image< PixelType, InputImageType::ImageDimension> MovingImageType;
+ const unsigned int SpaceDimension = InputImageType::ImageDimension;
+ //Whatever the pixel type, internally we work with an image represented in float
+ //Reading reference image
+ if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
+ //=======================================================
+ //Input
+ //=======================================================
+ typename FixedImageType::Pointer croppedFixedImage=fixedImage;
+ //=======================================================
+ // Regions
+ //=======================================================
+ // The original input region
+ typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
+
+ // The transform region with respect to the input region:
+ // where should the transform be DEFINED (depends on mask)
+ typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
+ typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
+ typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
+ typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
+
+ // The metric region with respect to the extracted transform region:
+ // where should the metric be CALCULATED (depends on transform)
+ typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
+ typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
+ typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
+
+
+ //===========================================================================
+ // If given, we connect a mask to reference or target
+ //============================================================================
+ typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
+ typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageLabelType;
+ typename MaskType::Pointer fixedMask = NULL;
+ typename ImageLabelType::Pointer labels = NULL;
+ if (m_ArgsInfo.referenceMask_given)
+ {
+ fixedMask = MaskType::New();
+ labels = ImageLabelType::New();
+ typedef itk::ImageFileReader< ImageLabelType > LabelReaderType;
+ typename LabelReaderType::Pointer labelReader = LabelReaderType::New();
+ labelReader->SetFileName(m_ArgsInfo.referenceMask_arg);
+ try
+ {
+ labelReader->Update();
+ }
+ catch( itk::ExceptionObject & err )
+ {
+ std::cerr << "ExceptionObject caught while reading mask or labels !" << std::endl;
+ std::cerr << err << std::endl;
+ return;
+ }
+ if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
+
+ // Resample labels
+ typedef itk::ResampleImageFilter<ImageLabelType, ImageLabelType> ResamplerType;
+ typename ResamplerType::Pointer resampler = ResamplerType::New();
+ typedef itk::NearestNeighborInterpolateImageFunction<ImageLabelType, TCoordRep> InterpolatorType;
+ typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
+ resampler->SetOutputParametersFromImage(fixedImage);
+ resampler->SetInterpolator(interpolator);
+ resampler->SetInput(labelReader->GetOutput());
+ resampler->Update();
+ labels = resampler->GetOutput();
+
+ // Set the image to the spatialObject
+ fixedMask->SetImage(labels);
+
+ // Find the bounding box of the "inside" label
+ typedef itk::LabelGeometryImageFilter<ImageLabelType> GeometryImageFilterType;
+ typename GeometryImageFilterType::Pointer geometryImageFilter=GeometryImageFilterType::New();
+ geometryImageFilter->SetInput(labels);
+ geometryImageFilter->Update();
+ typename GeometryImageFilterType::BoundingBoxType boundingBox = geometryImageFilter->GetBoundingBox(1);
+
+ // Limit the transform region to the mask
+ for (unsigned int i=0; i<InputImageType::ImageDimension; i++)
+ {
+ transformRegionIndex[i]=boundingBox[2*i];
+ transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
+ }
+ transformRegion.SetSize(transformRegionSize);
+ transformRegion.SetIndex(transformRegionIndex);
+ fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
+
+ // Crop the fixedImage to the bounding box to facilitate multi-resolution
+ typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
+ typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
+#if ITK_VERSION_MAJOR == 4
+ extractImageFilter->SetDirectionCollapseToSubmatrix();
+#endif
+ extractImageFilter->SetInput(fixedImage);
+ extractImageFilter->SetExtractionRegion(transformRegion);
+ extractImageFilter->Update();
+ croppedFixedImage=extractImageFilter->GetOutput();
+
+ // Update the metric region
+ metricRegion = croppedFixedImage->GetLargestPossibleRegion();
+ metricRegionIndex=metricRegion.GetIndex();
+ croppedFixedImage->TransformIndexToPhysicalPoint(metricRegionIndex, metricRegionOrigin);
+
+ // Set start index to zero (with respect to croppedFixedImage/transform region)
+ metricRegionIndex.Fill(0);
+ metricRegion.SetIndex(metricRegionIndex);
+ croppedFixedImage->SetRegions(metricRegion);
+ croppedFixedImage->SetOrigin(metricRegionOrigin);
+ }