mitkImage.h

Go to the documentation of this file.

/*============================================================================
 
The Medical Imaging Interaction Toolkit (MITK)
 
Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.
 
Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.
 
============================================================================*/
 
#ifndef mitkImage_h
#define mitkImage_h
 
#include <mitkBaseData.h>
#include <mitkImageAccessorBase.h>
#include <mitkImageDataItem.h>
#include <mitkImageDescriptor.h>
#include <mitkImageVtkWriteAccessor.h>
#include <mitkLevelWindow.h>
#include <mitkPlaneGeometry.h>
#include <mitkSlicedData.h>
#include <MitkCoreExports.h>
#include <mitkProportionalTimeGeometry.h>
 
#ifndef __itkHistogram_h
#include <itkHistogram.h>
#endif
 
class vtkImageData;
 
namespace itk
{
  template <class T>
  class MutexLockHolder;
}
 
namespace mitk
{
  class SubImageSelector;
  class ImageTimeSelector;
 
  class ImageStatisticsHolder;
 
  class MITKCORE_EXPORT Image : public SlicedData
  {
    friend class SubImageSelector;
 
    friend class ImageAccessorBase;
    friend class ImageVtkReadAccessor;
    friend class ImageVtkWriteAccessor;
    friend class ImageReadAccessor;
    friend class ImageWriteAccessor;
 
  public:
    mitkClassMacro(Image, SlicedData);
 
    itkFactorylessNewMacro(Self);
 
    itkCloneMacro(Self);
 
    typedef itk::SmartPointer<ImageDataItem> ImageDataItemPointer;
 
    typedef itk::Statistics::Histogram<double> HistogramType;
 
    typedef mitk::ImageStatisticsHolder *StatisticsHolderPointer;
 
    enum ImportMemoryManagementType
    {
      CopyMemory, 
      ManageMemory, 
      ReferenceMemory, 
      DontManageMemory = ReferenceMemory 
    };
 
    typedef std::vector<ImageDataItemPointer> ImageDataItemPointerArray;
 
  public:
    const mitk::PixelType GetPixelType(int n = 0) const;
 
    unsigned int GetDimension() const;
 
    unsigned int GetDimension(int i) const;
 
  public:
    virtual vtkImageData *GetVtkImageData(int t = 0, int n = 0);
 
    virtual const vtkImageData *GetVtkImageData(int t = 0, int n = 0) const;
 
    bool IsSliceSet(int s = 0, int t = 0, int n = 0) const override;
 
    bool IsVolumeSet(int t = 0, int n = 0) const override;
 
    bool IsChannelSet(int n = 0) const override;
 
    virtual bool SetSlice(const void *data, int s = 0, int t = 0, int n = 0);
 
    virtual bool SetVolume(const void *data, int t = 0, int n = 0);
 
    virtual bool SetChannel(const void *data, int n = 0);
 
    virtual bool SetImportSlice(
      void *data, int s = 0, int t = 0, int n = 0, ImportMemoryManagementType importMemoryManagement = CopyMemory);
 
    virtual bool SetImportVolume(void *data,
                                 int t = 0,
                                 int n = 0,
                                 ImportMemoryManagementType importMemoryManagement = CopyMemory);
 
    virtual bool SetImportVolume(const void *const_data, int t = 0, int n = 0);
 
    virtual void AllocateZeroedVolume(int t = 0, int n = 0);
 
    virtual bool SetImportChannel(void *data,
                                  int n = 0,
                                  ImportMemoryManagementType importMemoryManagement = CopyMemory);
 
    virtual void Initialize(const mitk::PixelType &type,
                            unsigned int dimension,
                            const unsigned int *dimensions,
                            unsigned int channels = 1);
 
    virtual void Initialize(const mitk::PixelType &type,
                            const mitk::BaseGeometry &geometry,
                            unsigned int channels = 1,
                            int tDim = 1);
 
    virtual void Initialize(const mitk::PixelType &type,
                            const mitk::TimeGeometry &geometry,
                            unsigned int channels = 1,
                            int tDim = -1);
 
