diff --git a/Modules/Core/include/mitkImage.h b/Modules/Core/include/mitkImage.h index fccd562f34..aeba9becfd 100644 --- a/Modules/Core/include/mitkImage.h +++ b/Modules/Core/include/mitkImage.h @@ -1,784 +1,792 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKIMAGE_H_HEADER_INCLUDED_C1C2FCD2 #define MITKIMAGE_H_HEADER_INCLUDED_C1C2FCD2 #include "mitkBaseData.h" #include "mitkImageAccessorBase.h" #include "mitkImageDataItem.h" #include "mitkImageDescriptor.h" #include "mitkImageVtkAccessor.h" #include "mitkLevelWindow.h" #include "mitkPlaneGeometry.h" #include "mitkSlicedData.h" #include #include // DEPRECATED #include #ifndef __itkHistogram_h #include #endif class vtkImageData; namespace itk { template class MutexLockHolder; } namespace mitk { class SubImageSelector; class ImageTimeSelector; class ImageStatisticsHolder; //##Documentation //## @brief Image class for storing images //## //## Can be asked for header information, the data vector, //## the mitkIpPicDescriptor struct or vtkImageData objects. If not the complete //## data is required, the appropriate SubImageSelector class should be used //## for access. //## Image organizes sets of slices (s x 2D), volumes (t x 3D) and channels (n //## x ND). Channels are for different kind of data, e.g., morphology in //## channel 0, velocities in channel 1. All channels must have the same Geometry! In //## particular, the dimensions of all channels are the same, only the pixel-type //## may differ between channels. //## //## For importing ITK images use of mitk::ITKImageImport is recommended, see //## \ref Adaptor. //## //## For ITK v3.8 and older: Converting coordinates from the ITK physical //## coordinate system (which does not support rotated images) to the MITK world //## coordinate system should be performed via the BaseGeometry of the Image, see //## BaseGeometry::WorldToItkPhysicalPoint. //## //## For more information, see \ref MitkImagePage . //## @ingroup Data class MITKCORE_EXPORT Image : public SlicedData { friend class SubImageSelector; friend class ImageAccessorBase; friend class ImageVtkAccessor; friend class ImageVtkReadAccessor; friend class ImageVtkWriteAccessor; friend class ImageReadAccessor; friend class ImageWriteAccessor; public: mitkClassMacro(Image, SlicedData); itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Smart Pointer type to a ImageDataItem. */ typedef itk::SmartPointer ImageDataItemPointer; typedef itk::Statistics::Histogram HistogramType; typedef mitk::ImageStatisticsHolder *StatisticsHolderPointer; //## @param ImportMemoryManagementType This parameter is evaluated when setting new data to an image. //## The different options are: //## CopyMemory: Data to be set is copied and assigned to a new memory block. Data memory block will be freed on // deletion of mitk::Image. //## MamageMemory: Data to be set will be referenced, and Data memory block will be freed on deletion of //mitk::Image. //## Reference Memory: Data to be set will be referenced, but Data memory block will not be freed on deletion of // mitk::Image. //## DontManageMemory = ReferenceMemory. enum ImportMemoryManagementType { CopyMemory, ManageMemory, ReferenceMemory, DontManageMemory = ReferenceMemory }; //##Documentation //## @brief Vector container of SmartPointers to ImageDataItems; //## Class is only for internal usage to allow convenient access to all slices over iterators; //## See documentation of ImageDataItem for details. typedef std::vector ImageDataItemPointerArray; public: //##Documentation //## @brief Returns the PixelType of channel @a n. const mitk::PixelType GetPixelType(int n = 0) const; //##Documentation //## @brief Get dimension of the image //## unsigned int GetDimension() const; //##Documentation //## @brief Get the size of dimension @a i (e.g., i=0 results in the number of pixels in x-direction). //## //## @sa GetDimensions() unsigned int GetDimension(int i) const; /** @brief Get the data vector of the complete image, i.e., of all channels linked together. If you only want to access a slice, volume at a specific time or single channel use one of the SubImageSelector classes. \deprecatedSince{2012_09} Please use image accessors instead: See Doxygen/Related-Pages/Concepts/Image. This method can be replaced by ImageWriteAccessor::GetData() or ImageReadAccessor::GetData() */ DEPRECATED(virtual void *GetData()); public: /** @brief Get the pixel value at one specific index position. The pixel type is always being converted to double. \deprecatedSince{2012_09} Please use image accessors instead: See Doxygen/Related-Pages/Concepts/Image. This method can be replaced by a method from ImagePixelWriteAccessor or ImagePixelReadAccessor */ DEPRECATED(double GetPixelValueByIndex(const itk::Index<3> &position, unsigned int timestep = 0, unsigned int component = 0)); /** @brief Get the pixel value at one specific world position. The pixel type is always being converted to double. \deprecatedSince{2012_09} Please use image accessors instead: See Doxygen/Related-Pages/Concepts/Image. This method can be replaced by a method from ImagePixelWriteAccessor or ImagePixelReadAccessor */ DEPRECATED(double GetPixelValueByWorldCoordinate(const mitk::Point3D &position, unsigned int timestep = 0, unsigned int component = 0)); //##Documentation //## @brief Get a volume at a specific time @a t of channel @a n as a vtkImageData. virtual vtkImageData *GetVtkImageData(int t = 0, int n = 0); virtual const vtkImageData *GetVtkImageData(int t = 0, int n = 0) const; //##Documentation //## @brief Get the complete image, i.e., all channels linked together, as a @a mitkIpPicDescriptor. //## //## If you only want to access a slice, volume at a specific time or single channel //## use one of the SubImageSelector classes. // virtual mitkIpPicDescriptor* GetPic(); //##Documentation //## @brief Check whether slice @a s at time @a t in channel @a n is set virtual bool IsSliceSet(int s = 0, int t = 0, int n = 0) const override; //##Documentation //## @brief Check whether volume at time @a t in channel @a n is set virtual bool IsVolumeSet(int t = 0, int n = 0) const override; //##Documentation //## @brief Check whether the channel @a n is set virtual bool IsChannelSet(int n = 0) const override; //##Documentation //## @brief Set @a data as slice @a s at time @a t in channel @a n. It is in //## the responsibility of the caller to ensure that the data vector @a data //## is really a slice (at least is not smaller than a slice), since there is //## no chance to check this. //## //## The data is copied to an array managed by the image. If the image shall //## reference the data, use SetImportSlice with ImportMemoryManagementType //## set to ReferenceMemory. For importing ITK images use of mitk:: //## ITKImageImport is recommended. //## @sa SetPicSlice, SetImportSlice, SetImportVolume virtual bool SetSlice(const void *data, int s = 0, int t = 0, int n = 0); //##Documentation //## @brief Set @a data as volume at time @a t in channel @a n. It is in //## the responsibility of the caller to ensure that the data vector @a data //## is really a volume (at least is not smaller than a volume), since there is //## no chance to check this. //## //## The data is copied to an array managed by the image. If the image shall //## reference the data, use SetImportVolume with ImportMemoryManagementType //## set to ReferenceMemory. For importing ITK images use of mitk:: //## ITKImageImport is recommended. //## @sa SetPicVolume, SetImportVolume virtual bool SetVolume(const void *data, int t = 0, int n = 0); //##Documentation //## @brief Set @a data in channel @a n. It is in //## the responsibility of the caller to ensure that the data vector @a data //## is really a channel (at least is not smaller than a channel), since there is //## no chance to check this. //## //## The data is copied to an array managed by the image. If the image shall //## reference the data, use SetImportChannel with ImportMemoryManagementType //## set to ReferenceMemory. For importing ITK images use of mitk:: //## ITKImageImport is recommended. //## @sa SetPicChannel, SetImportChannel virtual bool SetChannel(const void *data, int n = 0); //##Documentation //## @brief Set @a data as slice @a s at time @a t in channel @a n. It is in //## the responsibility of the caller to ensure that the data vector @a data //## is really a slice (at least is not smaller than a slice), since there is //## no chance to check this. //## //## The data is managed according to the parameter \a importMemoryManagement. //## @sa SetPicSlice virtual bool SetImportSlice( void *data, int s = 0, int t = 0, int n = 0, ImportMemoryManagementType importMemoryManagement = CopyMemory); //##Documentation //## @brief Set @a data as volume at time @a t in channel @a n. It is in //## the responsibility of the caller to ensure that the data vector @a data //## is really a volume (at least is not smaller than a volume), since there is //## no chance to check this. //## //## The data is managed according to the parameter \a importMemoryManagement. //## @sa SetPicVolume virtual bool SetImportVolume(void *data, int t = 0, int n = 0, ImportMemoryManagementType importMemoryManagement = CopyMemory); //##Documentation //## @brief Set @a data in channel @a n. It is in //## the responsibility of the caller to ensure that the data vector @a data //## is really a channel (at least is not smaller than a channel), since there is //## no chance to check this. //## //## The data is managed according to the parameter \a importMemoryManagement. //## @sa SetPicChannel virtual bool SetImportChannel(void *data, int n = 0, ImportMemoryManagementType importMemoryManagement = CopyMemory); //##Documentation //## initialize new (or re-initialize) image information //## @warning Initialize() by pic assumes a plane, evenly spaced geometry starting at (0,0,0). virtual void Initialize(const mitk::PixelType &type, unsigned int dimension, const unsigned int *dimensions, unsigned int channels = 1); //##Documentation //## initialize new (or re-initialize) image information by a BaseGeometry //## //## @param tDim defines the number of time steps for which the Image should be initialized virtual void Initialize(const mitk::PixelType &type, const mitk::BaseGeometry &geometry, unsigned int channels = 1, int tDim = 1); /** * initialize new (or re-initialize) image information by a TimeGeometry * * @param tDim defines the number of time steps for which the Image should be initialized * \deprecatedSince{2013_09} Please use TimeGeometry instead of TimeSlicedGeometry. For more information see * http://www.mitk.org/Development/Refactoring%20of%20the%20Geometry%20Classes%20-%20Part%201 */ DEPRECATED(virtual void Initialize(const mitk::PixelType & /*type*/, const mitk::TimeSlicedGeometry * /*geometry*/, unsigned int /*channels = 1*/, int /*tDim=1*/)) { } /** * \brief Initialize new (or re-initialize) image information by a TimeGeometry * * \param tDim override time dimension if the value is bigger than 0 (Default -1) */ virtual void Initialize(const mitk::PixelType &type, const mitk::TimeGeometry &geometry, unsigned int channels = 1, int tDim = -1); //##Documentation //## initialize new (or re-initialize) image information by a PlaneGeometry and number of slices //## //## Initializes the bounding box according to the width/height of the //## PlaneGeometry and @a sDim via SlicedGeometry3D::InitializeEvenlySpaced. //## The spacing is calculated from the PlaneGeometry. //## \sa SlicedGeometry3D::InitializeEvenlySpaced + //## \deprecatedSince{2016_11} Use a left-handed or right-handed PlaneGeometry to define the + //## direction of the image stack instead of the flipped parameter + DEPRECATED(virtual void Initialize(const mitk::PixelType &type, + int sDim, + const mitk::PlaneGeometry &geometry2d, + bool flipped, + unsigned int channels = 1, + int tDim = 1)); + virtual void Initialize(const mitk::PixelType &type, int sDim, const mitk::PlaneGeometry &geometry2d, - bool flipped = false, unsigned int channels = 1, int tDim = 1); //##Documentation //## initialize new (or re-initialize) image information by another //## mitk-image. //## Only the header is used, not the data vector! //## virtual void Initialize(const mitk::Image *image); virtual void Initialize(const mitk::ImageDescriptor::Pointer inDesc); //##Documentation //## initialize new (or re-initialize) image information by @a pic. //## Dimensions and @a Geometry3D /@a PlaneGeometry are set according //## to the tags in @a pic. //## Only the header is used, not the data vector! Use SetPicVolume(pic) //## to set the data vector. //## //## @param tDim override time dimension (@a n[3]) in @a pic (if >0) //## @param sDim override z-space dimension (@a n[2]) in @a pic (if >0) //## @warning Initialize() by pic assumes a plane, evenly spaced geometry starting at (0,0,0). // virtual void Initialize(const mitkIpPicDescriptor* pic, int channels = 1, int tDim = -1, int sDim = -1); //##Documentation //## initialize new (or re-initialize) image information by @a vtkimagedata, //## a vtk-image. //## Only the header is used, not the data vector! Use //## SetVolume(vtkimage->GetScalarPointer()) to set the data vector. //## //## @param tDim override time dimension in @a vtkimagedata (if >0 and <) //## @param sDim override z-space dimension in @a vtkimagedata (if >0 and <) //## @param pDim override y-space dimension in @a vtkimagedata (if >0 and <) virtual void Initialize(vtkImageData *vtkimagedata, int channels = 1, int tDim = -1, int sDim = -1, int pDim = -1); //##Documentation //## initialize new (or re-initialize) image information by @a itkimage, //## a templated itk-image. //## Only the header is used, not the data vector! Use //## SetVolume(itkimage->GetBufferPointer()) to set the data vector. //## //## @param tDim override time dimension in @a itkimage (if >0 and <) //## @param sDim override z-space dimension in @a itkimage (if >0 and <) template 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(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 determind 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 is 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(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(); } //##Documentation //## @brief Check whether slice @a s at time @a t in channel @a n is valid, i.e., //## is (or can be) inside of the image virtual bool IsValidSlice(int s = 0, int t = 0, int n = 0) const; //##Documentation //## @brief Check whether volume at time @a t in channel @a n is valid, i.e., //## is (or can be) inside of the image virtual bool IsValidVolume(int t = 0, int n = 0) const; //##Documentation //## @brief Check whether the channel @a n is valid, i.e., //## is (or can be) inside of the image virtual bool IsValidChannel(int n = 0) const; //##Documentation //## @brief Returns true if an image is rotated, i.e. its geometry's //## transformation matrix has nonzero elements besides the diagonal. //## Non-diagonal elements are checked if larger then 1/1000 of the matrix' trace. bool IsRotated() const; //##Documentation //## @brief Get the sizes of all dimensions as an integer-array. //## //## @sa GetDimension(int i); unsigned int *GetDimensions() const; ImageDescriptor::Pointer GetImageDescriptor() const { return m_ImageDescriptor; } ChannelDescriptor GetChannelDescriptor(int id = 0) const { return m_ImageDescriptor->GetChannelDescriptor(id); } /** \brief Sets a geometry to an image. */ virtual void SetGeometry(BaseGeometry *aGeometry3D) override; /** * @warning for internal use only */ virtual ImageDataItemPointer GetSliceData(int s = 0, int t = 0, int n = 0, void *data = nullptr, ImportMemoryManagementType importMemoryManagement = CopyMemory) const; /** * @warning for internal use only */ virtual ImageDataItemPointer GetVolumeData(int t = 0, int n = 0, void *data = nullptr, ImportMemoryManagementType importMemoryManagement = CopyMemory) const; /** * @warning for internal use only */ virtual ImageDataItemPointer GetChannelData(int n = 0, void *data = nullptr, ImportMemoryManagementType importMemoryManagement = CopyMemory) const; /** \brief (DEPRECATED) Get the minimum for scalar images */ DEPRECATED(ScalarType GetScalarValueMin(int t = 0) const); /** \brief (DEPRECATED) Get the maximum for scalar images \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetScalarValueMax(int t = 0) const); /** \brief (DEPRECATED) Get the second smallest value for scalar images \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetScalarValue2ndMin(int t = 0) const); /** \brief (DEPRECATED) Get the smallest value for scalar images, but do not recompute it first \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetScalarValueMinNoRecompute(unsigned int t = 0) const); /** \brief (DEPRECATED) Get the second smallest value for scalar images, but do not recompute it first \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetScalarValue2ndMinNoRecompute(unsigned int t = 0) const); /** \brief (DEPRECATED) Get the second largest value for scalar images \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetScalarValue2ndMax(int t = 0) const); /** \brief (DEPRECATED) Get the largest value for scalar images, but do not recompute it first \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetScalarValueMaxNoRecompute(unsigned int t = 0) const); /** \brief (DEPRECATED) Get the second largest value for scalar images, but do not recompute it first \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetScalarValue2ndMaxNoRecompute(unsigned int t = 0) const); /** \brief (DEPRECATED) Get the count of voxels with the smallest scalar value in the dataset \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetCountOfMinValuedVoxels(int t = 0) const); /** \brief (DEPRECATED) Get the count of voxels with the largest scalar value in the dataset \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(ScalarType GetCountOfMaxValuedVoxels(int t = 0) const); /** \brief (DEPRECATED) Get the count of voxels with the largest scalar value in the dataset \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(unsigned int GetCountOfMaxValuedVoxelsNoRecompute(unsigned int t = 0) const); /** \brief (DEPRECATED) Get the count of voxels with the smallest scalar value in the dataset \warning This method is deprecated and will not be available in the future. Use the \a GetStatistics instead */ DEPRECATED(unsigned int GetCountOfMinValuedVoxelsNoRecompute(unsigned int t = 0) const); /** \brief Returns a pointer to the ImageStatisticsHolder object that holds all statistics information for the image. All Get-methods for statistics properties formerly accessible directly from an Image object are now moved to the new \a ImageStatisticsHolder object. */ StatisticsHolderPointer GetStatistics() const { return m_ImageStatistics; } protected: mitkCloneMacro(Self); typedef itk::MutexLockHolder 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; virtual 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); virtual ~Image(); virtual void Clear() override; //## @warning Has to be called by every Initialize method! virtual void Initialize() override; virtual void PrintSelf(std::ostream &os, itk::Indent indent) const override; mutable ImageDataItemPointerArray m_Channels; mutable ImageDataItemPointerArray m_Volumes; mutable ImageDataItemPointerArray m_Slices; mutable itk::SimpleFastMutexLock 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; /** Stores all existing ImageReadAccessors */ mutable std::vector m_Readers; /** Stores all existing ImageWriteAccessors */ mutable std::vector m_Writers; /** Stores all existing ImageVtkAccessors */ mutable std::vector m_VtkReaders; /** A mutex, which needs to be locked to manage m_Readers and m_Writers */ itk::SimpleFastMutexLock m_ReadWriteLock; /** A mutex, which needs to be locked to manage m_VtkReaders */ itk::SimpleFastMutexLock m_VtkReadersLock; }; /** * @brief Equal A function comparing two images for beeing equal in meta- and imagedata * @warning This method is deprecated and will not be available in the future. Use the \a bool mitk::Equal(const * mitk::Image& i1, const mitk::Image& i2) instead. * * @ingroup MITKTestingAPI * * Following aspects are tested for equality: * - dimension of the images * - size of the images * - pixel type * - pixel values : pixel values are expected to be identical at each position ( for other options see * mitk::CompareImageFilter ) * * @param rightHandSide An image to be compared * @param leftHandSide An image to be compared * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return true, if all subsequent comparisons are true, false otherwise */ DEPRECATED(MITKCORE_EXPORT bool Equal( const mitk::Image *leftHandSide, const mitk::Image *rightHandSide, ScalarType eps, bool verbose)); /** * @brief Equal A function comparing two images for beeing equal in meta- and imagedata * * @ingroup MITKTestingAPI * * Following aspects are tested for equality: * - dimension of the images * - size of the images * - pixel type * - pixel values : pixel values are expected to be identical at each position ( for other options see * mitk::CompareImageFilter ) * * @param rightHandSide An image to be compared * @param leftHandSide An image to be compared * @param eps Tolarence for comparison. You can use mitk::eps in most cases. * @param verbose Flag indicating if the user wants detailed console output or not. * @return true, if all subsequent comparisons are true, false otherwise */ MITKCORE_EXPORT bool Equal(const mitk::Image &leftHandSide, const mitk::Image &rightHandSide, ScalarType eps, bool verbose); } // namespace mitk #endif /* MITKIMAGE_H_HEADER_INCLUDED_C1C2FCD2 */ diff --git a/Modules/Core/src/DataManagement/mitkImage.cpp b/Modules/Core/src/DataManagement/mitkImage.cpp index 69d74a2204..a49682db41 100644 --- a/Modules/Core/src/DataManagement/mitkImage.cpp +++ b/Modules/Core/src/DataManagement/mitkImage.cpp @@ -1,1537 +1,1546 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // MITK #include "mitkImage.h" #include "mitkCompareImageDataFilter.h" #include "mitkImageStatisticsHolder.h" #include "mitkImageVtkReadAccessor.h" #include "mitkImageVtkWriteAccessor.h" #include "mitkPixelTypeMultiplex.h" #include // VTK #include // ITK #include // Other #include #define FILL_C_ARRAY(_arr, _size, _value) \ for (unsigned int i = 0u; i < _size; i++) \ \ { \ _arr[i] = _value; \ } mitk::Image::Image() : m_Dimension(0), m_Dimensions(nullptr), m_ImageDescriptor(nullptr), m_OffsetTable(nullptr), m_CompleteData(nullptr), m_ImageStatistics(nullptr) { m_Dimensions = new unsigned int[MAX_IMAGE_DIMENSIONS]; FILL_C_ARRAY(m_Dimensions, MAX_IMAGE_DIMENSIONS, 0u); m_Initialized = false; } mitk::Image::Image(const Image &other) : SlicedData(other), m_Dimension(0), m_Dimensions(nullptr), m_ImageDescriptor(nullptr), m_OffsetTable(nullptr), m_CompleteData(nullptr), m_ImageStatistics(nullptr) { m_Dimensions = new unsigned int[MAX_IMAGE_DIMENSIONS]; FILL_C_ARRAY(m_Dimensions, MAX_IMAGE_DIMENSIONS, 0u); this->Initialize(other.GetPixelType(), other.GetDimension(), other.GetDimensions()); // Since the above called "Initialize" method doesn't take the geometry into account we need to set it // here manually TimeGeometry::Pointer cloned = other.GetTimeGeometry()->Clone(); this->SetTimeGeometry(cloned.GetPointer()); if (this->GetDimension() > 3) { const unsigned int time_steps = this->GetDimension(3); for (unsigned int i = 0u; i < time_steps; ++i) { ImageDataItemPointer volume = const_cast(other).GetVolumeData(i); this->SetVolume(volume->GetData(), i); } } else { ImageDataItemPointer volume = const_cast(other).GetVolumeData(0); this->SetVolume(volume->GetData(), 0); } } mitk::Image::~Image() { this->Clear(); m_ReferenceCount = 3; m_ReferenceCount = 0; delete[] m_OffsetTable; delete m_ImageStatistics; } const mitk::PixelType mitk::Image::GetPixelType(int n) const { return this->m_ImageDescriptor->GetChannelTypeById(n); } unsigned int mitk::Image::GetDimension() const { return m_Dimension; } unsigned int mitk::Image::GetDimension(int i) const { if ((i >= 0) && (i < (int)m_Dimension)) return m_Dimensions[i]; return 1; } void *mitk::Image::GetData() { if (m_Initialized == false) { if (GetSource().IsNull()) return nullptr; if (GetSource()->Updating() == false) GetSource()->UpdateOutputInformation(); } m_CompleteData = GetChannelData(); // update channel's data // if data was not available at creation point, the m_Data of channel descriptor is NULL // if data present, it won't be overwritten m_ImageDescriptor->GetChannelDescriptor(0).SetData(m_CompleteData->GetData()); return m_CompleteData->GetData(); } template void AccessPixel(const mitk::PixelType ptype, void *data, const unsigned int offset, double &value) { value = 0.0; if (data == nullptr) return; if (ptype.GetBpe() != 24) { value = (double)(((T *)data)[offset]); } else { const unsigned int rgboffset = offset; double returnvalue = (((T *)data)[rgboffset]); returnvalue += (((T *)data)[rgboffset + 1]); returnvalue += (((T *)data)[rgboffset + 2]); value = returnvalue; } } double mitk::Image::GetPixelValueByIndex(const itk::Index<3> &position, unsigned int timestep, unsigned int component) { double value = 0; if (this->GetTimeSteps() < timestep) { timestep = this->GetTimeSteps(); } value = 0.0; const unsigned int *imageDims = this->m_ImageDescriptor->GetDimensions(); const mitk::PixelType ptype = this->m_ImageDescriptor->GetChannelTypeById(0); // Comparison ?>=0 not needed since all position[i] and timestep are unsigned int // (position[0]>=0 && position[1] >=0 && position[2]>=0 && timestep>=0) // bug-11978 : we still need to catch index with negative values if (position[0] < 0 || position[1] < 0 || position[2] < 0) { MITK_WARN << "Given position (" << position << ") is out of image range, returning 0."; } // check if the given position is inside the index range of the image, the 3rd dimension needs to be compared only if // the dimension is not 0 else if ((unsigned int)position[0] >= imageDims[0] || (unsigned int)position[1] >= imageDims[1] || (imageDims[2] && (unsigned int)position[2] >= imageDims[2])) { MITK_WARN << "Given position (" << position << ") is out of image range, returning 0."; } else { const unsigned int offset = component + ptype.GetNumberOfComponents() * (position[0] + position[1] * imageDims[0] + position[2] * imageDims[0] * imageDims[1] + timestep * imageDims[0] * imageDims[1] * imageDims[2]); mitkPixelTypeMultiplex3(AccessPixel, ptype, this->GetData(), offset, value); } return value; } double mitk::Image::GetPixelValueByWorldCoordinate(const mitk::Point3D &position, unsigned int timestep, unsigned int component) { double value = 0.0; if (this->GetTimeSteps() < timestep) { timestep = this->GetTimeSteps(); } itk::Index<3> itkIndex; this->GetGeometry()->WorldToIndex(position, itkIndex); value = this->GetPixelValueByIndex(itkIndex, timestep, component); return value; } vtkImageData *mitk::Image::GetVtkImageData(int t, int n) { if (m_Initialized == false) { if (GetSource().IsNull()) return nullptr; if (GetSource()->Updating() == false) GetSource()->UpdateOutputInformation(); } ImageDataItemPointer volume = GetVolumeData(t, n); return volume.GetPointer() == nullptr ? nullptr : volume->GetVtkImageAccessor(this)->GetVtkImageData(); } const vtkImageData *mitk::Image::GetVtkImageData(int t, int n) const { if (m_Initialized == false) { if (GetSource().IsNull()) return nullptr; if (GetSource()->Updating() == false) GetSource()->UpdateOutputInformation(); } ImageDataItemPointer volume = GetVolumeData(t, n); return volume.GetPointer() == nullptr ? nullptr : volume->GetVtkImageAccessor(this)->GetVtkImageData(); } mitk::Image::ImageDataItemPointer mitk::Image::GetSliceData( int s, int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const { MutexHolder lock(m_ImageDataArraysLock); return GetSliceData_unlocked(s, t, n, data, importMemoryManagement); } mitk::Image::ImageDataItemPointer mitk::Image::GetSliceData_unlocked( int s, int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const { if (IsValidSlice(s, t, n) == false) return nullptr; const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); // slice directly available? int pos = GetSliceIndex(s, t, n); if (m_Slices[pos].GetPointer() != nullptr) { return m_Slices[pos]; } // is slice available as part of a volume that is available? ImageDataItemPointer sl, ch, vol; vol = m_Volumes[GetVolumeIndex(t, n)]; if ((vol.GetPointer() != nullptr) && (vol->IsComplete())) { sl = new ImageDataItem(*vol, m_ImageDescriptor, t, 2, data, importMemoryManagement == ManageMemory, ((size_t)s) * m_OffsetTable[2] * (ptypeSize)); sl->SetComplete(true); return m_Slices[pos] = sl; } // is slice available as part of a channel that is available? ch = m_Channels[n]; if ((ch.GetPointer() != nullptr) && (ch->IsComplete())) { sl = new ImageDataItem(*ch, m_ImageDescriptor, t, 2, data, importMemoryManagement == ManageMemory, (((size_t)s) * m_OffsetTable[2] + ((size_t)t) * m_OffsetTable[3]) * (ptypeSize)); sl->SetComplete(true); return m_Slices[pos] = sl; } // slice is unavailable. Can we calculate it? if ((GetSource().IsNotNull()) && (GetSource()->Updating() == false)) { // ... wir mussen rechnen!!! .... m_RequestedRegion.SetIndex(0, 0); m_RequestedRegion.SetIndex(1, 0); m_RequestedRegion.SetIndex(2, s); m_RequestedRegion.SetIndex(3, t); m_RequestedRegion.SetIndex(4, n); m_RequestedRegion.SetSize(0, m_Dimensions[0]); m_RequestedRegion.SetSize(1, m_Dimensions[1]); m_RequestedRegion.SetSize(2, 1); m_RequestedRegion.SetSize(3, 1); m_RequestedRegion.SetSize(4, 1); m_RequestedRegionInitialized = true; GetSource()->Update(); if (IsSliceSet_unlocked(s, t, n)) // yes: now we can call ourselves without the risk of a endless loop (see "if" above) return GetSliceData_unlocked(s, t, n, data, importMemoryManagement); else return nullptr; } else { ImageDataItemPointer item = AllocateSliceData_unlocked(s, t, n, data, importMemoryManagement); item->SetComplete(true); return item; } } mitk::Image::ImageDataItemPointer mitk::Image::GetVolumeData(int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const { MutexHolder lock(m_ImageDataArraysLock); return GetVolumeData_unlocked(t, n, data, importMemoryManagement); } mitk::Image::ImageDataItemPointer mitk::Image::GetVolumeData_unlocked( int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const { if (IsValidVolume(t, n) == false) return nullptr; ImageDataItemPointer ch, vol; // volume directly available? int pos = GetVolumeIndex(t, n); vol = m_Volumes[pos]; if ((vol.GetPointer() != nullptr) && (vol->IsComplete())) return vol; const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); // is volume available as part of a channel that is available? ch = m_Channels[n]; if ((ch.GetPointer() != nullptr) && (ch->IsComplete())) { vol = new ImageDataItem(*ch, m_ImageDescriptor, t, 3, data, importMemoryManagement == ManageMemory, (((size_t)t) * m_OffsetTable[3]) * (ptypeSize)); vol->SetComplete(true); return m_Volumes[pos] = vol; } // let's see if all slices of the volume are set, so that we can (could) combine them to a volume bool complete = true; unsigned int s; for (s = 0; s < m_Dimensions[2]; ++s) { if (m_Slices[GetSliceIndex(s, t, n)].GetPointer() == nullptr) { complete = false; break; } } if (complete) { // if there is only single slice we do not need to combine anything if (m_Dimensions[2] <= 1) { ImageDataItemPointer sl; sl = GetSliceData_unlocked(0, t, n, data, importMemoryManagement); vol = new ImageDataItem(*sl, m_ImageDescriptor, t, 3, data, importMemoryManagement == ManageMemory); vol->SetComplete(true); } else { mitk::PixelType chPixelType = this->m_ImageDescriptor->GetChannelTypeById(n); vol = m_Volumes[pos]; // ok, let's combine the slices! if (vol.GetPointer() == nullptr) { vol = new ImageDataItem(chPixelType, t, 3, m_Dimensions, nullptr, true); } vol->SetComplete(true); size_t size = m_OffsetTable[2] * (ptypeSize); for (s = 0; s < m_Dimensions[2]; ++s) { int posSl; ImageDataItemPointer sl; posSl = GetSliceIndex(s, t, n); sl = m_Slices[posSl]; if (sl->GetParent() != vol) { // copy data of slices in volume size_t offset = ((size_t)s) * size; std::memcpy(static_cast(vol->GetData()) + offset, sl->GetData(), size); // FIXME mitkIpPicDescriptor * pic = sl->GetPicDescriptor(); // replace old slice with reference to volume sl = new ImageDataItem( *vol, m_ImageDescriptor, t, 2, data, importMemoryManagement == ManageMemory, ((size_t)s) * size); sl->SetComplete(true); // mitkIpFuncCopyTags(sl->GetPicDescriptor(), pic); m_Slices[posSl] = sl; } } // if(vol->GetPicDescriptor()->info->tags_head==NULL) // mitkIpFuncCopyTags(vol->GetPicDescriptor(), m_Slices[GetSliceIndex(0,t,n)]->GetPicDescriptor()); } return m_Volumes[pos] = vol; } // volume is unavailable. Can we calculate it? if ((GetSource().IsNotNull()) && (GetSource()->Updating() == false)) { // ... wir muessen rechnen!!! .... m_RequestedRegion.SetIndex(0, 0); m_RequestedRegion.SetIndex(1, 0); m_RequestedRegion.SetIndex(2, 0); m_RequestedRegion.SetIndex(3, t); m_RequestedRegion.SetIndex(4, n); m_RequestedRegion.SetSize(0, m_Dimensions[0]); m_RequestedRegion.SetSize(1, m_Dimensions[1]); m_RequestedRegion.SetSize(2, m_Dimensions[2]); m_RequestedRegion.SetSize(3, 1); m_RequestedRegion.SetSize(4, 1); m_RequestedRegionInitialized = true; GetSource()->Update(); if (IsVolumeSet_unlocked(t, n)) // yes: now we can call ourselves without the risk of a endless loop (see "if" above) return GetVolumeData_unlocked(t, n, data, importMemoryManagement); else return nullptr; } else { ImageDataItemPointer item = AllocateVolumeData_unlocked(t, n, data, importMemoryManagement); item->SetComplete(true); return item; } } mitk::Image::ImageDataItemPointer mitk::Image::GetChannelData(int n, void *data, ImportMemoryManagementType importMemoryManagement) const { MutexHolder lock(m_ImageDataArraysLock); return GetChannelData_unlocked(n, data, importMemoryManagement); } mitk::Image::ImageDataItemPointer mitk::Image::GetChannelData_unlocked( int n, void *data, ImportMemoryManagementType importMemoryManagement) const { if (IsValidChannel(n) == false) return nullptr; ImageDataItemPointer ch, vol; ch = m_Channels[n]; if ((ch.GetPointer() != nullptr) && (ch->IsComplete())) return ch; // let's see if all volumes are set, so that we can (could) combine them to a channel if (IsChannelSet_unlocked(n)) { // if there is only one time frame we do not need to combine anything if (m_Dimensions[3] <= 1) { vol = GetVolumeData_unlocked(0, n, data, importMemoryManagement); ch = new ImageDataItem(*vol, m_ImageDescriptor, 0, m_ImageDescriptor->GetNumberOfDimensions(), data, importMemoryManagement == ManageMemory); ch->SetComplete(true); } else { const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); ch = m_Channels[n]; // ok, let's combine the volumes! if (ch.GetPointer() == nullptr) ch = new ImageDataItem(this->m_ImageDescriptor, -1, nullptr, true); ch->SetComplete(true); size_t size = m_OffsetTable[m_Dimension - 1] * (ptypeSize); unsigned int t; auto slicesIt = m_Slices.begin() + n * m_Dimensions[2] * m_Dimensions[3]; for (t = 0; t < m_Dimensions[3]; ++t) { int posVol; ImageDataItemPointer vol; posVol = GetVolumeIndex(t, n); vol = GetVolumeData_unlocked(t, n, data, importMemoryManagement); if (vol->GetParent() != ch) { // copy data of volume in channel size_t offset = ((size_t)t) * m_OffsetTable[3] * (ptypeSize); std::memcpy(static_cast(ch->GetData()) + offset, vol->GetData(), size); // REVEIW FIX mitkIpPicDescriptor * pic = vol->GetPicDescriptor(); // replace old volume with reference to channel vol = new ImageDataItem(*ch, m_ImageDescriptor, t, 3, data, importMemoryManagement == ManageMemory, offset); vol->SetComplete(true); // mitkIpFuncCopyTags(vol->GetPicDescriptor(), pic); m_Volumes[posVol] = vol; // get rid of slices - they may point to old volume ImageDataItemPointer dnull = nullptr; for (unsigned int i = 0; i < m_Dimensions[2]; ++i, ++slicesIt) { assert(slicesIt != m_Slices.end()); *slicesIt = dnull; } } } // REVIEW FIX // if(ch->GetPicDescriptor()->info->tags_head==NULL) // mitkIpFuncCopyTags(ch->GetPicDescriptor(), m_Volumes[GetVolumeIndex(0,n)]->GetPicDescriptor()); } return m_Channels[n] = ch; } // channel is unavailable. Can we calculate it? if ((GetSource().IsNotNull()) && (GetSource()->Updating() == false)) { // ... wir muessen rechnen!!! .... m_RequestedRegion.SetIndex(0, 0); m_RequestedRegion.SetIndex(1, 0); m_RequestedRegion.SetIndex(2, 0); m_RequestedRegion.SetIndex(3, 0); m_RequestedRegion.SetIndex(4, n); m_RequestedRegion.SetSize(0, m_Dimensions[0]); m_RequestedRegion.SetSize(1, m_Dimensions[1]); m_RequestedRegion.SetSize(2, m_Dimensions[2]); m_RequestedRegion.SetSize(3, m_Dimensions[3]); m_RequestedRegion.SetSize(4, 1); m_RequestedRegionInitialized = true; GetSource()->Update(); // did it work? if (IsChannelSet_unlocked(n)) // yes: now we can call ourselves without the risk of a endless loop (see "if" above) return GetChannelData_unlocked(n, data, importMemoryManagement); else return nullptr; } else { ImageDataItemPointer item = AllocateChannelData_unlocked(n, data, importMemoryManagement); item->SetComplete(true); return item; } } bool mitk::Image::IsSliceSet(int s, int t, int n) const { MutexHolder lock(m_ImageDataArraysLock); return IsSliceSet_unlocked(s, t, n); } bool mitk::Image::IsSliceSet_unlocked(int s, int t, int n) const { if (IsValidSlice(s, t, n) == false) return false; if (m_Slices[GetSliceIndex(s, t, n)].GetPointer() != nullptr) { return true; } ImageDataItemPointer ch, vol; vol = m_Volumes[GetVolumeIndex(t, n)]; if ((vol.GetPointer() != nullptr) && (vol->IsComplete())) { return true; } ch = m_Channels[n]; if ((ch.GetPointer() != nullptr) && (ch->IsComplete())) { return true; } return false; } bool mitk::Image::IsVolumeSet(int t, int n) const { MutexHolder lock(m_ImageDataArraysLock); return IsVolumeSet_unlocked(t, n); } bool mitk::Image::IsVolumeSet_unlocked(int t, int n) const { if (IsValidVolume(t, n) == false) return false; ImageDataItemPointer ch, vol; // volume directly available? vol = m_Volumes[GetVolumeIndex(t, n)]; if ((vol.GetPointer() != nullptr) && (vol->IsComplete())) return true; // is volume available as part of a channel that is available? ch = m_Channels[n]; if ((ch.GetPointer() != nullptr) && (ch->IsComplete())) return true; // let's see if all slices of the volume are set, so that we can (could) combine them to a volume unsigned int s; for (s = 0; s < m_Dimensions[2]; ++s) { if (m_Slices[GetSliceIndex(s, t, n)].GetPointer() == nullptr) { return false; } } return true; } bool mitk::Image::IsChannelSet(int n) const { MutexHolder lock(m_ImageDataArraysLock); return IsChannelSet_unlocked(n); } bool mitk::Image::IsChannelSet_unlocked(int n) const { if (IsValidChannel(n) == false) return false; ImageDataItemPointer ch, vol; ch = m_Channels[n]; if ((ch.GetPointer() != nullptr) && (ch->IsComplete())) return true; // let's see if all volumes are set, so that we can (could) combine them to a channel unsigned int t; for (t = 0; t < m_Dimensions[3]; ++t) { if (IsVolumeSet_unlocked(t, n) == false) { return false; } } return true; } bool mitk::Image::SetSlice(const void *data, int s, int t, int n) { // const_cast is no risk for ImportMemoryManagementType == CopyMemory return SetImportSlice(const_cast(data), s, t, n, CopyMemory); } bool mitk::Image::SetVolume(const void *data, int t, int n) { // const_cast is no risk for ImportMemoryManagementType == CopyMemory return SetImportVolume(const_cast(data), t, n, CopyMemory); } bool mitk::Image::SetChannel(const void *data, int n) { // const_cast is no risk for ImportMemoryManagementType == CopyMemory return SetImportChannel(const_cast(data), n, CopyMemory); } bool mitk::Image::SetImportSlice(void *data, int s, int t, int n, ImportMemoryManagementType importMemoryManagement) { if (IsValidSlice(s, t, n) == false) return false; ImageDataItemPointer sl; const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); if (IsSliceSet(s, t, n)) { sl = GetSliceData(s, t, n, data, importMemoryManagement); if (sl->GetManageMemory() == false) { sl = AllocateSliceData(s, t, n, data, importMemoryManagement); if (sl.GetPointer() == nullptr) return false; } if (sl->GetData() != data) std::memcpy(sl->GetData(), data, m_OffsetTable[2] * (ptypeSize)); sl->Modified(); // we have changed the data: call Modified()! Modified(); } else { sl = AllocateSliceData(s, t, n, data, importMemoryManagement); if (sl.GetPointer() == nullptr) return false; if (sl->GetData() != data) std::memcpy(sl->GetData(), data, m_OffsetTable[2] * (ptypeSize)); // we just added a missing slice, which is not regarded as modification. // Therefore, we do not call Modified()! } return true; } bool mitk::Image::SetImportVolume(void *data, int t, int n, ImportMemoryManagementType importMemoryManagement) { if (IsValidVolume(t, n) == false) return false; const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); ImageDataItemPointer vol; if (IsVolumeSet(t, n)) { vol = GetVolumeData(t, n, data, importMemoryManagement); if (vol->GetManageMemory() == false) { vol = AllocateVolumeData(t, n, data, importMemoryManagement); if (vol.GetPointer() == nullptr) return false; } if (vol->GetData() != data) std::memcpy(vol->GetData(), data, m_OffsetTable[3] * (ptypeSize)); vol->Modified(); vol->SetComplete(true); // we have changed the data: call Modified()! Modified(); } else { vol = AllocateVolumeData(t, n, data, importMemoryManagement); if (vol.GetPointer() == nullptr) return false; if (vol->GetData() != data) { std::memcpy(vol->GetData(), data, m_OffsetTable[3] * (ptypeSize)); } vol->SetComplete(true); this->m_ImageDescriptor->GetChannelDescriptor(n).SetData(vol->GetData()); // we just added a missing Volume, which is not regarded as modification. // Therefore, we do not call Modified()! } return true; } bool mitk::Image::SetImportChannel(void *data, int n, ImportMemoryManagementType importMemoryManagement) { if (IsValidChannel(n) == false) return false; // channel descriptor const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); ImageDataItemPointer ch; if (IsChannelSet(n)) { ch = GetChannelData(n, data, importMemoryManagement); if (ch->GetManageMemory() == false) { ch = AllocateChannelData(n, data, importMemoryManagement); if (ch.GetPointer() == nullptr) return false; } if (ch->GetData() != data) std::memcpy(ch->GetData(), data, m_OffsetTable[4] * (ptypeSize)); ch->Modified(); ch->SetComplete(true); // we have changed the data: call Modified()! Modified(); } else { ch = AllocateChannelData(n, data, importMemoryManagement); if (ch.GetPointer() == nullptr) return false; if (ch->GetData() != data) std::memcpy(ch->GetData(), data, m_OffsetTable[4] * (ptypeSize)); ch->SetComplete(true); this->m_ImageDescriptor->GetChannelDescriptor(n).SetData(ch->GetData()); // we just added a missing Channel, which is not regarded as modification. // Therefore, we do not call Modified()! } return true; } void mitk::Image::Initialize() { ImageDataItemPointerArray::iterator it, end; for (it = m_Slices.begin(), end = m_Slices.end(); it != end; ++it) { (*it) = nullptr; } for (it = m_Volumes.begin(), end = m_Volumes.end(); it != end; ++it) { (*it) = nullptr; } for (it = m_Channels.begin(), end = m_Channels.end(); it != end; ++it) { (*it) = nullptr; } m_CompleteData = nullptr; if (m_ImageStatistics == nullptr) { m_ImageStatistics = new mitk::ImageStatisticsHolder(this); } SetRequestedRegionToLargestPossibleRegion(); } void mitk::Image::Initialize(const mitk::ImageDescriptor::Pointer inDesc) { // store the descriptor this->m_ImageDescriptor = inDesc; // initialize image this->Initialize( inDesc->GetChannelDescriptor(0).GetPixelType(), inDesc->GetNumberOfDimensions(), inDesc->GetDimensions(), 1); } void mitk::Image::Initialize(const mitk::PixelType &type, unsigned int dimension, const unsigned int *dimensions, unsigned int channels) { Clear(); m_Dimension = dimension; if (!dimensions) itkExceptionMacro(<< "invalid zero dimension image"); unsigned int i; for (i = 0; i < dimension; ++i) { if (dimensions[i] < 1) itkExceptionMacro(<< "invalid dimension[" << i << "]: " << dimensions[i]); } // create new array since the old was deleted m_Dimensions = new unsigned int[MAX_IMAGE_DIMENSIONS]; // initialize the first four dimensions to 1, the remaining 4 to 0 FILL_C_ARRAY(m_Dimensions, 4, 1u); FILL_C_ARRAY((m_Dimensions + 4), 4, 0u); // copy in the passed dimension information std::memcpy(m_Dimensions, dimensions, sizeof(unsigned int) * m_Dimension); this->m_ImageDescriptor = mitk::ImageDescriptor::New(); this->m_ImageDescriptor->Initialize(this->m_Dimensions, this->m_Dimension); for (i = 0; i < 4; ++i) { m_LargestPossibleRegion.SetIndex(i, 0); m_LargestPossibleRegion.SetSize(i, m_Dimensions[i]); } m_LargestPossibleRegion.SetIndex(i, 0); m_LargestPossibleRegion.SetSize(i, channels); if (m_LargestPossibleRegion.GetNumberOfPixels() == 0) { delete[] m_Dimensions; m_Dimensions = nullptr; return; } for (unsigned int i = 0u; i < channels; i++) { this->m_ImageDescriptor->AddNewChannel(type); } PlaneGeometry::Pointer planegeometry = PlaneGeometry::New(); planegeometry->InitializeStandardPlane(m_Dimensions[0], m_Dimensions[1]); SlicedGeometry3D::Pointer slicedGeometry = SlicedGeometry3D::New(); slicedGeometry->InitializeEvenlySpaced(planegeometry, m_Dimensions[2]); ProportionalTimeGeometry::Pointer timeGeometry = ProportionalTimeGeometry::New(); timeGeometry->Initialize(slicedGeometry, m_Dimensions[3]); for (TimeStepType step = 0; step < timeGeometry->CountTimeSteps(); ++step) { timeGeometry->GetGeometryForTimeStep(step)->ImageGeometryOn(); } SetTimeGeometry(timeGeometry); ImageDataItemPointer dnull = nullptr; m_Channels.assign(GetNumberOfChannels(), dnull); m_Volumes.assign(GetNumberOfChannels() * m_Dimensions[3], dnull); m_Slices.assign(GetNumberOfChannels() * m_Dimensions[3] * m_Dimensions[2], dnull); ComputeOffsetTable(); Initialize(); m_Initialized = true; } void mitk::Image::Initialize(const mitk::PixelType &type, const mitk::BaseGeometry &geometry, unsigned int channels, int tDim) { mitk::ProportionalTimeGeometry::Pointer timeGeometry = ProportionalTimeGeometry::New(); timeGeometry->Initialize(geometry.Clone(), tDim); this->Initialize(type, *timeGeometry, channels, tDim); } void mitk::Image::Initialize(const mitk::PixelType &type, const mitk::TimeGeometry &geometry, unsigned int channels, int tDim) { unsigned int dimensions[5]; dimensions[0] = (unsigned int)(geometry.GetGeometryForTimeStep(0)->GetExtent(0) + 0.5); dimensions[1] = (unsigned int)(geometry.GetGeometryForTimeStep(0)->GetExtent(1) + 0.5); dimensions[2] = (unsigned int)(geometry.GetGeometryForTimeStep(0)->GetExtent(2) + 0.5); dimensions[3] = (tDim > 0) ? tDim : geometry.CountTimeSteps(); dimensions[4] = 0; unsigned int dimension = 2; if (dimensions[2] > 1) dimension = 3; if (dimensions[3] > 1) dimension = 4; Initialize(type, dimension, dimensions, channels); if (geometry.CountTimeSteps() > 1) { TimeGeometry::Pointer cloned = geometry.Clone(); SetTimeGeometry(cloned.GetPointer()); // make sure the image geometry flag is properly set for all time steps for (TimeStepType step = 0; step < cloned->CountTimeSteps(); ++step) { if (!cloned->GetGeometryCloneForTimeStep(step)->GetImageGeometry()) { MITK_WARN("Image.3DnT.Initialize") << " Attempt to initialize an image with a non-image geometry. " "Re-interpretting the initialization geometry for timestep " << step << " as image geometry, the original geometry remains unchanged."; cloned->GetGeometryForTimeStep(step)->ImageGeometryOn(); } } } else { // make sure the image geometry coming from outside has proper value of the image geometry flag BaseGeometry::Pointer cloned = geometry.GetGeometryCloneForTimeStep(0)->Clone(); if (!cloned->GetImageGeometry()) { MITK_WARN("Image.Initialize") << " Attempt to initialize an image with a non-image geometry. Re-interpretting " "the initialization geometry as image geometry, the original geometry remains " "unchanged."; cloned->ImageGeometryOn(); } Superclass::SetGeometry(cloned); } /* //Old //TODO_GOETZ Really necessary? mitk::BoundingBox::BoundsArrayType bounds = geometry.GetBoundingBoxInWorld()->GetBounds(); if( (bounds[0] != 0.0) || (bounds[2] != 0.0) || (bounds[4] != 0.0) ) { SlicedGeometry3D* slicedGeometry = GetSlicedGeometry(0); mitk::Point3D origin; origin.Fill(0.0); slicedGeometry->IndexToWorld(origin, origin); bounds[1]-=bounds[0]; bounds[3]-=bounds[2]; bounds[5]-=bounds[4]; bounds[0] = 0.0; bounds[2] = 0.0; bounds[4] = 0.0; this->m_ImageDescriptor->Initialize( this->m_Dimensions, this->m_Dimension ); slicedGeometry->SetBounds(bounds); slicedGeometry->GetIndexToWorldTransform()->SetOffset(origin.GetVnlVector().data_block()); ProportionalTimeGeometry::Pointer timeGeometry = ProportionalTimeGeometry::New(); timeGeometry->Initialize(slicedGeometry, m_Dimensions[3]); SetTimeGeometry(timeGeometry); }*/ } +void mitk::Image::Initialize(const mitk::PixelType &, + int, + const mitk::PlaneGeometry &, + bool, + unsigned int, + int) +{ + mitkThrow() << "Use this method without the flipped parameter (direction is specified by the handedness of the PlaneGeometry instead)."; +} + void mitk::Image::Initialize(const mitk::PixelType &type, int sDim, const mitk::PlaneGeometry &geometry2d, - bool flipped, unsigned int channels, int tDim) { SlicedGeometry3D::Pointer slicedGeometry = SlicedGeometry3D::New(); - slicedGeometry->InitializeEvenlySpaced(static_cast(geometry2d.Clone().GetPointer()), sDim, flipped); + slicedGeometry->InitializeEvenlySpaced(geometry2d.Clone(), sDim); Initialize(type, *slicedGeometry, channels, tDim); } void mitk::Image::Initialize(const mitk::Image *image) { Initialize(image->GetPixelType(), *image->GetTimeGeometry()); } void mitk::Image::Initialize(vtkImageData *vtkimagedata, int channels, int tDim, int sDim, int pDim) { if (vtkimagedata == nullptr) return; m_Dimension = vtkimagedata->GetDataDimension(); unsigned int i, *tmpDimensions = new unsigned int[m_Dimension > 4 ? m_Dimension : 4]; for (i = 0; i < m_Dimension; ++i) tmpDimensions[i] = vtkimagedata->GetDimensions()[i]; if (m_Dimension < 4) { unsigned int *p; for (i = 0, p = tmpDimensions + m_Dimension; i < 4 - m_Dimension; ++i, ++p) *p = 1; } if (pDim >= 0) { tmpDimensions[1] = pDim; if (m_Dimension < 2) m_Dimension = 2; } if (sDim >= 0) { tmpDimensions[2] = sDim; if (m_Dimension < 3) m_Dimension = 3; } if (tDim >= 0) { tmpDimensions[3] = tDim; if (m_Dimension < 4) m_Dimension = 4; } mitk::PixelType pixelType(MakePixelType(vtkimagedata)); Initialize(pixelType, m_Dimension, tmpDimensions, channels); const double *spacinglist = vtkimagedata->GetSpacing(); Vector3D spacing; FillVector3D(spacing, spacinglist[0], 1.0, 1.0); if (m_Dimension >= 2) spacing[1] = spacinglist[1]; if (m_Dimension >= 3) spacing[2] = spacinglist[2]; // access origin of vtkImage Point3D origin; double vtkorigin[3]; vtkimagedata->GetOrigin(vtkorigin); FillVector3D(origin, vtkorigin[0], 0.0, 0.0); if (m_Dimension >= 2) origin[1] = vtkorigin[1]; if (m_Dimension >= 3) origin[2] = vtkorigin[2]; SlicedGeometry3D *slicedGeometry = GetSlicedGeometry(0); // re-initialize PlaneGeometry with origin and direction PlaneGeometry *planeGeometry = static_cast(slicedGeometry->GetPlaneGeometry(0)); planeGeometry->SetOrigin(origin); // re-initialize SlicedGeometry3D slicedGeometry->SetOrigin(origin); slicedGeometry->SetSpacing(spacing); ProportionalTimeGeometry::Pointer timeGeometry = ProportionalTimeGeometry::New(); timeGeometry->Initialize(slicedGeometry, m_Dimensions[3]); SetTimeGeometry(timeGeometry); delete[] tmpDimensions; } bool mitk::Image::IsValidSlice(int s, int t, int n) const { if (m_Initialized) return ((s >= 0) && (s < (int)m_Dimensions[2]) && (t >= 0) && (t < (int)m_Dimensions[3]) && (n >= 0) && (n < (int)GetNumberOfChannels())); else return false; } bool mitk::Image::IsValidVolume(int t, int n) const { if (m_Initialized) return IsValidSlice(0, t, n); else return false; } bool mitk::Image::IsValidChannel(int n) const { if (m_Initialized) return IsValidSlice(0, 0, n); else return false; } void mitk::Image::ComputeOffsetTable() { if (m_OffsetTable != nullptr) delete[] m_OffsetTable; m_OffsetTable = new size_t[m_Dimension > 4 ? m_Dimension + 1 : 4 + 1]; unsigned int i; size_t num = 1; m_OffsetTable[0] = 1; for (i = 0; i < m_Dimension; ++i) { num *= m_Dimensions[i]; m_OffsetTable[i + 1] = num; } for (; i < 4; ++i) m_OffsetTable[i + 1] = num; } bool mitk::Image::IsValidTimeStep(int t) const { return ((m_Dimension >= 4 && t <= (int)m_Dimensions[3] && t > 0) || (t == 0)); } void mitk::Image::Expand(unsigned int timeSteps) { if (timeSteps < 1) itkExceptionMacro(<< "Invalid timestep in Image!"); Superclass::Expand(timeSteps); } int mitk::Image::GetSliceIndex(int s, int t, int n) const { if (IsValidSlice(s, t, n) == false) return false; return ((size_t)s) + ((size_t)t) * m_Dimensions[2] + ((size_t)n) * m_Dimensions[3] * m_Dimensions[2]; //?? } int mitk::Image::GetVolumeIndex(int t, int n) const { if (IsValidVolume(t, n) == false) return false; return ((size_t)t) + ((size_t)n) * m_Dimensions[3]; //?? } mitk::Image::ImageDataItemPointer mitk::Image::AllocateSliceData( int s, int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const { MutexHolder lock(m_ImageDataArraysLock); return AllocateSliceData_unlocked(s, t, n, data, importMemoryManagement); } mitk::Image::ImageDataItemPointer mitk::Image::AllocateSliceData_unlocked( int s, int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const { int pos; pos = GetSliceIndex(s, t, n); const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); // is slice available as part of a volume that is available? ImageDataItemPointer sl, ch, vol; vol = m_Volumes[GetVolumeIndex(t, n)]; if (vol.GetPointer() != nullptr) { sl = new ImageDataItem(*vol, m_ImageDescriptor, t, 2, data, importMemoryManagement == ManageMemory, ((size_t)s) * m_OffsetTable[2] * (ptypeSize)); sl->SetComplete(true); return m_Slices[pos] = sl; } // is slice available as part of a channel that is available? ch = m_Channels[n]; if (ch.GetPointer() != nullptr) { sl = new ImageDataItem(*ch, m_ImageDescriptor, t, 2, data, importMemoryManagement == ManageMemory, (((size_t)s) * m_OffsetTable[2] + ((size_t)t) * m_OffsetTable[3]) * (ptypeSize)); sl->SetComplete(true); return m_Slices[pos] = sl; } // allocate new volume (instead of a single slice to keep data together!) m_Volumes[GetVolumeIndex(t, n)] = vol = AllocateVolumeData_unlocked(t, n, nullptr, importMemoryManagement); sl = new ImageDataItem(*vol, m_ImageDescriptor, t, 2, data, importMemoryManagement == ManageMemory, ((size_t)s) * m_OffsetTable[2] * (ptypeSize)); sl->SetComplete(true); return m_Slices[pos] = sl; ////ALTERNATIVE: //// allocate new slice // sl=new ImageDataItem(*m_PixelType, 2, m_Dimensions); // m_Slices[pos]=sl; // return vol; } mitk::Image::ImageDataItemPointer mitk::Image::AllocateVolumeData( int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const { MutexHolder lock(m_ImageDataArraysLock); return AllocateVolumeData_unlocked(t, n, data, importMemoryManagement); } mitk::Image::ImageDataItemPointer mitk::Image::AllocateVolumeData_unlocked( int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) const { int pos; pos = GetVolumeIndex(t, n); const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); // is volume available as part of a channel that is available? ImageDataItemPointer ch, vol; ch = m_Channels[n]; if (ch.GetPointer() != nullptr) { vol = new ImageDataItem(*ch, m_ImageDescriptor, t, 3, data, importMemoryManagement == ManageMemory, (((size_t)t) * m_OffsetTable[3]) * (ptypeSize)); return m_Volumes[pos] = vol; } mitk::PixelType chPixelType = this->m_ImageDescriptor->GetChannelTypeById(n); // allocate new volume if (importMemoryManagement == CopyMemory) { vol = new ImageDataItem(chPixelType, t, 3, m_Dimensions, nullptr, true); if (data != nullptr) std::memcpy(vol->GetData(), data, m_OffsetTable[3] * (ptypeSize)); } else { vol = new ImageDataItem(chPixelType, t, 3, m_Dimensions, data, importMemoryManagement == ManageMemory); } m_Volumes[pos] = vol; return vol; } mitk::Image::ImageDataItemPointer mitk::Image::AllocateChannelData( int n, void *data, ImportMemoryManagementType importMemoryManagement) const { MutexHolder lock(m_ImageDataArraysLock); return AllocateChannelData_unlocked(n, data, importMemoryManagement); } mitk::Image::ImageDataItemPointer mitk::Image::AllocateChannelData_unlocked( int n, void *data, ImportMemoryManagementType importMemoryManagement) const { ImageDataItemPointer ch; // allocate new channel if (importMemoryManagement == CopyMemory) { const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); ch = new ImageDataItem(this->m_ImageDescriptor, -1, nullptr, true); if (data != nullptr) std::memcpy(ch->GetData(), data, m_OffsetTable[4] * (ptypeSize)); } else { ch = new ImageDataItem(this->m_ImageDescriptor, -1, data, importMemoryManagement == ManageMemory); } m_Channels[n] = ch; return ch; } unsigned int *mitk::Image::GetDimensions() const { return m_Dimensions; } void mitk::Image::Clear() { Superclass::Clear(); delete[] m_Dimensions; m_Dimensions = nullptr; } void mitk::Image::SetGeometry(BaseGeometry *aGeometry3D) { // Please be aware of the 0.5 offset/pixel-center issue! See Geometry documentation for further information if (aGeometry3D->GetImageGeometry() == false) { MITK_INFO << "WARNING: Applied a non-image geometry onto an image. Please be SURE that this geometry is " "pixel-center-based! If it is not, you need to call " "Geometry3D->ChangeImageGeometryConsideringOriginOffset(true) before calling image->setGeometry(..)\n"; } Superclass::SetGeometry(aGeometry3D); for (TimeStepType step = 0; step < GetTimeGeometry()->CountTimeSteps(); ++step) GetTimeGeometry()->GetGeometryForTimeStep(step)->ImageGeometryOn(); } void mitk::Image::PrintSelf(std::ostream &os, itk::Indent indent) const { if (m_Initialized) { unsigned char i; os << indent << " Dimension: " << m_Dimension << std::endl; os << indent << " Dimensions: "; for (i = 0; i < m_Dimension; ++i) os << GetDimension(i) << " "; os << std::endl; for (unsigned int ch = 0; ch < this->m_ImageDescriptor->GetNumberOfChannels(); ch++) { mitk::PixelType chPixelType = this->m_ImageDescriptor->GetChannelTypeById(ch); os << indent << " Channel: " << this->m_ImageDescriptor->GetChannelName(ch) << std::endl; os << indent << " PixelType: " << chPixelType.GetPixelTypeAsString() << std::endl; os << indent << " BytesPerElement: " << chPixelType.GetSize() << std::endl; os << indent << " ComponentType: " << chPixelType.GetComponentTypeAsString() << std::endl; os << indent << " NumberOfComponents: " << chPixelType.GetNumberOfComponents() << std::endl; os << indent << " BitsPerComponent: " << chPixelType.GetBitsPerComponent() << std::endl; } } else { os << indent << " Image not initialized: m_Initialized: false" << std::endl; } Superclass::PrintSelf(os, indent); } bool mitk::Image::IsRotated() const { const mitk::BaseGeometry *geo = this->GetGeometry(); bool ret = false; if (geo) { const vnl_matrix_fixed &mx = geo->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix(); mitk::ScalarType ref = 0; for (short k = 0; k < 3; ++k) ref += mx[k][k]; ref /= 1000; // Arbitrary value; if a non-diagonal (nd) element is bigger then this, matrix is considered nd. for (short i = 0; i < 3; ++i) { for (short j = 0; j < 3; ++j) { if (i != j) { if (std::abs(mx[i][j]) > ref) // matrix is nd ret = true; } } } } return ret; } mitk::ScalarType mitk::Image::GetScalarValueMin(int t) const { return m_ImageStatistics->GetScalarValueMin(t); } //## \brief Get the maximum for scalar images mitk::ScalarType mitk::Image::GetScalarValueMax(int t) const { return m_ImageStatistics->GetScalarValueMax(t); } //## \brief Get the second smallest value for scalar images mitk::ScalarType mitk::Image::GetScalarValue2ndMin(int t) const { return m_ImageStatistics->GetScalarValue2ndMin(t); } mitk::ScalarType mitk::Image::GetScalarValueMinNoRecompute(unsigned int t) const { return m_ImageStatistics->GetScalarValueMinNoRecompute(t); } mitk::ScalarType mitk::Image::GetScalarValue2ndMinNoRecompute(unsigned int t) const { return m_ImageStatistics->GetScalarValue2ndMinNoRecompute(t); } mitk::ScalarType mitk::Image::GetScalarValue2ndMax(int t) const { return m_ImageStatistics->GetScalarValue2ndMax(t); } mitk::ScalarType mitk::Image::GetScalarValueMaxNoRecompute(unsigned int t) const { return m_ImageStatistics->GetScalarValueMaxNoRecompute(t); } mitk::ScalarType mitk::Image::GetScalarValue2ndMaxNoRecompute(unsigned int t) const { return m_ImageStatistics->GetScalarValue2ndMaxNoRecompute(t); } mitk::ScalarType mitk::Image::GetCountOfMinValuedVoxels(int t) const { return m_ImageStatistics->GetCountOfMinValuedVoxels(t); } mitk::ScalarType mitk::Image::GetCountOfMaxValuedVoxels(int t) const { return m_ImageStatistics->GetCountOfMaxValuedVoxels(t); } unsigned int mitk::Image::GetCountOfMaxValuedVoxelsNoRecompute(unsigned int t) const { return m_ImageStatistics->GetCountOfMaxValuedVoxelsNoRecompute(t); } unsigned int mitk::Image::GetCountOfMinValuedVoxelsNoRecompute(unsigned int t) const { return m_ImageStatistics->GetCountOfMinValuedVoxelsNoRecompute(t); } bool mitk::Equal(const mitk::Image *leftHandSide, const mitk::Image *rightHandSide, ScalarType eps, bool verbose) { if ((leftHandSide == nullptr) || (rightHandSide == nullptr)) { MITK_ERROR << "mitk::Equal(const mitk::Image* leftHandSide, const mitk::Image* rightHandSide, ScalarType eps, bool " "verbose) does not work with NULL pointer input."; return false; } return mitk::Equal(*leftHandSide, *rightHandSide, eps, verbose); } bool mitk::Equal(const mitk::Image &leftHandSide, const mitk::Image &rightHandSide, ScalarType eps, bool verbose) { bool returnValue = true; // Dimensionality if (rightHandSide.GetDimension() != leftHandSide.GetDimension()) { if (verbose) { MITK_INFO << "[( Image )] Dimensionality differs."; MITK_INFO << "leftHandSide is " << leftHandSide.GetDimension() << "rightHandSide is " << rightHandSide.GetDimension(); } returnValue = false; } // Pair-wise dimension (size) comparison unsigned int minDimensionality = std::min(rightHandSide.GetDimension(), leftHandSide.GetDimension()); for (unsigned int i = 0; i < minDimensionality; ++i) { if (rightHandSide.GetDimension(i) != leftHandSide.GetDimension(i)) { returnValue = false; if (verbose) { MITK_INFO << "[( Image )] dimension differs."; MITK_INFO << "leftHandSide->GetDimension(" << i << ") is " << leftHandSide.GetDimension(i) << "rightHandSide->GetDimension(" << i << ") is " << rightHandSide.GetDimension(i); } } } // Pixeltype mitk::PixelType pixelTypeRightHandSide = rightHandSide.GetPixelType(); mitk::PixelType pixelTypeLeftHandSide = leftHandSide.GetPixelType(); if (!(pixelTypeRightHandSide == pixelTypeLeftHandSide)) { if (verbose) { MITK_INFO << "[( Image )] PixelType differs."; MITK_INFO << "leftHandSide is " << pixelTypeLeftHandSide.GetTypeAsString() << "rightHandSide is " << pixelTypeRightHandSide.GetTypeAsString(); } returnValue = false; } // Geometries if (!mitk::Equal(*leftHandSide.GetGeometry(), *rightHandSide.GetGeometry(), eps, verbose)) { if (verbose) { MITK_INFO << "[( Image )] Geometries differ."; } returnValue = false; } // Pixel values - default mode [ 0 threshold in difference ] // compare only if all previous checks were successfull, otherwise the ITK filter will throw an exception if (returnValue) { mitk::CompareImageDataFilter::Pointer compareFilter = mitk::CompareImageDataFilter::New(); compareFilter->SetInput(0, &rightHandSide); compareFilter->SetInput(1, &leftHandSide); compareFilter->SetTolerance(eps); compareFilter->Update(); if ((!compareFilter->GetResult())) { returnValue = false; if (verbose) { MITK_INFO << "[(Image)] Pixel values differ: "; compareFilter->GetCompareResults().PrintSelf(); } } } return returnValue; } diff --git a/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp b/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp index 60a3be37dc..a5ec10e50f 100644 --- a/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp +++ b/Modules/Core/src/DataManagement/mitkSlicedGeometry3D.cpp @@ -1,961 +1,962 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include "mitkSlicedGeometry3D.h" #include "mitkAbstractTransformGeometry.h" #include "mitkApplyTransformMatrixOperation.h" #include "mitkInteractionConst.h" #include "mitkPlaneGeometry.h" #include "mitkPlaneOperation.h" #include "mitkRestorePlanePositionOperation.h" #include "mitkRotationOperation.h" #include "mitkSliceNavigationController.h" const mitk::ScalarType PI = 3.14159265359; mitk::SlicedGeometry3D::SlicedGeometry3D() : m_EvenlySpaced(true), m_Slices(0), m_ReferenceGeometry(nullptr), m_SliceNavigationController(nullptr) { m_DirectionVector.Fill(0); this->InitializeSlicedGeometry(m_Slices); } mitk::SlicedGeometry3D::SlicedGeometry3D(const SlicedGeometry3D &other) : Superclass(other), m_EvenlySpaced(other.m_EvenlySpaced), m_Slices(other.m_Slices), m_ReferenceGeometry(other.m_ReferenceGeometry), m_SliceNavigationController(other.m_SliceNavigationController) { m_DirectionVector.Fill(0); SetSpacing(other.GetSpacing()); SetDirectionVector(other.GetDirectionVector()); if (m_EvenlySpaced) { + assert(!other.m_PlaneGeometries.empty() && "This may happen when you use one of the old Initialize methods, which had a bool parameter that is implicitly casted to the number of slices now."); PlaneGeometry::Pointer geometry = other.m_PlaneGeometries[0]->Clone(); assert(geometry.IsNotNull()); SetPlaneGeometry(geometry, 0); } else { unsigned int s; for (s = 0; s < other.m_Slices; ++s) { if (other.m_PlaneGeometries[s].IsNull()) { assert(other.m_EvenlySpaced); m_PlaneGeometries[s] = nullptr; } else { PlaneGeometry *geometry2D = other.m_PlaneGeometries[s]->Clone(); assert(geometry2D != nullptr); SetPlaneGeometry(geometry2D, s); } } } } mitk::SlicedGeometry3D::~SlicedGeometry3D() { } mitk::PlaneGeometry *mitk::SlicedGeometry3D::GetPlaneGeometry(int s) const { mitk::PlaneGeometry::Pointer geometry2D = nullptr; if (this->IsValidSlice(s)) { geometry2D = m_PlaneGeometries[s]; // If (a) m_EvenlySpaced==true, (b) we don't have a PlaneGeometry stored // for the requested slice, and (c) the first slice (s=0) // is a PlaneGeometry instance, then we calculate the geometry of the // requested as the plane of the first slice shifted by m_Spacing[2]*s // in the direction of m_DirectionVector. if ((m_EvenlySpaced) && (geometry2D.IsNull())) { PlaneGeometry *firstSlice = m_PlaneGeometries[0]; if (firstSlice != nullptr && dynamic_cast(m_PlaneGeometries[0].GetPointer()) == nullptr) { if ((m_DirectionVector[0] == 0.0) && (m_DirectionVector[1] == 0.0) && (m_DirectionVector[2] == 0.0)) { m_DirectionVector = firstSlice->GetNormal(); m_DirectionVector.Normalize(); } Vector3D direction; direction = m_DirectionVector * this->GetSpacing()[2]; mitk::PlaneGeometry::Pointer requestedslice; requestedslice = static_cast(firstSlice->Clone().GetPointer()); requestedslice->SetOrigin(requestedslice->GetOrigin() + direction * s); geometry2D = requestedslice; m_PlaneGeometries[s] = geometry2D; } } return geometry2D; } else { return nullptr; } } const mitk::BoundingBox *mitk::SlicedGeometry3D::GetBoundingBox() const { assert(this->IsBoundingBoxNull() == false); return Superclass::GetBoundingBox(); } bool mitk::SlicedGeometry3D::SetPlaneGeometry(mitk::PlaneGeometry *geometry2D, int s) { if (this->IsValidSlice(s)) { m_PlaneGeometries[s] = geometry2D; m_PlaneGeometries[s]->SetReferenceGeometry(m_ReferenceGeometry); return true; } return false; } void mitk::SlicedGeometry3D::InitializeSlicedGeometry(unsigned int slices) { Superclass::Initialize(); m_Slices = slices; PlaneGeometry::Pointer gnull = nullptr; m_PlaneGeometries.assign(m_Slices, gnull); Vector3D spacing; spacing.Fill(1.0); this->SetSpacing(spacing); m_DirectionVector.Fill(0); } void mitk::SlicedGeometry3D::InitializeEvenlySpaced(mitk::PlaneGeometry *geometry2D, unsigned int slices) { assert(geometry2D != nullptr); this->InitializeEvenlySpaced(geometry2D, geometry2D->GetExtentInMM(2) / geometry2D->GetExtent(2), slices); } void mitk::SlicedGeometry3D::InitializeEvenlySpaced(mitk::PlaneGeometry *geometry2D, mitk::ScalarType zSpacing, unsigned int slices) { assert(geometry2D != nullptr); assert(geometry2D->GetExtent(0) > 0); assert(geometry2D->GetExtent(1) > 0); geometry2D->Register(); Superclass::Initialize(); m_Slices = slices; BoundingBox::BoundsArrayType bounds = geometry2D->GetBounds(); bounds[4] = 0; bounds[5] = slices; // clear and reserve PlaneGeometry::Pointer gnull = nullptr; m_PlaneGeometries.assign(m_Slices, gnull); Vector3D directionVector = geometry2D->GetAxisVector(2); directionVector.Normalize(); directionVector *= zSpacing; // Normally we should use the following four lines to create a copy of // the transform contrained in geometry2D, because it may not be changed // by us. But we know that SetSpacing creates a new transform without // changing the old (coming from geometry2D), so we can use the fifth // line instead. We check this at (**). // // AffineTransform3D::Pointer transform = AffineTransform3D::New(); // transform->SetMatrix(geometry2D->GetIndexToWorldTransform()->GetMatrix()); // transform->SetOffset(geometry2D->GetIndexToWorldTransform()->GetOffset()); // SetIndexToWorldTransform(transform); this->SetIndexToWorldTransform(const_cast(geometry2D->GetIndexToWorldTransform())); mitk::Vector3D spacing; FillVector3D(spacing, geometry2D->GetExtentInMM(0) / bounds[1], geometry2D->GetExtentInMM(1) / bounds[3], zSpacing); this->SetDirectionVector(directionVector); this->SetBounds(bounds); this->SetPlaneGeometry(geometry2D, 0); this->SetSpacing(spacing, true); this->SetEvenlySpaced(); // this->SetTimeBounds( geometry2D->GetTimeBounds() ); assert(this->GetIndexToWorldTransform() != geometry2D->GetIndexToWorldTransform()); // (**) see above. this->SetFrameOfReferenceID(geometry2D->GetFrameOfReferenceID()); this->SetImageGeometry(geometry2D->GetImageGeometry()); geometry2D->UnRegister(); } void mitk::SlicedGeometry3D::InitializePlanes(const mitk::BaseGeometry *geometry3D, mitk::PlaneGeometry::PlaneOrientation planeorientation, bool top, bool frontside, bool rotated) { m_ReferenceGeometry = geometry3D; PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New(); planeGeometry->InitializeStandardPlane(geometry3D, top, planeorientation, frontside, rotated); int worldAxis = planeorientation == PlaneGeometry::Sagittal ? 0 : planeorientation == PlaneGeometry::Frontal ? 1 : 2; // Inspired by: // http://www.na-mic.org/Wiki/index.php/Coordinate_System_Conversion_Between_ITK_and_Slicer3 mitk::AffineTransform3D::MatrixType matrix = geometry3D->GetIndexToWorldTransform()->GetMatrix(); matrix.GetVnlMatrix().normalize_columns(); mitk::AffineTransform3D::MatrixType::InternalMatrixType inverseMatrix = matrix.GetInverse(); int dominantAxis = itk::Function::Max3( inverseMatrix[0][worldAxis], inverseMatrix[1][worldAxis], inverseMatrix[2][worldAxis]); ScalarType viewSpacing = geometry3D->GetSpacing()[dominantAxis]; unsigned int slices = static_cast(geometry3D->GetExtent(dominantAxis)); #ifndef NDEBUG int upDirection = itk::Function::Sign(inverseMatrix[dominantAxis][worldAxis]); /// The normal vector of an imaginary plane that points from the world origin (bottom left back /// corner or the world, with the lowest physical coordinates) towards the inside of the volume, /// along the renderer axis. Length is the slice thickness. Vector3D worldPlaneNormal = inverseMatrix.get_row(dominantAxis) * (upDirection * viewSpacing); /// The normal of the standard plane geometry just created. Vector3D standardPlaneNormal = planeGeometry->GetNormal(); /// The standard plane must be parallel to the 'world plane'. The normal of the standard plane /// must point against the world plane if and only if 'top' is 'false'. The length of the /// standard plane normal must be equal to the slice thickness. assert((standardPlaneNormal - (top ? 1.0 : -1.0) * worldPlaneNormal).GetSquaredNorm() < 0.000001); #endif this->InitializeEvenlySpaced(planeGeometry, viewSpacing, slices); #ifndef NDEBUG /// The standard plane normal and the z axis vector of the sliced geometry must point in /// the same direction. Vector3D zAxisVector = this->GetAxisVector(2); Vector3D upscaledStandardPlaneNormal = standardPlaneNormal; upscaledStandardPlaneNormal *= slices; assert((zAxisVector - upscaledStandardPlaneNormal).GetSquaredNorm() < 0.000001); /// You can use this test is to check the handedness of the coordinate system of the current /// geometry. In principle, you can use either left- or right-handed coordinate systems, but /// you normally want it to be consistent, that is the handedness should be the same across /// the renderers of the same viewer. // ScalarType det = vnl_det(this->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix()); // MITK_DEBUG << "world axis: " << worldAxis << (det > 0 ? " ; right-handed" : " ; left-handed"); #endif } void mitk::SlicedGeometry3D::ReinitializePlanes(const Point3D ¢er, const Point3D &referencePoint) { // Need a reference frame to align the rotated planes if (!m_ReferenceGeometry) { return; } // Get first plane of plane stack PlaneGeometry *firstPlane = m_PlaneGeometries[0]; // If plane stack is empty, exit if (!firstPlane || dynamic_cast(firstPlane)) { return; } // Calculate the "directed" spacing when taking the plane (defined by its axes // vectors and normal) as the reference coordinate frame. // // This is done by calculating the radius of the ellipsoid defined by the // original volume spacing axes, in the direction of the respective axis of the // reference frame. mitk::Vector3D axis0 = firstPlane->GetAxisVector(0); mitk::Vector3D axis1 = firstPlane->GetAxisVector(1); mitk::Vector3D normal = firstPlane->GetNormal(); normal.Normalize(); Vector3D spacing; spacing[0] = this->CalculateSpacing(axis0); spacing[1] = this->CalculateSpacing(axis1); spacing[2] = this->CalculateSpacing(normal); Superclass::SetSpacing(spacing); // Now we need to calculate the number of slices in the plane's normal // direction, so that the entire volume is covered. This is done by first // calculating the dot product between the volume diagonal (the maximum // distance inside the volume) and the normal, and dividing this value by // the directed spacing calculated above. ScalarType directedExtent = std::abs(m_ReferenceGeometry->GetExtentInMM(0) * normal[0]) + std::abs(m_ReferenceGeometry->GetExtentInMM(1) * normal[1]) + std::abs(m_ReferenceGeometry->GetExtentInMM(2) * normal[2]); if (directedExtent >= spacing[2]) { m_Slices = static_cast(directedExtent / spacing[2] + 0.5); } else { m_Slices = 1; } // The origin of our "first plane" needs to be adapted to this new extent. // To achieve this, we first calculate the current distance to the volume's // center, and then shift the origin in the direction of the normal by the // difference between this distance and half of the new extent. double centerOfRotationDistance = firstPlane->SignedDistanceFromPlane(center); if (centerOfRotationDistance > 0) { firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * (centerOfRotationDistance - directedExtent / 2.0)); m_DirectionVector = normal; } else { firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * (directedExtent / 2.0 + centerOfRotationDistance)); m_DirectionVector = -normal; } // Now we adjust this distance according with respect to the given reference // point: we need to make sure that the point is touched by one slice of the // new slice stack. double referencePointDistance = firstPlane->SignedDistanceFromPlane(referencePoint); int referencePointSlice = static_cast(referencePointDistance / spacing[2]); double alignmentValue = referencePointDistance / spacing[2] - referencePointSlice; firstPlane->SetOrigin(firstPlane->GetOrigin() + normal * alignmentValue * spacing[2]); // Finally, we can clear the previous geometry stack and initialize it with // our re-initialized "first plane". m_PlaneGeometries.assign(m_Slices, PlaneGeometry::Pointer(nullptr)); if (m_Slices > 0) { m_PlaneGeometries[0] = firstPlane; } // Reinitialize SNC with new number of slices m_SliceNavigationController->GetSlice()->SetSteps(m_Slices); this->Modified(); } double mitk::SlicedGeometry3D::CalculateSpacing(const mitk::Vector3D &d) const { // Need the spacing of the underlying dataset / geometry if (!m_ReferenceGeometry) { return 1.0; } const mitk::Vector3D &spacing = m_ReferenceGeometry->GetSpacing(); return SlicedGeometry3D::CalculateSpacing(spacing, d); } double mitk::SlicedGeometry3D::CalculateSpacing(const mitk::Vector3D &spacing, const mitk::Vector3D &d) { // The following can be derived from the ellipsoid equation // // 1 = x^2/a^2 + y^2/b^2 + z^2/c^2 // // where (a,b,c) = spacing of original volume (ellipsoid radii) // and (x,y,z) = scaled coordinates of vector d (according to ellipsoid) // double scaling = d[0] * d[0] / (spacing[0] * spacing[0]) + d[1] * d[1] / (spacing[1] * spacing[1]) + d[2] * d[2] / (spacing[2] * spacing[2]); scaling = sqrt(scaling); return (sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]) / scaling); } mitk::Vector3D mitk::SlicedGeometry3D::AdjustNormal(const mitk::Vector3D &normal) const { TransformType::Pointer inverse = TransformType::New(); m_ReferenceGeometry->GetIndexToWorldTransform()->GetInverse(inverse); Vector3D transformedNormal = inverse->TransformVector(normal); transformedNormal.Normalize(); return transformedNormal; } void mitk::SlicedGeometry3D::SetImageGeometry(const bool isAnImageGeometry) { Superclass::SetImageGeometry(isAnImageGeometry); unsigned int s; for (s = 0; s < m_Slices; ++s) { mitk::BaseGeometry *geometry = m_PlaneGeometries[s]; if (geometry != nullptr) { geometry->SetImageGeometry(isAnImageGeometry); } } } void mitk::SlicedGeometry3D::ChangeImageGeometryConsideringOriginOffset(const bool isAnImageGeometry) { unsigned int s; for (s = 0; s < m_Slices; ++s) { mitk::BaseGeometry *geometry = m_PlaneGeometries[s]; if (geometry != nullptr) { geometry->ChangeImageGeometryConsideringOriginOffset(isAnImageGeometry); } } Superclass::ChangeImageGeometryConsideringOriginOffset(isAnImageGeometry); } bool mitk::SlicedGeometry3D::IsValidSlice(int s) const { return ((s >= 0) && (s < (int)m_Slices)); } const mitk::BaseGeometry *mitk::SlicedGeometry3D::GetReferenceGeometry() const { return m_ReferenceGeometry; } void mitk::SlicedGeometry3D::SetReferenceGeometry(const BaseGeometry *referenceGeometry) { m_ReferenceGeometry = referenceGeometry; std::vector::iterator it; for (it = m_PlaneGeometries.begin(); it != m_PlaneGeometries.end(); ++it) { (*it)->SetReferenceGeometry(referenceGeometry); } } bool mitk::SlicedGeometry3D::HasReferenceGeometry() const { return ( m_ReferenceGeometry != nullptr ); } void mitk::SlicedGeometry3D::PreSetSpacing(const mitk::Vector3D &aSpacing) { bool hasEvenlySpacedPlaneGeometry = false; mitk::Point3D origin; mitk::Vector3D rightDV, bottomDV; BoundingBox::BoundsArrayType bounds; // Check for valid spacing if (!(aSpacing[0] > 0 && aSpacing[1] > 0 && aSpacing[2] > 0)) { mitkThrow() << "You try to set a spacing with at least one element equal or " "smaller to \"0\". This might lead to a crash during rendering. Please double" " check your data!"; } // In case of evenly-spaced data: re-initialize instances of PlaneGeometry, // since the spacing influences them if ((m_EvenlySpaced) && (m_PlaneGeometries.size() > 0)) { const PlaneGeometry *planeGeometry = m_PlaneGeometries[0]; if (planeGeometry && !dynamic_cast(planeGeometry)) { this->WorldToIndex(planeGeometry->GetOrigin(), origin); this->WorldToIndex(planeGeometry->GetAxisVector(0), rightDV); this->WorldToIndex(planeGeometry->GetAxisVector(1), bottomDV); bounds = planeGeometry->GetBounds(); hasEvenlySpacedPlaneGeometry = true; } } BaseGeometry::_SetSpacing(aSpacing); mitk::PlaneGeometry::Pointer firstGeometry; // In case of evenly-spaced data: re-initialize instances of PlaneGeometry, // since the spacing influences them if (hasEvenlySpacedPlaneGeometry) { // create planeGeometry according to new spacing this->IndexToWorld(origin, origin); this->IndexToWorld(rightDV, rightDV); this->IndexToWorld(bottomDV, bottomDV); mitk::PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New(); planeGeometry->SetImageGeometry(this->GetImageGeometry()); planeGeometry->SetReferenceGeometry(m_ReferenceGeometry); // Store spacing, as Initialize... needs a pointer mitk::Vector3D lokalSpacing = this->GetSpacing(); planeGeometry->InitializeStandardPlane(rightDV.GetVnlVector(), bottomDV.GetVnlVector(), &lokalSpacing); planeGeometry->SetOrigin(origin); planeGeometry->SetBounds(bounds); firstGeometry = planeGeometry; } else if ((m_EvenlySpaced) && (m_PlaneGeometries.size() > 0)) { firstGeometry = m_PlaneGeometries[0].GetPointer(); } // clear and reserve PlaneGeometry::Pointer gnull = nullptr; m_PlaneGeometries.assign(m_Slices, gnull); if (m_Slices > 0) { m_PlaneGeometries[0] = firstGeometry; } this->Modified(); } void mitk::SlicedGeometry3D::SetSliceNavigationController(SliceNavigationController *snc) { m_SliceNavigationController = snc; } mitk::SliceNavigationController *mitk::SlicedGeometry3D::GetSliceNavigationController() { return m_SliceNavigationController; } void mitk::SlicedGeometry3D::SetEvenlySpaced(bool on) { if (m_EvenlySpaced != on) { m_EvenlySpaced = on; this->Modified(); } } void mitk::SlicedGeometry3D::SetDirectionVector(const mitk::Vector3D &directionVector) { Vector3D newDir = directionVector; newDir.Normalize(); if (newDir != m_DirectionVector) { m_DirectionVector = newDir; this->Modified(); } } // void // mitk::SlicedGeometry3D::SetTimeBounds( const mitk::TimeBounds& timebounds ) //{ // Superclass::SetTimeBounds( timebounds ); // // unsigned int s; // for ( s = 0; s < m_Slices; ++s ) // { // if(m_Geometry2Ds[s].IsNotNull()) // { // m_Geometry2Ds[s]->SetTimeBounds( timebounds ); // } // } // m_TimeBounds = timebounds; //} itk::LightObject::Pointer mitk::SlicedGeometry3D::InternalClone() const { Self::Pointer newGeometry = new SlicedGeometry3D(*this); newGeometry->UnRegister(); return newGeometry.GetPointer(); } void mitk::SlicedGeometry3D::PrintSelf(std::ostream &os, itk::Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << " EvenlySpaced: " << m_EvenlySpaced << std::endl; if (m_EvenlySpaced) { os << indent << " DirectionVector: " << m_DirectionVector << std::endl; } os << indent << " Slices: " << m_Slices << std::endl; os << std::endl; os << indent << " GetPlaneGeometry(0): "; if (this->GetPlaneGeometry(0) == nullptr) { os << "NULL" << std::endl; } else { this->GetPlaneGeometry(0)->Print(os, indent); } } void mitk::SlicedGeometry3D::ExecuteOperation(Operation *operation) { PlaneGeometry::Pointer geometry2D; ApplyTransformMatrixOperation *applyMatrixOp; Point3D center; switch (operation->GetOperationType()) { case OpNOTHING: break; case OpROTATE: if (m_EvenlySpaced) { // Need a reference frame to align the rotation if (m_ReferenceGeometry) { // Clear all generated geometries and then rotate only the first slice. // The other slices will be re-generated on demand // Save first slice PlaneGeometry::Pointer geometry2D = m_PlaneGeometries[0]; RotationOperation *rotOp = dynamic_cast(operation); // Generate a RotationOperation using the dataset center instead of // the supplied rotation center. This is necessary so that the rotated // zero-plane does not shift away. The supplied center is instead used // to adjust the slice stack afterwards. Point3D center = m_ReferenceGeometry->GetCenter(); RotationOperation centeredRotation( rotOp->GetOperationType(), center, rotOp->GetVectorOfRotation(), rotOp->GetAngleOfRotation()); // Rotate first slice geometry2D->ExecuteOperation(¢eredRotation); // Clear the slice stack and adjust it according to the center of // the dataset and the supplied rotation center (see documentation of // ReinitializePlanes) this->ReinitializePlanes(center, rotOp->GetCenterOfRotation()); geometry2D->SetSpacing(this->GetSpacing()); if (m_SliceNavigationController) { m_SliceNavigationController->SelectSliceByPoint(rotOp->GetCenterOfRotation()); m_SliceNavigationController->AdjustSliceStepperRange(); } BaseGeometry::ExecuteOperation(¢eredRotation); } else { // we also have to consider the case, that there is no reference geometry available. if (m_PlaneGeometries.size() > 0) { // Reach through to all slices in my container for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter) { // Test for empty slices, which can happen if evenly spaced geometry if ((*iter).IsNotNull()) { (*iter)->ExecuteOperation(operation); } } // rotate overall geometry RotationOperation *rotOp = dynamic_cast(operation); BaseGeometry::ExecuteOperation(rotOp); } } } else { // Reach through to all slices for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter) { (*iter)->ExecuteOperation(operation); } } break; case OpORIENT: if (m_EvenlySpaced) { // get operation data PlaneOperation *planeOp = dynamic_cast(operation); // Get first slice PlaneGeometry::Pointer planeGeometry = m_PlaneGeometries[0]; // Need a PlaneGeometry, a PlaneOperation and a reference frame to // carry out the re-orientation. If not all avaialble, stop here if (!m_ReferenceGeometry || (!planeGeometry || dynamic_cast(planeGeometry.GetPointer())) || !planeOp) { break; } // General Behavior: // Clear all generated geometries and then rotate only the first slice. // The other slices will be re-generated on demand // // 1st Step: Reorient Normal Vector of first plane // Point3D center = planeOp->GetPoint(); // m_ReferenceGeometry->GetCenter(); mitk::Vector3D currentNormal = planeGeometry->GetNormal(); mitk::Vector3D newNormal; if (planeOp->AreAxisDefined()) { // If planeOp was defined by one centerpoint and two axis vectors newNormal = CrossProduct(planeOp->GetAxisVec0(), planeOp->GetAxisVec1()); } else { // If planeOp was defined by one centerpoint and one normal vector newNormal = planeOp->GetNormal(); } // Get Rotation axis und angle currentNormal.Normalize(); newNormal.Normalize(); ScalarType rotationAngle = angle(currentNormal.GetVnlVector(), newNormal.GetVnlVector()); rotationAngle *= 180.0 / vnl_math::pi; // from rad to deg Vector3D rotationAxis = itk::CrossProduct(currentNormal, newNormal); if (std::abs(rotationAngle - 180) < mitk::eps) { // current Normal and desired normal are not linear independent!!(e.g 1,0,0 and -1,0,0). // Rotation Axis should be ANY vector that is 90� to current Normal mitk::Vector3D helpNormal; helpNormal = currentNormal; helpNormal[0] += 1; helpNormal[1] -= 1; helpNormal[2] += 1; helpNormal.Normalize(); rotationAxis = itk::CrossProduct(helpNormal, currentNormal); } RotationOperation centeredRotation(mitk::OpROTATE, center, rotationAxis, rotationAngle); // Rotate first slice planeGeometry->ExecuteOperation(¢eredRotation); // Reinitialize planes and select slice, if my rotations are all done. if (!planeOp->AreAxisDefined()) { // Clear the slice stack and adjust it according to the center of // rotation and plane position (see documentation of ReinitializePlanes) this->ReinitializePlanes(center, planeOp->GetPoint()); planeGeometry->SetSpacing(this->GetSpacing()); if (m_SliceNavigationController) { m_SliceNavigationController->SelectSliceByPoint(planeOp->GetPoint()); m_SliceNavigationController->AdjustSliceStepperRange(); } } // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry) BaseGeometry::ExecuteOperation(¢eredRotation); // // 2nd step. If axis vectors were defined, rotate the plane around its normal to fit these // if (planeOp->AreAxisDefined()) { mitk::Vector3D vecAxixNew = planeOp->GetAxisVec0(); vecAxixNew.Normalize(); mitk::Vector3D VecAxisCurr = planeGeometry->GetAxisVector(0); VecAxisCurr.Normalize(); ScalarType rotationAngle = angle(VecAxisCurr.GetVnlVector(), vecAxixNew.GetVnlVector()); rotationAngle = rotationAngle * 180 / PI; // Rad to Deg // we rotate around the normal of the plane, but we do not know, if we need to rotate clockwise // or anti-clockwise. So we rotate around the crossproduct of old and new Axisvector. // Since both axis vectors lie in the plane, the crossproduct is the planes normal or the negative planes // normal rotationAxis = itk::CrossProduct(VecAxisCurr, vecAxixNew); if (std::abs(rotationAngle - 180) < mitk::eps) { // current axisVec and desired axisVec are not linear independent!!(e.g 1,0,0 and -1,0,0). // Rotation Axis can be just plane Normal. (have to rotate by 180�) rotationAxis = newNormal; } // Perfom Rotation mitk::RotationOperation op(mitk::OpROTATE, center, rotationAxis, rotationAngle); planeGeometry->ExecuteOperation(&op); // Apply changes on first slice to whole slice stack this->ReinitializePlanes(center, planeOp->GetPoint()); planeGeometry->SetSpacing(this->GetSpacing()); if (m_SliceNavigationController) { m_SliceNavigationController->SelectSliceByPoint(planeOp->GetPoint()); m_SliceNavigationController->AdjustSliceStepperRange(); } // Also apply rotation on the slicedGeometry - Geometry3D (Bounding geometry) BaseGeometry::ExecuteOperation(&op); } } else { // Reach through to all slices for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter) { (*iter)->ExecuteOperation(operation); } } break; case OpRESTOREPLANEPOSITION: if (m_EvenlySpaced) { // Save first slice PlaneGeometry::Pointer planeGeometry = m_PlaneGeometries[0]; RestorePlanePositionOperation *restorePlaneOp = dynamic_cast(operation); // Need a PlaneGeometry, a PlaneOperation and a reference frame to // carry out the re-orientation if (m_ReferenceGeometry && (planeGeometry && dynamic_cast(planeGeometry.GetPointer()) == nullptr) && restorePlaneOp) { // Clear all generated geometries and then rotate only the first slice. // The other slices will be re-generated on demand // Rotate first slice planeGeometry->ExecuteOperation(restorePlaneOp); m_DirectionVector = restorePlaneOp->GetDirectionVector(); double centerOfRotationDistance = planeGeometry->SignedDistanceFromPlane(m_ReferenceGeometry->GetCenter()); if (centerOfRotationDistance <= 0) { m_DirectionVector = -m_DirectionVector; } Vector3D spacing = restorePlaneOp->GetSpacing(); Superclass::SetSpacing(spacing); // /*Now we need to calculate the number of slices in the plane's normal // direction, so that the entire volume is covered. This is done by first // calculating the dot product between the volume diagonal (the maximum // distance inside the volume) and the normal, and dividing this value by // the directed spacing calculated above.*/ ScalarType directedExtent = std::abs(m_ReferenceGeometry->GetExtentInMM(0) * m_DirectionVector[0]) + std::abs(m_ReferenceGeometry->GetExtentInMM(1) * m_DirectionVector[1]) + std::abs(m_ReferenceGeometry->GetExtentInMM(2) * m_DirectionVector[2]); if (directedExtent >= spacing[2]) { m_Slices = static_cast(directedExtent / spacing[2] + 0.5); } else { m_Slices = 1; } m_PlaneGeometries.assign(m_Slices, PlaneGeometry::Pointer(nullptr)); if (m_Slices > 0) { m_PlaneGeometries[0] = planeGeometry; } m_SliceNavigationController->GetSlice()->SetSteps(m_Slices); this->Modified(); // End Reinitialization if (m_SliceNavigationController) { m_SliceNavigationController->GetSlice()->SetPos(restorePlaneOp->GetPos()); m_SliceNavigationController->AdjustSliceStepperRange(); } BaseGeometry::ExecuteOperation(restorePlaneOp); } } else { // Reach through to all slices for (auto iter = m_PlaneGeometries.begin(); iter != m_PlaneGeometries.end(); ++iter) { (*iter)->ExecuteOperation(operation); } } break; case OpAPPLYTRANSFORMMATRIX: // Clear all generated geometries and then transform only the first slice. // The other slices will be re-generated on demand // Save first slice geometry2D = m_PlaneGeometries[0]; applyMatrixOp = dynamic_cast(operation); // Apply transformation to first plane geometry2D->ExecuteOperation(applyMatrixOp); // Generate a ApplyTransformMatrixOperation using the dataset center instead of // the supplied rotation center. The supplied center is instead used to adjust the // slice stack afterwards (see OpROTATE). center = m_ReferenceGeometry->GetCenter(); // Clear the slice stack and adjust it according to the center of // the dataset and the supplied rotation center (see documentation of // ReinitializePlanes) this->ReinitializePlanes(center, applyMatrixOp->GetReferencePoint()); BaseGeometry::ExecuteOperation(applyMatrixOp); break; default: // let handle by base class if we don't do anything BaseGeometry::ExecuteOperation(operation); } this->Modified(); } diff --git a/Modules/Core/test/mitkImageTest.cpp b/Modules/Core/test/mitkImageTest.cpp index bb6f1b6da4..ee4f316c67 100644 --- a/Modules/Core/test/mitkImageTest.cpp +++ b/Modules/Core/test/mitkImageTest.cpp @@ -1,555 +1,556 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // mitk includes #include "mitkException.h" #include "mitkIOUtil.h" #include "mitkImageGenerator.h" #include "mitkImagePixelReadAccessor.h" #include "mitkImageReadAccessor.h" #include "mitkPixelTypeMultiplex.h" #include #include #include #include #include #include "mitkImageSliceSelector.h" // itk includes #include #include // stl includes #include // vtk includes #include // Checks if reference count is correct after using GetVtkImageData() bool ImageVtkDataReferenceCheck(const char *fname) { const std::string filename = std::string(fname); try { mitk::Image::Pointer image = mitk::IOUtil::LoadImage(filename); MITK_TEST_CONDITION_REQUIRED(image.IsNotNull(), "Non-NULL image") vtkImageData *vtk = image->GetVtkImageData(); if (vtk == NULL) return false; } catch (...) { MITK_TEST_FAILED_MSG(<< "Could not read file for testing: " << filename); return false; } return true; } template void TestRandomPixelAccess(const mitk::PixelType /*ptype*/, mitk::Image::Pointer image, mitk::Point3D &point, mitk::ScalarType &value) { // generate a random point in world coordinates mitk::Point3D xMax, yMax, zMax, xMaxIndex, yMaxIndex, zMaxIndex; xMaxIndex.Fill(0.0f); yMaxIndex.Fill(0.0f); zMaxIndex.Fill(0.0f); xMaxIndex[0] = image->GetLargestPossibleRegion().GetSize()[0]; yMaxIndex[1] = image->GetLargestPossibleRegion().GetSize()[1]; zMaxIndex[2] = image->GetLargestPossibleRegion().GetSize()[2]; image->GetGeometry()->IndexToWorld(xMaxIndex, xMax); image->GetGeometry()->IndexToWorld(yMaxIndex, yMax); image->GetGeometry()->IndexToWorld(zMaxIndex, zMax); MITK_INFO << "Origin " << image->GetGeometry()->GetOrigin()[0] << " " << image->GetGeometry()->GetOrigin()[1] << " " << image->GetGeometry()->GetOrigin()[2] << ""; MITK_INFO << "MaxExtend " << xMax[0] << " " << yMax[1] << " " << zMax[2] << ""; itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randomGenerator = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); randomGenerator->Initialize(std::rand()); // initialize with random value, to get sensible random points for the image point[0] = randomGenerator->GetUniformVariate(image->GetGeometry()->GetOrigin()[0], xMax[0]); point[1] = randomGenerator->GetUniformVariate(image->GetGeometry()->GetOrigin()[1], yMax[1]); point[2] = randomGenerator->GetUniformVariate(image->GetGeometry()->GetOrigin()[2], zMax[2]); MITK_INFO << "RandomPoint " << point[0] << " " << point[1] << " " << point[2] << ""; // test values and max/min mitk::ScalarType imageMin = image->GetStatistics()->GetScalarValueMin(); mitk::ScalarType imageMax = image->GetStatistics()->GetScalarValueMax(); // test accessing PixelValue with coordinate leading to a negative index const mitk::Point3D geom_origin = image->GetGeometry()->GetOrigin(); const mitk::Point3D geom_center = image->GetGeometry()->GetCenter(); // shift position from origin outside of the image ( in the opposite direction to [center-origin] vector which points // in the inside) mitk::Point3D position = geom_origin + (geom_origin - geom_center); MITK_INFO << "Testing access outside of the image"; unsigned int dim = image->GetDimension(); if (dim == 3 || dim == 4) { mitk::ImagePixelReadAccessor imAccess3(image, image->GetVolumeData(0)); // Comparison ?>=0 not needed since all position[i] and timestep are unsigned int // (position[0]>=0 && position[1] >=0 && position[2]>=0 && timestep>=0) // bug-11978 : we still need to catch index with negative values if (point[0] < 0 || point[1] < 0 || point[2] < 0) { MITK_WARN << "Given position (" << point << ") is out of image range, returning 0."; } else { value = static_cast(imAccess3.GetPixelByWorldCoordinates(point)); MITK_TEST_CONDITION((value >= imageMin && value <= imageMax), "Value returned is between max/min"); } itk::Index<3> itkIndex; image->GetGeometry()->WorldToIndex(position, itkIndex); MITK_TEST_FOR_EXCEPTION_BEGIN(mitk::Exception); imAccess3.GetPixelByIndexSafe(itkIndex); MITK_TEST_FOR_EXCEPTION_END(mitk::Exception); } MITK_INFO << imageMin << " " << imageMax << " " << value << ""; } class mitkImageTestClass { public: void SetClonedGeometry_None_ClonedEqualInput() { mitk::Image::Pointer image = mitk::ImageGenerator::GenerateRandomImage(100, 100, 100, 1, 0.2, 0.3, 0.4); //----------------- // geometry information for image mitk::Point3D origin; mitk::Vector3D right, bottom; mitk::Vector3D spacing; mitk::FillVector3D(origin, 17.0, 19.92, 7.83); mitk::FillVector3D(right, 1.0, 2.0, 3.0); mitk::FillVector3D(bottom, 0.0, -3.0, 2.0); mitk::FillVector3D(spacing, 0.78, 0.91, 2.23); // InitializeStandardPlane(rightVector, downVector, spacing) mitk::PlaneGeometry::Pointer planegeometry = mitk::PlaneGeometry::New(); planegeometry->InitializeStandardPlane(100, 100, right, bottom, &spacing); planegeometry->SetOrigin(origin); planegeometry->ChangeImageGeometryConsideringOriginOffset(true); image->SetClonedGeometry(planegeometry); mitk::BaseGeometry::Pointer imageGeometry = image->GetGeometry(); MITK_ASSERT_EQUAL(imageGeometry, planegeometry, "Matrix elements of cloned matrix equal original matrix"); } }; int mitkImageTest(int argc, char *argv[]) { MITK_TEST_BEGIN(mitkImageTest); mitkImageTestClass tester; tester.SetClonedGeometry_None_ClonedEqualInput(); // Create Image out of nowhere mitk::Image::Pointer imgMem = mitk::Image::New(); mitk::PixelType pt = mitk::MakeScalarPixelType(); unsigned int dim[] = {100, 100, 20}; MITK_TEST_CONDITION_REQUIRED(imgMem.IsNotNull(), "An image was created. "); // Initialize image imgMem->Initialize(pt, 3, dim); MITK_TEST_CONDITION_REQUIRED(imgMem->IsInitialized(), "Image::IsInitialized() ?"); MITK_TEST_CONDITION_REQUIRED(imgMem->GetPixelType() == pt, "PixelType was set correctly."); int *p = NULL; int *p2 = NULL; try { mitk::ImageReadAccessor imgMemAcc(imgMem); p = (int *)imgMemAcc.GetData(); } catch (mitk::Exception &e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION(p != NULL, "GetData() returned not-NULL pointer."); // filling image const unsigned int size = dim[0] * dim[1] * dim[2]; for (unsigned int i = 0; i < size; ++i, ++p) *p = (signed int)i; // Getting it again and compare with filled values: try { mitk::ImageReadAccessor imgMemAcc(imgMem); p2 = (int *)imgMemAcc.GetData(); } catch (mitk::Exception &e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION(p2 != NULL, "GetData() returned not-NULL pointer."); bool isEqual = true; for (unsigned int i = 0; i < size; ++i, ++p2) { if (*p2 != (signed int)i) { isEqual = false; } } MITK_TEST_CONDITION(isEqual, "The values previously set as data are correct [pixelwise comparison]."); // Testing GetSliceData() and compare with filled values: try { mitk::ImageReadAccessor imgMemAcc(imgMem, imgMem->GetSliceData(dim[2] / 2)); p2 = (int *)imgMemAcc.GetData(); } catch (mitk::Exception &e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION_REQUIRED(p2 != NULL, "Valid slice data returned"); unsigned int xy_size = dim[0] * dim[1]; unsigned int start_mid_slice = (dim[2] / 2) * xy_size; isEqual = true; for (unsigned int i = 0; i < xy_size; ++i, ++p2) { if (*p2 != (signed int)(i + start_mid_slice)) { isEqual = false; } } MITK_TEST_CONDITION(isEqual, "The SliceData are correct [pixelwise comparison]. "); imgMem = mitk::Image::New(); // testing re-initialization of test image mitk::PixelType pType = mitk::MakePixelType(); imgMem->Initialize(pType, 3, dim); MITK_TEST_CONDITION_REQUIRED(imgMem->GetDimension() == 3, "Testing initialization parameter dimension!"); MITK_TEST_CONDITION_REQUIRED(imgMem->GetPixelType() == pType, "Testing initialization parameter pixeltype!"); MITK_TEST_CONDITION_REQUIRED( imgMem->GetDimension(0) == dim[0] && imgMem->GetDimension(1) == dim[1] && imgMem->GetDimension(2) == dim[2], "Testing initialization of dimensions!"); MITK_TEST_CONDITION(imgMem->IsInitialized(), "Image is initialized."); // Setting volume again: try { mitk::ImageReadAccessor imgMemAcc(imgMem); imgMem->SetVolume(imgMemAcc.GetData()); } catch (mitk::Exception &e) { MITK_ERROR << e.what(); } //----------------- // geometry information for image mitk::Point3D origin; mitk::Vector3D right, bottom; mitk::Vector3D spacing; mitk::FillVector3D(origin, 17.0, 19.92, 7.83); mitk::FillVector3D(right, 1.0, 2.0, 3.0); mitk::FillVector3D(bottom, 0.0, -3.0, 2.0); mitk::FillVector3D(spacing, 0.78, 0.91, 2.23); // InitializeStandardPlane(rightVector, downVector, spacing) mitk::PlaneGeometry::Pointer planegeometry = mitk::PlaneGeometry::New(); planegeometry->InitializeStandardPlane(100, 100, right, bottom, &spacing); planegeometry->SetOrigin(origin); + planegeometry->SetImageGeometry(true); // Testing Initialize(const mitk::PixelType& type, const mitk::Geometry3D& geometry, unsigned int slices) with // PlaneGeometry and GetData(): "; imgMem->Initialize(mitk::MakePixelType(), *planegeometry); MITK_TEST_CONDITION_REQUIRED( imgMem->GetGeometry()->GetOrigin() == static_cast(planegeometry)->GetOrigin(), "Testing correct setting of geometry via initialize!"); try { mitk::ImageReadAccessor imgMemAcc(imgMem); p = (int *)imgMemAcc.GetData(); } catch (mitk::Exception &e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION_REQUIRED(p != NULL, "GetData() returned valid pointer."); // Testing Initialize(const mitk::PixelType& type, int sDim, const mitk::PlaneGeometry& geometry) and GetData(): "; imgMem->Initialize(mitk::MakePixelType(), 40, *planegeometry); try { mitk::ImageReadAccessor imgMemAcc(imgMem); p = (int *)imgMemAcc.GetData(); } catch (mitk::Exception &e) { MITK_ERROR << e.what(); } MITK_TEST_CONDITION_REQUIRED(p != NULL, "GetData() returned valid pointer."); //----------------- // testing origin information and methods MITK_TEST_CONDITION_REQUIRED(mitk::Equal(imgMem->GetGeometry()->GetOrigin(), origin), "Testing correctness of origin via GetGeometry()->GetOrigin(): "); // Setting origin via SetOrigin(origin): "; mitk::FillVector3D(origin, 37.0, 17.92, 27.83); imgMem->SetOrigin(origin); // Test origin MITK_TEST_CONDITION_REQUIRED(mitk::Equal(imgMem->GetGeometry()->GetOrigin(), origin), "Testing correctness of changed origin via GetGeometry()->GetOrigin(): "); MITK_TEST_CONDITION_REQUIRED( mitk::Equal(imgMem->GetSlicedGeometry()->GetPlaneGeometry(0)->GetOrigin(), origin), "Testing correctness of changed origin via GetSlicedGeometry()->GetPlaneGeometry(0)->GetOrigin(): "); //----------------- // testing spacing information and methodsunsigned int dim[]={100,100,20}; MITK_TEST_CONDITION_REQUIRED(mitk::Equal(imgMem->GetGeometry()->GetSpacing(), spacing), "Testing correct spacing from Geometry3D!"); mitk::FillVector3D(spacing, 7.0, 0.92, 1.83); imgMem->SetSpacing(spacing); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(imgMem->GetGeometry()->GetSpacing(), spacing), "Testing correctness of changed spacing via GetGeometry()->GetSpacing(): "); MITK_TEST_CONDITION_REQUIRED( mitk::Equal(imgMem->GetSlicedGeometry()->GetPlaneGeometry(0)->GetSpacing(), spacing), "Testing correctness of changed spacing via GetSlicedGeometry()->GetPlaneGeometry(0)->GetSpacing(): "); mitk::Image::Pointer vecImg = mitk::Image::New(); try { mitk::ImageReadAccessor imgMemAcc(imgMem); vecImg->Initialize(imgMem->GetPixelType(), *imgMem->GetGeometry(), 2 /* #channels */, 0 /*tDim*/); vecImg->SetImportChannel(const_cast(imgMemAcc.GetData()), 0, mitk::Image::CopyMemory); vecImg->SetImportChannel(const_cast(imgMemAcc.GetData()), 1, mitk::Image::CopyMemory); mitk::ImageReadAccessor vecImgAcc(vecImg); mitk::ImageReadAccessor vecImgAcc0(vecImg, vecImg->GetChannelData(0)); mitk::ImageReadAccessor vecImgAcc1(vecImg, vecImg->GetChannelData(1)); MITK_TEST_CONDITION_REQUIRED(vecImgAcc0.GetData() != NULL && vecImgAcc1.GetData() != NULL, "Testing set and return of channel data!"); MITK_TEST_CONDITION_REQUIRED(vecImg->IsValidSlice(0, 0, 1), ""); MITK_TEST_OUTPUT(<< " Testing whether CopyMemory worked"); MITK_TEST_CONDITION_REQUIRED(imgMemAcc.GetData() != vecImgAcc.GetData(), ""); MITK_TEST_OUTPUT(<< " Testing destruction after SetImportChannel"); vecImg = NULL; MITK_TEST_CONDITION_REQUIRED(vecImg.IsNull(), "testing destruction!"); } catch (mitk::Exception &e) { MITK_ERROR << e.what(); } //----------------- MITK_TEST_OUTPUT(<< "Testing initialization via vtkImageData"); MITK_TEST_OUTPUT(<< " Setting up vtkImageData"); vtkImageData *vtkimage = vtkImageData::New(); vtkimage->Initialize(); vtkimage->SetDimensions(2, 3, 4); double vtkorigin[] = {-350, -358.203, -1363.5}; vtkimage->SetOrigin(vtkorigin); mitk::Point3D vtkoriginAsMitkPoint; mitk::vtk2itk(vtkorigin, vtkoriginAsMitkPoint); double vtkspacing[] = {1.367, 1.367, 2}; vtkimage->SetSpacing(vtkspacing); vtkimage->AllocateScalars(VTK_SHORT, 1); std::cout << "[PASSED]" << std::endl; MITK_TEST_OUTPUT(<< " Testing mitk::Image::Initialize(vtkImageData*, ...)"); mitk::Image::Pointer mitkByVtkImage = mitk::Image::New(); mitkByVtkImage->Initialize(vtkimage); MITK_TEST_CONDITION_REQUIRED(mitkByVtkImage->IsInitialized(), ""); vtkimage->Delete(); MITK_TEST_OUTPUT(<< " Testing whether spacing has been correctly initialized from vtkImageData"); mitk::Vector3D spacing2 = mitkByVtkImage->GetGeometry()->GetSpacing(); mitk::Vector3D vtkspacingAsMitkVector; mitk::vtk2itk(vtkspacing, vtkspacingAsMitkVector); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(spacing2, vtkspacingAsMitkVector), ""); MITK_TEST_OUTPUT( << " Testing whether GetSlicedGeometry(0)->GetOrigin() has been correctly initialized from vtkImageData"); mitk::Point3D origin2 = mitkByVtkImage->GetSlicedGeometry(0)->GetOrigin(); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(origin2, vtkoriginAsMitkPoint), ""); MITK_TEST_OUTPUT(<< " Testing whether GetGeometry()->GetOrigin() has been correctly initialized from vtkImageData"); origin2 = mitkByVtkImage->GetGeometry()->GetOrigin(); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(origin2, vtkoriginAsMitkPoint), ""); MITK_TEST_OUTPUT(<< " Testing if vtkOrigin is (0, 0, 0). This behaviour is due to historical development of MITK. " "Aslo see bug 5050!"); vtkImageData *vtkImage = imgMem->GetVtkImageData(); auto vtkOrigin = vtkImage->GetOrigin(); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(vtkOrigin[0], 0), "testing vtkOrigin[0] to be 0"); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(vtkOrigin[1], 0), "testing vtkOrigin[1] to be 0"); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(vtkOrigin[2], 0), "testing vtkOrigin[2] to be 0"); // TODO test the following initializers on channel-incorporation // void mitk::Image::Initialize(const mitk::PixelType& type, unsigned int dimension, unsigned int *dimensions, // unsigned int channels) // void mitk::Image::Initialize(const mitk::PixelType& type, int sDim, const mitk::PlaneGeometry& geometry2d, bool // flipped, unsigned int channels, int tDim ) // void mitk::Image::Initialize(const mitk::Image* image) // void mitk::Image::Initialize(const mitkIpPicDescriptor* pic, int channels, int tDim, int sDim) // mitk::Image::Pointer vecImg = mitk::Image::New(); // vecImg->Initialize(PixelType(typeid(float), 6, itk::ImageIOBase::SYMMETRICSECONDRANKTENSOR), // *imgMem->GetGeometry(), 2 /* #channels */, 0 /*tDim*/, false /*shiftBoundingBoxMinimumToZero*/ ); // vecImg->Initialize(PixelType(typeid(itk::Vector)), *imgMem->GetGeometry(), 2 /* #channels */, 0 /*tDim*/, // false /*shiftBoundingBoxMinimumToZero*/ ); // testing access by index coordinates and by world coordinates MITK_TEST_CONDITION_REQUIRED(argc == 2, "Check if test image is accessible!"); const std::string filename = std::string(argv[1]); mitk::Image::Pointer image; try { image = mitk::IOUtil::LoadImage(filename); MITK_TEST_CONDITION_REQUIRED(image.IsNotNull(), "Non-NULL image") } catch (...) { MITK_TEST_FAILED_MSG(<< "Could not read file for testing: " << filename); return 0; } mitk::Point3D point; mitk::ScalarType value = -1.; mitkPixelTypeMultiplex3( TestRandomPixelAccess, image->GetImageDescriptor()->GetChannelTypeById(0), image, point, value) { // testing the clone method of mitk::Image mitk::Image::Pointer cloneImage = image->Clone(); MITK_TEST_CONDITION_REQUIRED(cloneImage->GetDimension() == image->GetDimension(), "Clone (testing dimension)"); MITK_TEST_CONDITION_REQUIRED(cloneImage->GetPixelType() == image->GetPixelType(), "Clone (testing pixel type)"); // After cloning an image the geometry of both images should be equal too MITK_TEST_CONDITION_REQUIRED(cloneImage->GetGeometry()->GetOrigin() == image->GetGeometry()->GetOrigin(), "Clone (testing origin)"); MITK_TEST_CONDITION_REQUIRED(cloneImage->GetGeometry()->GetSpacing() == image->GetGeometry()->GetSpacing(), "Clone (testing spacing)"); MITK_TEST_CONDITION_REQUIRED( mitk::MatrixEqualElementWise(cloneImage->GetGeometry()->GetIndexToWorldTransform()->GetMatrix(), image->GetGeometry()->GetIndexToWorldTransform()->GetMatrix()), "Clone (testing transformation matrix)"); MITK_TEST_CONDITION_REQUIRED( mitk::MatrixEqualElementWise(cloneImage->GetTimeGeometry() ->GetGeometryForTimeStep(cloneImage->GetDimension(3) - 1) ->GetIndexToWorldTransform() ->GetMatrix(), cloneImage->GetTimeGeometry() ->GetGeometryForTimeStep(image->GetDimension(3) - 1) ->GetIndexToWorldTransform() ->GetMatrix()), "Clone(testing time sliced geometry)"); for (unsigned int i = 0u; i < cloneImage->GetDimension(); ++i) { MITK_TEST_CONDITION_REQUIRED(cloneImage->GetDimension(i) == image->GetDimension(i), "Clone (testing dimension " << i << ")"); } } // access via itk if (image->GetDimension() > 3) // CastToItk only works with 3d images so we need to check for 4d images { mitk::ImageTimeSelector::Pointer selector = mitk::ImageTimeSelector::New(); selector->SetTimeNr(0); selector->SetInput(image); selector->Update(); image = selector->GetOutput(); } if (image->GetDimension() == 3) { typedef itk::Image ItkFloatImage3D; ItkFloatImage3D::Pointer itkimage; try { mitk::CastToItkImage(image, itkimage); MITK_TEST_CONDITION_REQUIRED(itkimage.IsNotNull(), "Test conversion to itk::Image!"); } catch (std::exception &e) { MITK_INFO << e.what(); } mitk::Point3D itkPhysicalPoint; image->GetGeometry()->WorldToItkPhysicalPoint(point, itkPhysicalPoint); MITK_INFO << "ITKPoint " << itkPhysicalPoint[0] << " " << itkPhysicalPoint[1] << " " << itkPhysicalPoint[2] << ""; mitk::Point3D backTransformedPoint; image->GetGeometry()->ItkPhysicalPointToWorld(itkPhysicalPoint, backTransformedPoint); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(point, backTransformedPoint), "Testing world->itk-physical->world consistency"); itk::Index<3> idx; bool status = itkimage->TransformPhysicalPointToIndex(itkPhysicalPoint, idx); MITK_INFO << "ITK Index " << idx[0] << " " << idx[1] << " " << idx[2] << ""; if (status && value != -1.) { float valByItk = itkimage->GetPixel(idx); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(valByItk, value), "Compare value of pixel returned by mitk in comparison to itk"); } else { MITK_WARN << "Index is out buffered region!"; } } else { MITK_INFO << "Image does not contain three dimensions, some test cases are skipped!"; } // clone generated 3D image with one slice in z direction (cf. bug 11058) unsigned int *threeDdim = new unsigned int[3]; threeDdim[0] = 100; threeDdim[1] = 200; threeDdim[2] = 1; mitk::Image::Pointer threeDImage = mitk::Image::New(); threeDImage->Initialize(mitk::MakeScalarPixelType(), 3, threeDdim); mitk::Image::Pointer cloneThreeDImage = threeDImage->Clone(); // check that the clone image has the same dimensionality as the source image MITK_TEST_CONDITION_REQUIRED(cloneThreeDImage->GetDimension() == 3, "Testing if the clone image initializes with 3D!"); MITK_TEST_CONDITION_REQUIRED(ImageVtkDataReferenceCheck(argv[1]), "Checking reference count of Image after using GetVtkImageData()"); MITK_TEST_END(); }