diff --git a/Modules/Core/include/mitkExtractSliceFilter.h b/Modules/Core/include/mitkExtractSliceFilter.h index c9a44b701e..3c1c4af499 100644 --- a/Modules/Core/include/mitkExtractSliceFilter.h +++ b/Modules/Core/include/mitkExtractSliceFilter.h @@ -1,198 +1,204 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkExtractSliceFilter_h_Included #define mitkExtractSliceFilter_h_Included #include "MitkCoreExports.h" #include "mitkImageToImageFilter.h" #include #include #include #include #include #include #include namespace mitk { /** \brief ExtractSliceFilter extracts a 2D abitrary oriented slice from a 3D volume. The filter can reslice in all orthogonal planes such as sagittal, coronal and axial, and is also able to reslice an arbitrary oriented oblique plane. Curved planes are specified via an AbstractTransformGeometry as the input worldgeometry. Additionally the filter extracts the specified component of a multi-component input image. This is done only if the caller requests an mitk::Image output ('m_VtkOutputRequested' set to false). The default component to be extracted is '0'. The convenient workflow is: 1. Set an image as input. 2. Set the worldPlaneGeometry. This defines a grid where the slice is being extracted 3. And then start the pipeline. There are a few more properties that can be set to modify the behavior of the slicing. The properties are: - interpolation mode either Nearestneighbor, Linear or Cubic. - a transform this is a convenient way to adapt the reslice axis for the case that the image is transformed e.g. rotated. - time step the time step in a times volume. - the component to extract from a multi-component input image - vtkoutputrequested, to define whether an mitk::image should be initialized - resample by geometry whether the resampling grid corresponds to the specs of the worldgeometry or is directly derived from the input image By default the properties are set to: - interpolation mode Nearestneighbor. - a transform nullptr (No transform is set). - time step 0. - component 0. - resample by geometry false (Corresponds to input image). */ class MITKCORE_EXPORT ExtractSliceFilter : public ImageToImageFilter { public: mitkClassMacro(ExtractSliceFilter, ImageToImageFilter); itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(Self, vtkImageReslice *); /** \brief Set the axis where to reslice at.*/ void SetWorldGeometry(const PlaneGeometry *geometry) { this->m_WorldGeometry = geometry; this->Modified(); } /** \brief Set the time step in the 4D volume */ void SetTimeStep(unsigned int timestep) { m_TimeStep = timestep; } unsigned int GetTimeStep() { return m_TimeStep; } /** \brief Set the component of an image to be extracted */ void SetComponent(unsigned int component) { m_Component = component; } /** \brief Set a transform for the reslice axes. * This transform is needed if the image volume itself is transformed. (Effects the reslice axis) */ void SetResliceTransformByGeometry(const BaseGeometry *transform) { this->m_ResliceTransform = transform; } /** \brief Resampling grid corresponds to: false->image true->worldgeometry*/ void SetInPlaneResampleExtentByGeometry(bool inPlaneResampleExtentByGeometry) { this->m_InPlaneResampleExtentByGeometry = inPlaneResampleExtentByGeometry; } /** \brief Sets the output dimension of the slice*/ void SetOutputDimensionality(unsigned int dimension) { this->m_OutputDimension = dimension; } /** \brief Set the spacing in z direction manually. * Required if the outputDimension is > 2. */ void SetOutputSpacingZDirection(double zSpacing) { this->m_ZSpacing = zSpacing; } /** \brief Set the extent in pixel for direction z manualy. Required if the output dimension is > 2. */ void SetOutputExtentZDirection(int zMin, int zMax) { this->m_ZMin = zMin; this->m_ZMax = zMax; } /** \brief Get the bounding box of the slice [xMin, xMax, yMin, yMax, zMin, zMax] * The method uses the input of the filter to calculate the bounds. * It is recommended to use * GetClippedPlaneBounds(const BaseGeometry*, const PlaneGeometry*, double*) * if you are not sure about the input. */ bool GetClippedPlaneBounds(double bounds[6]); /** \brief Get the bounding box of the slice [xMin, xMax, yMin, yMax, zMin, zMax]*/ bool GetClippedPlaneBounds(const BaseGeometry *boundingGeometry, const PlaneGeometry *planeGeometry, double *bounds); /** \brief Get the spacing of the slice. returns mitk::ScalarType[2] */ mitk::ScalarType *GetOutputSpacing(); /** \brief Get Output as vtkImageData. * Note: * SetVtkOutputRequest(true) has to be called at least once before * GetVtkOutput(). Otherwise the output is empty for the first update step. */ vtkImageData *GetVtkOutput() { m_VtkOutputRequested = true; return m_Reslicer->GetOutput(); } /** Set VtkOutPutRequest to suppress the convertion of the image. * It is suggested to use this with GetVtkOutput(). * Note: * SetVtkOutputRequest(true) has to be called at least once before * GetVtkOutput(). Otherwise the output is empty for the first update step. */ void SetVtkOutputRequest(bool isRequested) { m_VtkOutputRequested = isRequested; } /** \brief Get the reslices axis matrix. * Note: the axis are recalculated when calling SetResliceTransformByGeometry. */ vtkMatrix4x4 *GetResliceAxes() { return this->m_Reslicer->GetResliceAxes(); } void SetBackgroundLevel(double backgroundLevel) { m_BackgroundLevel = backgroundLevel; } enum ResliceInterpolation { RESLICE_NEAREST = 0, RESLICE_LINEAR = 1, RESLICE_CUBIC = 3 }; void SetInterpolationMode(ExtractSliceFilter::ResliceInterpolation interpolation) { this->m_InterpolationMode = interpolation; } protected: ExtractSliceFilter(vtkImageReslice *reslicer = nullptr); ~ExtractSliceFilter() override; void GenerateData() override; void GenerateOutputInformation() override; void GenerateInputRequestedRegion() override; const PlaneGeometry *m_WorldGeometry; vtkSmartPointer m_Reslicer; unsigned int m_TimeStep; unsigned int m_OutputDimension; double m_ZSpacing; int m_ZMin; int m_ZMax; ResliceInterpolation m_InterpolationMode; - BaseGeometry::ConstPointer m_ResliceTransform; - bool m_InPlaneResampleExtentByGeometry; // Resampling grid corresponds to: false->image true->worldgeometry mitk::ScalarType *m_OutPutSpacing; bool m_VtkOutputRequested; double m_BackgroundLevel; unsigned int m_Component; + + private: + BaseGeometry::ConstPointer m_ResliceTransform; + /* Axis vectors of the relevant geometry. Set in GenerateOutputInformation() and also used in GenerateData().*/ + Vector3D m_Right, m_Bottom; + /* Bounds of the relevant plane. Set in GenerateOutputInformation() and also used in GenerateData().*/ + int m_XMin, m_XMax, m_YMin, m_YMax; + }; } #endif // mitkExtractSliceFilter_h_Included diff --git a/Modules/Core/src/Algorithms/mitkExtractSliceFilter.cpp b/Modules/Core/src/Algorithms/mitkExtractSliceFilter.cpp index 0cf7695e50..f7c9c34e0d 100644 --- a/Modules/Core/src/Algorithms/mitkExtractSliceFilter.cpp +++ b/Modules/Core/src/Algorithms/mitkExtractSliceFilter.cpp @@ -1,476 +1,510 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkExtractSliceFilter.h" #include #include #include #include #include #include #include -mitk::ExtractSliceFilter::ExtractSliceFilter(vtkImageReslice *reslicer) +mitk::ExtractSliceFilter::ExtractSliceFilter(vtkImageReslice *reslicer): m_XMin(0), m_XMax(0), m_YMin(0), m_YMax(0) { if (reslicer == nullptr) { m_Reslicer = vtkSmartPointer::New(); } else { m_Reslicer = reslicer; } m_TimeStep = 0; m_Reslicer->ReleaseDataFlagOn(); m_InterpolationMode = ExtractSliceFilter::RESLICE_NEAREST; m_ResliceTransform = nullptr; m_InPlaneResampleExtentByGeometry = false; m_OutPutSpacing = new mitk::ScalarType[2]; m_OutputDimension = 2; m_ZSpacing = 1.0; m_ZMin = 0; m_ZMax = 0; m_VtkOutputRequested = false; m_BackgroundLevel = -32768.0; m_Component = 0; } mitk::ExtractSliceFilter::~ExtractSliceFilter() { m_ResliceTransform = nullptr; m_WorldGeometry = nullptr; delete[] m_OutPutSpacing; } void mitk::ExtractSliceFilter::GenerateOutputInformation() { - Superclass::GenerateOutputInformation(); - // TODO try figure out how to set the specs of the slice before it is actually extracted - /*Image::Pointer output = this->GetOutput(); Image::ConstPointer input = this->GetInput(); - if (input.IsNull()) return; - unsigned int dimensions[2]; - dimensions[0] = m_WorldGeometry->GetExtent(0); - dimensions[1] = m_WorldGeometry->GetExtent(1); - output->Initialize(input->GetPixelType(), 2, dimensions, 1);*/ + + if (input.IsNull()) + { + return; + } + + if (nullptr == m_WorldGeometry) + { + return; + } + + Vector3D right, bottom; + double widthInMM, heightInMM; + Vector2D extent; + + // set the geometry from current worldgeometry for the reusultimage + // this is needed that the image has the correct mitk geometry + // the sliceGeometry is the Geometry of the result slice + PlaneGeometry::Pointer sliceGeometry = m_WorldGeometry->Clone(); + + sliceGeometry->GetIndexToWorldTransform()->SetMatrix(m_WorldGeometry->GetIndexToWorldTransform()->GetMatrix()); + + // the origin of the worldGeometry is transformed to center based coordinates to be an imageGeometry + Point3D sliceOrigin = sliceGeometry->GetOrigin(); + + auto abstractGeometry = + dynamic_cast(m_WorldGeometry); + + auto planeGeometry = dynamic_cast(m_WorldGeometry); + + if (abstractGeometry != nullptr) + { + extent[0] = abstractGeometry->GetParametricExtent(0); + extent[1] = abstractGeometry->GetParametricExtent(1); + + widthInMM = abstractGeometry->GetParametricExtentInMM(0); + heightInMM = abstractGeometry->GetParametricExtentInMM(1); + + right = abstractGeometry->GetPlane()->GetAxisVector(0); + bottom = abstractGeometry->GetPlane()->GetAxisVector(1); + } + else + { + if (planeGeometry != nullptr) + { + // if the worldGeomatry is a PlaneGeometry everything is straight forward + right = planeGeometry->GetAxisVector(0); + bottom = planeGeometry->GetAxisVector(1); + + if (m_InPlaneResampleExtentByGeometry) + { + // Resampling grid corresponds to the current world geometry. This + // means that the spacing of the output 2D image depends on the + // currently selected world geometry, and *not* on the image itself. + extent[0] = m_WorldGeometry->GetExtent(0); + extent[1] = m_WorldGeometry->GetExtent(1); + } + else + { + const TimeGeometry *inputTimeGeometry = input->GetTimeGeometry(); + if ((inputTimeGeometry == nullptr) || (inputTimeGeometry->CountTimeSteps() <= 0)) + { + itkWarningMacro(<< "Error reading input image TimeGeometry."); + return; + } + + // Resampling grid corresponds to the input geometry. This means that + // the spacing of the output 2D image is directly derived from the + // associated input image, regardless of the currently selected world + // geometry. + Vector3D rightInIndex, bottomInIndex; + inputTimeGeometry->GetGeometryForTimeStep(m_TimeStep)->WorldToIndex(right, rightInIndex); + inputTimeGeometry->GetGeometryForTimeStep(m_TimeStep)->WorldToIndex(bottom, bottomInIndex); + extent[0] = rightInIndex.GetNorm(); + extent[1] = bottomInIndex.GetNorm(); + } + + // Get the extent of the current world geometry and calculate resampling + // spacing therefrom. + widthInMM = m_WorldGeometry->GetExtentInMM(0); + heightInMM = m_WorldGeometry->GetExtentInMM(1); + } + else + { + mitkThrow()<<"mitk::ExtractSliceFilter: No fitting geometry for reslice axis!"; + return; + } + } + + right.Normalize(); + bottom.Normalize(); + + m_OutPutSpacing[0] = widthInMM / extent[0]; + m_OutPutSpacing[1] = heightInMM / extent[1]; + + /*========== BEGIN setup extent of the slice ==========*/ + int xMin, xMax, yMin, yMax; + + xMin = yMin = 0; + xMax = static_cast(extent[0]); + yMax = static_cast(extent[1]); + + if (m_WorldGeometry->GetReferenceGeometry()) + { + double sliceBounds[6]; + for (auto &sliceBound : sliceBounds) + { + sliceBound = 0.0; + } + + if (this->GetClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds)) + { + // Calculate output extent (integer values) + xMin = static_cast(sliceBounds[0] / m_OutPutSpacing[0] + 0.5); + xMax = static_cast(sliceBounds[1] / m_OutPutSpacing[0] + 0.5); + yMin = static_cast(sliceBounds[2] / m_OutPutSpacing[1] + 0.5); + yMax = static_cast(sliceBounds[3] / m_OutPutSpacing[1] + 0.5); + } // ELSE we use the default values + } + + + sliceOrigin += right * (m_OutPutSpacing[0] * 0.5); + sliceOrigin += bottom * (m_OutPutSpacing[1] * 0.5); + + // a worldGeometry is no imageGeometry, thus it is manually set to true + sliceGeometry->ImageGeometryOn(); + + /*At this point we have to adjust the geometry because the origin isn't correct. + The wrong origin is related to the rotation of the current world geometry plane. + This causes errors on transferring world to index coordinates. We just shift the + origin in each direction about the amount of the expanding (needed while rotating + the plane). + */ + Vector3D axis0 = sliceGeometry->GetAxisVector(0); + Vector3D axis1 = sliceGeometry->GetAxisVector(1); + axis0.Normalize(); + axis1.Normalize(); + + // adapt the origin. Note that for orthogonal planes the minima are '0' and thus the origin stays the same. + sliceOrigin += (axis0 * (xMin * m_OutPutSpacing[0])) + (axis1 * (yMin * m_OutPutSpacing[1])); + + sliceGeometry->SetOrigin(sliceOrigin); + + /*the bounds as well as the extent of the worldGeometry are not adapted correctly during crosshair rotation. + This is only a quick fix and has to be evaluated. + The new bounds are set via the max values of the calculated slice extent. It will look like [ 0, x, 0, y, 0, 1]. + */ + mitk::BoundingBox::BoundsArrayType boundsCopy; + boundsCopy[0] = boundsCopy[2] = boundsCopy[4] = 0; + boundsCopy[5] = 1; + boundsCopy[1] = xMax - xMin; + boundsCopy[3] = yMax - yMin; + sliceGeometry->SetBounds(boundsCopy); + + sliceGeometry->Modified(); + + Image::Pointer output = this->GetOutput(); + output->Initialize(input->GetPixelType(), 2, *sliceGeometry); + + m_XMin = xMin; + m_XMax = xMax; + m_YMin = yMin; + m_YMax = yMax; + m_Right = right; + m_Bottom = bottom; } void mitk::ExtractSliceFilter::GenerateInputRequestedRegion() { // As we want all pixel information fo the image in our plane, the requested region // is set to the largest possible region in the image. // This is needed because an oblique plane has a larger extent then the image // and the in pipeline it is checked via PropagateResquestedRegion(). But the // extent of the slice is actually fitting because it is oblique within the image. ImageToImageFilter::InputImagePointer input = this->GetInput(); input->SetRequestedRegionToLargestPossibleRegion(); } mitk::ScalarType *mitk::ExtractSliceFilter::GetOutputSpacing() { return m_OutPutSpacing; } void mitk::ExtractSliceFilter::GenerateData() { mitk::Image *input = this->GetInput(); if (!input) { MITK_ERROR << "mitk::ExtractSliceFilter: No input image available. Please set the input!" << std::endl; itkExceptionMacro("mitk::ExtractSliceFilter: No input image available. Please set the input!"); return; } if (!m_WorldGeometry) { MITK_ERROR << "mitk::ExtractSliceFilter: No Geometry for reslicing available." << std::endl; itkExceptionMacro("mitk::ExtractSliceFilter: No Geometry for reslicing available."); return; } const TimeGeometry *inputTimeGeometry = this->GetInput()->GetTimeGeometry(); if ((inputTimeGeometry == nullptr) || (inputTimeGeometry->CountTimeSteps() <= 0)) { itkWarningMacro(<< "Error reading input image TimeGeometry."); return; } // is it a valid timeStep? if (inputTimeGeometry->IsValidTimeStep(m_TimeStep) == false) { itkWarningMacro(<< "This is not a valid timestep: " << m_TimeStep); return; } // check if there is something to display. if (!input->IsVolumeSet(m_TimeStep)) { itkWarningMacro(<< "No volume data existent at given timestep " << m_TimeStep); return; } /*================#BEGIN setup vtkImageReslice properties================*/ Point3D origin; - Vector3D right, bottom, normal; - double widthInMM, heightInMM; - Vector2D extent; + Vector3D normal; const auto *planeGeometry = dynamic_cast(m_WorldGeometry); // Code for curved planes, mostly taken 1:1 from imageVtkMapper2D and not tested yet. // Do we have an AbstractTransformGeometry? // This is the case for AbstractTransformGeometry's (e.g. a ThinPlateSplineCurvedGeometry ) const auto *abstractGeometry = dynamic_cast(m_WorldGeometry); if (abstractGeometry != nullptr) { m_ResliceTransform = abstractGeometry; - extent[0] = abstractGeometry->GetParametricExtent(0); - extent[1] = abstractGeometry->GetParametricExtent(1); - - widthInMM = abstractGeometry->GetParametricExtentInMM(0); - heightInMM = abstractGeometry->GetParametricExtentInMM(1); - - m_OutPutSpacing[0] = widthInMM / extent[0]; - m_OutPutSpacing[1] = heightInMM / extent[1]; - origin = abstractGeometry->GetPlane()->GetOrigin(); - right = abstractGeometry->GetPlane()->GetAxisVector(0); - right.Normalize(); - - bottom = abstractGeometry->GetPlane()->GetAxisVector(1); - bottom.Normalize(); - normal = abstractGeometry->GetPlane()->GetNormal(); normal.Normalize(); // Use a combination of the InputGeometry *and* the possible non-rigid // AbstractTransformGeometry for reslicing the 3D Image vtkSmartPointer composedResliceTransform = vtkSmartPointer::New(); composedResliceTransform->Identity(); composedResliceTransform->Concatenate( inputTimeGeometry->GetGeometryForTimeStep(m_TimeStep)->GetVtkTransform()->GetLinearInverse()); composedResliceTransform->Concatenate(abstractGeometry->GetVtkAbstractTransform()); m_Reslicer->SetResliceTransform(composedResliceTransform); // Set background level to BLACK instead of translucent, to avoid // boundary artifacts (see PlaneGeometryDataVtkMapper3D) // Note: Backgroundlevel was hardcoded before to -1023 m_Reslicer->SetBackgroundLevel(m_BackgroundLevel); } else { if (planeGeometry != nullptr) { // if the worldGeomatry is a PlaneGeometry everything is straight forward origin = planeGeometry->GetOrigin(); - right = planeGeometry->GetAxisVector(0); - bottom = planeGeometry->GetAxisVector(1); normal = planeGeometry->GetNormal(); - - if (m_InPlaneResampleExtentByGeometry) - { - // Resampling grid corresponds to the current world geometry. This - // means that the spacing of the output 2D image depends on the - // currently selected world geometry, and *not* on the image itself. - extent[0] = m_WorldGeometry->GetExtent(0); - extent[1] = m_WorldGeometry->GetExtent(1); - } - else - { - // Resampling grid corresponds to the input geometry. This means that - // the spacing of the output 2D image is directly derived from the - // associated input image, regardless of the currently selected world - // geometry. - Vector3D rightInIndex, bottomInIndex; - inputTimeGeometry->GetGeometryForTimeStep(m_TimeStep)->WorldToIndex(right, rightInIndex); - inputTimeGeometry->GetGeometryForTimeStep(m_TimeStep)->WorldToIndex(bottom, bottomInIndex); - extent[0] = rightInIndex.GetNorm(); - extent[1] = bottomInIndex.GetNorm(); - } - - // Get the extent of the current world geometry and calculate resampling - // spacing therefrom. - widthInMM = m_WorldGeometry->GetExtentInMM(0); - heightInMM = m_WorldGeometry->GetExtentInMM(1); - - m_OutPutSpacing[0] = widthInMM / extent[0]; - m_OutPutSpacing[1] = heightInMM / extent[1]; - - right.Normalize(); - bottom.Normalize(); normal.Normalize(); /* * Transform the origin to center based coordinates. * Note: * This is needed besause vtk's origin is center based too (!!!) ( see 'The VTK book' page 88 ) * and the worldGeometry surrouding the image is no imageGeometry. So the worldGeometry * has its origin at the corner of the voxel and needs to be transformed. */ - origin += right * (m_OutPutSpacing[0] * 0.5); - origin += bottom * (m_OutPutSpacing[1] * 0.5); + origin += m_Right * (m_OutPutSpacing[0] * 0.5); + origin += m_Bottom * (m_OutPutSpacing[1] * 0.5); // set the tranform for reslicing. // Use inverse transform of the input geometry for reslicing the 3D image. // This is needed if the image volume already transformed if (m_ResliceTransform.IsNotNull()) m_Reslicer->SetResliceTransform(m_ResliceTransform->GetVtkTransform()->GetLinearInverse()); // Set background level to TRANSLUCENT (see PlaneGeometryDataVtkMapper3D), // else the background of the image turns out gray // Note: Backgroundlevel was hardcoded to -32768 m_Reslicer->SetBackgroundLevel(m_BackgroundLevel); } else { itkExceptionMacro("mitk::ExtractSliceFilter: No fitting geometry for reslice axis!"); return; } } if (m_ResliceTransform.IsNotNull()) { // if the resliceTransform is set the reslice axis are recalculated. // Thus the geometry information is not fitting. Therefor a unitSpacingFilter // is used to set up a global spacing of 1 and compensate the transform. vtkSmartPointer unitSpacingImageFilter = vtkSmartPointer::New(); unitSpacingImageFilter->ReleaseDataFlagOn(); unitSpacingImageFilter->SetOutputSpacing(1.0, 1.0, 1.0); unitSpacingImageFilter->SetInputData(input->GetVtkImageData(m_TimeStep)); m_Reslicer->SetInputConnection(unitSpacingImageFilter->GetOutputPort()); } else { // if no transform is set the image can be used directly m_Reslicer->SetInputData(input->GetVtkImageData(m_TimeStep)); } /*setup the plane where vktImageReslice extracts the slice*/ // ResliceAxesOrigin is the anchor point of the plane double originInVtk[3]; itk2vtk(origin, originInVtk); m_Reslicer->SetResliceAxesOrigin(originInVtk); // the cosines define the plane: x and y are the direction vectors, n is the planes normal // this specifies a matrix 3x3 // x1 y1 n1 // x2 y2 n2 // x3 y3 n3 double cosines[9]; - vnl2vtk(right.GetVnlVector(), cosines); // x + vnl2vtk(m_Right.GetVnlVector(), cosines); // x - vnl2vtk(bottom.GetVnlVector(), cosines + 3); // y + vnl2vtk(m_Bottom.GetVnlVector(), cosines + 3); // y vnl2vtk(normal.GetVnlVector(), cosines + 6); // n m_Reslicer->SetResliceAxesDirectionCosines(cosines); // we only have one slice, not a volume m_Reslicer->SetOutputDimensionality(m_OutputDimension); // set the interpolation mode for slicing switch (this->m_InterpolationMode) { case RESLICE_NEAREST: m_Reslicer->SetInterpolationModeToNearestNeighbor(); break; case RESLICE_LINEAR: m_Reslicer->SetInterpolationModeToLinear(); break; case RESLICE_CUBIC: m_Reslicer->SetInterpolationModeToCubic(); break; default: // the default interpolation used by mitk m_Reslicer->SetInterpolationModeToNearestNeighbor(); } /*========== BEGIN setup extent of the slice ==========*/ - int xMin, xMax, yMin, yMax; - - xMin = yMin = 0; - xMax = static_cast(extent[0]); - yMax = static_cast(extent[1]); - - if (m_WorldGeometry->GetReferenceGeometry()) - { - double sliceBounds[6]; - for (auto &sliceBound : sliceBounds) - { - sliceBound = 0.0; - } - - if (this->GetClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds)) - { - // Calculate output extent (integer values) - xMin = static_cast(sliceBounds[0] / m_OutPutSpacing[0] + 0.5); - xMax = static_cast(sliceBounds[1] / m_OutPutSpacing[0] + 0.5); - yMin = static_cast(sliceBounds[2] / m_OutPutSpacing[1] + 0.5); - yMax = static_cast(sliceBounds[3] / m_OutPutSpacing[1] + 0.5); - } // ELSE we use the default values - } - // Set the output extents! First included pixel index and last included pixel index // xMax and yMax are one after the last pixel. so they have to be decremented by 1. // In case we have a 2D image, xMax or yMax might be 0. in this case, do not decrement, but take 0. - m_Reslicer->SetOutputExtent(xMin, std::max(0, xMax - 1), yMin, std::max(0, yMax - 1), m_ZMin, m_ZMax); + m_Reslicer->SetOutputExtent(m_XMin, std::max(0, m_XMax - 1), m_YMin, std::max(0, m_YMax - 1), m_ZMin, m_ZMax); /*========== END setup extent of the slice ==========*/ m_Reslicer->SetOutputOrigin(0.0, 0.0, 0.0); m_Reslicer->SetOutputSpacing(m_OutPutSpacing[0], m_OutPutSpacing[1], m_ZSpacing); // TODO check the following lines, they are responsible whether vtk error outputs appear or not m_Reslicer->UpdateWholeExtent(); // this produces a bad allocation error for 2D images // m_Reslicer->GetOutput()->UpdateInformation(); // m_Reslicer->GetOutput()->SetUpdateExtentToWholeExtent(); // start the pipeline m_Reslicer->Update(); /*================ #END setup vtkImageReslice properties================*/ if (m_VtkOutputRequested) { // no conversion to mitk // no mitk geometry will be set, as the output is vtkImageData only!!! // no image component will be extracted, as the caller might need the whole multi-component image as vtk output return; } else { auto reslicedImage = vtkSmartPointer::New(); reslicedImage = m_Reslicer->GetOutput(); if (nullptr == reslicedImage) { itkWarningMacro(<< "Reslicer returned empty image"); return; } /*================ #BEGIN Extract component from image slice ================*/ int numberOfScalarComponent = reslicedImage->GetNumberOfScalarComponents(); if (numberOfScalarComponent > 1 && static_cast(numberOfScalarComponent) >= m_Component) { // image has more than one component, extract the correct component information with the given 'component' parameter auto vectorComponentExtractor = vtkSmartPointer::New(); vectorComponentExtractor->SetInputData(reslicedImage); vectorComponentExtractor->SetComponents(m_Component); vectorComponentExtractor->Update(); reslicedImage = vectorComponentExtractor->GetOutput(); } /*================ #END Extract component from image slice ================*/ /*================ #BEGIN Convert the slice to an mitk::Image ================*/ mitk::Image::Pointer resultImage = GetOutput(); + /*Temporary store the geometry that is already correct (set in GeneratOutputInformation()) + but will be reset due to initialize.*/ + mitk::BaseGeometry::Pointer resultGeometry = resultImage->GetGeometry(); + // initialize resultimage with the specs of the vtkImageData object returned from vtkImageReslice if (reslicedImage->GetDataDimension() == 1) { // If original image was 2D, the slice might have an y extent of 0. // Still i want to ensure here that Image is 2D resultImage->Initialize(reslicedImage, 1, -1, -1, 1); } else { resultImage->Initialize(reslicedImage); } - - // transfer the voxel data - resultImage->SetVolume(reslicedImage->GetScalarPointer()); - - // set the geometry from current worldgeometry for the reusultimage - // this is needed that the image has the correct mitk geometry - // the originalGeometry is the Geometry of the result slice - - // mitk::AffineGeometryFrame3D::Pointer originalGeometryAGF = m_WorldGeometry->Clone(); - // PlaneGeometry::Pointer originalGeometry = dynamic_cast( originalGeometryAGF.GetPointer() ); - PlaneGeometry::Pointer originalGeometry = m_WorldGeometry->Clone(); - - originalGeometry->GetIndexToWorldTransform()->SetMatrix(m_WorldGeometry->GetIndexToWorldTransform()->GetMatrix()); - - // the origin of the worldGeometry is transformed to center based coordinates to be an imageGeometry - Point3D sliceOrigin = originalGeometry->GetOrigin(); - - sliceOrigin += right * (m_OutPutSpacing[0] * 0.5); - sliceOrigin += bottom * (m_OutPutSpacing[1] * 0.5); - - // a worldGeometry is no imageGeometry, thus it is manually set to true - originalGeometry->ImageGeometryOn(); - - /*At this point we have to adjust the geometry because the origin isn't correct. - The wrong origin is related to the rotation of the current world geometry plane. - This causes errors on transferring world to index coordinates. We just shift the - origin in each direction about the amount of the expanding (needed while rotating - the plane). - */ - Vector3D axis0 = originalGeometry->GetAxisVector(0); - Vector3D axis1 = originalGeometry->GetAxisVector(1); - axis0.Normalize(); - axis1.Normalize(); - - // adapt the origin. Note that for orthogonal planes the minima are '0' and thus the origin stays the same. - sliceOrigin += (axis0 * (xMin * m_OutPutSpacing[0])) + (axis1 * (yMin * m_OutPutSpacing[1])); - - originalGeometry->SetOrigin(sliceOrigin); - - originalGeometry->Modified(); - - resultImage->SetGeometry(originalGeometry); - - /*the bounds as well as the extent of the worldGeometry are not adapted correctly during crosshair rotation. - This is only a quick fix and has to be evaluated. - The new bounds are set via the max values of the calculated slice extent. It will look like [ 0, x, 0, y, 0, 1]. - */ - mitk::BoundingBox::BoundsArrayType boundsCopy; - boundsCopy[0] = boundsCopy[2] = boundsCopy[4] = 0; - boundsCopy[5] = 1; - boundsCopy[1] = xMax - xMin; - boundsCopy[3] = yMax - yMin; - resultImage->GetGeometry()->SetBounds(boundsCopy); - /*================ #END Convert the slice to an mitk::Image ================*/ + + resultImage->SetGeometry(resultGeometry); } } bool mitk::ExtractSliceFilter::GetClippedPlaneBounds(double bounds[6]) { if (!m_WorldGeometry || !this->GetInput()) return false; return this->GetClippedPlaneBounds( m_WorldGeometry->GetReferenceGeometry(), dynamic_cast(m_WorldGeometry), bounds); } bool mitk::ExtractSliceFilter::GetClippedPlaneBounds(const BaseGeometry *boundingGeometry, const PlaneGeometry *planeGeometry, double *bounds) { bool b = mitk::PlaneClipping::CalculateClippedPlaneBounds(boundingGeometry, planeGeometry, bounds); return b; } diff --git a/Modules/Core/src/DataManagement/mitkTimeGeometry.cpp b/Modules/Core/src/DataManagement/mitkTimeGeometry.cpp index 48edadedec..f253f8c02e 100644 --- a/Modules/Core/src/DataManagement/mitkTimeGeometry.cpp +++ b/Modules/Core/src/DataManagement/mitkTimeGeometry.cpp @@ -1,296 +1,296 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include mitk::TimeGeometry::TimeGeometry() : m_BoundingBox(BoundingBox::New()) { typedef BoundingBox::PointsContainer ContainerType; ContainerType::Pointer points = ContainerType::New(); m_BoundingBox->SetPoints(points.GetPointer()); } mitk::TimeGeometry::~TimeGeometry() = default; void mitk::TimeGeometry::Initialize() { } /* \brief short description * parameters * */ mitk::Point3D mitk::TimeGeometry::GetCornerPointInWorld(int id) const { assert(id >= 0); assert(m_BoundingBox.IsNotNull()); - BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); + auto& bounds = m_BoundingBox->GetBounds(); Point3D cornerpoint; switch (id) { case 0: FillVector3D(cornerpoint, bounds[0], bounds[2], bounds[4]); break; case 1: FillVector3D(cornerpoint, bounds[0], bounds[2], bounds[5]); break; case 2: FillVector3D(cornerpoint, bounds[0], bounds[3], bounds[4]); break; case 3: FillVector3D(cornerpoint, bounds[0], bounds[3], bounds[5]); break; case 4: FillVector3D(cornerpoint, bounds[1], bounds[2], bounds[4]); break; case 5: FillVector3D(cornerpoint, bounds[1], bounds[2], bounds[5]); break; case 6: FillVector3D(cornerpoint, bounds[1], bounds[3], bounds[4]); break; case 7: FillVector3D(cornerpoint, bounds[1], bounds[3], bounds[5]); break; default: { itkExceptionMacro(<< "A cube only has 8 corners. These are labeled 0-7."); return Point3D(); } } // TimeGeometry has no Transformation. Therefore the bounding box // contains all data in world coordinates return cornerpoint; } mitk::Point3D mitk::TimeGeometry::GetCornerPointInWorld(bool xFront, bool yFront, bool zFront) const { assert(m_BoundingBox.IsNotNull()); - BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); + auto& bounds = m_BoundingBox->GetBounds(); Point3D cornerpoint; cornerpoint[0] = (xFront ? bounds[0] : bounds[1]); cornerpoint[1] = (yFront ? bounds[2] : bounds[3]); cornerpoint[2] = (zFront ? bounds[4] : bounds[5]); return cornerpoint; } mitk::Point3D mitk::TimeGeometry::GetCenterInWorld() const { assert(m_BoundingBox.IsNotNull()); return m_BoundingBox->GetCenter(); } double mitk::TimeGeometry::GetDiagonalLength2InWorld() const { const Vector3D diagonalvector = GetCornerPointInWorld()-GetCornerPointInWorld(false, false, false); return diagonalvector.GetSquaredNorm(); } double mitk::TimeGeometry::GetDiagonalLengthInWorld() const { return sqrt(GetDiagonalLength2InWorld()); } bool mitk::TimeGeometry::IsWorldPointInside(const mitk::Point3D &p) const { return m_BoundingBox->IsInside(p); } void mitk::TimeGeometry::UpdateBoundingBox() { assert(m_BoundingBox.IsNotNull()); typedef BoundingBox::PointsContainer ContainerType; unsigned long lastModifiedTime = 0; unsigned long currentModifiedTime = 0; ContainerType::Pointer points = ContainerType::New(); const TimeStepType numberOfTimesteps = CountTimeSteps(); points->reserve(2*numberOfTimesteps); for (TimeStepType step = 0; step GetMTime(); if (currentModifiedTime > lastModifiedTime) lastModifiedTime = currentModifiedTime; for (int i = 0; i < 8; ++i) { Point3D cornerPoint = GetGeometryForTimeStep(step)->GetCornerPoint(i); points->push_back(cornerPoint); } } m_BoundingBox->SetPoints(points); m_BoundingBox->ComputeBoundingBox(); if (this->GetMTime() < lastModifiedTime) this->Modified(); } mitk::ScalarType mitk::TimeGeometry::GetExtentInWorld(unsigned int direction) const { assert(direction < 3); assert(m_BoundingBox.IsNotNull()); - BoundingBox::BoundsArrayType bounds = m_BoundingBox->GetBounds(); + auto& bounds = m_BoundingBox->GetBounds(); return bounds[direction * 2 + 1] - bounds[direction * 2]; } void mitk::TimeGeometry::Update() { this->UpdateBoundingBox(); this->UpdateWithoutBoundingBox(); } void mitk::TimeGeometry::ExecuteOperation(mitk::Operation *op) { for (TimeStepType step = 0; step < CountTimeSteps(); ++step) { GetGeometryForTimeStep(step)->ExecuteOperation(op); } } void mitk::TimeGeometry::PrintSelf(std::ostream &os, itk::Indent indent) const { // Superclass::PrintSelf(os,indent); os << indent << " TimeSteps: " << this->CountTimeSteps() << std::endl; os << std::endl; os << indent << " GetGeometryForTimeStep(0): "; if (GetGeometryForTimeStep(0).IsNull()) os << "nullptr" << std::endl; else GetGeometryForTimeStep(0)->Print(os, indent); } itk::LightObject::Pointer mitk::TimeGeometry::InternalClone() const { itk::LightObject::Pointer parent = Superclass::InternalClone(); Self::Pointer rval = dynamic_cast(parent.GetPointer()); if (rval.IsNull()) { mitkThrow() << " Downcast to type " << this->GetNameOfClass() << " failed."; } rval->m_BoundingBox = m_BoundingBox->DeepCopy(); return parent; } bool mitk::Equal(const TimeGeometry &leftHandSide, const TimeGeometry &rightHandSide, ScalarType eps, bool verbose) { bool result = true; // Compare BoundingBoxInWorld if (!mitk::Equal(*(leftHandSide.GetBoundingBoxInWorld()), *(rightHandSide.GetBoundingBoxInWorld()), eps, verbose)) { if (verbose) { MITK_INFO << "[( TimeGeometry )] BoundingBoxInWorld differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetBoundingBoxInWorld() << " : leftHandSide is " << leftHandSide.GetBoundsInWorld() << " and tolerance is " << eps; } result = false; } if (!mitk::Equal(leftHandSide.CountTimeSteps(), rightHandSide.CountTimeSteps(), eps, verbose)) { if (verbose) { MITK_INFO << "[( TimeGeometry )] CountTimeSteps differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.CountTimeSteps() << " : leftHandSide is " << leftHandSide.CountTimeSteps() << " and tolerance is " << eps; } result = false; } if (!mitk::Equal(leftHandSide.GetMinimumTimePoint(), rightHandSide.GetMinimumTimePoint(), eps, verbose)) { if (verbose) { MITK_INFO << "[( TimeGeometry )] MinimumTimePoint differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetMinimumTimePoint() << " : leftHandSide is " << leftHandSide.GetMinimumTimePoint() << " and tolerance is " << eps; } result = false; } if (!mitk::Equal(leftHandSide.GetMaximumTimePoint(), rightHandSide.GetMaximumTimePoint(), eps, verbose)) { if (verbose) { MITK_INFO << "[( TimeGeometry )] MaximumTimePoint differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.GetMaximumTimePoint() << " : leftHandSide is " << leftHandSide.GetMaximumTimePoint() << " and tolerance is " << eps; } result = false; } if (!result) return false; // further tests require that both parts have identical number of time steps for (mitk::TimeStepType t = 0; t < leftHandSide.CountTimeSteps(); ++t) { if (!mitk::Equal(leftHandSide.TimeStepToTimePoint(t), rightHandSide.TimeStepToTimePoint(t), eps, verbose)) { if (verbose) { MITK_INFO << "[( TimeGeometry )] TimeStepToTimePoint(" << t << ") differs."; MITK_INFO << "rightHandSide is " << setprecision(12) << rightHandSide.TimeStepToTimePoint(t) << " : leftHandSide is " << leftHandSide.TimeStepToTimePoint(t) << " and tolerance is " << eps; } result = false; } BaseGeometry::Pointer leftGeometry = leftHandSide.GetGeometryForTimeStep(t); BaseGeometry::Pointer rightGeometry = rightHandSide.GetGeometryForTimeStep(t); if (leftGeometry.IsNotNull() && rightGeometry.IsNull()) continue; // identical if (leftGeometry.IsNull()) { if (verbose) { MITK_INFO << "[( TimeGeometry )] TimeStepToTimePoint(" << t << ") differs."; MITK_INFO << "rightHandSide is an object : leftHandSide is nullptr"; } result = false; continue; // next geometry } if (rightGeometry.IsNull()) { if (verbose) { MITK_INFO << "[( TimeGeometry )] TimeStepToTimePoint(" << t << ") differs."; MITK_INFO << "rightHandSide is nullptr : leftHandSide is an object"; } result = false; continue; // next geometry } if (!mitk::Equal(*leftGeometry, *rightGeometry, eps, verbose)) { if (verbose) { MITK_INFO << "[( TimeGeometry )] GetGeometryForTimeStep(" << t << ") differs."; } result = false; } } // end for each t return result; }