    virtual void Initialize(const mitk::PixelType &type,
                            int sDim,
                            const mitk::PlaneGeometry &geometry2d,
                            unsigned int channels = 1,
                            int tDim = 1);
 
    virtual void Initialize(const mitk::Image *image);
 
    virtual void Initialize(const mitk::ImageDescriptor::Pointer inDesc);
 
    virtual void Initialize(vtkImageData *vtkimagedata, int channels = 1, int tDim = -1, int sDim = -1, int pDim = -1);
 
    template <typename itkImageType>
    void InitializeByItk(const itkImageType *itkimage, int channels = 1, int tDim = -1, int sDim = -1)
    {
      if (itkimage == nullptr)
        return;
 
      MITK_DEBUG << "Initializing MITK image from ITK image.";
      // build array with dimensions in each direction with at least 4 entries
      m_Dimension = itkimage->GetImageDimension();
      unsigned int i, *tmpDimensions = new unsigned int[m_Dimension > 4 ? m_Dimension : 4];
      for (i = 0; i < m_Dimension; ++i)
        tmpDimensions[i] = itkimage->GetLargestPossibleRegion().GetSize().GetSize()[i];
      if (m_Dimension < 4)
      {
        unsigned int *p;
        for (i = 0, p = tmpDimensions + m_Dimension; i < 4 - m_Dimension; ++i, ++p)
          *p = 1;
      }
 
      // overwrite number of slices if sDim is set
      if ((m_Dimension > 2) && (sDim >= 0))
        tmpDimensions[2] = sDim;
      // overwrite number of time points if tDim is set
      if ((m_Dimension > 3) && (tDim >= 0))
        tmpDimensions[3] = tDim;
 
      // rough initialization of Image
      // mitk::PixelType importType = ImportItkPixelType( itkimage::PixelType );
 
      Initialize(
        MakePixelType<itkImageType>(itkimage->GetNumberOfComponentsPerPixel()), m_Dimension, tmpDimensions, channels);
      const typename itkImageType::SpacingType &itkspacing = itkimage->GetSpacing();
 
      MITK_DEBUG << "ITK spacing " << itkspacing;
      // access spacing of itk::Image
      Vector3D spacing;
      FillVector3D(spacing, itkspacing[0], 1.0, 1.0);
      if (m_Dimension >= 2)
        spacing[1] = itkspacing[1];
      if (m_Dimension >= 3)
        spacing[2] = itkspacing[2];
 
      // access origin of itk::Image
      Point3D origin;
      const typename itkImageType::PointType &itkorigin = itkimage->GetOrigin();
      MITK_DEBUG << "ITK origin " << itkorigin;
      FillVector3D(origin, itkorigin[0], 0.0, 0.0);
      if (m_Dimension >= 2)
        origin[1] = itkorigin[1];
      if (m_Dimension >= 3)
        origin[2] = itkorigin[2];
 
      // access direction of itk::Imagm_PixelType = new mitk::PixelType(type);e and include spacing
      const typename itkImageType::DirectionType &itkdirection = itkimage->GetDirection();
      MITK_DEBUG << "ITK direction " << itkdirection;
      mitk::Matrix3D matrix;
      matrix.SetIdentity();
      unsigned int j, itkDimMax3 = (m_Dimension >= 3 ? 3 : m_Dimension);
      // check if spacing has no zero entry and itkdirection has no zero columns
      bool itkdirectionOk = true;
      mitk::ScalarType columnSum;
      for (j = 0; j < itkDimMax3; ++j)
      {
        columnSum = 0.0;
        for (i = 0; i < itkDimMax3; ++i)
        {
          columnSum += fabs(itkdirection[i][j]);
        }
        if (columnSum < mitk::eps)
        {
          itkdirectionOk = false;
        }
        if ((spacing[j] < -mitk::eps) // (normally sized) negative value
            &&
            (j == 2) && (m_Dimensions[2] == 1))
        {
          // Negative spacings can occur when reading single DICOM slices with ITK via GDCMIO
          // In these cases spacing is not determined by ITK correctly (because it distinguishes correctly
          // between slice thickness and inter slice distance -- slice distance is meaningless for
          // single slices).
          // I experienced that ITK produced something meaningful nonetheless because it is
          // evaluating the tag "(0018,0088) Spacing between slices" as a fallback. This tag is not
          // reliable (http://www.itk.org/pipermail/insight-users/2005-September/014711.html)
          // but gives at least a hint.
          // In real world cases I experienced that this tag contained the correct inter slice distance
          // with a negative sign, so we just invert such negative spacings.
          MITK_WARN << "Illegal value of itk::Image::GetSpacing()[" << j << "]=" << spacing[j]
                    << ". Using inverted value " << -spacing[j];
          spacing[j] = -spacing[j];
        }
        else if (spacing[j] < mitk::eps) // value near zero
        {
          MITK_ERROR << "Illegal value of itk::Image::GetSpacing()[" << j << "]=" << spacing[j]
                     << ". Using 1.0 instead.";
          spacing[j] = 1.0;
        }
      }
      if (itkdirectionOk == false)
      {
        MITK_ERROR << "Illegal matrix returned by itk::Image::GetDirection():" << itkdirection
                   << " Using identity instead.";
        for (i = 0; i < itkDimMax3; ++i)
          for (j = 0; j < itkDimMax3; ++j)
            if (i == j)
              matrix[i][j] = spacing[j];
            else
              matrix[i][j] = 0.0;
      }
      else
      {
        for (i = 0; i < itkDimMax3; ++i)
          for (j = 0; j < itkDimMax3; ++j)
            matrix[i][j] = itkdirection[i][j] * spacing[j];
      }
 
      // re-initialize PlaneGeometry with origin and direction
      PlaneGeometry *planeGeometry = static_cast<PlaneGeometry *>(GetSlicedGeometry(0)->GetPlaneGeometry(0));
      planeGeometry->SetOrigin(origin);
      planeGeometry->GetIndexToWorldTransform()->SetMatrix(matrix);
 
      // re-initialize SlicedGeometry3D
      SlicedGeometry3D *slicedGeometry = GetSlicedGeometry(0);
      slicedGeometry->InitializeEvenlySpaced(planeGeometry, m_Dimensions[2]);
      slicedGeometry->SetSpacing(spacing);
 
      // re-initialize TimeGeometry
      ProportionalTimeGeometry::Pointer timeGeometry = ProportionalTimeGeometry::New();
      timeGeometry->Initialize(slicedGeometry, m_Dimensions[3]);
      SetTimeGeometry(timeGeometry);
 
      // clean-up
      delete[] tmpDimensions;
 
      this->Initialize();
    }
 
