SMILX  1.01
milxImage.h
1 /*=========================================================================
2  The Software is copyright (c) Commonwealth Scientific and Industrial Research Organisation (CSIRO)
3  ABN 41 687 119 230.
4  All rights reserved.
5 
6  Licensed under the CSIRO BSD 3-Clause License
7  You may not use this file except in compliance with the License.
8  You may obtain a copy of the License in the file LICENSE.md or at
9 
10  https://stash.csiro.au/projects/SMILI/repos/smili/browse/license.txt
11 
12  Unless required by applicable law or agreed to in writing, software
13  distributed under the License is distributed on an "AS IS" BASIS,
14  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15  See the License for the specific language governing permissions and
16  limitations under the License.
17 =========================================================================*/
18 #ifndef __MILXIMAGE_H
19 #define __MILXIMAGE_H
20 
21 //ITK
22 #include <itkVectorImage.h>
23 #include <itkImageDuplicator.h>
24 #include <itkSpatialOrientationAdapter.h>
25 #include <itkRescaleIntensityImageFilter.h>
26 #include <itkAdaptiveHistogramEqualizationImageFilter.h>
27 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3) //Review only members
28  #include <itkLabelImageToLabelMapFilter.h>
29  #include <itkLabelMapToLabelImageFilter.h>
30  #include <itkAutoCropLabelMapFilter.h>
31  #include <itkBinaryContourImageFilter.h>
32  #include <itkLabelContourImageFilter.h>
33  #include <itkLabelOverlayImageFilter.h>
34  #include <itkMergeLabelMapFilter.h>
35 #endif // (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
36 #include <itkImportImageFilter.h>
37 #include <itkCastImageFilter.h>
38 #include <itkGradientMagnitudeRecursiveGaussianImageFilter.h>
39 #include <itkNormalizeImageFilter.h>
40 #include <itkInvertIntensityImageFilter.h>
41 #include <itkChangeInformationImageFilter.h>
42 #include <itkHistogramMatchingImageFilter.h>
43 #include <itkCheckerBoardImageFilter.h>
44 #include <itkDanielssonDistanceMapImageFilter.h>
45 #include <itkApproximateSignedDistanceMapImageFilter.h>
46 #include <itkSignedMaurerDistanceMapImageFilter.h>
47 #include <itkThresholdImageFilter.h>
48 #include <itkBinaryThresholdImageFilter.h>
49 #include <itkOtsuThresholdImageFilter.h>
50 #include <itkOtsuMultipleThresholdsImageFilter.h>
51 #include <itkMinimumMaximumImageCalculator.h>
52 #include <itkFlipImageFilter.h>
53 #include <itkConstantPadImageFilter.h>
54 #include <itkMaskImageFilter.h>
55 #include <itkGradientAnisotropicDiffusionImageFilter.h>
56 #include <itkBilateralImageFilter.h>
57 #include <itkSmoothingRecursiveGaussianImageFilter.h>
58 #include <itkSobelEdgeDetectionImageFilter.h>
59 #include <itkMedianImageFilter.h>
60 #include <itkLaplacianRecursiveGaussianImageFilter.h>
61 #include <itkCannyEdgeDetectionImageFilter.h>
62 #include <itkIdentityTransform.h>
63 #include <itkExtractImageFilter.h>
64 #include <itkVectorIndexSelectionCastImageFilter.h>
65 #include <itkResampleImageFilter.h>
66 #include <itkAffineTransform.h>
67 #include <itkVersorRigid3DTransform.h>
68 #include <itkCenteredEuler3DTransform.h>
69 #include <itkEuler3DTransform.h>
70 #include <itkTransformFileWriter.h>
71 #include <itkShrinkImageFilter.h>
72 #include <itkAddImageFilter.h>
73 #include <itkSubtractImageFilter.h>
74 #include <itkMultiplyImageFilter.h>
75 #include <itkRegionOfInterestImageFilter.h>
76 #include <itkNearestNeighborInterpolateImageFunction.h>
77 #include <itkBSplineInterpolateImageFunction.h>
78 #include <itkJoinSeriesImageFilter.h>
79 #include <itkConnectedComponentImageFilter.h>
80 #include <itkRelabelComponentImageFilter.h>
81 //ITK4 Compatibility
82 #if (ITK_VERSION_MAJOR > 3)
83  #include <itkv3Rigid3DTransform.h> //ITK4 (4.3.1) version has New() member issues
84  #include <itkFFTConvolutionImageFilter.h>
85  #include <itkVectorMagnitudeImageFilter.h> //vector images
86 #else
87  #include <itkRigid3DTransform.h>
88  #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3) //Review only members
89  #include <itkConvolutionImageFilter.h>
90  #endif
91  #include <itkGradientToMagnitudeImageFilter.h> //vector images
92 #endif
93 //Exporter/Importer
94 #include "itkImageToVTKImageFilter.h"
95 #include "itkVTKImageToImageFilter.h"
96 #ifndef ITK_ONLY
97  //VTK
98  #include <vtkMatrix4x4.h>
99  #include <vtkTransform.h>
100  #include <vtkImageReslice.h>
101 #endif
102 //ITK Extensions
103 #include "itkVectorShiftScaleImageFilter.h"
104 //SMILI
105 #include "milxGlobal.h"
106 
107 const double CoordinateTolerance = 1e-3;
108 const double DirectionTolerance = 1e-3;
109 
110 namespace milx
111 {
112 
123 {
124 public:
129  ImageBase();
134  virtual ~ImageBase() {};
135 };
136 
173 template<class TImage>
175 {
176 public:
181  Image();
186 // Image(itk::SmartPointer<TImage> image);
191  virtual ~Image() {};
192 
197 // void SetInput(itk::SmartPointer<TImage> image);
204  inline itk::SmartPointer<TImage> Result()
205  {
206  return CurrentImage;
207  }
214  inline itk::SmartPointer<TImage> PreviousResult()
215  {
216  return PreviousImage;
217  }
222  inline itk::SmartPointer<TImage> GetOutput()
223  {
224  return Result();
225  }
226 
232  //Conversions
237  static vtkSmartPointer<vtkImageData> ConvertITKImageToVTKImage(itk::SmartPointer<TImage> img);
242  static vtkSmartPointer<vtkImageData> ConvertITKVectorImageToVTKImage(itk::SmartPointer<TImage> img);
243 
244  template<typename TRefImage, class TPrecision>
245  static itk::SmartPointer<TImage> ApplyOrientationToITKImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TRefImage> refImage, const bool labelledImage, const bool flipY = true, const bool ignoreDirection = false);
246 #ifndef ITK_ONLY //Requires VTK
247 
259  static vtkSmartPointer<vtkImageData> ApplyOrientationToVTKImage(vtkSmartPointer<vtkImageData> img, itk::SmartPointer<TImage> refImage, vtkSmartPointer<vtkMatrix4x4> &transformMatrix, const bool labelledImage, const bool flipY = true);
260 #endif
261 
265  static itk::SmartPointer<TImage> ConvertVTKImageToITKImage(vtkSmartPointer<vtkImageData> img);
270  template<typename TOutImage>
271  static itk::SmartPointer<TOutImage> CastImage(itk::SmartPointer<TImage> img);
276  //static void DeepCopy(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToCopyFrom);
281  static itk::SmartPointer<TImage> DuplicateImage(itk::SmartPointer<TImage> img);
286  static itk::SmartPointer<TImage> DuplicateImage(const TImage *img);
293  template<typename TVector>
294  static itk::SmartPointer<TImage> ImportVectorToImage(vnl_vector<TVector> &vec, typename TImage::SizeType size, itk::SmartPointer<TImage> image = NULL);
301  template<typename TMatrix>
302  static itk::SmartPointer<TImage> ImportMatrixToImage(vnl_matrix<TMatrix> &matrix, itk::SmartPointer<TImage> image = NULL);
303 
308  static itk::SmartPointer<TImage> BlankImage(typename TImage::PixelType value, typename TImage::SizeType imgSize);
309 
310  //Geometry
315  template<typename TSubImage>
316  static itk::SmartPointer<TSubImage> ExtractSubImage(itk::SmartPointer<TImage> img, typename TImage::RegionType imgRegion);
323  template<typename TImageSlice>
324  static itk::SmartPointer<TImageSlice> ExtractSlice(itk::SmartPointer<TImage> img, int *extent);
331  template<typename TImageComponent>
332  static itk::SmartPointer<TImageComponent> ExtractComponent(itk::SmartPointer<TImage> img, int component);
337  static itk::SmartPointer<TImage> ResizeImage(itk::SmartPointer<TImage> img, typename TImage::SizeType imgSize);
342  static itk::SmartPointer<TImage> ResizeImage(itk::SmartPointer<TImage> img, typename TImage::SizeType imgSize, typename TImage::SpacingType outputSpacing);
347  static itk::SmartPointer<TImage> ResizeImage(itk::SmartPointer<TImage> img, typename TImage::SizeType imgSize, typename TImage::SpacingType outputSpacing, typename TImage::PointType outputOrigin, typename TImage::DirectionType outputDirection);
355  static itk::SmartPointer<TImage> SubsampleImage(itk::SmartPointer<TImage> img, typename TImage::SizeType factors);
356 
357  //Transform
365  template<typename TOutImage, typename TTransform, typename TPrecision>
366  static itk::SmartPointer<TOutImage> TransformImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TOutImage> refImg, itk::SmartPointer<TTransform> transf, const bool inverse, const int interp = 1);
367  template<typename TOutImage, typename TTransform, typename TPrecision>
368  static itk::SmartPointer<TOutImage> TransformImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TTransform> transf, const bool inverse, const int interp = 1);
373  template<typename TOutImage>
374  static itk::SmartPointer<TOutImage> ResampleImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TOutImage> refImg, const bool linearInterp = false);
379  template<typename TOutImage>
380  static itk::SmartPointer<TOutImage> ResampleLabel(itk::SmartPointer<TImage> img, itk::SmartPointer<TOutImage> refImg);
381 
382  //Arithmetic
387  static itk::SmartPointer<TImage> AddImages(itk::SmartPointer<TImage> img1, itk::SmartPointer<TImage> img2);
392  static itk::SmartPointer<TImage> SubtractImages(itk::SmartPointer<TImage> img1, itk::SmartPointer<TImage> img2);
397  static inline itk::SmartPointer<TImage> DifferenceImages(itk::SmartPointer<TImage> img1, itk::SmartPointer<TImage> img2)
398  {
399  return SubtractImages(img1, img2);
400  }
405  static itk::SmartPointer<TImage> MultiplyImages(itk::SmartPointer<TImage> img1, itk::SmartPointer<TImage> img2);
410  template<class TOutImage>
411  static itk::SmartPointer<TOutImage> ScaleImage(itk::SmartPointer<TImage> img, float scaling);
416  static itk::SmartPointer<TImage> ScaleVectorImage(itk::SmartPointer<TImage> img, float scaling, int numberOfComponents);
417 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
418 
424  static itk::SmartPointer<TImage> ConvolveImages(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> kernelImg);
425 #endif // (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
426 
427  //Info
432  static void Information(itk::SmartPointer<TImage> img);
437  static double ImageMaximum(itk::SmartPointer<TImage> img);
442  static double ImageMinimum(itk::SmartPointer<TImage> img);
447  static std::string ImageOrientation(itk::SmartPointer<TImage> img);
448 
449  //Filters
454  static itk::SmartPointer<TImage> CheckerBoard(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToCheckerBoard, const int numberOfSquares = 0);
459  static itk::SmartPointer<TImage> RescaleIntensity(itk::SmartPointer<TImage> img, float minValue, float maxValue);
464  static itk::SmartPointer<TImage> InvertIntensity(itk::SmartPointer<TImage> img, float maxValue);
473  static itk::SmartPointer<TImage> HistogramEqualisation(itk::SmartPointer<TImage> img, float alpha = 0.3, float beta = 0.3, float radius = 5);
478  static itk::SmartPointer<TImage> GradientMagnitude(itk::SmartPointer<TImage> img);
483  static itk::SmartPointer<TImage> SobelEdges(itk::SmartPointer<TImage> img);
488  static itk::SmartPointer<TImage> Laplacian(itk::SmartPointer<TImage> img);
493  static itk::SmartPointer<TImage> CannyEdges(itk::SmartPointer<TImage> img, float variance, float lowerThreshold, float upperThreshold);
498  static itk::SmartPointer< itk::Image<float, TImage::ImageDimension> > Normalization(itk::SmartPointer<TImage> img);
499 
506  static itk::SmartPointer<TImage> MatchInformation(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToMatch, bool originOnly = false);
513  static itk::SmartPointer<TImage> CopyInformation(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToMatch, bool originOnly = false)
514  {
515  return MatchInformation(img, imgToMatch, originOnly);
516  }
517 
522  static itk::SmartPointer<TImage> MatchHistogram(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToMatch, const int bins = 128);
531  template<typename TOutImage>
532  static itk::SmartPointer<TOutImage> DistanceMap(itk::SmartPointer<TImage> img, const bool binaryImage = false, const bool signedDistances = false, const bool computeInsideObject = false, const bool squaredDistance = false);
537  static itk::SmartPointer<TImage> FlipImage(itk::SmartPointer<TImage> img, bool xAxis = false, bool yAxis = true, bool zAxis = false, bool aboutOrigin = true);
542  static itk::SmartPointer<TImage> PadImageByConstant(itk::SmartPointer<TImage> img, size_t xAxis, size_t yAxis, size_t zAxis, typename TImage::PixelType value);
547  template<typename TOutImage>
548  static itk::SmartPointer<TOutImage> AnisotropicDiffusion(itk::SmartPointer<TImage> img, const int iterations, const float timestep);
553  static itk::SmartPointer<TImage> Bilateral(itk::SmartPointer<TImage> img, const float sigmaRange, const float sigmaSpatial);
558  static itk::SmartPointer<TImage> GaussianSmooth(itk::SmartPointer<TImage> img, const float variance);
563  static itk::SmartPointer<TImage> Median(itk::SmartPointer<TImage> img, const int radius);
564 
565  //Labelling
570  template<typename TMaskImage>
571  static itk::SmartPointer<TImage> MaskImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TMaskImage> maskImg);
576  static itk::SmartPointer<TImage> RelabelImage(itk::SmartPointer<TImage> labelledImg);
577 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
578 
582  static itk::SmartPointer<TImage> BinaryContour(itk::SmartPointer<TImage> img, const float foregroundValue, const float backgroundValue);
587  static itk::SmartPointer<TImage> LabelContour(itk::SmartPointer<TImage> img, const bool fullyConnected, const float backgroundValue);
592  template<typename TMaskImage>
593  static itk::SmartPointer<TImage> MaskAndCropImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TMaskImage> maskImg, const size_t pixelPadding = 1);
598  static std::vector<unsigned char> LabelValues(itk::SmartPointer<TImage> labelledImg);
603  template<typename TLabelImage>
604  static itk::SmartPointer< itk::Image<itk::RGBPixel<unsigned char>, TImage::ImageDimension> > Overlay(itk::SmartPointer<TImage> img, itk::SmartPointer<TLabelImage> overlayImg);
609  template<typename TLabelImage>
610  static itk::SmartPointer< itk::Image<itk::RGBPixel<unsigned char>, TImage::ImageDimension> > OverlayContour(itk::SmartPointer<TImage> img, itk::SmartPointer<TLabelImage> overlayImg);
611 #endif // (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
612 
613  //Thresholding
618  static itk::SmartPointer<TImage> ThresholdAboveImage(itk::SmartPointer<TImage> img, float outsideValue, float aboveValue);
623  static itk::SmartPointer<TImage> ThresholdBelowImage(itk::SmartPointer<TImage> img, float outsideValue, float aboveValue);
628  static itk::SmartPointer<TImage> ThresholdImage(itk::SmartPointer<TImage> img, float outsideValue, float belowValue, float aboveValue);
633  template<typename TOutImage>
634  static itk::SmartPointer<TOutImage> BinaryThresholdImage(itk::SmartPointer<TImage> img, float outsideValue, float insideValue, float belowValue, float aboveValue);
639  static double OtsuThreshold(itk::SmartPointer<TImage> img, const int bins);
646  template<typename TOutImage>
647  static itk::SmartPointer<TOutImage> OtsuThresholdImage(itk::SmartPointer<TImage> img, const int bins);
654  template<typename TOutImage>
655  static itk::SmartPointer<TOutImage> OtsuMultipleThresholdImage(itk::SmartPointer<TImage> img, const int bins, const int noOfLabels = 1);
656 
657  //Vector operations
662  template<typename TScalarImage>
663  static itk::SmartPointer<TScalarImage> VectorMagnitude(itk::SmartPointer<TImage> img);
665 
666  //Collection operations
676  static void InformationCollection(std::vector< typename itk::SmartPointer<TImage> > &images);
681  static void MatchHistogramCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TImage> refImage, const int bins = 128);
686  static void CheckerboardCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TImage> refImage, const int squares = 10);
691  static void RescaleIntensityCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float belowValue, float aboveValue);
696  static void InvertIntensityCollection(std::vector< typename itk::SmartPointer<TImage> > &images);
701  static void FlipCollection(std::vector< typename itk::SmartPointer<TImage> > &images, bool xAxis = false, bool yAxis = true, bool zAxis = false, bool aboutOrigin = true);
708  template<typename TOutImage>
709  static std::vector< typename itk::SmartPointer<TOutImage> > AnisotropicDiffusionCollection(const std::vector< typename itk::SmartPointer<TImage> > &images, const int iterations, float timestep = -1.0);
714  static void BilateralCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float sigmaRange, float sigmaSpatial);
719  template<typename TOutImage>
720  static std::vector< typename itk::SmartPointer<TOutImage> > CastCollection(const std::vector< typename itk::SmartPointer<TImage> > &images);
725  static void MedianCollection(std::vector< typename itk::SmartPointer<TImage> > &images, const int radius);
730  static void GradientMagnitudeCollection(std::vector< typename itk::SmartPointer<TImage> > &images);
735  static void LaplacianCollection(std::vector< typename itk::SmartPointer<TImage> > &images);
740  template<typename TOutImage>
741  static std::vector< typename itk::SmartPointer<TOutImage> > DistanceMapCollection(const std::vector< typename itk::SmartPointer<TImage> > &images, const bool binaryImage = false, const bool signedDistances = false, const bool computeInsideObject = false, const bool squaredDistance = false);
742 
743  //Thresholding
748  static void ThresholdAboveCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float outsideValue, float aboveValue);
753  static void ThresholdBelowCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float outsideValue, float belowValue);
758  static void ThresholdCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float outsideValue, float belowValue, float aboveValue);
763  template<typename TOutImage>
764  static std::vector< typename itk::SmartPointer<TOutImage> > BinaryThresholdCollection(const std::vector< typename itk::SmartPointer<TImage> > &images, float outsideValue, float insideValue, float belowValue, float aboveValue);
769  template<typename TOutImage>
770  static std::vector< typename itk::SmartPointer<TOutImage> > OtsuMultipleThresholdCollection(const std::vector< typename itk::SmartPointer<TImage> > &images, const int bins, const int noOfLabels);
771 
772 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
773 
777  template<class TMaskImage>
778  static void MaskAndCropCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TMaskImage> maskImage, const size_t pixelPadding = 1);
779 #endif // (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
780 
784  template<class TMaskImage>
785  static void MaskCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TMaskImage> maskImage);
790  static void RelabelCollection(std::vector< typename itk::SmartPointer<TImage> > &images);
795  template<class TRefImage>
796  static void ResampleCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TRefImage> refImage);
801  template<class TRefImage>
802  static void ResampleLabelCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TRefImage> refImage);
807  static itk::SmartPointer<TImage> AddCollection(std::vector< typename itk::SmartPointer<TImage> > &images);
812  static itk::SmartPointer<TImage> SubtractCollection(std::vector< typename itk::SmartPointer<TImage> > &images);
817  static inline itk::SmartPointer<TImage> DifferenceCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
818  {
819  return SubtractCollection(images);
820  }
825  template<class TOutImage>
826  static itk::SmartPointer<TOutImage> AverageCollection(std::vector< typename itk::SmartPointer<TImage> > &images);
831  static itk::SmartPointer<TImage> AverageVectorCollection(std::vector< typename itk::SmartPointer<TImage> > &images, int numberOfComponents);
832  template<class TJoinedImage>
833  itk::SmartPointer<TJoinedImage> JoinCollection(std::vector< typename itk::SmartPointer<TImage> > &images, const double origin, const double spacing);
834 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
835 
845  static itk::SmartPointer<TImage> MergeLabelledImages(std::vector< typename itk::SmartPointer<TImage> > &images, const size_t mergeType = 0);
846 #endif // (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
847 
848 
849 protected:
850  itk::SmartPointer<TImage> CurrentImage;
851  itk::SmartPointer<TImage> PreviousImage;
852 
853 private:
854 
855 };
856 
857 template<class TImage>
859 {
860 
861 }
862 
863 //template<class TImage>
864 //Image<TImage>::Image(itk::SmartPointer<TImage> image)
865 //{
866 // SetInput(image);
867 //}
868 
869 //template<class TImage>
870 //void Image<TImage>::SetInput(itk::SmartPointer<TImage> image)
871 //{
872 // if(!CurrentImage)
873 // CurrentImage = itk::SmartPointer<TImage>::New();
874 // CurrentImage->DuplicateImage(image);
875 //}
876 
877 //Arithmetic
878 template<class TImage>
879 vtkSmartPointer<vtkImageData> Image<TImage>::ConvertITKImageToVTKImage(itk::SmartPointer<TImage> img)
880 {
881  typedef itk::ImageToVTKImageFilter<TImage> ConvertImageType;
882 
883  typename ConvertImageType::Pointer convertFilter = ConvertImageType::New();
884  convertFilter->SetInput(img);
885  convertFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
886  try
887  {
888  convertFilter->Update();
889  }
890  catch (itk::ExceptionObject & ex )
891  {
892  PrintError("Failed Converting ITK Image to VTK Image");
893  PrintError(ex.GetDescription());
894  }
895 
896  return convertFilter->GetOutput();
897 }
898 
899 template<class TImage>
900 vtkSmartPointer<vtkImageData> Image<TImage>::ConvertITKVectorImageToVTKImage(itk::SmartPointer<TImage> img)
901 {
902  typedef itk::Image<typename TImage::InternalPixelType, TImage::ImageDimension> TImageComponent;
903 
904  //Get the size etc. right of result image
905  typename TImageComponent::Pointer componentImg = milx::Image<TImage>::ExtractComponent<TImageComponent>(img, 0);
906  vtkSmartPointer<vtkImageData> fieldComponent = milx::Image<TImageComponent>::ConvertITKImageToVTKImage(componentImg);
907 
908  vtkSmartPointer<vtkImageData> field = vtkSmartPointer<vtkImageData>::New();
909  field->SetDimensions(fieldComponent->GetDimensions());
910  field->SetOrigin(fieldComponent->GetOrigin());
911  field->SetSpacing(fieldComponent->GetSpacing());
912  #if VTK_MAJOR_VERSION <= 5
913  field->SetNumberOfScalarComponents(img->GetNumberOfComponentsPerPixel());
914  field->SetScalarTypeToFloat();
915  field->AllocateScalars();
916  #else
917  field->AllocateScalars(VTK_FLOAT, 3);
918  #endif
919  /*vtkSmartPointer<vtkImageAppendComponents> appendImage = vtkSmartPointer<vtkImageAppendComponents>::New(); //doesnt work!
920  appendImage->DebugOn();
921  appendImage->AddObserver(vtkCommand::ProgressEvent, VTKProgressUpdates);
922  appendImage->SetInput(fieldComponent);*/
923 
924  for(vtkIdType j = 0; j < img->GetNumberOfComponentsPerPixel(); j ++)
925  {
926  typename TImageComponent::Pointer componentImgTmp = milx::Image<TImage>::ExtractComponent<TImageComponent>(img, j);
927  vtkSmartPointer<vtkImageData> fieldComponentTmp = milx::Image<TImageComponent>::ConvertITKImageToVTKImage(componentImgTmp);
928 
929  for(int k = 0; k < fieldComponent->GetDimensions()[0]; k ++)
930  for(int l = 0; l < fieldComponent->GetDimensions()[1]; l ++)
931  for(int m = 0; m < fieldComponent->GetDimensions()[2]; m ++)
932  field->SetScalarComponentFromFloat(k, l, m, j, fieldComponentTmp->GetScalarComponentAsFloat(k, l, m, 0));
933  }
934 
935  return field;
936 }
937 
938 template<class TImage>
939 template<class TRefImage, class TPrecision>
940 itk::SmartPointer<TImage> Image<TImage>::ApplyOrientationToITKImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TRefImage> refImage, const bool labelledImage, const bool flipY, const bool ignoreDirection)
941 {
942  typename TImage::DirectionType direction = refImage->GetDirection();
943  typename TImage::PointType origin = refImage->GetOrigin();
944 
945  typedef itk::InterpolateImageFunction<TImage, TPrecision> InterpolatorType;
946  typename InterpolatorType::Pointer interpolator;
947  if(labelledImage)
948  interpolator = itk::NearestNeighborInterpolateImageFunction<TImage, TPrecision>::New();
949  else
950  interpolator = itk::LinearInterpolateImageFunction<TImage, TPrecision>::New();
951 
952  typedef itk::IdentityTransform<TPrecision, 3> TransformType;
953  typename TransformType::Pointer identityTransform = TransformType::New();
954 
955  /*
956  typename TImage::DirectionType identityMatrix;
957  identityMatrix.SetIdentity();
958 // typename TImage::PointType axesOrigin;
959 // axesOrigin.Fill(0.0);
960 
961  //Remove origin and direction to get raw volume
962  typedef itk::ChangeInformationImageFilter<TImage> ChangeFilterType;
963  typename ChangeFilterType::Pointer stripInfo = ChangeFilterType::New();
964  stripInfo->SetInput( img );
965  stripInfo->ChangeDirectionOn();
966  stripInfo->SetOutputDirection(identityMatrix);
967 // stripInfo->ChangeOriginOn();
968 // stripInfo->SetOutputOrigin(axesOrigin);
969 
970  vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New(); //start with identity matrix
971  matrix->Identity();
972  for (int i = 0; i < 3; i ++)
973  for (int k = 0; k < 3; k ++)
974  matrix->SetElement(i,k, direction(i,k));
975  matrix->Transpose(); //img volume to correct space, comment to go from space to img volume
976 
977  vtkSmartPointer<vtkMatrix4x4> transformMatrix = vtkSmartPointer<vtkMatrix4x4>::New(); //start with identity matrix
978  transformMatrix->Identity();
979  //Flip the image for VTK coordinate system
980 // if(flipY)
981 // transformMatrix->SetElement(1,1,-1); //flip
982 
983  vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
984  transform->Identity();
985  transform->PostMultiply();
986  transform->Concatenate(transformMatrix); //flip
987  transform->Translate(-origin[0], -origin[1], -origin[2]); //remove origin
988  transform->Concatenate(matrix); //direction
989  transform->Translate(origin.GetDataPointer()); //add origin displacement
990  transform->Update();
991  vtkSmartPointer<vtkMatrix4x4> matrix2 = transform->GetMatrix();
992 
993  typename TImage::DirectionType orientMatrix;
994  for (int i = 0; i < 3; i ++)
995  for (int k = 0; k < 3; k ++)
996  orientMatrix(i,k) = matrix2->GetElement(i,k);
997  typename itk::Vector<TPrecision, milx::imgDimension> offset;
998 // typename TImage::PointType new_origin;
999  for (int i = 0; i < 3; i ++)
1000  {
1001 // new_origin[i] = matrix2->GetElement(i,3);
1002  offset[i] = matrix2->GetElement(i,3);
1003  }
1004  matrix2->Print(cout);*/
1005 
1006  typedef itk::ResampleImageFilter<TImage, TImage, TPrecision> ResampleImageFilterType;
1007  typename ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
1008  resample->SetInput(img);
1009  resample->SetDefaultPixelValue(0.0);
1010  resample->SetTransform(identityTransform);
1011  resample->SetSize(refImage->GetLargestPossibleRegion().GetSize());
1012 // resample->SetOutputOrigin( axesOrigin );
1013  resample->SetOutputSpacing(refImage->GetSpacing());
1014 // resample->SetOutputDirection( identityMatrix ); // 1 1 1 direction
1015  resample->SetInterpolator(interpolator);
1016  resample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1017 
1018  try
1019  {
1020  resample->UpdateLargestPossibleRegion();
1021  resample->Update();
1022  }
1023  catch (itk::ExceptionObject & ex )
1024  {
1025  PrintError("Failed applying orientation to Image");
1026  PrintError(ex.GetDescription());
1027  }
1028 
1029  return resample->GetOutput();
1030 }
1031 
1032 #ifndef ITK_ONLY //Requires VTK
1033 template<class TImage>
1034 vtkSmartPointer<vtkImageData> Image<TImage>::ApplyOrientationToVTKImage(vtkSmartPointer<vtkImageData> img, itk::SmartPointer<TImage> refImage, vtkSmartPointer<vtkMatrix4x4> &transformMatrix, const bool labelledImage, const bool flipY)
1035 {
1036  typename TImage::DirectionType direction = refImage->GetDirection();
1037  typename TImage::PointType origin = refImage->GetOrigin();
1038 
1039  vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New(); //start with identity matrix
1040  matrix->Identity();
1041  for (int i = 0; i < 3; i ++)
1042  for (int k = 0; k < 3; k ++)
1043  matrix->SetElement(i,k, direction(i,k));
1044  matrix->Transpose(); //img volume to correct space, comment to go from space to img volume
1045 
1046  if(!transformMatrix)
1047  transformMatrix = vtkSmartPointer<vtkMatrix4x4>::New(); //start with identity matrix
1048  transformMatrix->Identity();
1049  //Flip the image for VTK coordinate system
1050  if(flipY)
1051  transformMatrix->SetElement(1,1,-1); //flip
1052 
1053  vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
1054  transform->Identity();
1055  transform->PostMultiply();
1056  transform->Concatenate(transformMatrix); //flip
1057  transform->Translate(-origin[0], -origin[1], -origin[2]); //remove origin
1058  transform->Concatenate(matrix); //direction
1059  transform->Translate(origin.GetDataPointer()); //add origin displacement
1060 
1061  vtkSmartPointer<vtkImageReslice> reslice = vtkSmartPointer<vtkImageReslice>::New();
1062 #if VTK_MAJOR_VERSION <= 5
1063  reslice->SetInput(img);
1064 #else
1065  reslice->SetInputData(img);
1066 #endif
1067  reslice->AddObserver(vtkCommand::ProgressEvent, VTKProgressUpdates);
1068  reslice->AutoCropOutputOn();
1069  reslice->TransformInputSamplingOn();
1070  reslice->SetOutputDimensionality(milx::imgDimension);
1071  reslice->SetResliceAxes(transform->GetMatrix());
1072 
1073  if(!labelledImage)
1074  reslice->SetInterpolationModeToLinear(); //reduce artefacts for normal images
1075  else
1076  {
1077  reslice->SetInterpolationModeToNearestNeighbor(); //reduce artefacts for labels
1078  PrintDebug("Using NN Interpolator for Image Orientation.");
1079 // #if VTK_MAJOR_VERSION > 5
1080 // reslice->SetOutputScalarType(VTK_CHAR);
1081 // #endif
1082  }
1083  reslice->Update();
1084 
1085  std::string resliceType = reslice->GetOutput()->GetScalarTypeAsString();
1086  PrintDebug("Orientated Image Type as " + resliceType);
1087 
1088  return reslice->GetOutput();
1089 }
1090 #endif
1091 
1092 template<class TImage>
1093 itk::SmartPointer<TImage> Image<TImage>::ConvertVTKImageToITKImage(vtkSmartPointer<vtkImageData> img)
1094 {
1095  typedef itk::VTKImageToImageFilter<TImage> ConvertImageType;
1096 
1097  typename ConvertImageType::Pointer convertFilter = ConvertImageType::New();
1098  convertFilter->SetInput(img);
1099  convertFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1100  try
1101  {
1102  convertFilter->Update();
1103  }
1104  catch (itk::ExceptionObject & ex )
1105  {
1106  PrintError("Failed Converting VTK Image to ITK Image");
1107  PrintError(ex.GetDescription());
1108  }
1109 
1110  itk::SmartPointer<TImage> tmpImage = DuplicateImage(convertFilter->GetOutput());
1111 
1112  return tmpImage;
1113 }
1114 
1115 template<class TImage>
1116 template<typename TOutImage>
1117 itk::SmartPointer<TOutImage> Image<TImage>::CastImage(itk::SmartPointer<TImage> img)
1118 {
1119  typedef itk::CastImageFilter<TImage, TOutImage> CastImageType;
1120 
1121  typename CastImageType::Pointer castFilter = CastImageType::New();
1122  castFilter->SetInput(img);
1123  castFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1124  try
1125  {
1126  castFilter->Update();
1127  }
1128  catch (itk::ExceptionObject & ex )
1129  {
1130  PrintError("Failed Casting");
1131  PrintError(ex.GetDescription());
1132  }
1133 
1134  return castFilter->GetOutput();
1135 }
1136 
1137 /*template<class TImage>
1138 template<typename TOutImage>
1139 void Image<TImage>::DeepCopy(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToCopyFrom, TImage::RegionType region)
1140 {
1141  itk::ImageRegionConstIterator<TImage> inputIterator(imgToCopyFrom, region);
1142  itk::ImageRegionIterator<TImage> outputIterator(img, region);
1143  while(!inputIterator.IsAtEnd())
1144  {
1145  outputIterator.Set(inputIterator.Get());
1146  ++inputIterator;
1147  ++outputIterator;
1148  }
1149 }*/
1150 
1151 template<class TImage>
1152 itk::SmartPointer<TImage> Image<TImage>::DuplicateImage(itk::SmartPointer<TImage> img)
1153 {
1154  typedef itk::ImageDuplicator<TImage> DuplicateType;
1155 
1156  typename DuplicateType::Pointer duplicator = DuplicateType::New();
1157  duplicator->SetInputImage(img);
1158  duplicator->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1159  try
1160  {
1161  duplicator->Update();
1162  }
1163  catch (itk::ExceptionObject & ex )
1164  {
1165  PrintError("Failed Generating Duplicate Image");
1166  PrintError(ex.GetDescription());
1167  }
1168 
1169  return duplicator->GetOutput();
1170 }
1171 
1172 template<class TImage>
1173 itk::SmartPointer<TImage> Image<TImage>::DuplicateImage(const TImage *img)
1174 {
1175  typedef itk::ImageDuplicator<TImage> DuplicateType;
1176 
1177  typename DuplicateType::Pointer duplicator = DuplicateType::New();
1178  duplicator->SetInputImage(img);
1179  duplicator->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1180  try
1181  {
1182  duplicator->Update();
1183  }
1184  catch (itk::ExceptionObject & ex )
1185  {
1186  PrintError("Failed Generating Duplicate Image");
1187  PrintError(ex.GetDescription());
1188  }
1189 
1190  return duplicator->GetOutput();
1191 }
1192 
1193 template<class TImage>
1194 template<typename TVector>
1195 itk::SmartPointer<TImage> Image<TImage>::ImportVectorToImage(vnl_vector<TVector> &vec, typename TImage::SizeType size, itk::SmartPointer<TImage> image)
1196 {
1197  //Pass by reference is necessary here.
1198  typedef itk::ImportImageFilter<typename TImage::PixelType, TImage::ImageDimension> ImportFilterType;
1199  typename ImportFilterType::Pointer importFilter = ImportFilterType::New();
1200 
1201  typename ImportFilterType::IndexType start;
1202  start.Fill( 0 );
1203 
1204  typename ImportFilterType::RegionType region;
1205  region.SetIndex(start);
1206  region.SetSize(size);
1207 
1208  importFilter->SetRegion( region );
1209  if(image)
1210  {
1211  importFilter->SetOrigin( image->GetOrigin() );
1212  importFilter->SetSpacing( image->GetSpacing() );
1213  importFilter->SetDirection( image->GetDirection() );
1214  }
1215 
1216  const bool importImageFilterWillOwnTheBuffer = false;
1217  importFilter->SetImportPointer( vec.data_block(), vec.size(), importImageFilterWillOwnTheBuffer );
1218  importFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1219  try
1220  {
1221  importFilter->Update();
1222  }
1223  catch (itk::ExceptionObject & ex )
1224  {
1225  PrintError("Failed Converting vector to ITK Image");
1226  PrintError(ex.GetDescription());
1227  }
1228 
1229  return importFilter->GetOutput();
1230 }
1231 
1232 template<class TImage>
1233 template<typename TMatrix>
1234 itk::SmartPointer<TImage> Image<TImage>::ImportMatrixToImage(vnl_matrix<TMatrix> &matrix, itk::SmartPointer<TImage> image)
1235 {
1236  //Pass by reference is necessary here.
1237  typedef itk::ImportImageFilter<typename TImage::PixelType, TImage::ImageDimension> ImportFilterType;
1238  typename ImportFilterType::Pointer importFilter = ImportFilterType::New();
1239 
1240  typename ImportFilterType::SizeType size;
1241  size[0] = matrix.rows(); // size along X
1242  size[1] = matrix.cols(); // size along Y
1243  size[2] = 1; // size along Z
1244 
1245  typename ImportFilterType::IndexType start;
1246  start.Fill( 0 );
1247 
1248  typename ImportFilterType::RegionType region;
1249  region.SetIndex(start);
1250  region.SetSize(size);
1251 
1252  importFilter->SetRegion( region );
1253  if(image)
1254  {
1255  importFilter->SetOrigin( image->GetOrigin() );
1256  importFilter->SetSpacing( image->GetSpacing() );
1257  importFilter->SetDirection( image->GetDirection() );
1258  }
1259 
1260  const bool importImageFilterWillOwnTheBuffer = false;
1261  importFilter->SetImportPointer( matrix.data_block(), matrix.size(), importImageFilterWillOwnTheBuffer );
1262  importFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1263  try
1264  {
1265  importFilter->Update();
1266  }
1267  catch (itk::ExceptionObject & ex )
1268  {
1269  PrintError("Failed Converting matrix to ITK Image");
1270  PrintError(ex.GetDescription());
1271  }
1272 
1273  return importFilter->GetOutput();
1274 }
1275 
1276 template<class TImage>
1277 itk::SmartPointer<TImage> Image<TImage>::BlankImage(typename TImage::PixelType value, typename TImage::SizeType imgSize)
1278 {
1279  typename TImage::IndexType blankStart;
1280  blankStart.Fill(0);
1281 
1282  typename TImage::RegionType blankRegion;
1283  blankRegion.SetIndex(blankStart);
1284  blankRegion.SetSize(imgSize);
1285 
1286  //create image
1287  typename TImage::Pointer img = TImage::New();
1288  img->SetRegions(blankRegion);
1289  img->Allocate();
1290  img->FillBuffer(value);
1291 
1292  return img;
1293 }
1294 
1295 //Geometry
1296 template<class TImage>
1297 template<typename TSubImage>
1298 itk::SmartPointer<TSubImage> Image<TImage>::ExtractSubImage(itk::SmartPointer<TImage> img, typename TImage::RegionType imgRegion)
1299 {
1300  typedef itk::ExtractImageFilter<TImage, TSubImage> FilterType;
1301  typename FilterType::Pointer extractor = FilterType::New();
1302  extractor->SetExtractionRegion(imgRegion);
1303  extractor->SetInput(img);
1304 #if ITK_VERSION_MAJOR >= 4
1305  extractor->SetDirectionCollapseToIdentity(); // This is required.
1306 #endif
1307  extractor->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1308 
1309  try
1310  {
1311  extractor->Update();
1312  }
1313  catch (itk::ExceptionObject & ex )
1314  {
1315  PrintError("Failed Extracting Sub-Image");
1316  PrintError(ex.GetDescription());
1317  }
1318 
1319  return extractor->GetOutput();
1320 }
1321 
1322 template<class TImage>
1323 template<typename TImageSlice>
1324 itk::SmartPointer<TImageSlice> Image<TImage>::ExtractSlice(itk::SmartPointer<TImage> img, int *extent)
1325 {
1326  typename TImage::IndexType desiredStart;
1327  desiredStart[0] = extent[0];
1328  desiredStart[1] = extent[2];
1329  desiredStart[2] = extent[4];
1330  milx::PrintDebug("Slice Start Indices: " + milx::NumberToString(desiredStart[0]) + ", " + milx::NumberToString(desiredStart[1]) + ", " + milx::NumberToString(desiredStart[2]));
1331 
1332  typename TImage::SizeType desiredSize;
1333  desiredSize[0] = extent[1]-extent[0];
1334  desiredSize[1] = extent[3]-extent[2];
1335  desiredSize[2] = extent[5]-extent[4];
1336 
1337  const size_t sliceDim = TImageSlice::ImageDimension;
1338  const size_t imgDim = TImage::ImageDimension;
1339 
1340  if(sliceDim == imgDim) //then just slicing data and not compacting a diemnsion
1341  {
1342  for(size_t j = 0; j < imgDim; j ++)
1343  {
1344  if(desiredSize[j] == 0) //dont compact
1345  desiredSize[j] += 1;
1346  }
1347  }
1348 
1349  milx::PrintDebug("Slice Size: " + milx::NumberToString(desiredSize[0]) + ", " + milx::NumberToString(desiredSize[1]) + ", " + milx::NumberToString(desiredSize[2]));
1350 
1351  typename TImage::RegionType desiredRegion(desiredStart, desiredSize);
1352 
1353  return ExtractSubImage<TImageSlice>(img, desiredRegion);
1354 }
1355 
1356 template<class TImage>
1357 template<typename TImageComponent>
1358 itk::SmartPointer<TImageComponent> Image<TImage>::ExtractComponent(itk::SmartPointer<TImage> img, int component)
1359 {
1360  typedef itk::VectorIndexSelectionCastImageFilter<TImage, TImageComponent> IndexSelectionType;
1361  typename IndexSelectionType::Pointer indexSelectionFilter = IndexSelectionType::New();
1362  indexSelectionFilter->SetIndex(component);
1363  indexSelectionFilter->SetInput(img);
1364 
1365  try
1366  {
1367  indexSelectionFilter->Update();
1368  }
1369  catch (itk::ExceptionObject & ex )
1370  {
1371  PrintError("Failed Extracting Component");
1372  PrintError(ex.GetDescription());
1373  }
1374 
1375  return indexSelectionFilter->GetOutput();
1376 }
1377 
1378 template<class TImage>
1379 itk::SmartPointer<TImage> Image<TImage>::ResizeImage(itk::SmartPointer<TImage> img, typename TImage::SizeType imgSize)
1380 {
1381  typedef itk::IdentityTransform<double, imgDimension> TransformType;
1382  typedef itk::ResampleImageFilter<TImage, TImage> ResampleImageFilterType;
1383  typename ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
1384  resample->SetInput(img);
1385  resample->SetSize(imgSize);
1386  resample->SetOutputSpacing(img->GetSpacing());
1387  resample->SetTransform(TransformType::New());
1388  resample->UpdateLargestPossibleRegion();
1389  resample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1390 
1391  try
1392  {
1393  resample->Update();
1394  }
1395  catch (itk::ExceptionObject & ex )
1396  {
1397  PrintError("Failed Resizing");
1398  PrintError(ex.GetDescription());
1399  }
1400 
1401  return resample->GetOutput();
1402 }
1403 
1404 template<class TImage>
1405 itk::SmartPointer<TImage> Image<TImage>::ResizeImage(itk::SmartPointer<TImage> img, typename TImage::SizeType imgSize, typename TImage::SpacingType outputSpacing)
1406 {
1407  typedef itk::IdentityTransform<double, imgDimension> TransformType;
1408  typedef itk::ResampleImageFilter<TImage, TImage> ResampleImageFilterType;
1409  typename ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
1410  resample->SetInput(img);
1411  resample->SetSize(imgSize);
1412  resample->SetOutputSpacing(outputSpacing);
1413  resample->SetTransform(TransformType::New());
1414  resample->UpdateLargestPossibleRegion();
1415  resample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1416 
1417  try
1418  {
1419  resample->Update();
1420  }
1421  catch (itk::ExceptionObject & ex )
1422  {
1423  PrintError("Failed Resizing");
1424  PrintError(ex.GetDescription());
1425  }
1426 
1427  return resample->GetOutput();
1428 }
1429 
1430 template<class TImage>
1431 itk::SmartPointer<TImage> Image<TImage>::ResizeImage(itk::SmartPointer<TImage> img, typename TImage::SizeType imgSize, typename TImage::SpacingType outputSpacing, typename TImage::PointType outputOrigin, typename TImage::DirectionType outputDirection)
1432 {
1433  typedef itk::IdentityTransform<double, imgDimension> TransformType;
1434  typedef itk::ResampleImageFilter<TImage, TImage> ResampleImageFilterType;
1435  typename ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
1436  resample->SetInput(img);
1437  resample->SetSize(imgSize);
1438  resample->SetOutputSpacing(outputSpacing);
1439  resample->SetOutputOrigin(outputOrigin);
1440  resample->SetOutputDirection(outputDirection);
1441  resample->SetTransform(TransformType::New());
1442  resample->UpdateLargestPossibleRegion();
1443  resample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1444 
1445  try
1446  {
1447  resample->Update();
1448  }
1449  catch (itk::ExceptionObject & ex )
1450  {
1451  PrintError("Failed Resizing");
1452  PrintError(ex.GetDescription());
1453  }
1454 
1455  return resample->GetOutput();
1456 }
1457 
1458 template<class TImage>
1459 itk::SmartPointer<TImage> Image<TImage>::SubsampleImage(itk::SmartPointer<TImage> img, typename TImage::SizeType factors)
1460 {
1461  typedef itk::ShrinkImageFilter<TImage, TImage> SubSampleImageFilterType;
1462  typename SubSampleImageFilterType::Pointer subsample = SubSampleImageFilterType::New();
1463  subsample->SetInput(img);
1464  for(size_t j = 0; j < TImage::ImageDimension; j ++)
1465  subsample->SetShrinkFactor(j, factors[j]); // shrink the first dimension by a factor of 2 (i.e. 100 gets changed to 50)
1466  subsample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1467 
1468  try
1469  {
1470  subsample->Update();
1471  }
1472  catch (itk::ExceptionObject & ex )
1473  {
1474  PrintError("Failed Subsampling");
1475  PrintError(ex.GetDescription());
1476  }
1477 
1478  return subsample->GetOutput();
1479 }
1480 
1481 //Transform
1482 template<class TImage>
1483 template<typename TOutImage, typename TTransform, typename TPrecision>
1484 itk::SmartPointer<TOutImage> Image<TImage>::TransformImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TOutImage> refImg, itk::SmartPointer<TTransform> transf, const bool inverse, const int interp)
1485 {
1486  typedef itk::InterpolateImageFunction<TImage, TPrecision> InterpolatorType;
1487  typename InterpolatorType::Pointer interpolator;
1488  if(interp == 0)
1489  {
1490  interpolator = itk::NearestNeighborInterpolateImageFunction<TImage, TPrecision>::New();
1491  PrintDebug("Using nearest neighbour interpolation for resampling.");
1492  }
1493  else if(interp == 1)
1494  {
1495  interpolator = itk::LinearInterpolateImageFunction<TImage, TPrecision>::New();
1496  PrintDebug("Using linear interpolation for resampling.");
1497  }
1498  else
1499  {
1500  typedef itk::BSplineInterpolateImageFunction<TImage, TPrecision, TPrecision> BSplineInterpolatorType;
1501  typename BSplineInterpolatorType::Pointer bsplineInterpolator = BSplineInterpolatorType::New();
1502  bsplineInterpolator->SetSplineOrder(3);
1503  interpolator = bsplineInterpolator;
1504  PrintDebug("Using BSpline interpolation for resampling.");
1505  }
1506 
1507  typedef itk::ResampleImageFilter<TImage, TOutImage, TPrecision> ResampleImageFilterType;
1508  typename ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
1509  resample->SetInput(img);
1510  resample->SetDefaultPixelValue(0.0);
1511  resample->UseReferenceImageOn();
1512  resample->SetReferenceImage(refImg);
1513 // resample->SetSize( refImg->GetLargestPossibleRegion().GetSize() );
1514 // resample->SetOutputOrigin( refImg->GetOrigin() );
1515 // resample->SetOutputSpacing( refImg->GetSpacing() );
1516 // resample->SetOutputDirection( refImg->GetDirection() );
1517  resample->SetInterpolator(interpolator);
1518  resample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1519 
1521  typedef itk::AffineTransform< TPrecision, milx::imgDimension > AffineTransformType;
1522  typedef typename AffineTransformType::Pointer AffineTransformPointer;
1523  typedef itk::VersorRigid3DTransform<TPrecision> VersorRigidTransformType;
1524  typedef typename VersorRigidTransformType::Pointer VersorRigidTransformPointer;
1525 #if ITK_VERSION_MAJOR > 3 //New() member issues, see ITK4 tests for detail
1526  typedef itkv3::Rigid3DTransform<TPrecision> RigidTransformType;
1527 #else
1528  typedef itk::Rigid3DTransform<TPrecision> RigidTransformType;
1529 #endif
1530  typedef typename RigidTransformType::Pointer RigidTransformPointer;
1531  typedef itk::CenteredEuler3DTransform<TPrecision> CenteredEuler3DTransformType;
1532  typedef typename CenteredEuler3DTransformType::Pointer CenteredEuler3DTransformPointer;
1533  typedef itk::Euler3DTransform<TPrecision> Euler3DTransformType;
1534  typedef typename Euler3DTransformType::Pointer Euler3DTransformPointer;
1535 
1536  AffineTransformPointer affineTransform;
1537  VersorRigidTransformPointer versorRigidTransform;
1538  RigidTransformPointer rigidTransform;
1539  CenteredEuler3DTransformPointer centeredEulerTransform;
1540  Euler3DTransformPointer euler3DTransform;
1541 
1542  if( !strcmp(transf->GetNameOfClass(),"AffineTransform") )
1543  {
1544  AffineTransformPointer affine_read = static_cast<AffineTransformType*>(transf.GetPointer());
1545  affineTransform = dynamic_cast< AffineTransformType * >( affine_read.GetPointer() );
1546  if( affineTransform )
1547  {
1548  PrintInfo("Successfully Read Affine Transform file.");
1549  AffineTransformPointer affineInvTransform;
1550  if(inverse)
1551  {
1552  affineInvTransform = AffineTransformType::New();
1553  affineTransform->GetInverse(affineInvTransform);
1554  }
1555  else
1556  affineInvTransform = affineTransform;
1557  resample->SetTransform(affineInvTransform);
1558  }
1559  }
1560  else if( !strcmp(transf->GetNameOfClass(),"VersorRigid3DTransform") )
1561  {
1562  VersorRigidTransformPointer rigid_read = static_cast<VersorRigidTransformType*>(transf.GetPointer());
1563  versorRigidTransform = dynamic_cast< VersorRigidTransformType * >( rigid_read.GetPointer() );
1564  if( versorRigidTransform )
1565  {
1566  PrintInfo("Successfully Read Versor Rigid Transform file.");
1567  VersorRigidTransformPointer versorRigidInvTransform;
1568  if(inverse)
1569  {
1570  versorRigidInvTransform = VersorRigidTransformType::New();
1571  versorRigidTransform->GetInverse(versorRigidInvTransform);
1572  }
1573  else
1574  versorRigidInvTransform = versorRigidTransform;
1575  resample->SetTransform(versorRigidInvTransform);
1576  }
1577  }
1578  else if( !strcmp(transf->GetNameOfClass(),"Rigid3DTransform") )
1579  {
1580  RigidTransformPointer rigid_read = static_cast<RigidTransformType*>(transf.GetPointer());
1581  rigidTransform = dynamic_cast< RigidTransformType * >( rigid_read.GetPointer() );
1582  if( rigidTransform )
1583  {
1584  PrintInfo("Successfully Read Rigid Transform file.");
1585  RigidTransformPointer rigidInvTransform;
1586  if(inverse)
1587  {
1588  rigidInvTransform = RigidTransformType::New();
1589  rigidTransform->GetInverse(rigidInvTransform);
1590  }
1591  else
1592  rigidInvTransform = rigidTransform;
1593  resample->SetTransform(rigidInvTransform);
1594  }
1595  }
1596  else if( !strcmp(transf->GetNameOfClass(),"CenteredEuler3DTransform") )
1597  {
1598  CenteredEuler3DTransformPointer rigid_read = static_cast<CenteredEuler3DTransformType*>(transf.GetPointer());
1599  centeredEulerTransform = dynamic_cast< CenteredEuler3DTransformType * >( rigid_read.GetPointer() );
1600  if( centeredEulerTransform )
1601  {
1602  PrintInfo("Successfully Read Centered Euler 3D Transform file.");
1603  CenteredEuler3DTransformPointer centeredEulerInvTransform;
1604  if(inverse)
1605  {
1606  centeredEulerInvTransform = CenteredEuler3DTransformType::New();
1607  centeredEulerTransform->GetInverse(centeredEulerInvTransform);
1608  }
1609  else
1610  centeredEulerInvTransform = centeredEulerTransform;
1611  resample->SetTransform(centeredEulerInvTransform);
1612  }
1613  }
1614  else if( !strcmp(transf->GetNameOfClass(),"Euler3DTransform") )
1615  {
1616  Euler3DTransformPointer rigid_read = static_cast<Euler3DTransformType*>(transf.GetPointer());
1617  euler3DTransform = dynamic_cast< Euler3DTransformType * >( rigid_read.GetPointer() );
1618  if( euler3DTransform )
1619  {
1620  PrintInfo("Successfully Read Euler 3D Transform file.");
1621  Euler3DTransformPointer euler3DInvTransform;
1622  if(inverse)
1623  {
1624  euler3DInvTransform = Euler3DTransformType::New();
1625  euler3DTransform->GetInverse(euler3DInvTransform);
1626  }
1627  else
1628  euler3DInvTransform = euler3DTransform;
1629  resample->SetTransform(euler3DInvTransform);
1630  }
1631  }
1632  else
1633  {
1634  PrintError("milxImage: Transform Input is not of known type");
1635  }
1636 
1637  try
1638  {
1639  resample->Update();
1640  }
1641  catch (itk::ExceptionObject & ex )
1642  {
1643  PrintError("Failed Transforming Image");
1644  PrintError(ex.GetDescription());
1645  }
1646 
1647  return resample->GetOutput();
1648 }
1649 
1650 template<class TImage>
1651 template<typename TOutImage, typename TTransform, typename TPrecision>
1652 itk::SmartPointer<TOutImage> Image<TImage>::TransformImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TTransform> transf, const bool inverse, const int interp)
1653 {
1654  typedef itk::InterpolateImageFunction<TImage, TPrecision> InterpolatorType;
1655  typename InterpolatorType::Pointer interpolator;
1656  if(interp == 0)
1657  {
1658  interpolator = itk::NearestNeighborInterpolateImageFunction<TImage, TPrecision>::New();
1659  PrintDebug("Using nearest neighbour interpolation for resampling.");
1660  }
1661  else if(interp == 1)
1662  {
1663  interpolator = itk::LinearInterpolateImageFunction<TImage, TPrecision>::New();
1664  PrintDebug("Using linear interpolation for resampling.");
1665  }
1666  else
1667  {
1668  typedef itk::BSplineInterpolateImageFunction<TImage, TPrecision, TPrecision> BSplineInterpolatorType;
1669  typename BSplineInterpolatorType::Pointer bsplineInterpolator = BSplineInterpolatorType::New();
1670  bsplineInterpolator->SetSplineOrder(3);
1671  interpolator = bsplineInterpolator;
1672  PrintDebug("Using BSpline interpolation for resampling.");
1673  }
1674 
1675  //transform origin
1676  typename TImage::PointType transformedOrigin = transf->TransformPoint(img->GetOrigin());
1677 
1678  typedef itk::ResampleImageFilter<TImage, TOutImage, TPrecision> ResampleImageFilterType;
1679  typename ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
1680 
1682  typedef itk::AffineTransform< TPrecision, milx::imgDimension > AffineTransformType;
1683  typedef typename AffineTransformType::Pointer AffineTransformPointer;
1684  typedef itk::VersorRigid3DTransform<TPrecision> VersorRigidTransformType;
1685  typedef typename VersorRigidTransformType::Pointer VersorRigidTransformPointer;
1686 #if ITK_VERSION_MAJOR > 3 //New() member issues, see ITK4 tests for detail
1687  typedef itkv3::Rigid3DTransform<TPrecision> RigidTransformType;
1688 #else
1689  typedef itk::Rigid3DTransform<TPrecision> RigidTransformType;
1690 #endif
1691  typedef typename RigidTransformType::Pointer RigidTransformPointer;
1692  typedef itk::CenteredEuler3DTransform<TPrecision> CenteredEuler3DTransformType;
1693  typedef typename CenteredEuler3DTransformType::Pointer CenteredEuler3DTransformPointer;
1694  typedef itk::Euler3DTransform<TPrecision> Euler3DTransformType;
1695  typedef typename Euler3DTransformType::Pointer Euler3DTransformPointer;
1696 
1697  AffineTransformPointer affineTransform;
1698  VersorRigidTransformPointer versorRigidTransform;
1699  RigidTransformPointer rigidTransform;
1700  CenteredEuler3DTransformPointer centeredEulerTransform;
1701  Euler3DTransformPointer euler3DTransform;
1702 
1703  if( !strcmp(transf->GetNameOfClass(),"AffineTransform") )
1704  {
1705  AffineTransformPointer affine_read = static_cast<AffineTransformType*>(transf.GetPointer());
1706  affineTransform = dynamic_cast< AffineTransformType * >( affine_read.GetPointer() );
1707  if( affineTransform )
1708  {
1709  PrintInfo("Successfully Read Affine Transform file.");
1710  AffineTransformPointer affineInvTransform;
1711  if(inverse)
1712  {
1713  affineInvTransform = AffineTransformType::New();
1714  affineTransform->GetInverse(affineInvTransform);
1715  }
1716  else
1717  affineInvTransform = affineTransform;
1718  resample->SetTransform(affineInvTransform);
1719  transformedOrigin = affineInvTransform->TransformPoint(img->GetOrigin());
1720  }
1721  }
1722  else if( !strcmp(transf->GetNameOfClass(),"VersorRigid3DTransform") )
1723  {
1724  VersorRigidTransformPointer rigid_read = static_cast<VersorRigidTransformType*>(transf.GetPointer());
1725  versorRigidTransform = dynamic_cast< VersorRigidTransformType * >( rigid_read.GetPointer() );
1726  if( versorRigidTransform )
1727  {
1728  PrintInfo("Successfully Read Versor Rigid Transform file.");
1729  VersorRigidTransformPointer versorRigidInvTransform;
1730  if(inverse)
1731  {
1732  versorRigidInvTransform = VersorRigidTransformType::New();
1733  versorRigidTransform->GetInverse(versorRigidInvTransform);
1734  }
1735  else
1736  versorRigidInvTransform = versorRigidTransform;
1737  resample->SetTransform(versorRigidInvTransform);
1738  transformedOrigin = versorRigidInvTransform->TransformPoint(img->GetOrigin());
1739  }
1740  }
1741  else if( !strcmp(transf->GetNameOfClass(),"Rigid3DTransform") )
1742  {
1743  RigidTransformPointer rigid_read = static_cast<RigidTransformType*>(transf.GetPointer());
1744  rigidTransform = dynamic_cast< RigidTransformType * >( rigid_read.GetPointer() );
1745  if( rigidTransform )
1746  {
1747  PrintInfo("Successfully Read Rigid Transform file.");
1748  RigidTransformPointer rigidInvTransform;
1749  if(inverse)
1750  {
1751  rigidInvTransform = RigidTransformType::New();
1752  rigidTransform->GetInverse(rigidInvTransform);
1753  }
1754  else
1755  rigidInvTransform = rigidTransform;
1756  resample->SetTransform(rigidInvTransform);
1757  transformedOrigin = rigidInvTransform->TransformPoint(img->GetOrigin());
1758  }
1759  }
1760  else if( !strcmp(transf->GetNameOfClass(),"CenteredEuler3DTransform") )
1761  {
1762  CenteredEuler3DTransformPointer rigid_read = static_cast<CenteredEuler3DTransformType*>(transf.GetPointer());
1763  centeredEulerTransform = dynamic_cast< CenteredEuler3DTransformType * >( rigid_read.GetPointer() );
1764  if( centeredEulerTransform )
1765  {
1766  PrintInfo("Successfully Read Centered Euler 3D Transform file.");
1767  CenteredEuler3DTransformPointer centeredEulerInvTransform;
1768  if(inverse)
1769  {
1770  centeredEulerInvTransform = CenteredEuler3DTransformType::New();
1771  centeredEulerTransform->GetInverse(centeredEulerInvTransform);
1772  }
1773  else
1774  centeredEulerInvTransform = centeredEulerTransform;
1775  resample->SetTransform(centeredEulerInvTransform);
1776  transformedOrigin = centeredEulerInvTransform->TransformPoint(img->GetOrigin());
1777  }
1778  }
1779  else if( !strcmp(transf->GetNameOfClass(),"Euler3DTransform") )
1780  {
1781  Euler3DTransformPointer rigid_read = static_cast<Euler3DTransformType*>(transf.GetPointer());
1782  euler3DTransform = dynamic_cast< Euler3DTransformType * >( rigid_read.GetPointer() );
1783  if( euler3DTransform )
1784  {
1785  PrintInfo("Successfully Read Euler 3D Transform file.");
1786  Euler3DTransformPointer euler3DInvTransform;
1787  if(inverse)
1788  {
1789  euler3DInvTransform = Euler3DTransformType::New();
1790  euler3DTransform->GetInverse(euler3DInvTransform);
1791  }
1792  else
1793  euler3DInvTransform = euler3DTransform;
1794  resample->SetTransform(euler3DInvTransform);
1795  transformedOrigin = euler3DInvTransform->TransformPoint(img->GetOrigin());
1796  }
1797  }
1798  else
1799  {
1800  PrintError("milxImage: Transform Input is not of known type");
1801  }
1802 
1803  resample->SetInput(img);
1804  resample->SetDefaultPixelValue(0.0);
1805  resample->SetSize( img->GetLargestPossibleRegion().GetSize() );
1806  resample->SetOutputOrigin( transformedOrigin );
1807  resample->SetOutputSpacing( img->GetSpacing() );
1808  resample->SetOutputDirection( img->GetDirection() );
1809  resample->SetInterpolator(interpolator);
1810  resample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1811  try
1812  {
1813  resample->Update();
1814  }
1815  catch (itk::ExceptionObject & ex )
1816  {
1817  PrintError("Failed Transforming Image");
1818  PrintError(ex.GetDescription());
1819  }
1820 
1821  return resample->GetOutput();
1822 }
1823 
1824 template<class TImage>
1825 template<typename TOutImage>
1826 itk::SmartPointer<TOutImage> Image<TImage>::ResampleImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TOutImage> refImg, const bool linearInterp)
1827 {
1828 // typedef itk::IdentityTransform<float, imgDimension> TransformType;
1829  typedef itk::InterpolateImageFunction<TImage, float> InterpolatorType;
1830  typename InterpolatorType::Pointer interpolator;
1831  if(linearInterp)
1832  {
1833  interpolator = itk::LinearInterpolateImageFunction<TImage, float>::New();
1834  PrintDebug("Using linear interpolation for resampling.");
1835  }
1836  else
1837  {
1838  typedef itk::BSplineInterpolateImageFunction<TImage, float, float> BSplineInterpolatorType;
1839  typename BSplineInterpolatorType::Pointer bsplineInterpolator = BSplineInterpolatorType::New();
1840  bsplineInterpolator->SetSplineOrder(3);
1841  interpolator = bsplineInterpolator;
1842  PrintDebug("Using BSpline interpolation for resampling.");
1843  }
1844 
1845  typedef itk::ResampleImageFilter<TImage, TOutImage, float> ResampleImageFilterType;
1846  typename ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
1847  resample->SetInput(img);
1848  resample->SetDefaultPixelValue(0.0);
1849  resample->UseReferenceImageOn();
1850  resample->SetReferenceImage(refImg);
1851 // resample->SetOutputSpacing(refImg->GetSpacing());
1852 // resample->SetOutputOrigin(refImg->GetOrigin());
1853 // resample->SetOutputDirection(refImg->GetDirection());
1854 // resample->SetTransform(TransformType::New());
1855  resample->SetInterpolator(interpolator);
1856  resample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1857 // if(matchSize)
1858 // {
1859 // PrintDebug("Resampling with size changes.");
1860 // resample->SetSize(refImg->GetLargestPossibleRegion().GetSize());
1861 // resample->SetOutputStartIndex(refImg->GetLargestPossibleRegion().GetIndex());
1862 // }
1863 
1864  try
1865  {
1866  resample->UpdateLargestPossibleRegion();
1867  resample->Update();
1868  }
1869  catch (itk::ExceptionObject & ex )
1870  {
1871  PrintError("Failed Resampling");
1872  PrintError(ex.GetDescription());
1873  }
1874 
1875  return resample->GetOutput();
1876 }
1877 
1878 template<class TImage>
1879 template<typename TOutImage>
1880 itk::SmartPointer<TOutImage> Image<TImage>::ResampleLabel(itk::SmartPointer<TImage> img, itk::SmartPointer<TOutImage> refImg)
1881 {
1882 // typedef itk::IdentityTransform<float, imgDimension> TransformType;
1883  typedef itk::NearestNeighborInterpolateImageFunction<TImage, float> NNInterpolatorType;
1884  typename NNInterpolatorType::Pointer interpolator = NNInterpolatorType::New();
1885  PrintDebug("Using Nearest Neighbour interpolation for resampling.");
1886 
1887  typedef itk::ResampleImageFilter<TImage, TOutImage, float> ResampleImageFilterType;
1888  typename ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();
1889  resample->SetInput(img);
1890  resample->SetDefaultPixelValue(0.0);
1891  resample->UseReferenceImageOn();
1892  resample->SetReferenceImage(refImg);
1893 // resample->SetOutputSpacing(refImg->GetSpacing());
1894 // resample->SetOutputOrigin(refImg->GetOrigin());
1895 // resample->SetOutputDirection(refImg->GetDirection());
1896 // resample->SetTransform(TransformType::New());
1897  resample->SetInterpolator(interpolator);
1898  resample->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1899 // if(matchSize)
1900 // {
1901 // PrintDebug("Resampling with size changes.");
1902 // resample->SetSize(refImg->GetLargestPossibleRegion().GetSize());
1903 // resample->SetOutputStartIndex(refImg->GetLargestPossibleRegion().GetIndex());
1904 // }
1905 
1906  try
1907  {
1908  resample->UpdateLargestPossibleRegion();
1909  resample->Update();
1910  }
1911  catch (itk::ExceptionObject & ex )
1912  {
1913  PrintError("Failed Resampling");
1914  PrintError(ex.GetDescription());
1915  }
1916 
1917  return resample->GetOutput();
1918 }
1919 
1920 //Arithmetic
1921 template<class TImage>
1922 itk::SmartPointer<TImage> Image<TImage>::AddImages(itk::SmartPointer<TImage> img1, itk::SmartPointer<TImage> img2)
1923 {
1924  typedef itk::AddImageFilter<TImage, TImage> AddImageType;
1925 
1926  typename AddImageType::Pointer addFilter = AddImageType::New();
1927  addFilter->SetInput1(img1);
1928  addFilter->SetInput2(img2);
1929  addFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1930  try
1931  {
1932  addFilter->Update();
1933  }
1934  catch (itk::ExceptionObject & ex )
1935  {
1936  PrintError("Failed Addition");
1937  PrintError(ex.GetDescription());
1938  }
1939 
1940  return addFilter->GetOutput();
1941 }
1942 
1943 template<class TImage>
1944 itk::SmartPointer<TImage> Image<TImage>::SubtractImages(itk::SmartPointer<TImage> img1, itk::SmartPointer<TImage> img2)
1945 {
1946  typedef itk::SubtractImageFilter<TImage, TImage> SubtractImageType;
1947 
1948  typename SubtractImageType::Pointer subtractFilter = SubtractImageType::New();
1949  subtractFilter->SetInput1(img1);
1950  subtractFilter->SetInput2(img2);
1951  subtractFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1952  try
1953  {
1954  subtractFilter->Update();
1955  }
1956  catch (itk::ExceptionObject & ex )
1957  {
1958  PrintError("Failed Subtraction");
1959  PrintError(ex.GetDescription());
1960  }
1961 
1962  return subtractFilter->GetOutput();
1963 }
1964 
1965 template<class TImage>
1966 itk::SmartPointer<TImage> Image<TImage>::MultiplyImages(itk::SmartPointer<TImage> img1, itk::SmartPointer<TImage> img2)
1967 {
1968  typedef itk::MultiplyImageFilter<TImage, TImage> MultiplyImageType;
1969 
1970  typename MultiplyImageType::Pointer multiplyFilter = MultiplyImageType::New();
1971  multiplyFilter->SetInput1(img1);
1972  multiplyFilter->SetInput2(img2);
1973  multiplyFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1974  try
1975  {
1976  multiplyFilter->Update();
1977  }
1978  catch (itk::ExceptionObject & ex)
1979  {
1980  PrintError("Failed Multiplication");
1981  PrintError(ex.GetDescription());
1982  }
1983 
1984  return multiplyFilter->GetOutput();
1985 }
1986 
1987 template<class TImage>
1988 template<class TOutImage>
1989 itk::SmartPointer<TOutImage> Image<TImage>::ScaleImage(itk::SmartPointer<TImage> img, float scaling)
1990 {
1991  typedef itk::ShiftScaleImageFilter<TImage, TOutImage> ScaleImageType;
1992 
1993  typename ScaleImageType::Pointer scaleFilter = ScaleImageType::New();
1994  scaleFilter->SetInput(img);
1995  scaleFilter->SetScale(scaling);
1996  scaleFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
1997  try
1998  {
1999  scaleFilter->Update();
2000  }
2001  catch (itk::ExceptionObject & ex )
2002  {
2003  PrintError("Failed Scaling");
2004  PrintError(ex.GetDescription());
2005  }
2006 
2007  return scaleFilter->GetOutput();
2008 }
2009 
2010 template<class TImage>
2011 itk::SmartPointer<TImage> Image<TImage>::ScaleVectorImage(itk::SmartPointer<TImage> img, float scaling, int numberOfComponents)
2012 {
2013  typedef itk::VectorShiftScaleImageFilter<TImage, TImage> ScaleImageType;
2014 
2015  typename ScaleImageType::Pointer scaleFilter = ScaleImageType::New();
2016  scaleFilter->SetInput(img);
2017  scaleFilter->SetScale(scaling);
2018  scaleFilter->SetDimension(numberOfComponents);
2019  scaleFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2020  try
2021  {
2022  scaleFilter->Update();
2023  }
2024  catch (itk::ExceptionObject & ex )
2025  {
2026  PrintError("Failed Scaling");
2027  PrintError(ex.GetDescription());
2028  }
2029 
2030  return scaleFilter->GetOutput();
2031 }
2032 
2033 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
2034 template<class TImage>
2035 itk::SmartPointer<TImage> Image<TImage>::ConvolveImages(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> kernelImg)
2036 {
2037 #if ITK_VERSION_MAJOR > 3
2038  typedef itk::FFTConvolutionImageFilter<TImage, TImage> ConvolveImageType;
2039 #else
2040  typedef itk::ConvolutionImageFilter<TImage, TImage> ConvolveImageType;
2041 #endif
2042 
2043  typename ConvolveImageType::Pointer convolveFilter = ConvolveImageType::New();
2044  convolveFilter->SetInput(img);
2045 #if ITK_VERSION_MAJOR > 3
2046  convolveFilter->SetKernelImage(kernelImg);
2047 #else
2048  convolveFilter->SetImageKernelInput(kernelImg);
2049 #endif
2050  convolveFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2051  try
2052  {
2053  convolveFilter->Update();
2054  }
2055  catch (itk::ExceptionObject & ex )
2056  {
2057  PrintError("Failed Convolution");
2058  PrintError(ex.GetDescription());
2059  }
2060 
2061  return convolveFilter->GetOutput();
2062 }
2063 #endif //(ITK_REVIEW || ITK_VERSION_MAJOR > 3)
2064 
2065 //Info
2066 template<class TImage>
2067 void Image<TImage>::Information(itk::SmartPointer<TImage> img)
2068 {
2069  typename TImage::DirectionType direction = img->GetDirection();
2070  typename TImage::PointType origin = img->GetOrigin();
2071  typename TImage::SpacingType spacing = img->GetSpacing();
2072  typename TImage::SizeType imageSize = img->GetLargestPossibleRegion().GetSize();
2073 
2074  std::cout << "Origin: ";
2075  for (typename TImage::SizeValueType j = 0; j < TImage::ImageDimension; j++)
2076  std::cout << milx::NumberToString(origin[j]) << " ";
2077  std::cout << "\nSpacing: ";
2078  for (typename TImage::SizeValueType j = 0; j < TImage::ImageDimension; j++)
2079  std::cout << milx::NumberToString(spacing[j]) << " ";
2080  std::cout << "\nSize: ";
2081  for (typename TImage::SizeValueType j = 0; j < TImage::ImageDimension; j++)
2082  std::cout << milx::NumberToString(imageSize[j]) << " ";
2083  std::cout << "\nReal Size: ";
2084  for (typename TImage::SizeValueType j = 0; j < TImage::ImageDimension; j++)
2085  std::cout << milx::NumberToString(spacing[j] * imageSize[j]) << " ";
2086  std::cout << "\nDirection/Orientation: " << std::endl;
2087  for (typename TImage::SizeValueType j = 0; j < TImage::ImageDimension; j++)
2088  {
2089  std::cout << "| ";
2090  for (typename TImage::SizeValueType k = 0; k < TImage::ImageDimension; k++)
2091  {
2092  std::cout << milx::NumberToString(direction(j, k));
2093  if (k < TImage::ImageDimension - 1)
2094  std::cout << ",\t";
2095  }
2096  std::cout << " |" << std::endl;
2097  }
2098  std::cout << "\nOrientation Flag: " << ImageOrientation(img) << std::endl;
2099 }
2100 
2101 template<class TImage>
2102 double Image<TImage>::ImageMaximum(itk::SmartPointer<TImage> img)
2103 {
2104  typedef itk::MinimumMaximumImageCalculator<TImage> ImageCalculatorFilterType;
2105  typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New();
2106  imageCalculatorFilter->SetImage(img);
2107  imageCalculatorFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2108  imageCalculatorFilter->Compute();
2109 
2110  return imageCalculatorFilter->GetMaximum();
2111 }
2112 
2113 template<class TImage>
2114 double Image<TImage>::ImageMinimum(itk::SmartPointer<TImage> img)
2115 {
2116  typedef itk::MinimumMaximumImageCalculator<TImage> ImageCalculatorFilterType;
2117  typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New();
2118  imageCalculatorFilter->SetImage(img);
2119  imageCalculatorFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2120  imageCalculatorFilter->Compute();
2121 
2122  return imageCalculatorFilter->GetMinimum();
2123 }
2124 
2125 template<class TImage>
2126 std::string Image<TImage>::ImageOrientation(itk::SmartPointer<TImage> img)
2127 {
2128  std::map<itk::SpatialOrientation::ValidCoordinateOrientationFlags, std::string> codeToString;
2129 
2130  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIL] = "AIL";
2131  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ASL] = "ASL";
2132  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAI] = "RAI";
2133  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LAI] = "LAI";
2134  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPS] = "RPS";
2135  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPS] = "LPS";
2136  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIP] = "RIP";
2137  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LIP] = "LIP";
2138  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSP] = "RSP";
2139  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LSP] = "LSP";
2140  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIA] = "RIA";
2141  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LIA] = "LIA";
2142  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RSA] = "RSA";
2143  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LSA] = "LSA";
2144  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IRP] = "IRP";
2145  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ILP] = "ILP";
2146  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SRP] = "SRP";
2147  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SLP] = "SLP";
2148  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IRA] = "IRA";
2149  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ILA] = "ILA";
2150  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SRA] = "SRA";
2151  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SLA] = "SLA";
2152  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPI] = "RPI";
2153  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPI] = "LPI";
2154  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RAS] = "RAS";
2155  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LAS] = "LAS";
2156  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PRI] = "PRI";
2157  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PLI] = "PLI";
2158  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ARI] = "ARI";
2159  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ALI] = "ALI";
2160  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PRS] = "PRS";
2161  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PLS] = "PLS";
2162  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ARS] = "ARS";
2163  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ALS] = "ALS";
2164  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IPR] = "IPR";
2165  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SPR] = "SPR";
2166  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IAR] = "IAR";
2167  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SAR] = "SAR";
2168  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IPL] = "IPL";
2169  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SPL] = "SPL";
2170  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_IAL] = "IAL";
2171  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_SAL] = "SAL";
2172  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PIR] = "PIR";
2173  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PSR] = "PSR";
2174  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_AIR] = "AIR";
2175  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_ASR] = "ASR";
2176  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PIL] = "PIL";
2177  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_PSL] = "PSL";
2178  codeToString[itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_INVALID] = "Unknown";
2179 
2180  itk::SpatialOrientation::ValidCoordinateOrientationFlags orientFlag = itk::SpatialOrientationAdapter().FromDirectionCosines(img->GetDirection());
2181  std::string orientFlagStr = codeToString[orientFlag];
2182 
2183  return orientFlagStr;
2184 }
2185 
2186 //Filters
2187 template<class TImage>
2188 itk::SmartPointer<TImage> Image<TImage>::CheckerBoard(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToCheckerBoard, const int numberOfSquares)
2189 {
2190  typedef itk::CheckerBoardImageFilter<TImage> CheckerBoardType;
2191  typename CheckerBoardType::Pointer checkerBoardFilter = CheckerBoardType::New();
2192  checkerBoardFilter->SetInput1(img);
2193  checkerBoardFilter->SetInput2(imgToCheckerBoard);
2194  if(numberOfSquares != 0)
2195  {
2196  typename CheckerBoardType::PatternArrayType squares = checkerBoardFilter->GetCheckerPattern();
2197  squares.Fill(numberOfSquares);
2198  checkerBoardFilter->SetCheckerPattern(squares);
2199  }
2200  checkerBoardFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2201  try
2202  {
2203  checkerBoardFilter->Update();
2204  }
2205  catch (itk::ExceptionObject & ex )
2206  {
2207  PrintError("Failed Generating Checker Board");
2208  PrintError(ex.GetDescription());
2209  }
2210 
2211  return checkerBoardFilter->GetOutput();
2212 }
2213 
2214 template<class TImage>
2215 itk::SmartPointer<TImage> Image<TImage>::RescaleIntensity(itk::SmartPointer<TImage> img, float minValue, float maxValue)
2216 {
2217  typedef itk::RescaleIntensityImageFilter< TImage, TImage > RescaleFilterType;
2218  typename RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
2219  rescaleFilter->SetInput(img);
2220  rescaleFilter->SetOutputMinimum(minValue);
2221  rescaleFilter->SetOutputMaximum(maxValue);
2222  rescaleFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2223  try
2224  {
2225  rescaleFilter->Update();
2226  }
2227  catch (itk::ExceptionObject & ex )
2228  {
2229  PrintError("Failed Rescaling Intensities");
2230  PrintError(ex.GetDescription());
2231  }
2232 
2233  return rescaleFilter->GetOutput();
2234 }
2235 
2236 template<class TImage>
2237 itk::SmartPointer<TImage> Image<TImage>::InvertIntensity(itk::SmartPointer<TImage> img, float maxValue)
2238 {
2239  typedef itk::InvertIntensityImageFilter<TImage, TImage> InvertIntensityImageFilterType;
2240  typename InvertIntensityImageFilterType::Pointer invert = InvertIntensityImageFilterType::New();
2241  invert->SetInput(img);
2242  invert->SetMaximum(maxValue);
2243  invert->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2244  try
2245  {
2246  invert->Update();
2247  }
2248  catch (itk::ExceptionObject & ex )
2249  {
2250  PrintError("Failed Computing Invert Intensity");
2251  PrintError(ex.GetDescription());
2252  }
2253 
2254  return invert->GetOutput();
2255 }
2256 
2257 template<class TImage>
2258 itk::SmartPointer<TImage> Image<TImage>::HistogramEqualisation(itk::SmartPointer<TImage> img, float alpha, float beta, float radius)
2259 {
2260  typename TImage::SizeType radii;
2261  radii.Fill(radius);
2262  typedef itk::AdaptiveHistogramEqualizationImageFilter<TImage> EqualiseFilterType;
2263  typename EqualiseFilterType::Pointer equaliseFilter = EqualiseFilterType::New();
2264  equaliseFilter->SetInput(img);
2265  equaliseFilter->SetAlpha(alpha);
2266  equaliseFilter->SetBeta(beta);
2267  equaliseFilter->SetRadius(radii);
2268  equaliseFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2269  try
2270  {
2271  equaliseFilter->Update();
2272  }
2273  catch (itk::ExceptionObject & ex )
2274  {
2275  PrintError("Failed Histogram Equalisation");
2276  PrintError(ex.GetDescription());
2277  }
2278 
2279  return equaliseFilter->GetOutput();
2280 }
2281 
2282 template<class TImage>
2283 itk::SmartPointer<TImage> Image<TImage>::GradientMagnitude(itk::SmartPointer<TImage> img)
2284 {
2285  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<TImage, TImage> GradImageFilterType;
2286  typename GradImageFilterType::Pointer gradientMagnitudeFilter = GradImageFilterType::New();
2287  gradientMagnitudeFilter->SetInput(img);
2288  gradientMagnitudeFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2289  try
2290  {
2291  gradientMagnitudeFilter->Update();
2292  }
2293  catch (itk::ExceptionObject & ex )
2294  {
2295  PrintError("Failed Computing Gradient Magnitude");
2296  PrintError(ex.GetDescription());
2297  }
2298 
2299  return gradientMagnitudeFilter->GetOutput();
2300 }
2301 
2302 template<class TImage>
2303 itk::SmartPointer<TImage> Image<TImage>::SobelEdges(itk::SmartPointer<TImage> img)
2304 {
2305  typedef itk::SobelEdgeDetectionImageFilter<TImage, TImage> SobelImageFilterType;
2306  typename SobelImageFilterType::Pointer sobelEdgeFilter = SobelImageFilterType::New();
2307  sobelEdgeFilter->SetInput(img);
2308  sobelEdgeFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2309  try
2310  {
2311  sobelEdgeFilter->Update();
2312  }
2313  catch (itk::ExceptionObject & ex )
2314  {
2315  PrintError("Failed Computing Sobel Edges");
2316  PrintError(ex.GetDescription());
2317  }
2318 
2319  return sobelEdgeFilter->GetOutput();
2320 }
2321 
2322 template<class TImage>
2323 itk::SmartPointer<TImage> Image<TImage>::Laplacian(itk::SmartPointer<TImage> img)
2324 {
2325  typedef itk::LaplacianRecursiveGaussianImageFilter<TImage, TImage> LaplacianImageFilterType;
2326  typename LaplacianImageFilterType::Pointer laplacianEdgeFilter = LaplacianImageFilterType::New();
2327  laplacianEdgeFilter->SetInput(img);
2328  laplacianEdgeFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2329  try
2330  {
2331  laplacianEdgeFilter->Update();
2332  }
2333  catch (itk::ExceptionObject & ex )
2334  {
2335  PrintError("Failed Computing Laplacian");
2336  PrintError(ex.GetDescription());
2337  }
2338 
2339  return laplacianEdgeFilter->GetOutput();
2340 }
2341 
2342 template<class TImage>
2343 itk::SmartPointer<TImage> Image<TImage>::CannyEdges(itk::SmartPointer<TImage> img, float variance, float lowerThreshold, float upperThreshold)
2344 {
2345  typedef itk::CannyEdgeDetectionImageFilter<TImage, TImage> CannyImageFilterType;
2346  typename CannyImageFilterType::Pointer cannyEdgeFilter = CannyImageFilterType::New();
2347  cannyEdgeFilter->SetInput(img);
2348  cannyEdgeFilter->SetVariance( variance );
2349  cannyEdgeFilter->SetUpperThreshold( upperThreshold );
2350  cannyEdgeFilter->SetLowerThreshold( lowerThreshold );
2351  cannyEdgeFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2352  try
2353  {
2354  cannyEdgeFilter->Update();
2355  }
2356  catch (itk::ExceptionObject & ex )
2357  {
2358  PrintError("Failed Computing Canny Edges");
2359  PrintError(ex.GetDescription());
2360  }
2361 
2362  return cannyEdgeFilter->GetOutput();
2363 }
2364 
2365 template<class TImage>
2366 itk::SmartPointer< itk::Image<float, TImage::ImageDimension> > Image<TImage>::Normalization(itk::SmartPointer<TImage> img)
2367 {
2368  typedef itk::NormalizeImageFilter< TImage, itk::Image<float, TImage::ImageDimension> > NormalizeImageFilterType;
2369  typename NormalizeImageFilterType::Pointer normalizeFilter = NormalizeImageFilterType::New();
2370  normalizeFilter->SetInput(img);
2371  normalizeFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2372  try
2373  {
2374  normalizeFilter->Update();
2375  }
2376  catch (itk::ExceptionObject & ex )
2377  {
2378  PrintError("Failed Computing Normalization");
2379  PrintError(ex.GetDescription());
2380  }
2381 
2382  return normalizeFilter->GetOutput();
2383 }
2384 
2385 template<class TImage>
2386 itk::SmartPointer<TImage> Image<TImage>::MatchInformation(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToMatch, bool originOnly)
2387 {
2388  typedef itk::ChangeInformationImageFilter<TImage> ChangeInfoImageFilterType;
2389  typename ChangeInfoImageFilterType::Pointer change = ChangeInfoImageFilterType::New();
2390  change->SetInput(img);
2391  change->SetReferenceImage(imgToMatch);
2392  change->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2393  change->UseReferenceImageOn();
2394  if(!originOnly)
2395  {
2396  change->ChangeDirectionOn();
2397  change->ChangeSpacingOn();
2398  change->ChangeRegionOn();
2399  }
2400  else
2401  {
2402  change->ChangeDirectionOff();
2403  change->ChangeSpacingOff();
2404  change->ChangeRegionOff();
2405  }
2406  change->ChangeOriginOn();
2407  try
2408  {
2409  change->Update();
2410  }
2411  catch (itk::ExceptionObject & ex )
2412  {
2413  PrintError("Failed Matching Info");
2414  PrintError(ex.GetDescription());
2415  }
2416 
2417  return change->GetOutput();
2418 }
2419 
2420 template<class TImage>
2421 itk::SmartPointer<TImage> Image<TImage>::MatchHistogram(itk::SmartPointer<TImage> img, itk::SmartPointer<TImage> imgToMatch, const int bins)
2422 {
2423  typedef itk::HistogramMatchingImageFilter <TImage, TImage> MatchingFilterType;
2424  typename MatchingFilterType::Pointer matcher = MatchingFilterType::New();
2425  matcher->SetInput(img);
2426  matcher->SetReferenceImage(imgToMatch);
2427  matcher->SetNumberOfHistogramLevels(bins);
2428  matcher->SetNumberOfMatchPoints(7);
2429  matcher->ThresholdAtMeanIntensityOn();
2430  matcher->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2431  try
2432  {
2433  matcher->Update();
2434  }
2435  catch (itk::ExceptionObject & ex )
2436  {
2437  PrintError("Failed Matching Histogram");
2438  PrintError(ex.GetDescription());
2439  }
2440 
2441  return matcher->GetOutput();
2442 }
2443 
2444 template<class TImage>
2445 template<typename TOutImage>
2446 itk::SmartPointer<TOutImage> Image<TImage>::DistanceMap(itk::SmartPointer<TImage> img, const bool binaryImage, const bool signedDistances, const bool computeInsideObject, const bool squaredDistance)
2447 {
2448  typedef itk::ImageToImageFilter<TImage, TOutImage> ImageFilterType;
2449 
2450  typename ImageFilterType::Pointer filter;
2451  if(computeInsideObject)
2452  {
2453  typedef itk::ApproximateSignedDistanceMapImageFilter<TImage, TOutImage> DistanceMapImageFilterType;
2454  typename DistanceMapImageFilterType::Pointer distanceMapFilter = DistanceMapImageFilterType::New();
2455  filter = distanceMapFilter;
2456  }
2457  else if(signedDistances)
2458  {
2459  typedef itk::SignedMaurerDistanceMapImageFilter<TImage, TOutImage> DistanceMapImageFilterType;
2460  typename DistanceMapImageFilterType::Pointer distanceMapFilter = DistanceMapImageFilterType::New();
2461  distanceMapFilter->UseImageSpacingOn();
2462  if(!squaredDistance)
2463  distanceMapFilter->SquaredDistanceOff();
2464  else
2465  distanceMapFilter->SquaredDistanceOn();
2466  filter = distanceMapFilter;
2467  }
2468  else
2469  {
2470  typedef itk::DanielssonDistanceMapImageFilter<TImage, TOutImage> DistanceMapImageFilterType;
2471  typename DistanceMapImageFilterType::Pointer distanceMapFilter = DistanceMapImageFilterType::New();
2472  distanceMapFilter->UseImageSpacingOn();
2473  if(!binaryImage)
2474  distanceMapFilter->InputIsBinaryOff();
2475  else
2476  distanceMapFilter->InputIsBinaryOn();
2477  if(!squaredDistance)
2478  distanceMapFilter->SquaredDistanceOff();
2479  else
2480  distanceMapFilter->SquaredDistanceOn();
2481  filter = distanceMapFilter;
2482  }
2483 
2484 // std::cout << filter->GetNameOfClass() << std::endl;
2485  filter->SetInput(img);
2486  filter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2487  try
2488  {
2489  filter->Update();
2490  }
2491  catch (itk::ExceptionObject & ex )
2492  {
2493  PrintError("Failed Computing Distance Map");
2494  PrintError(ex.GetDescription());
2495  }
2496 
2497  return filter->GetOutput();
2498 }
2499 
2500 template<class TImage>
2501 itk::SmartPointer<TImage> Image<TImage>::FlipImage(itk::SmartPointer<TImage> img, bool xAxis, bool yAxis, bool zAxis, bool aboutOrigin)
2502 {
2503  itk::FixedArray<bool, imgDimension> flipAxes;
2504  flipAxes[0] = xAxis;
2505  flipAxes[1] = yAxis;
2506  flipAxes[2] = zAxis;
2507 // flipAxes[3] = false;
2508 
2509  typedef itk::FlipImageFilter<TImage> FlipImageFilterType;
2510  typename FlipImageFilterType::Pointer flipFilter = FlipImageFilterType::New();
2511  flipFilter->SetInput(img);
2512  flipFilter->SetFlipAxes(flipAxes);
2513  if(aboutOrigin)
2514  flipFilter->FlipAboutOriginOn();
2515  else
2516  flipFilter->FlipAboutOriginOff();
2517  flipFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2518  try
2519  {
2520  flipFilter->Update();
2521  }
2522  catch (itk::ExceptionObject & ex )
2523  {
2524  PrintError("Failed Computing Flip Image");
2525  PrintError(ex.GetDescription());
2526  }
2527 
2528  return flipFilter->GetOutput();
2529 }
2530 
2531 template<class TImage>
2532 itk::SmartPointer<TImage> Image<TImage>::PadImageByConstant(itk::SmartPointer<TImage> img, size_t xAxis, size_t yAxis, size_t zAxis, typename TImage::PixelType value)
2533 {
2534  typename TImage::SizeType lowerExtendRegion;
2535  lowerExtendRegion[0] = xAxis;
2536  lowerExtendRegion[1] = yAxis;
2537  lowerExtendRegion[2] = zAxis;
2538 
2539  typename TImage::SizeType upperExtendRegion;
2540  upperExtendRegion[0] = xAxis;
2541  upperExtendRegion[1] = yAxis;
2542  upperExtendRegion[2] = zAxis;
2543 
2544  typename TImage::PixelType constantPixel = value;
2545 
2546  typedef itk::ConstantPadImageFilter<TImage, TImage> ConstantPadImageFilterType;
2547  typename ConstantPadImageFilterType::Pointer padFilter = ConstantPadImageFilterType::New();
2548  padFilter->SetInput(img);
2549  padFilter->SetPadLowerBound(lowerExtendRegion);
2550  padFilter->SetPadUpperBound(upperExtendRegion);
2551  padFilter->SetConstant(constantPixel);
2552  padFilter->Update();
2553  padFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2554  try
2555  {
2556  padFilter->Update();
2557  }
2558  catch (itk::ExceptionObject & ex )
2559  {
2560  PrintError("Failed Padding Image");
2561  PrintError(ex.GetDescription());
2562  }
2563 
2564  return padFilter->GetOutput();
2565 }
2566 
2567 template<class TImage>
2568 template<typename TOutImage>
2569 itk::SmartPointer<TOutImage> Image<TImage>::AnisotropicDiffusion(itk::SmartPointer<TImage> img, const int iterations, const float timestep)
2570 {
2571  typedef itk::GradientAnisotropicDiffusionImageFilter<TImage, TOutImage> GradientAnisotropicDiffusionImageFilterType;
2572  typename GradientAnisotropicDiffusionImageFilterType::Pointer gradAniDiff = GradientAnisotropicDiffusionImageFilterType::New();
2573  gradAniDiff->SetInput(img);
2574  gradAniDiff->SetUseImageSpacing(true);
2575  gradAniDiff->SetTimeStep(timestep);
2576  PrintInfo("Time step " + milx::NumberToString(gradAniDiff->GetTimeStep()));
2577  gradAniDiff->SetNumberOfIterations(iterations);
2578  gradAniDiff->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2579  PrintInfo("Conductance " + milx::NumberToString(gradAniDiff->GetConductanceParameter()));
2580  try
2581  {
2582  gradAniDiff->Update();
2583  }
2584  catch (itk::ExceptionObject & ex )
2585  {
2586  PrintError("Failed Computing Invert Intensity");
2587  PrintError(ex.GetDescription());
2588  }
2589 
2590  return gradAniDiff->GetOutput();
2591 }
2592 
2593 template<class TImage>
2594 itk::SmartPointer<TImage> Image<TImage>::Bilateral(itk::SmartPointer<TImage> img, const float sigmaRange, const float sigmaSpatial)
2595 {
2596  typedef itk::BilateralImageFilter<TImage, TImage> BilateralImageFilterType;
2597  typename BilateralImageFilterType::Pointer bilateralFilter = BilateralImageFilterType::New();
2598  bilateralFilter->SetInput(img);
2599  bilateralFilter->SetRangeSigma(sigmaRange);
2600  bilateralFilter->SetDomainSigma(sigmaSpatial);
2601  bilateralFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2602  try
2603  {
2604  bilateralFilter->Update();
2605  }
2606  catch (itk::ExceptionObject & ex )
2607  {
2608  PrintError("Failed Computing Bilateral Smoothing");
2609  PrintError(ex.GetDescription());
2610  }
2611 
2612  return bilateralFilter->GetOutput();
2613 }
2614 
2615 template<class TImage>
2616 itk::SmartPointer<TImage> Image<TImage>::GaussianSmooth(itk::SmartPointer<TImage> img, const float variance)
2617 {
2618  typedef itk::SmoothingRecursiveGaussianImageFilter<TImage, TImage> GuassianImageFilterType;
2619  typename GuassianImageFilterType::Pointer smoothGaussian = GuassianImageFilterType::New();
2620  smoothGaussian->SetInput(img);
2621  smoothGaussian->SetSigma(variance);
2622  smoothGaussian->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2623  try
2624  {
2625  smoothGaussian->Update();
2626  }
2627  catch (itk::ExceptionObject & ex )
2628  {
2629  PrintError("Failed Computing Gaussian Smoothing");
2630  PrintError(ex.GetDescription());
2631  }
2632 
2633  return smoothGaussian->GetOutput();
2634 }
2635 
2636 template<class TImage>
2637 itk::SmartPointer<TImage> Image<TImage>::Median(itk::SmartPointer<TImage> img, const int radius)
2638 {
2639  typedef itk::MedianImageFilter<TImage, TImage> MedianImageFilterType;
2640  typename MedianImageFilterType::Pointer medianFilter = MedianImageFilterType::New();
2641  medianFilter->SetInput(img);
2642  typename MedianImageFilterType::InputSizeType radiusInput;
2643  radiusInput.Fill(radius);
2644  medianFilter->SetRadius(radiusInput);
2645  medianFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2646  try
2647  {
2648  medianFilter->Update();
2649  }
2650  catch (itk::ExceptionObject & ex )
2651  {
2652  PrintError("Failed Computing Median Image");
2653  PrintError(ex.GetDescription());
2654  }
2655 
2656  return medianFilter->GetOutput();
2657 }
2658 
2659 //Labelling
2660 template<class TImage>
2661 template<typename TMaskImage>
2662 itk::SmartPointer<TImage> Image<TImage>::MaskImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TMaskImage> maskImg)
2663 {
2664  typedef itk::MaskImageFilter<TImage, TMaskImage> MaskFilterType;
2665  typename MaskFilterType::Pointer maskFilter = MaskFilterType::New();
2666  maskFilter->SetInput(img);
2667  maskFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2668 #if ITK_VERSION_MAJOR > 3
2669  maskFilter->SetCoordinateTolerance(CoordinateTolerance);
2670  maskFilter->SetDirectionTolerance(DirectionTolerance);
2671  maskFilter->SetMaskImage(maskImg);
2672 #else
2673  maskFilter->SetInput2(maskImg);
2674 #endif
2675 
2676  try
2677  {
2678  maskFilter->Update();
2679  }
2680  catch (itk::ExceptionObject & ex )
2681  {
2682  PrintError("Failed Masking Image");
2683  PrintError(ex.GetDescription());
2684  }
2685 
2686  return maskFilter->GetOutput();
2687 }
2688 
2689 template<class TImage>
2690 itk::SmartPointer<TImage> Image<TImage>::RelabelImage(itk::SmartPointer<TImage> labelledImg)
2691 {
2692  typedef itk::ConnectedComponentImageFilter<TImage, TImage> ConnectedComponentImageFilterType;
2693  typedef itk::RelabelComponentImageFilter<TImage, TImage> RelabelImageType;
2694 
2695  typename ConnectedComponentImageFilterType::Pointer connected = ConnectedComponentImageFilterType::New ();
2696  connected->SetInput(labelledImg);
2697 
2698  typename RelabelImageType::Pointer relabelImageFilter = RelabelImageType::New();
2699  relabelImageFilter->SetInput(connected->GetOutput());
2700 
2701  try
2702  {
2703  relabelImageFilter->Update();
2704  }
2705  catch (itk::ExceptionObject & ex )
2706  {
2707  PrintError("Failed Relabelling Image");
2708  PrintError(ex.GetDescription());
2709  }
2710 
2711  return relabelImageFilter->GetOutput();
2712 }
2713 
2714 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
2715 template<class TImage>
2716 itk::SmartPointer<TImage> Image<TImage>::BinaryContour(itk::SmartPointer<TImage> img, const float foregroundValue, const float backgroundValue)
2717 {
2718  typedef itk::BinaryContourImageFilter<TImage, TImage> BinaryContourType;
2719 
2720  typename BinaryContourType::Pointer contourFilter = BinaryContourType::New();
2721  contourFilter->SetInput(img);
2722  contourFilter->SetForegroundValue(foregroundValue);
2723  contourFilter->SetBackgroundValue(backgroundValue);
2724  contourFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2725  try
2726  {
2727  contourFilter->Update();
2728  }
2729  catch (itk::ExceptionObject & ex )
2730  {
2731  PrintError("Failed Generating Contour");
2732  PrintError(ex.GetDescription());
2733  }
2734 
2735  return contourFilter->GetOutput();
2736 }
2737 
2738 template<class TImage>
2739 itk::SmartPointer<TImage> Image<TImage>::LabelContour(itk::SmartPointer<TImage> img, const bool fullyConnected, const float backgroundValue)
2740 {
2741  typedef itk::LabelContourImageFilter<TImage, TImage> LabelContourType;
2742 
2743  typename LabelContourType::Pointer contourFilter = LabelContourType::New();
2744  contourFilter->SetInput(img);
2745  contourFilter->SetFullyConnected(fullyConnected);
2746  contourFilter->SetBackgroundValue(backgroundValue);
2747  contourFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2748  try
2749  {
2750  contourFilter->Update();
2751  }
2752  catch (itk::ExceptionObject & ex )
2753  {
2754  PrintError("Failed Generating Label Contours");
2755  PrintError(ex.GetDescription());
2756  }
2757 
2758  return contourFilter->GetOutput();
2759 }
2760 
2761 template<class TImage>
2762 template<typename TMaskImage>
2763 itk::SmartPointer<TImage> Image<TImage>::MaskAndCropImage(itk::SmartPointer<TImage> img, itk::SmartPointer<TMaskImage> maskImg, const size_t pixelPadding)
2764 {
2765  typedef itk::LabelObject<unsigned char, TImage::ImageDimension> LabelObjectType;
2766  typedef itk::LabelMap<LabelObjectType> LabelMapType;
2767 
2768  typedef itk::LabelImageToLabelMapFilter<TMaskImage, LabelMapType> Image2LabelMapType;
2769  typename Image2LabelMapType::Pointer img2label = Image2LabelMapType::New();
2770  img2label->SetInput(maskImg);
2771  img2label->SetBackgroundValue(0);
2772 
2773  typedef itk::AutoCropLabelMapFilter<LabelMapType> AutoCropType;
2774  typename AutoCropType::Pointer autoCropper = AutoCropType::New();
2775  autoCropper->SetInput(img2label->GetOutput());
2776  typename AutoCropType::SizeType padSize = {{pixelPadding, pixelPadding, pixelPadding}};
2777  autoCropper->SetCropBorder(padSize);
2778  autoCropper->Update();
2779 
2780  typename AutoCropType::SizeType actualPadSize = autoCropper->GetCropBorder();
2781  milx::PrintDebug("Auto cropping with " + milx::NumberToString(actualPadSize[0]) + " pixel padding.");
2782 
2783  return ExtractSubImage<TImage>(img, autoCropper->GetRegion());
2784 }
2785 
2786 template<class TImage>
2787 std::vector<unsigned char> Image<TImage>::LabelValues(itk::SmartPointer<TImage> labelledImg)
2788 {
2789  typedef itk::LabelObject<unsigned char, milx::imgDimension> LabelObjectType;
2790  typedef itk::LabelMap<LabelObjectType> LabelMapType;
2791  typedef itk::LabelImageToLabelMapFilter<TImage, LabelMapType> LabelImageToLabelMapType;
2792 
2793  typename LabelImageToLabelMapType::Pointer labelImageToLabelMapFilter = LabelImageToLabelMapType::New();
2794  labelImageToLabelMapFilter->SetInput(labelledImg);
2795  labelImageToLabelMapFilter->Update();
2796  typename LabelMapType::Pointer labelMap = labelImageToLabelMapFilter->GetOutput();
2797 
2798  std::vector<unsigned char> labelValues = labelMap->GetLabels();
2799  std::cout << "Found " << labelMap->GetNumberOfLabelObjects() << std::endl;
2800 
2801  return labelValues;
2802 }
2803 
2804 template<class TImage>
2805 template<typename TLabelImage>
2806 itk::SmartPointer< itk::Image<itk::RGBPixel<unsigned char>, TImage::ImageDimension> > Image<TImage>::Overlay(itk::SmartPointer<TImage> img, itk::SmartPointer<TLabelImage> overlayImg)
2807 {
2808  itk::SmartPointer<TImage> rescaledImg = RescaleIntensity(img, 0.0, 255.0);
2809 
2810  typedef unsigned char charPixelType;
2811  typedef itk::Image<charPixelType, TImage::ImageDimension> charImageType;
2812 
2813  typedef itk::CastImageFilter<TLabelImage, charImageType> CastImageType;
2814  typename CastImageType::Pointer castFilter = CastImageType::New();
2815  castFilter->SetInput(overlayImg);
2816  castFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2817  try
2818  {
2819  castFilter->Update();
2820  }
2821  catch (itk::ExceptionObject & ex )
2822  {
2823  PrintError("Failed Casting in Overlay");
2824  PrintError(ex.GetDescription());
2825  }
2826 
2828  typedef itk::LabelOverlayImageFilter< TImage, charImageType, itk::Image<itk::RGBPixel<unsigned char>, TImage::ImageDimension> > OverlayType;
2829  typename OverlayType::Pointer overlayFilter = OverlayType::New();
2830  overlayFilter->SetInput( rescaledImg );
2831  overlayFilter->SetLabelImage( castFilter->GetOutput() );
2832  overlayFilter->SetOpacity( 0.6 );
2833  overlayFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2834  try
2835  {
2836  overlayFilter->Update();
2837  }
2838  catch (itk::ExceptionObject & ex )
2839  {
2840  PrintError("Failed Generating Overlay");
2841  PrintError(ex.GetDescription());
2842  }
2843 
2844  return overlayFilter->GetOutput();
2845 }
2846 
2847 template<class TImage>
2848 template<typename TLabelImage>
2849 itk::SmartPointer< itk::Image<itk::RGBPixel<unsigned char>, TImage::ImageDimension> > Image<TImage>::OverlayContour(itk::SmartPointer<TImage> img, itk::SmartPointer<TLabelImage> overlayImg)
2850 {
2851  itk::SmartPointer<TImage> rescaledImg = RescaleIntensity(img, 0.0, 255.0);
2852 
2853  typedef unsigned char charPixelType;
2854  typedef itk::Image<charPixelType, TImage::ImageDimension> charImageType;
2855 
2856  typedef itk::LabelContourImageFilter<TLabelImage, charImageType> LabelContourType;
2857  typename LabelContourType::Pointer contourFilter = LabelContourType::New();
2858  contourFilter->SetInput(overlayImg);
2859  contourFilter->SetFullyConnected(true);
2860  contourFilter->SetBackgroundValue(0.0);
2861  contourFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2862  try
2863  {
2864  contourFilter->Update();
2865  }
2866  catch (itk::ExceptionObject & ex )
2867  {
2868  PrintError("Failed Generating Label Contours in Overlay");
2869  PrintError(ex.GetDescription());
2870  }
2871 
2873  typedef itk::LabelOverlayImageFilter< TImage, charImageType, itk::Image<itk::RGBPixel<unsigned char>, TImage::ImageDimension> > OverlayType;
2874  typename OverlayType::Pointer overlayFilter = OverlayType::New();
2875  overlayFilter->SetInput( rescaledImg );
2876  overlayFilter->SetLabelImage( contourFilter->GetOutput() );
2877  overlayFilter->SetOpacity( 0.6 );
2878  overlayFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2879  try
2880  {
2881  overlayFilter->Update();
2882  }
2883  catch (itk::ExceptionObject & ex )
2884  {
2885  PrintError("Failed Generating Overlay Contour");
2886  PrintError(ex.GetDescription());
2887  }
2888 
2889  return overlayFilter->GetOutput();
2890 }
2891 #endif // (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
2892 
2893 //Thresholding
2894 template<class TImage>
2895 itk::SmartPointer<TImage> Image<TImage>::ThresholdAboveImage(itk::SmartPointer<TImage> img, float outsideValue, float aboveValue)
2896 {
2897  typedef itk::ThresholdImageFilter<TImage> ThresholdImageFilterType;
2898  typename ThresholdImageFilterType::Pointer thresholding = ThresholdImageFilterType::New();
2899  thresholding->SetInput(img);
2900  thresholding->SetOutsideValue(outsideValue);
2901  thresholding->ThresholdAbove(aboveValue);
2902  thresholding->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2903  try
2904  {
2905  thresholding->Update();
2906  }
2907  catch (itk::ExceptionObject & ex )
2908  {
2909  PrintError("Failed Computing Threshold");
2910  PrintError(ex.GetDescription());
2911  }
2912 
2913  return thresholding->GetOutput();
2914 }
2915 
2916 template<class TImage>
2917 itk::SmartPointer<TImage> Image<TImage>::ThresholdBelowImage(itk::SmartPointer<TImage> img, float outsideValue, float belowValue)
2918 {
2919  typedef itk::ThresholdImageFilter<TImage> ThresholdImageFilterType;
2920  typename ThresholdImageFilterType::Pointer thresholding = ThresholdImageFilterType::New();
2921  thresholding->SetInput(img);
2922  thresholding->SetOutsideValue(outsideValue);
2923  thresholding->ThresholdBelow(belowValue);
2924  thresholding->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2925  try
2926  {
2927  thresholding->Update();
2928  }
2929  catch (itk::ExceptionObject & ex )
2930  {
2931  PrintError("Failed Computing Threshold");
2932  PrintError(ex.GetDescription());
2933  }
2934 
2935  return thresholding->GetOutput();
2936 }
2937 
2938 template<class TImage>
2939 itk::SmartPointer<TImage> Image<TImage>::ThresholdImage(itk::SmartPointer<TImage> img, float outsideValue, float belowValue, float aboveValue)
2940 {
2941  typedef itk::ThresholdImageFilter<TImage> ThresholdImageFilterType;
2942  typename ThresholdImageFilterType::Pointer thresholding = ThresholdImageFilterType::New();
2943  thresholding->SetInput(img);
2944  thresholding->SetOutsideValue(outsideValue);
2945  thresholding->ThresholdOutside(belowValue, aboveValue);
2946  thresholding->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2947  try
2948  {
2949  thresholding->Update();
2950  }
2951  catch (itk::ExceptionObject & ex )
2952  {
2953  PrintError("Failed Computing Threshold");
2954  PrintError(ex.GetDescription());
2955  }
2956 
2957  return thresholding->GetOutput();
2958 }
2959 
2960 template<class TImage>
2961 template<typename TOutImage>
2962 itk::SmartPointer<TOutImage> Image<TImage>::BinaryThresholdImage(itk::SmartPointer<TImage> img, float outsideValue, float insideValue, float belowValue, float aboveValue)
2963 {
2964  typedef itk::BinaryThresholdImageFilter<TImage, TOutImage> ThresholdImageFilterType;
2965  typename ThresholdImageFilterType::Pointer thresholding = ThresholdImageFilterType::New();
2966  thresholding->SetInput(img);
2967  thresholding->SetOutsideValue(outsideValue);
2968  thresholding->SetInsideValue(insideValue);
2969  thresholding->SetLowerThreshold(belowValue);
2970  thresholding->SetUpperThreshold(aboveValue);
2971  thresholding->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2972  try
2973  {
2974  thresholding->Update();
2975  }
2976  catch (itk::ExceptionObject & ex )
2977  {
2978  PrintError("Failed Computing Threshold");
2979  PrintError(ex.GetDescription());
2980  }
2981 
2982  return thresholding->GetOutput();
2983 }
2984 
2985 template<class TImage>
2986 double Image<TImage>::OtsuThreshold(itk::SmartPointer<TImage> img, const int bins)
2987 {
2988  typedef itk::OtsuThresholdImageFilter<TImage, TImage> OtsuThresholdImageCalculatorType;
2989  typename OtsuThresholdImageCalculatorType::Pointer otsuThresholdImageCalculator = OtsuThresholdImageCalculatorType::New();
2990  otsuThresholdImageCalculator->SetInput(img);
2991  otsuThresholdImageCalculator->SetNumberOfHistogramBins(bins);
2992  otsuThresholdImageCalculator->AddObserver(itk::ProgressEvent(), ProgressUpdates);
2993  try
2994  {
2995  otsuThresholdImageCalculator->Update();
2996  }
2997  catch (itk::ExceptionObject & ex )
2998  {
2999  PrintError("Failed Computing Otsu Threshold");
3000  PrintError(ex.GetDescription());
3001  }
3002 
3003  //create binary image
3004  return otsuThresholdImageCalculator->GetThreshold();
3005 }
3006 
3007 template<class TImage>
3008 template<typename TOutImage>
3009 itk::SmartPointer<TOutImage> Image<TImage>::OtsuThresholdImage(itk::SmartPointer<TImage> img, const int bins)
3010 {
3011  typedef itk::OtsuThresholdImageFilter<TImage, TOutImage> OtsuThresholdImageCalculatorType;
3012  typename OtsuThresholdImageCalculatorType::Pointer otsuThresholdImageCalculator = OtsuThresholdImageCalculatorType::New();
3013  otsuThresholdImageCalculator->SetInput(img);
3014  otsuThresholdImageCalculator->SetNumberOfHistogramBins(bins);
3015  otsuThresholdImageCalculator->AddObserver(itk::ProgressEvent(), ProgressUpdates);
3016  try
3017  {
3018  otsuThresholdImageCalculator->Update();
3019  }
3020  catch (itk::ExceptionObject & ex )
3021  {
3022  PrintError("Failed Computing Otsu Threshold of Image");
3023  PrintError(ex.GetDescription());
3024  }
3025  PrintDebug("Otsu Threshold - " + milx::NumberToString(otsuThresholdImageCalculator->GetThreshold()));
3026 
3027  //create binary image
3028  return otsuThresholdImageCalculator->GetOutput();
3029 }
3030 
3031 template<class TImage>
3032 template<typename TOutImage>
3033 itk::SmartPointer<TOutImage> Image<TImage>::OtsuMultipleThresholdImage(itk::SmartPointer<TImage> img, const int bins, const int noOfLabels)
3034 {
3035  typedef itk::OtsuMultipleThresholdsImageFilter<TImage, TOutImage> OtsuThresholdImageType;
3036  typename OtsuThresholdImageType::Pointer otsuThresholdImageCalculator = OtsuThresholdImageType::New();
3037  otsuThresholdImageCalculator->SetInput(img);
3038  otsuThresholdImageCalculator->SetNumberOfHistogramBins(bins);
3039  otsuThresholdImageCalculator->SetNumberOfThresholds(noOfLabels);
3040  otsuThresholdImageCalculator->SetLabelOffset(0);
3041  otsuThresholdImageCalculator->AddObserver(itk::ProgressEvent(), ProgressUpdates);
3042  try
3043  {
3044  otsuThresholdImageCalculator->Update();
3045  }
3046  catch (itk::ExceptionObject & ex )
3047  {
3048  PrintError("Failed Computing Otsu Multiple Threshold");
3049  PrintError(ex.GetDescription());
3050  }
3051  //print thresholds
3052  typename OtsuThresholdImageType::ThresholdVectorType thresholds = otsuThresholdImageCalculator->GetThresholds();
3053  std::string thresholdStrs;
3054  for(size_t j = 0; j < thresholds.size(); j ++)
3055  thresholdStrs += milx::NumberToString(thresholds[j]) + ", ";
3056  PrintDebug("Otsu Thresholds - " + thresholdStrs);
3057 
3058  return otsuThresholdImageCalculator->GetOutput();
3059 }
3060 
3061 template<class TImage>
3062 template<typename TScalarImage>
3063 itk::SmartPointer<TScalarImage> Image<TImage>::VectorMagnitude(itk::SmartPointer<TImage> img)
3064 {
3065 #if (ITK_VERSION_MAJOR > 3)
3066  typedef itk::VectorMagnitudeImageFilter<TImage, TScalarImage> MagnitudeImageFilterType;
3067 #else
3068  typedef itk::GradientToMagnitudeImageFilter<TImage, TScalarImage> MagnitudeImageFilterType;
3069 #endif
3070  typename MagnitudeImageFilterType::Pointer magnitudeFilter = MagnitudeImageFilterType::New();
3071  magnitudeFilter->SetInput(img);
3072  magnitudeFilter->AddObserver(itk::ProgressEvent(), ProgressUpdates);
3073  try
3074  {
3075  magnitudeFilter->Update();
3076  }
3077  catch (itk::ExceptionObject & ex )
3078  {
3079  PrintError("Failed Computing Magnitude of Image");
3080  PrintError(ex.GetDescription());
3081  }
3082 
3083  return magnitudeFilter->GetOutput();
3084 }
3085 
3086 //Collection members
3087 template<class TImage>
3088 void Image<TImage>::InformationCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
3089 {
3090  const size_t n = images.size();
3091 
3092  for(size_t j = 0; j < n; j ++)
3093  {
3094  Information(images[j]);
3095  std::cout << std::endl;
3096  }
3097 }
3098 
3099 template<class TImage>
3100 void Image<TImage>::MatchHistogramCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TImage> refImage, const int bins)
3101 {
3102  const size_t n = images.size();
3103 
3104  for(size_t j = 0; j < n; j ++)
3105  images[j] = MatchHistogram(images[j], refImage, bins);
3106 }
3107 
3108 template<class TImage>
3109 void Image<TImage>::CheckerboardCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TImage> refImage, const int squares)
3110 {
3111  const size_t n = images.size();
3112 
3113  for(size_t j = 0; j < n; j ++)
3114  images[j] = CheckerBoard(images[j], refImage, squares);
3115 }
3116 
3117 template<class TImage>
3118 void Image<TImage>::RescaleIntensityCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float belowValue, float aboveValue)
3119 {
3120  const size_t n = images.size();
3121 
3122  for(size_t j = 0; j < n; j ++)
3123  images[j] = RescaleIntensity(images[j], belowValue, aboveValue);
3124 }
3125 
3126 template<class TImage>
3127 void Image<TImage>::InvertIntensityCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
3128 {
3129  const size_t n = images.size();
3130 
3131  for(size_t j = 0; j < n; j ++)
3132  {
3133  float maxValue = ImageMaximum(images[j]);
3134  images[j] = InvertIntensity(images[j], maxValue);
3135  }
3136 }
3137 
3138 template<class TImage>
3139 void Image<TImage>::FlipCollection(std::vector< typename itk::SmartPointer<TImage> > &images, bool xAxis, bool yAxis, bool zAxis, bool aboutOrigin)
3140 {
3141  const size_t n = images.size();
3142 
3143  for(size_t j = 0; j < n; j ++)
3144  {
3145  images[j] = FlipImage(images[j], xAxis, yAxis, zAxis, aboutOrigin);
3146  }
3147 }
3148 
3149 template<class TImage>
3150 template<typename TOutImage>
3151 std::vector< typename itk::SmartPointer<TOutImage> > Image<TImage>::AnisotropicDiffusionCollection(const std::vector< typename itk::SmartPointer<TImage> > &images, const int iterations, float timestep)
3152 {
3153  const size_t n = images.size();
3154 
3155  std::vector< itk::SmartPointer<TOutImage> > collection;
3156  for(size_t j = 0; j < n; j ++)
3157  {
3158  typename TImage::SizeType imgSize = images[j]->GetLargestPossibleRegion().GetSize();
3159 
3160  if(timestep < 0.0) //auto set timestep
3161  {
3162  size_t maxDimension = 0;
3163  for(size_t k = 0; k < imgSize.GetSizeDimension(); k ++)
3164  {
3165  if(imgSize[k] > maxDimension)
3166  maxDimension = imgSize[k];
3167  }
3168  milx::PrintDebug("Timestep will be reciprocal of " + milx::NumberToString(maxDimension));
3169  timestep = 2.0/maxDimension;
3170  }
3171 
3172  itk::SmartPointer<TOutImage> newImage = AnisotropicDiffusion<TOutImage>(images[j], iterations, timestep);
3173  collection.push_back(newImage);
3174  }
3175 
3176  return collection;
3177 }
3178 
3179 template<class TImage>
3180 void Image<TImage>::BilateralCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float sigmaRange, float sigmaSpatial)
3181 {
3182  const size_t n = images.size();
3183 
3184  for(size_t j = 0; j < n; j ++)
3185  images[j] = Bilateral(images[j], sigmaRange, sigmaSpatial);
3186 }
3187 
3188 template<class TImage>
3189 template<typename TOutImage>
3190 std::vector< typename itk::SmartPointer<TOutImage> > Image<TImage>::CastCollection(const std::vector< typename itk::SmartPointer<TImage> > &images)
3191 {
3192  const size_t n = images.size();
3193 
3194  std::vector< itk::SmartPointer<TOutImage> > collection;
3195  for(size_t j = 0; j < n; j ++)
3196  {
3197  itk::SmartPointer<TOutImage> newImage = CastImage<TOutImage>(images[j]);
3198  collection.push_back(newImage);
3199  }
3200 
3201  return collection;
3202 }
3203 
3204 template<class TImage>
3205 void Image<TImage>::MedianCollection(std::vector< typename itk::SmartPointer<TImage> > &images, const int radius)
3206 {
3207  const size_t n = images.size();
3208 
3209  for(size_t j = 0; j < n; j ++)
3210  images[j] = Median(images[j], radius);
3211 }
3212 
3213 template<class TImage>
3214 void Image<TImage>::GradientMagnitudeCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
3215 {
3216  const size_t n = images.size();
3217 
3218  for(size_t j = 0; j < n; j ++)
3219  images[j] = GradientMagnitude(images[j]);
3220 }
3221 
3222 template<class TImage>
3223 void Image<TImage>::LaplacianCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
3224 {
3225  const size_t n = images.size();
3226 
3227  for(size_t j = 0; j < n; j ++)
3228  images[j] = Laplacian(images[j]);
3229 }
3230 
3231 template<class TImage>
3232 template<typename TOutImage>
3233 std::vector< typename itk::SmartPointer<TOutImage> > Image<TImage>::DistanceMapCollection(const std::vector< typename itk::SmartPointer<TImage> > &images, const bool binaryImage, const bool signedDistances, const bool computeInsideObject, const bool squaredDistance)
3234 {
3235  const size_t n = images.size();
3236 
3237  std::vector< itk::SmartPointer<TOutImage> > collection;
3238  for(size_t j = 0; j < n; j ++)
3239  {
3240  itk::SmartPointer<TOutImage> newImage = DistanceMap<TOutImage>(images[j], binaryImage, signedDistances, computeInsideObject, squaredDistance);
3241  collection.push_back(newImage);
3242  }
3243 
3244  return collection;
3245 }
3246 
3247 //Thresholding
3248 template<class TImage>
3249 void Image<TImage>::ThresholdAboveCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float outsideValue, float aboveValue)
3250 {
3251  const size_t n = images.size();
3252 
3253  for(size_t j = 0; j < n; j ++)
3254  images[j] = ThresholdAboveImage(images[j], outsideValue, aboveValue);
3255 }
3256 
3257 template<class TImage>
3258 void Image<TImage>::ThresholdBelowCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float outsideValue, float belowValue)
3259 {
3260  const size_t n = images.size();
3261 
3262  for(size_t j = 0; j < n; j ++)
3263  images[j] = ThresholdBelowImage(images[j], outsideValue, belowValue);
3264 }
3265 
3266 template<class TImage>
3267 void Image<TImage>::ThresholdCollection(std::vector< typename itk::SmartPointer<TImage> > &images, float outsideValue, float belowValue, float aboveValue)
3268 {
3269  const size_t n = images.size();
3270 
3271  for(size_t j = 0; j < n; j ++)
3272  images[j] = ThresholdImage(images[j], outsideValue, belowValue, aboveValue);
3273 }
3274 
3275 template<class TImage>
3276 template<typename TOutImage>
3277 std::vector< typename itk::SmartPointer<TOutImage> > Image<TImage>::BinaryThresholdCollection(const std::vector< typename itk::SmartPointer<TImage> > &images, float outsideValue, float insideValue, float belowValue, float aboveValue)
3278 {
3279  const size_t n = images.size();
3280 
3281  std::vector< itk::SmartPointer<TOutImage> > collection;
3282  for(size_t j = 0; j < n; j ++)
3283  {
3284  itk::SmartPointer<TOutImage> newImage = BinaryThresholdImage<TOutImage>(images[j], outsideValue, insideValue, belowValue, aboveValue);
3285  collection.push_back(newImage);
3286  }
3287 
3288  return collection;
3289 }
3290 
3291 template<class TImage>
3292 template<typename TOutImage>
3293 std::vector< typename itk::SmartPointer<TOutImage> > Image<TImage>::OtsuMultipleThresholdCollection(const std::vector< typename itk::SmartPointer<TImage> > &images, const int bins, const int noOfLabels)
3294 {
3295  const size_t n = images.size();
3296 
3297  std::vector< itk::SmartPointer<TOutImage> > collection;
3298  for(size_t j = 0; j < n; j ++)
3299  {
3300  itk::SmartPointer<TOutImage> newImage = OtsuMultipleThresholdImage<TOutImage>(images[j], bins, noOfLabels);
3301  collection.push_back(newImage);
3302  }
3303 
3304  return collection;
3305 }
3306 
3307 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
3308 template<class TImage>
3309 template<class TMaskImage>
3310 void Image<TImage>::MaskAndCropCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TMaskImage> maskImage, const size_t pixelPadding)
3311 {
3312  const size_t n = images.size();
3313 
3314  for(size_t j = 0; j < n; j ++)
3315  images[j] = MaskAndCropImage<TMaskImage>(images[j], maskImage, pixelPadding);
3316 }
3317 #endif // (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
3318 
3319 template<class TImage>
3320 template<class TMaskImage>
3321 void Image<TImage>::MaskCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TMaskImage> maskImage)
3322 {
3323  const size_t n = images.size();
3324 
3325  for(size_t j = 0; j < n; j ++)
3326  images[j] = MaskImage<TMaskImage>(images[j], maskImage);
3327 }
3328 
3329 template<class TImage>
3330 void Image<TImage>::RelabelCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
3331 {
3332  const size_t n = images.size();
3333 
3334  for(size_t j = 0; j < n; j ++)
3335  images[j] = RelabelImage(images[j]);
3336 }
3337 
3338 template<class TImage>
3339 template<class TRefImage>
3340 void Image<TImage>::ResampleCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TRefImage> refImage)
3341 {
3342  const size_t n = images.size();
3343 
3344  for(size_t j = 0; j < n; j ++)
3345  images[j] = ResampleImage<TRefImage>(images[j], refImage);
3346 }
3347 
3348 template<class TImage>
3349 template<class TRefImage>
3350 void Image<TImage>::ResampleLabelCollection(std::vector< typename itk::SmartPointer<TImage> > &images, itk::SmartPointer<TRefImage> refImage)
3351 {
3352  const size_t n = images.size();
3353 
3354  for(size_t j = 0; j < n; j ++)
3355  images[j] = ResampleLabel<TRefImage>(images[j], refImage);
3356 }
3357 
3358 template<class TImage>
3359 itk::SmartPointer<TImage> Image<TImage>::AddCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
3360 {
3361  const size_t n = images.size();
3362 
3363  itk::SmartPointer<TImage> result = DuplicateImage(images[0]);
3364  for(size_t j = 1; j < n; j ++)
3365  result = AddImages(result, images[j]);
3366 
3367  return result;
3368 }
3369 
3370 template<class TImage>
3371 itk::SmartPointer<TImage> Image<TImage>::SubtractCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
3372 {
3373  const size_t n = images.size();
3374 
3375  itk::SmartPointer<TImage> result = DuplicateImage(images[0]);
3376  for(size_t j = 1; j < n; j ++)
3377  result = SubtractImages(result, images[j]);
3378 
3379  return result;
3380 }
3381 
3382 template<class TImage>
3383 template<class TOutImage>
3384 itk::SmartPointer<TOutImage> Image<TImage>::AverageCollection(std::vector< typename itk::SmartPointer<TImage> > &images)
3385 {
3386  const size_t n = images.size();
3387 
3388  itk::SmartPointer<TImage> sumImg = AddCollection(images);
3389  itk::SmartPointer<TOutImage> averageImg = ScaleImage<TOutImage>(sumImg, 1.0/n);
3390 
3391  return averageImg;
3392 }
3393 
3394 template<class TImage>
3395 itk::SmartPointer<TImage> Image<TImage>::AverageVectorCollection(std::vector< typename itk::SmartPointer<TImage> > &images, int numberOfComponents)
3396 {
3397  const size_t n = images.size();
3398 
3399  itk::SmartPointer<TImage> sumImg = AddCollection(images);
3400  itk::SmartPointer<TImage> averageImg = ScaleVectorImage(sumImg, 1.0/n, numberOfComponents);
3401 
3402  return averageImg;
3403 }
3404 
3405 template<class TImage>
3406 template<class TJoinedImage>
3407 itk::SmartPointer<TJoinedImage> Image<TImage>::JoinCollection(std::vector< typename itk::SmartPointer<TImage> > &images, const double origin, const double spacing)
3408 {
3409  const size_t n = images.size();
3410 
3411  typedef itk::JoinSeriesImageFilter<TImage, TJoinedImage> JoinSeriesFilterType;
3412  typename JoinSeriesFilterType::Pointer joinSeries = JoinSeriesFilterType::New();
3413  joinSeries->SetOrigin(origin);
3414  joinSeries->SetSpacing(spacing);
3415  for(size_t j = 0; j < n; j ++)
3416  joinSeries->PushBackInput(images[j]);
3417  try
3418  {
3419  joinSeries->Update();
3420  }
3421  catch (itk::ExceptionObject & ex )
3422  {
3423  PrintError("Failed Joining Images");
3424  PrintError(ex.GetDescription());
3425  }
3426 
3427  return joinSeries->GetOutput();
3428 }
3429 
3430 #if (ITK_REVIEW || ITK_VERSION_MAJOR > 3)
3431 template<class TImage>
3432 itk::SmartPointer<TImage> Image<TImage>::MergeLabelledImages(std::vector< typename itk::SmartPointer<TImage> > &images, const size_t mergeType)
3433 {
3434  const size_t n = images.size();
3435 
3436  typedef itk::LabelObject<unsigned char, milx::imgDimension> LabelObjectType;
3437  typedef itk::LabelMap<LabelObjectType> LabelMapType;
3438  typedef itk::LabelImageToLabelMapFilter<TImage, LabelMapType> LabelImageToLabelMapType;
3439  typedef itk::MergeLabelMapFilter<LabelMapType> MergerType;
3440  typename MergerType::Pointer merger = MergerType::New();
3441  if(mergeType == 0)
3442  merger->SetMethod(MergerType::KEEP);
3443  else if(mergeType == 1)
3444  merger->SetMethod(MergerType::AGGREGATE);
3445  else if(mergeType == 2)
3446  merger->SetMethod(MergerType::PACK);
3447  else if(mergeType == 3)
3448  merger->SetMethod(MergerType::STRICT);
3449 
3450  for (size_t i = 0; i < n; i ++)
3451  {
3452  typename LabelImageToLabelMapType::Pointer labelImageToLabelMapFilter = LabelImageToLabelMapType::New();
3453  labelImageToLabelMapFilter->SetInput(images[i]);
3454  labelImageToLabelMapFilter->Update();
3455 
3456  std::cout << "Found " << labelImageToLabelMapFilter->GetOutput()->GetNumberOfLabelObjects() << " in the labelled image " << i << std::endl;
3457 
3458  merger->SetInput(i, labelImageToLabelMapFilter->GetOutput());
3459  }
3460 
3461  try
3462  {
3463  merger->Update();
3464  }
3465  catch (itk::ExceptionObject & ex )
3466  {
3467  PrintError("Failed Merging Labelled Images");
3468  PrintError(ex.GetDescription());
3469  }
3470 
3471  //Convert to image again
3472  typedef itk::LabelMapToLabelImageFilter<LabelMapType, TImage> LabelMapToLabelImageType;
3473  typename LabelMapToLabelImageType::Pointer labelMapToLabelImageFilter = LabelMapToLabelImageType::New();
3474  labelMapToLabelImageFilter->SetInput(merger->GetOutput());
3475  labelMapToLabelImageFilter->Update();
3476 
3477  return labelMapToLabelImageFilter->GetOutput();
3478 }
3479 #endif
3480 
3481 } //end namespace milx
3482 
3483 #endif //__MILXIMAGE_H
void PrintDebug(const std::string msg)
Displays a generic msg to standard error with carriage return if in Debug mode.
Definition: milxGlobal.h:192
static itk::SmartPointer< TOutImage > BinaryThresholdImage(itk::SmartPointer< TImage > img, float outsideValue, float insideValue, float belowValue, float aboveValue)
Generates a binary image thresholded at the intensities below and above a certain level to the outsid...
Definition: milxImage.h:2962
static void FlipCollection(std::vector< typename itk::SmartPointer< TImage > > &images, bool xAxis=false, bool yAxis=true, bool zAxis=false, bool aboutOrigin=true)
Batch process images by flipping each image along axis provided.
Definition: milxImage.h:3139
static itk::SmartPointer< TImage > SubtractCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by the subtracting the images from the first image (pixel-wise).
Definition: milxImage.h:3371
static itk::SmartPointer< TOutImage > AnisotropicDiffusion(itk::SmartPointer< TImage > img, const int iterations, const float timestep)
Generates the gradient anisotropic diffusion (smoothing) of an image using the number of iterations...
Definition: milxImage.h:2569
static itk::SmartPointer< TImage > MultiplyImages(itk::SmartPointer< TImage > img1, itk::SmartPointer< TImage > img2)
Multiplies (element-wise) image 2 from image 1, returning the result.
Definition: milxImage.h:1966
static double ImageMaximum(itk::SmartPointer< TImage > img)
Returns the maximum pixel value of an image.
Definition: milxImage.h:2102
static itk::SmartPointer< TImage > ThresholdImage(itk::SmartPointer< TImage > img, float outsideValue, float belowValue, float aboveValue)
Generates an image with the intensities below and above a certain level thresholded (capped) to the o...
Definition: milxImage.h:2939
static itk::SmartPointer< TOutImage > ScaleImage(itk::SmartPointer< TImage > img, float scaling)
Scales the image intensities by scaling factor and returns the result.
Definition: milxImage.h:1989
static itk::SmartPointer< TOutImage > OtsuMultipleThresholdImage(itk::SmartPointer< TImage > img, const int bins, const int noOfLabels=1)
Generates the multiple Otsu threshold of an image of given the number of histogram bins...
Definition: milxImage.h:3033
virtual ~Image()
Standard Destructor.
Definition: milxImage.h:191
static vtkSmartPointer< vtkImageData > ConvertITKImageToVTKImage(itk::SmartPointer< TImage > img)
Converts a ITK image object to an VTK image object. You MUST DeepCopy the result as the ITK smartpoin...
Definition: milxImage.h:879
static void MedianCollection(std::vector< typename itk::SmartPointer< TImage > > &images, const int radius)
Batch process images by smoothing each image using the Median of the neighbourhood.
Definition: milxImage.h:3205
static itk::SmartPointer< TOutImage > CastImage(itk::SmartPointer< TImage > img)
Casts an image from one type to another.
Definition: milxImage.h:1117
static itk::SmartPointer< TImage > ConvertVTKImageToITKImage(vtkSmartPointer< vtkImageData > img)
Converts a VTK image object to an ITK image object.
Definition: milxImage.h:1093
itk::SmartPointer< TImage > PreviousImage
Holds the previous image in the pipeline.
Definition: milxImage.h:851
static itk::SmartPointer< TImage > RelabelImage(itk::SmartPointer< TImage > labelledImg)
Returns a labelled image with labelled values relabelled consecutively based on connectivity.
Definition: milxImage.h:2690
static double ImageMinimum(itk::SmartPointer< TImage > img)
Returns the minimum pixel value of an image.
Definition: milxImage.h:2114
static std::vector< typename itk::SmartPointer< TOutImage > > AnisotropicDiffusionCollection(const std::vector< typename itk::SmartPointer< TImage > > &images, const int iterations, float timestep=-1.0)
Batch process images by smoothing each image using Gradient Anisotropic Diffusion.
Definition: milxImage.h:3151
static itk::SmartPointer< TImage > CheckerBoard(itk::SmartPointer< TImage > img, itk::SmartPointer< TImage > imgToCheckerBoard, const int numberOfSquares=0)
Generates a checker board pattern where each box alternates between two images. Ideal for comparison ...
Definition: milxImage.h:2188
Image()
Standard constructor.
Definition: milxImage.h:858
void PrintError(const std::string msg)
Displays a generic msg to standard error with carriage return.
Definition: milxGlobal.h:202
static itk::SmartPointer< TOutImage > ResampleLabel(itk::SmartPointer< TImage > img, itk::SmartPointer< TOutImage > refImg)
Resample the a labelled image into a new reference image space given ensuring no interpolator artefac...
Definition: milxImage.h:1880
static void ThresholdAboveCollection(std::vector< typename itk::SmartPointer< TImage > > &images, float outsideValue, float aboveValue)
Batch process images by thresholding each image from above.
Definition: milxImage.h:3249
static vtkSmartPointer< vtkImageData > ConvertITKVectorImageToVTKImage(itk::SmartPointer< TImage > img)
Converts a ITK Vector image object to an VTK image object. You MUST DeepCopy the result as the ITK sm...
Definition: milxImage.h:900
static itk::SmartPointer< TImage > CopyInformation(itk::SmartPointer< TImage > img, itk::SmartPointer< TImage > imgToMatch, bool originOnly=false)
Changes the image info to match that of provided image. By default, only the spacing, region and origin are changed.
Definition: milxImage.h:513
static itk::SmartPointer< TScalarImage > VectorMagnitude(itk::SmartPointer< TImage > img)
Generates an image comprised of the magnitude of the vectors in the image.
Definition: milxImage.h:3063
static itk::SmartPointer< TImage > SubtractImages(itk::SmartPointer< TImage > img1, itk::SmartPointer< TImage > img2)
Subtracts image 2 from image 1, returning the result.
Definition: milxImage.h:1944
static itk::SmartPointer< TImageSlice > ExtractSlice(itk::SmartPointer< TImage > img, int *extent)
Extract a sub image (such as a slice) from the current image given the extent (typically obtained fro...
Definition: milxImage.h:1324
static itk::SmartPointer< TImage > GradientMagnitude(itk::SmartPointer< TImage > img)
Generates an image with the gradient magnitude of the given image.
Definition: milxImage.h:2283
itk::SmartPointer< TImage > CurrentImage
Holds the current image in the pipeline.
Definition: milxImage.h:850
Base class for images that contains non-templated features so to allow dynamic template instantiation...
Definition: milxImage.h:122
static void ThresholdCollection(std::vector< typename itk::SmartPointer< TImage > > &images, float outsideValue, float belowValue, float aboveValue)
Batch process images by thresholding each image within band provided.
Definition: milxImage.h:3267
#define SMILI_EXPORT
DLL Function Symbol for Windows. It is empty for other OSes.
Definition: milxGlobal.h:227
static itk::SmartPointer< TImage > ThresholdAboveImage(itk::SmartPointer< TImage > img, float outsideValue, float aboveValue)
Generates an image with the intensities above a certain level thresholded (capped) to the outsideValu...
Definition: milxImage.h:2895
static itk::SmartPointer< TImage > ImportVectorToImage(vnl_vector< TVector > &vec, typename TImage::SizeType size, itk::SmartPointer< TImage > image=NULL)
Imports a VNL vector to an ITK image object.
Definition: milxImage.h:1195
static void BilateralCollection(std::vector< typename itk::SmartPointer< TImage > > &images, float sigmaRange, float sigmaSpatial)
Batch process images by bilateral smoothing each image.
Definition: milxImage.h:3180
static itk::SmartPointer< TOutImage > OtsuThresholdImage(itk::SmartPointer< TImage > img, const int bins)
Generates the Otsu threshold of an image of given the number of histogram bins.
Definition: milxImage.h:3009
static std::string ImageOrientation(itk::SmartPointer< TImage > img)
Returns the orientation flag of an image.
Definition: milxImage.h:2126
static itk::SmartPointer< TImage > GaussianSmooth(itk::SmartPointer< TImage > img, const float variance)
Generates the Gaussian smoothing (by convolution) of an image using the variance provided.
Definition: milxImage.h:2616
static itk::SmartPointer< TImage > ImportMatrixToImage(vnl_matrix< TMatrix > &matrix, itk::SmartPointer< TImage > image=NULL)
Imports a VNL matrix to an ITK image object.
Definition: milxImage.h:1234
static itk::SmartPointer< TImage > ThresholdBelowImage(itk::SmartPointer< TImage > img, float outsideValue, float aboveValue)
Generates an image with the intensities below a certain level thresholded (capped) to the outsideValu...
Definition: milxImage.h:2917
static itk::SmartPointer< TImage > Laplacian(itk::SmartPointer< TImage > img)
Generates an Laplacian image from the input image.
Definition: milxImage.h:2323
static itk::SmartPointer< TImage > ResizeImage(itk::SmartPointer< TImage > img, typename TImage::SizeType imgSize)
Resizes current image using current spacing.
Definition: milxImage.h:1379
static itk::SmartPointer< TImage > MatchHistogram(itk::SmartPointer< TImage > img, itk::SmartPointer< TImage > imgToMatch, const int bins=128)
Changes the image gray levels to match histogram of image provided.
Definition: milxImage.h:2421
static itk::SmartPointer< TOutImage > DistanceMap(itk::SmartPointer< TImage > img, const bool binaryImage=false, const bool signedDistances=false, const bool computeInsideObject=false, const bool squaredDistance=false)
Generates a distance map image, where the distances are from the object boundary to the outside or vi...
Definition: milxImage.h:2446
static itk::SmartPointer< TOutImage > TransformImage(itk::SmartPointer< TImage > img, itk::SmartPointer< TOutImage > refImg, itk::SmartPointer< TTransform > transf, const bool inverse, const int interp=1)
Transform the image into a new reference image space given the transform.
Definition: milxImage.h:1484
static void ResampleCollection(std::vector< typename itk::SmartPointer< TImage > > &images, itk::SmartPointer< TRefImage > refImage)
Batch process images by resampling images using the reference image provided.
Definition: milxImage.h:3340
static double OtsuThreshold(itk::SmartPointer< TImage > img, const int bins)
Returns the Otsu threshold of an image of given the number of histogram bins.
Definition: milxImage.h:2986
static std::vector< typename itk::SmartPointer< TOutImage > > DistanceMapCollection(const std::vector< typename itk::SmartPointer< TImage > > &images, const bool binaryImage=false, const bool signedDistances=false, const bool computeInsideObject=false, const bool squaredDistance=false)
Batch process images by computing the distance map of each image. Output should always be a float ima...
Definition: milxImage.h:3233
static itk::SmartPointer< TImage > Bilateral(itk::SmartPointer< TImage > img, const float sigmaRange, const float sigmaSpatial)
Generates the (non-linear) Bilateral smoothing of an image using the sigmas provided.
Definition: milxImage.h:2594
static itk::SmartPointer< TImage > SobelEdges(itk::SmartPointer< TImage > img)
Generates an image with Sobel edges, i.e. uses the Sobel edge detection on the input image...
Definition: milxImage.h:2303
Represents an image (i.e. an regular rectangular array with scalar values) and their common operation...
Definition: milxImage.h:174
static itk::SmartPointer< TImage > HistogramEqualisation(itk::SmartPointer< TImage > img, float alpha=0.3, float beta=0.3, float radius=5)
Generates an image with the intensities after histogram equalisation. Defaults to classic histogram e...
Definition: milxImage.h:2258
static itk::SmartPointer< TImage > DifferenceImages(itk::SmartPointer< TImage > img1, itk::SmartPointer< TImage > img2)
Same as SubtractImages().
Definition: milxImage.h:397
static itk::SmartPointer< TOutImage > ResampleImage(itk::SmartPointer< TImage > img, itk::SmartPointer< TOutImage > refImg, const bool linearInterp=false)
Resample the image into a new reference image space given.
Definition: milxImage.h:1826
static std::vector< typename itk::SmartPointer< TOutImage > > BinaryThresholdCollection(const std::vector< typename itk::SmartPointer< TImage > > &images, float outsideValue, float insideValue, float belowValue, float aboveValue)
Batch process images by binary thresholding each image within band and inside value provided...
Definition: milxImage.h:3277
static vtkSmartPointer< vtkImageData > ApplyOrientationToVTKImage(vtkSmartPointer< vtkImageData > img, itk::SmartPointer< TImage > refImage, vtkSmartPointer< vtkMatrix4x4 > &transformMatrix, const bool labelledImage, const bool flipY=true)
Applies orientation/direction and origin to a VTK image from a reference image.
Definition: milxImage.h:1034
static void ThresholdBelowCollection(std::vector< typename itk::SmartPointer< TImage > > &images, float outsideValue, float belowValue)
Batch process images by thresholding each image from below.
Definition: milxImage.h:3258
static void InvertIntensityCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by inverting the intensities for each image.
Definition: milxImage.h:3127
itk::SmartPointer< TImage > Result()
Returns the current image, i.e. the result of the latest operation.
Definition: milxImage.h:204
static itk::SmartPointer< TSubImage > ExtractSubImage(itk::SmartPointer< TImage > img, typename TImage::RegionType imgRegion)
Extract a sub image (such as a slice) from the current image given the region.
Definition: milxImage.h:1298
static itk::SmartPointer< TImage > AddCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by the adding the images together (pixel-wise).
Definition: milxImage.h:3359
void PrintInfo(const std::string msg)
Displays a generic msg to standard output with carriage return.
Definition: milxGlobal.h:174
std::string NumberToString(double num, unsigned zeroPad=0)
Number to string converter.
Definition: milxGlobal.h:112
itk::SmartPointer< TImage > PreviousResult()
Returns the previous image, i.e. the result of the penultimate operation.
Definition: milxImage.h:214
static itk::SmartPointer< TImage > MatchInformation(itk::SmartPointer< TImage > img, itk::SmartPointer< TImage > imgToMatch, bool originOnly=false)
Changes the image info to match that of provided image. By default, only the spacing, region and origin are changed.
Definition: milxImage.h:2386
static itk::SmartPointer< TImage > MaskImage(itk::SmartPointer< TImage > img, itk::SmartPointer< TMaskImage > maskImg)
Mask an image by given binary image.
Definition: milxImage.h:2662
static itk::SmartPointer< TImage > AverageVectorCollection(std::vector< typename itk::SmartPointer< TImage > > &images, int numberOfComponents)
Batch process images by the averaging the vector images (pixel-wise).
Definition: milxImage.h:3395
static itk::SmartPointer< TImage > DifferenceCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Same as SubtractCollection.
Definition: milxImage.h:817
static itk::SmartPointer< TImage > AddImages(itk::SmartPointer< TImage > img1, itk::SmartPointer< TImage > img2)
Adds image 1 to image 2, returning the result.
Definition: milxImage.h:1922
static void Information(itk::SmartPointer< TImage > img)
Prints the information of the image to standard output.
Definition: milxImage.h:2067
static itk::SmartPointer< TImage > FlipImage(itk::SmartPointer< TImage > img, bool xAxis=false, bool yAxis=true, bool zAxis=false, bool aboutOrigin=true)
Generates an image with axes flipped. Currently flips the y-axis as default, but can do other axes al...
Definition: milxImage.h:2501
static itk::SmartPointer< TImage > DuplicateImage(itk::SmartPointer< TImage > img)
Duplicates the image into a new image.
Definition: milxImage.h:1152
static void ResampleLabelCollection(std::vector< typename itk::SmartPointer< TImage > > &images, itk::SmartPointer< TRefImage > refImage)
Batch process labelled images by resampling images using the reference image provided. This uses the nearest neighbour interploator.
Definition: milxImage.h:3350
static void RelabelCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by relabelling unconnected regions consecutatively.
Definition: milxImage.h:3330
static itk::SmartPointer< TImage > ScaleVectorImage(itk::SmartPointer< TImage > img, float scaling, int numberOfComponents)
Scales each component of the vector image intensities by scaling factor and returns the result...
Definition: milxImage.h:2011
static itk::SmartPointer< TImage > CannyEdges(itk::SmartPointer< TImage > img, float variance, float lowerThreshold, float upperThreshold)
Generates an Canny edge image from the input image.
Definition: milxImage.h:2343
static itk::SmartPointer< itk::Image< float, TImage::ImageDimension > > Normalization(itk::SmartPointer< TImage > img)
Generates an image which is statistically normalised so having pixel values between -1 and 1...
Definition: milxImage.h:2366
static itk::SmartPointer< TImage > BlankImage(typename TImage::PixelType value, typename TImage::SizeType imgSize)
Creates an empty image filled with value given.
Definition: milxImage.h:1277
static std::vector< typename itk::SmartPointer< TOutImage > > CastCollection(const std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by casting each image to type provided in templates.
Definition: milxImage.h:3190
static void InformationCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by output the image information for each image.
Definition: milxImage.h:3088
static void GradientMagnitudeCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by hightlighting the edges of each image using Gradient magnitude information of...
Definition: milxImage.h:3214
static itk::SmartPointer< TImage > InvertIntensity(itk::SmartPointer< TImage > img, float maxValue)
Generates an image with the intensities reversed based on the max pixel value.
Definition: milxImage.h:2237
Applies a linear transformation to the magnitude of pixel vectors in a VectorImage.
static itk::SmartPointer< TOutImage > AverageCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by the averaging the images (pixel-wise).
Definition: milxImage.h:3384
static itk::SmartPointer< TImage > Median(itk::SmartPointer< TImage > img, const int radius)
Generates the (non-linear) median filtering (smoothing) of an image of given neighbourhood radius...
Definition: milxImage.h:2637
static std::vector< typename itk::SmartPointer< TOutImage > > OtsuMultipleThresholdCollection(const std::vector< typename itk::SmartPointer< TImage > > &images, const int bins, const int noOfLabels)
Batch process images by Otsu (histogram-based) thresholding each image having number of partitions/la...
Definition: milxImage.h:3293
static itk::SmartPointer< TImage > PadImageByConstant(itk::SmartPointer< TImage > img, size_t xAxis, size_t yAxis, size_t zAxis, typename TImage::PixelType value)
Generates an image padded by extending each axes by size given in both directions. The created areas have the value given by value.
Definition: milxImage.h:2532
itk::SmartPointer< TImage > GetOutput()
Returns the current image, i.e. the result of the latest operation ITK/VTK style. ...
Definition: milxImage.h:222
static void LaplacianCollection(std::vector< typename itk::SmartPointer< TImage > > &images)
Batch process images by hightlighting the edges of each image using the Laplacian of the image...
Definition: milxImage.h:3223
static itk::SmartPointer< TImageComponent > ExtractComponent(itk::SmartPointer< TImage > img, int component)
Extract a component from the current vector image given the component index.
Definition: milxImage.h:1358
static itk::SmartPointer< TImage > SubsampleImage(itk::SmartPointer< TImage > img, typename TImage::SizeType factors)
Subsamples or shrinks the current image using the factors provided.
Definition: milxImage.h:1459
static void CheckerboardCollection(std::vector< typename itk::SmartPointer< TImage > > &images, itk::SmartPointer< TImage > refImage, const int squares=10)
Batch process images by checkerboarding all of the images to the reference image provided.
Definition: milxImage.h:3109
static void MatchHistogramCollection(std::vector< typename itk::SmartPointer< TImage > > &images, itk::SmartPointer< TImage > refImage, const int bins=128)
Batch process images by matching the histograms of the images to the reference image provided...
Definition: milxImage.h:3100
static void MaskCollection(std::vector< typename itk::SmartPointer< TImage > > &images, itk::SmartPointer< TMaskImage > maskImage)
Batch process images by masking using the mask image provided.
Definition: milxImage.h:3321