    virtual bool IsValidSlice(int s = 0, int t = 0, int n = 0) const;
 
    virtual bool IsValidVolume(int t = 0, int n = 0) const;
 
    virtual bool IsValidChannel(int n = 0) const;
 
    bool IsRotated() const;
 
    unsigned int *GetDimensions() const;
 
    ImageDescriptor::Pointer GetImageDescriptor() const { return m_ImageDescriptor; }
 
    ChannelDescriptor GetChannelDescriptor(int id = 0) const { return m_ImageDescriptor->GetChannelDescriptor(id); }
 
    void SetGeometry(BaseGeometry *aGeometry3D) override;
 
    virtual ImageDataItemPointer GetSliceData(int s = 0,
                                              int t = 0,
                                              int n = 0,
                                              void *data = nullptr,
                                              ImportMemoryManagementType importMemoryManagement = CopyMemory) const;
 
    virtual ImageDataItemPointer GetVolumeData(int t = 0,
                                               int n = 0,
                                               void *data = nullptr,
                                               ImportMemoryManagementType importMemoryManagement = CopyMemory) const;
 
    virtual ImageDataItemPointer GetChannelData(int n = 0,
                                                void *data = nullptr,
                                                ImportMemoryManagementType importMemoryManagement = CopyMemory) const;
 
    void *GetData(int t = 0, int n = 0);
 
    const void *GetData(int t = 0, int n = 0) const;
 
    StatisticsHolderPointer GetStatistics() const { return m_ImageStatistics; }
 
  protected:
    mitkCloneMacro(Self);
 
    typedef std::lock_guard<std::mutex> MutexHolder;
 
    int GetSliceIndex(int s = 0, int t = 0, int n = 0) const;
 
    int GetVolumeIndex(int t = 0, int n = 0) const;
 
    void ComputeOffsetTable();
 
    virtual bool IsValidTimeStep(int t) const;
 
    void Expand(unsigned int timeSteps) override;
 
    virtual ImageDataItemPointer AllocateSliceData(
      int s = 0,
      int t = 0,
      int n = 0,
      void *data = nullptr,
      ImportMemoryManagementType importMemoryManagement = CopyMemory) const;
 
    virtual ImageDataItemPointer AllocateVolumeData(
      int t = 0, int n = 0, void *data = nullptr, ImportMemoryManagementType importMemoryManagement = CopyMemory) const;
 
    virtual ImageDataItemPointer AllocateChannelData(
      int n = 0, void *data = nullptr, ImportMemoryManagementType importMemoryManagement = CopyMemory) const;
 
    Image();
 
    Image(const Image &other);
 
    ~Image() override;
 
    void Clear() override;
 
    void Initialize() override;
 
    void PrintSelf(std::ostream &os, itk::Indent indent) const override;
 
    mutable ImageDataItemPointerArray m_Channels;
    mutable ImageDataItemPointerArray m_Volumes;
    mutable ImageDataItemPointerArray m_Slices;
    mutable std::mutex m_ImageDataArraysLock;
 
    unsigned int m_Dimension;
 
    unsigned int *m_Dimensions;
 
    ImageDescriptor::Pointer m_ImageDescriptor;
 
    size_t *m_OffsetTable;
    ImageDataItemPointer m_CompleteData;
 
    // Image statistics Holder replaces the former implementation directly inside this class
    friend class ImageStatisticsHolder;
    StatisticsHolderPointer m_ImageStatistics;
 
  private:
    ImageDataItemPointer GetSliceData_unlocked(
      int s, int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const;
    ImageDataItemPointer GetVolumeData_unlocked(int t,
                                                int n,
                                                void *data,
                                                ImportMemoryManagementType importMemoryManagement) const;
    ImageDataItemPointer GetChannelData_unlocked(int n,
                                                 void *data,
                                                 ImportMemoryManagementType importMemoryManagement) const;
 
    ImageDataItemPointer AllocateSliceData_unlocked(
      int s, int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const;
    ImageDataItemPointer AllocateVolumeData_unlocked(int t,
                                                     int n,
                                                     void *data,
                                                     ImportMemoryManagementType importMemoryManagement) const;
    ImageDataItemPointer AllocateChannelData_unlocked(int n,
                                                      void *data,
                                                      ImportMemoryManagementType importMemoryManagement) const;
 
    bool IsSliceSet_unlocked(int s, int t, int n) const;
    bool IsVolumeSet_unlocked(int t, int n) const;
    bool IsChannelSet_unlocked(int n) const;
 
    mutable std::vector<ImageAccessorBase *> m_Readers;
    mutable std::vector<ImageAccessorBase *> m_Writers;
    mutable std::vector<ImageAccessorBase *> m_VtkReaders;
 
    mutable std::mutex m_ReadWriteLock;
    mutable std::mutex m_VtkReadersLock;
  };
 
 
  using ImageDimensionVectorType = std::vector<unsigned int>;
  MITKCORE_EXPORT ImageDimensionVectorType DetermineImageDimensionsFromTimeGeometry(const TimeGeometry* timeGeometry);
 
  MITKCORE_EXPORT bool Equal(const mitk::Image &leftHandSide,
                             const mitk::Image &rightHandSide,
                             ScalarType eps,
                             bool verbose);
 
} // namespace mitk
 
#endif