diff --git a/Core/Code/Algorithms/mitkExtractSliceFilter.cpp b/Core/Code/Algorithms/mitkExtractSliceFilter.cpp index 6ad58c9103..4d33011ff4 100644 --- a/Core/Code/Algorithms/mitkExtractSliceFilter.cpp +++ b/Core/Code/Algorithms/mitkExtractSliceFilter.cpp @@ -1,587 +1,601 @@ /*=================================================================== 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 "mitkExtractSliceFilter.h" #include #include #include #include #include #include mitk::ExtractSliceFilter::ExtractSliceFilter(vtkImageReslice* reslicer ){ if(reslicer == NULL){ m_Reslicer = vtkSmartPointer::New(); } else { m_Reslicer = reslicer; } m_TimeStep = 0; m_Reslicer->ReleaseDataFlagOn(); m_InterpolationMode = ExtractSliceFilter::RESLICE_NEAREST; m_ResliceTransform = NULL; 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; } mitk::ExtractSliceFilter::~ExtractSliceFilter(){ m_ResliceTransform = NULL; m_WorldGeometry = NULL; delete [] m_OutPutSpacing; } void mitk::ExtractSliceFilter::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);*/ } 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 = const_cast< ImageToImageFilter::InputImageType* > ( this->GetInput() ); input->SetRequestedRegionToLargestPossibleRegion(); } mitk::ScalarType* mitk::ExtractSliceFilter::GetOutputSpacing(){ return m_OutPutSpacing; } void mitk::ExtractSliceFilter::GenerateData(){ mitk::Image *input = const_cast< mitk::Image * >( 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 TimeSlicedGeometry *inputTimeGeometry = this->GetInput()->GetTimeSlicedGeometry(); if ( ( inputTimeGeometry == NULL ) || ( inputTimeGeometry->GetTimeSteps() == 0 ) ) { itkWarningMacro(<<"Error reading input image TimeSlicedGeometry."); return; } // is it a valid timeStep? if ( inputTimeGeometry->IsValidTime( 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 vtkImageRslice properties================*/ Point3D origin; Vector3D right, bottom, normal; double widthInMM, heightInMM; Vector2D extent; const PlaneGeometry* planeGeometry = dynamic_cast(m_WorldGeometry); if ( planeGeometry != NULL ) { //if the worldGeomatry is a PlaneGeometry everthing 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->GetGeometry3D( m_TimeStep )->WorldToIndex( right, rightInIndex ); inputTimeGeometry->GetGeometry3D( 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 ); //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 Geometry2DDataVtkMapper3D), // else the background of the image turns out gray m_Reslicer->SetBackgroundLevel( -32768 ); } else{ //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 mitk::AbstractTransformGeometry* abstractGeometry = dynamic_cast< const AbstractTransformGeometry * >(m_WorldGeometry); if(abstractGeometry != NULL) { 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 = vtkGeneralTransform::New(); composedResliceTransform->Identity(); composedResliceTransform->Concatenate( inputTimeGeometry->GetGeometry3D( m_TimeStep )->GetVtkTransform()->GetLinearInverse() ); composedResliceTransform->Concatenate( abstractGeometry->GetVtkAbstractTransform() ); m_Reslicer->SetResliceTransform( composedResliceTransform ); composedResliceTransform->UnRegister( NULL ); // decrease RC // Set background level to BLACK instead of translucent, to avoid // boundary artifacts (see Geometry2DDataVtkMapper3D) m_Reslicer->SetBackgroundLevel( -1023 ); } 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->SetInput( input->GetVtkImageData(m_TimeStep) ); m_Reslicer->SetInput(unitSpacingImageFilter->GetOutput() ); } else { //if no tranform is set the image can be used directly m_Reslicer->SetInput(input->GetVtkImageData(m_TimeStep)); } /*setup the plane where vktImageReslice extracts the slice*/ //ResliceAxesOrigin is the ancor 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(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(); 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< int >( extent[0]); yMax = static_cast< int >( extent[1]); vtkFloatingPointType sliceBounds[6]; if (m_WorldGeometry->GetReferenceGeometry()) { for ( int i = 0; i < 6; ++i ) { sliceBounds[i] = 0.0; } if (this->GetClippedPlaneBounds( m_WorldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds )) { // Calculate output extent (integer values) xMin = static_cast< int >( sliceBounds[0] / m_OutPutSpacing[0] + 0.5 ); xMax = static_cast< int >( sliceBounds[1] / m_OutPutSpacing[0] + 0.5 ); yMin = static_cast< int >( sliceBounds[2] / m_OutPutSpacing[1] + 0.5 ); yMax = static_cast< int >( sliceBounds[3] / m_OutPutSpacing[1] + 0.5 ); } // ELSE we use the default values } - m_Reslicer->SetOutputExtent(xMin, xMax-1, yMin, yMax-1, m_ZMin, m_ZMax ); + // 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 ); /*========== 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 wether 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 vtkImageRslice properties================*/ if(m_VtkOutputRequested){ return; //no converting to mitk //no mitk geometry will be set, as the output is vtkImageData only!!! } else { /*================ #BEGIN Get the slice from vtkImageReslice and convert it to mit::Image================*/ vtkImageData* reslicedImage; reslicedImage = m_Reslicer->GetOutput(); if(!reslicedImage) { itkWarningMacro(<<"Reslicer returned empty image"); return; } mitk::Image::Pointer resultImage = this->GetOutput(); //initialize resultimage with the specs of the vtkImageData object returned from vtkImageReslice - resultImage->Initialize(reslicedImage); + 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 AffineGeometryFrame3D::Pointer originalGeometryAGF = m_WorldGeometry->Clone(); Geometry2D::Pointer originalGeometry = dynamic_cast( originalGeometryAGF.GetPointer() ); 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 transfering 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 calcuted 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 Get the slice from vtkImageReslice and convert it to mitk Image================*/ } } bool mitk::ExtractSliceFilter::GetClippedPlaneBounds(vtkFloatingPointType bounds[6]){ if(!m_WorldGeometry || !this->GetInput()) return false; return this->GetClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), dynamic_cast< const PlaneGeometry * >( m_WorldGeometry ), bounds); } bool mitk::ExtractSliceFilter::GetClippedPlaneBounds( const Geometry3D *boundingGeometry, const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ) { bool b = this->CalculateClippedPlaneBounds(boundingGeometry, planeGeometry, bounds); return b; } bool mitk::ExtractSliceFilter ::CalculateClippedPlaneBounds( const Geometry3D *boundingGeometry, const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ) { // Clip the plane with the bounding geometry. To do so, the corner points // of the bounding box are transformed by the inverse transformation // matrix, and the transformed bounding box edges derived therefrom are // clipped with the plane z=0. The resulting min/max values are taken as // bounds for the image reslicer. const BoundingBox *boundingBox = boundingGeometry->GetBoundingBox(); BoundingBox::PointType bbMin = boundingBox->GetMinimum(); BoundingBox::PointType bbMax = boundingBox->GetMaximum(); vtkPoints *points = vtkPoints::New(); if(boundingGeometry->GetImageGeometry()) { points->InsertPoint( 0, bbMin[0]-0.5, bbMin[1]-0.5, bbMin[2]-0.5 ); points->InsertPoint( 1, bbMin[0]-0.5, bbMin[1]-0.5, bbMax[2]-0.5 ); points->InsertPoint( 2, bbMin[0]-0.5, bbMax[1]-0.5, bbMax[2]-0.5 ); points->InsertPoint( 3, bbMin[0]-0.5, bbMax[1]-0.5, bbMin[2]-0.5 ); points->InsertPoint( 4, bbMax[0]-0.5, bbMin[1]-0.5, bbMin[2]-0.5 ); points->InsertPoint( 5, bbMax[0]-0.5, bbMin[1]-0.5, bbMax[2]-0.5 ); points->InsertPoint( 6, bbMax[0]-0.5, bbMax[1]-0.5, bbMax[2]-0.5 ); points->InsertPoint( 7, bbMax[0]-0.5, bbMax[1]-0.5, bbMin[2]-0.5 ); } else { points->InsertPoint( 0, bbMin[0], bbMin[1], bbMin[2] ); points->InsertPoint( 1, bbMin[0], bbMin[1], bbMax[2] ); points->InsertPoint( 2, bbMin[0], bbMax[1], bbMax[2] ); points->InsertPoint( 3, bbMin[0], bbMax[1], bbMin[2] ); points->InsertPoint( 4, bbMax[0], bbMin[1], bbMin[2] ); points->InsertPoint( 5, bbMax[0], bbMin[1], bbMax[2] ); points->InsertPoint( 6, bbMax[0], bbMax[1], bbMax[2] ); points->InsertPoint( 7, bbMax[0], bbMax[1], bbMin[2] ); } vtkPoints *newPoints = vtkPoints::New(); vtkTransform *transform = vtkTransform::New(); transform->Identity(); transform->Concatenate( planeGeometry->GetVtkTransform()->GetLinearInverse() ); transform->Concatenate( boundingGeometry->GetVtkTransform() ); transform->TransformPoints( points, newPoints ); transform->Delete(); bounds[0] = bounds[2] = 10000000.0; bounds[1] = bounds[3] = -10000000.0; bounds[4] = bounds[5] = 0.0; this->LineIntersectZero( newPoints, 0, 1, bounds ); this->LineIntersectZero( newPoints, 1, 2, bounds ); this->LineIntersectZero( newPoints, 2, 3, bounds ); this->LineIntersectZero( newPoints, 3, 0, bounds ); this->LineIntersectZero( newPoints, 0, 4, bounds ); this->LineIntersectZero( newPoints, 1, 5, bounds ); this->LineIntersectZero( newPoints, 2, 6, bounds ); this->LineIntersectZero( newPoints, 3, 7, bounds ); this->LineIntersectZero( newPoints, 4, 5, bounds ); this->LineIntersectZero( newPoints, 5, 6, bounds ); this->LineIntersectZero( newPoints, 6, 7, bounds ); this->LineIntersectZero( newPoints, 7, 4, bounds ); // clean up vtk data points->Delete(); newPoints->Delete(); if ( (bounds[0] > 9999999.0) || (bounds[2] > 9999999.0) || (bounds[1] < -9999999.0) || (bounds[3] < -9999999.0) ) { return false; } else { // The resulting bounds must be adjusted by the plane spacing, since we // we have so far dealt with index coordinates const float *planeSpacing = planeGeometry->GetFloatSpacing(); bounds[0] *= planeSpacing[0]; bounds[1] *= planeSpacing[0]; bounds[2] *= planeSpacing[1]; bounds[3] *= planeSpacing[1]; bounds[4] *= planeSpacing[2]; bounds[5] *= planeSpacing[2]; return true; } } bool mitk::ExtractSliceFilter ::LineIntersectZero( vtkPoints *points, int p1, int p2, vtkFloatingPointType *bounds ) { vtkFloatingPointType point1[3]; vtkFloatingPointType point2[3]; points->GetPoint( p1, point1 ); points->GetPoint( p2, point2 ); if ( (point1[2] * point2[2] <= 0.0) && (point1[2] != point2[2]) ) { double x, y; x = ( point1[0] * point2[2] - point1[2] * point2[0] ) / ( point2[2] - point1[2] ); y = ( point1[1] * point2[2] - point1[2] * point2[1] ) / ( point2[2] - point1[2] ); if ( x < bounds[0] ) { bounds[0] = x; } if ( x > bounds[1] ) { bounds[1] = x; } if ( y < bounds[2] ) { bounds[2] = y; } if ( y > bounds[3] ) { bounds[3] = y; } bounds[4] = bounds[5] = 0.0; return true; } return false; } diff --git a/Core/Code/DataManagement/mitkImage.cpp b/Core/Code/DataManagement/mitkImage.cpp index c25cb7443b..d06fed5c33 100644 --- a/Core/Code/DataManagement/mitkImage.cpp +++ b/Core/Code/DataManagement/mitkImage.cpp @@ -1,1278 +1,1284 @@ /*=================================================================== 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 "mitkImage.h" #include "mitkImageStatisticsHolder.h" #include "mitkPixelTypeMultiplex.h" #include #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(NULL), m_ImageDescriptor(NULL), m_OffsetTable(NULL), m_CompleteData(NULL), m_ImageStatistics(NULL) { 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(NULL), m_ImageDescriptor(NULL), m_OffsetTable(NULL), m_CompleteData(NULL), m_ImageStatistics(NULL) { 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 this->SetGeometry(dynamic_cast(other.GetTimeSlicedGeometry()->Clone().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() { Clear(); m_ReferenceCountLock.Lock(); m_ReferenceCount = 3; m_ReferenceCountLock.Unlock(); m_ReferenceCountLock.Lock(); m_ReferenceCount = 0; m_ReferenceCountLock.Unlock(); if(m_OffsetTable != NULL) delete [] m_OffsetTable; if(m_ImageStatistics != NULL) 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 NULL; 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 == NULL ) return; if(ptype.GetBpe() != 24) { value = (double) (((T*) data)[ offset ]); } else { const unsigned int rgboffset = 3 * offset; double returnvalue = (((T*) data)[rgboffset ]); returnvalue += (((T*) data)[rgboffset + 1]); returnvalue += (((T*) data)[rgboffset + 2]); value = returnvalue; } } double mitk::Image::GetPixelValueByIndex(const mitk::Index3D &position, unsigned int timestep) { 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 = 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) { double value = 0.0; if (this->GetTimeSteps() < timestep) { timestep = this->GetTimeSteps(); } Index3D itkIndex; this->GetGeometry()->WorldToIndex(position, itkIndex); value = this->GetPixelValueByIndex( itkIndex, timestep); return value; } mitk::ImageVtkAccessor* mitk::Image::GetVtkImageData(int t, int n) { if(m_Initialized==false) { if(GetSource().IsNull()) return NULL; if(GetSource()->Updating()==false) GetSource()->UpdateOutputInformation(); } ImageDataItemPointer volume=GetVolumeData(t, n); if(volume.GetPointer()==NULL || volume->GetVtkImageData(this) == NULL) return NULL; float *fspacing = const_cast(GetSlicedGeometry(t)->GetFloatSpacing()); double dspacing[3] = {fspacing[0],fspacing[1],fspacing[2]}; volume->GetVtkImageData(this)->SetSpacing( dspacing ); return volume->GetVtkImageData(this); } mitk::Image::ImageDataItemPointer mitk::Image::GetSliceData(int s, int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) { if(IsValidSlice(s,t,n)==false) return NULL; const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); // slice directly available? int pos=GetSliceIndex(s,t,n); if(m_Slices[pos].GetPointer()!=NULL) 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()!=NULL) && (vol->IsComplete())) { sl=new ImageDataItem(*vol, m_ImageDescriptor, 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()!=NULL) && (ch->IsComplete())) { sl=new ImageDataItem(*ch, m_ImageDescriptor, 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(s,t,n)) //yes: now we can call ourselves without the risk of a endless loop (see "if" above) return GetSliceData(s,t,n,data,importMemoryManagement); else return NULL; } else { ImageDataItemPointer item = AllocateSliceData(s,t,n,data,importMemoryManagement); item->SetComplete(true); return item; } } mitk::Image::ImageDataItemPointer mitk::Image::GetVolumeData(int t, int n, void *data, ImportMemoryManagementType importMemoryManagement) { if(IsValidVolume(t,n)==false) return NULL; ImageDataItemPointer ch, vol; // volume directly available? int pos=GetVolumeIndex(t,n); vol=m_Volumes[pos]; if((vol.GetPointer()!=NULL) && (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()!=NULL) && (ch->IsComplete())) { vol=new ImageDataItem(*ch, m_ImageDescriptor, 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;sSetComplete(true); } else { mitk::PixelType chPixelType = this->m_ImageDescriptor->GetChannelTypeById(n); vol=m_Volumes[pos]; // ok, let's combine the slices! if(vol.GetPointer()==NULL) vol=new ImageDataItem( chPixelType, 3, m_Dimensions, NULL, true); vol->SetComplete(true); size_t size=m_OffsetTable[2]*(ptypeSize); for(s=0;sGetParent()!=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, 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(t,n)) //yes: now we can call ourselves without the risk of a endless loop (see "if" above) return GetVolumeData(t,n,data,importMemoryManagement); else return NULL; } else { ImageDataItemPointer item = AllocateVolumeData(t,n,data,importMemoryManagement); item->SetComplete(true); return item; } } mitk::Image::ImageDataItemPointer mitk::Image::GetChannelData(int n, void *data, ImportMemoryManagementType importMemoryManagement) { if(IsValidChannel(n)==false) return NULL; ImageDataItemPointer ch, vol; ch=m_Channels[n]; if((ch.GetPointer()!=NULL) && (ch->IsComplete())) return ch; // let's see if all volumes are set, so that we can (could) combine them to a channel if(IsChannelSet(n)) { // if there is only one time frame we do not need to combine anything if(m_Dimensions[3]<=1) { vol=GetVolumeData(0,n,data,importMemoryManagement); ch=new ImageDataItem(*vol, m_ImageDescriptor, 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()==NULL) ch=new ImageDataItem(this->m_ImageDescriptor, NULL, true); ch->SetComplete(true); size_t size=m_OffsetTable[m_Dimension-1]*(ptypeSize); unsigned int t; ImageDataItemPointerArray::iterator slicesIt = m_Slices.begin()+n*m_Dimensions[2]*m_Dimensions[3]; for(t=0;tGetParent()!=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, 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=NULL; 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(n)) //yes: now we can call ourselves without the risk of a endless loop (see "if" above) return GetChannelData(n,data,importMemoryManagement); else return NULL; } else { ImageDataItemPointer item = AllocateChannelData(n,data,importMemoryManagement); item->SetComplete(true); return item; } } bool mitk::Image::IsSliceSet(int s, int t, int n) const { if(IsValidSlice(s,t,n)==false) return false; if(m_Slices[GetSliceIndex(s,t,n)].GetPointer()!=NULL) return true; ImageDataItemPointer ch, vol; vol=m_Volumes[GetVolumeIndex(t,n)]; if((vol.GetPointer()!=NULL) && (vol->IsComplete())) return true; ch=m_Channels[n]; if((ch.GetPointer()!=NULL) && (ch->IsComplete())) return true; return false; } bool mitk::Image::IsVolumeSet(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()!=NULL) && (vol->IsComplete())) return true; // is volume available as part of a channel that is available? ch=m_Channels[n]; if((ch.GetPointer()!=NULL) && (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;sIsComplete())) 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(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()==NULL) 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()==NULL) 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()==NULL) 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()==NULL) 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()==NULL) 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()==NULL) 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)=NULL; } for( it=m_Volumes.begin(), end=m_Volumes.end(); it!=end; ++it ) { (*it)=NULL; } for( it=m_Channels.begin(), end=m_Channels.end(); it!=end; ++it ) { (*it)=NULL; } m_CompleteData = NULL; if( m_ImageStatistics == NULL) { 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;im_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 = NULL; return; } for( unsigned int i=0u; im_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]); if(dimension>=4) { TimeBounds timebounds; timebounds[0] = 0.0; timebounds[1] = 1.0; slicedGeometry->SetTimeBounds(timebounds); } TimeSlicedGeometry::Pointer timeSliceGeometry = TimeSlicedGeometry::New(); timeSliceGeometry->InitializeEvenlyTimed(slicedGeometry, m_Dimensions[3]); timeSliceGeometry->ImageGeometryOn(); SetGeometry(timeSliceGeometry); ImageDataItemPointer dnull=NULL; 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::Geometry3D& geometry, unsigned int channels, int tDim ) { unsigned int dimensions[5]; dimensions[0] = (unsigned int)(geometry.GetExtent(0)+0.5); dimensions[1] = (unsigned int)(geometry.GetExtent(1)+0.5); dimensions[2] = (unsigned int)(geometry.GetExtent(2)+0.5); dimensions[3] = 0; dimensions[4] = 0; unsigned int dimension = 2; if ( dimensions[2] > 1 ) dimension = 3; if ( tDim > 0) { dimensions[3] = tDim; } else { const mitk::TimeSlicedGeometry* timeGeometry = dynamic_cast(&geometry); if ( timeGeometry != NULL ) { dimensions[3] = timeGeometry->GetTimeSteps(); } } if ( dimensions[3] > 1 ) dimension = 4; Initialize( type, dimension, dimensions, channels ); SetGeometry(static_cast(geometry.Clone().GetPointer())); mitk::BoundingBox::BoundsArrayType bounds = geometry.GetBoundingBox()->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.Get_vnl_vector().data_block()); GetTimeSlicedGeometry()->InitializeEvenlyTimed(slicedGeometry, m_Dimensions[3]); } } void mitk::Image::Initialize(const mitk::PixelType& type, int sDim, const mitk::Geometry2D& geometry2d, bool flipped, unsigned int channels, int tDim ) { SlicedGeometry3D::Pointer slicedGeometry = SlicedGeometry3D::New(); slicedGeometry->InitializeEvenlySpaced(static_cast(geometry2d.Clone().GetPointer()), sDim, flipped); Initialize(type, *slicedGeometry, channels, tDim); } void mitk::Image::Initialize(const mitk::Image* image) { Initialize(image->GetPixelType(), *image->GetTimeSlicedGeometry()); } -void mitk::Image::Initialize(vtkImageData* vtkimagedata, int channels, int tDim, int sDim) +void mitk::Image::Initialize(vtkImageData* vtkimagedata, int channels, int tDim, int sDim, int pDim) { if(vtkimagedata==NULL) return; m_Dimension=vtkimagedata->GetDataDimension(); unsigned int i, *tmpDimensions=new unsigned int[m_Dimension>4?m_Dimension:4]; for(i=0;iGetDimensions()[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; } switch ( vtkimagedata->GetScalarType() ) { case VTK_BIT: case VTK_CHAR: //pixelType.Initialize(typeid(char), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_UNSIGNED_CHAR: //pixelType.Initialize(typeid(unsigned char), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_SHORT: //pixelType.Initialize(typeid(short), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_UNSIGNED_SHORT: //pixelType.Initialize(typeid(unsigned short), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_INT: //pixelType.Initialize(typeid(int), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_UNSIGNED_INT: //pixelType.Initialize(typeid(unsigned int), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_LONG: //pixelType.Initialize(typeid(long), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_UNSIGNED_LONG: //pixelType.Initialize(typeid(unsigned long), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_FLOAT: //pixelType.Initialize(typeid(float), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; case VTK_DOUBLE: //pixelType.Initialize(typeid(double), vtkimagedata->GetNumberOfScalarComponents()); Initialize(mitk::MakeScalarPixelType(), m_Dimension, tmpDimensions, channels); break; default: break; } /* 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; vtkFloatingPointType 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->GetGeometry2D(0)); planeGeometry->SetOrigin(origin); // re-initialize SlicedGeometry3D slicedGeometry->SetOrigin(origin); slicedGeometry->SetSpacing(spacing); GetTimeSlicedGeometry()->InitializeEvenlyTimed(slicedGeometry, m_Dimensions[3]); 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!=NULL) 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) { 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()!=NULL) { sl=new ImageDataItem(*vol, m_ImageDescriptor, 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()!=NULL) { sl=new ImageDataItem(*ch, m_ImageDescriptor, 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(t,n,NULL,importMemoryManagement); sl=new ImageDataItem(*vol, m_ImageDescriptor, 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) { 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()!=NULL) { vol=new ImageDataItem(*ch, m_ImageDescriptor, 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, 3, m_Dimensions, NULL, true); if(data != NULL) std::memcpy(vol->GetData(), data, m_OffsetTable[3]*(ptypeSize)); } else { vol=new ImageDataItem( chPixelType, 3, m_Dimensions, data, importMemoryManagement == ManageMemory); } m_Volumes[pos]=vol; return vol; } mitk::Image::ImageDataItemPointer mitk::Image::AllocateChannelData(int n, void *data, ImportMemoryManagementType importMemoryManagement) { ImageDataItemPointer ch; // allocate new channel if(importMemoryManagement == CopyMemory) { const size_t ptypeSize = this->m_ImageDescriptor->GetChannelTypeById(n).GetSize(); ch=new ImageDataItem(this->m_ImageDescriptor, NULL, true); if(data != NULL) std::memcpy(ch->GetData(), data, m_OffsetTable[4]*(ptypeSize)); } else { ch=new ImageDataItem(this->m_ImageDescriptor, 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 = NULL; } void mitk::Image::SetGeometry(Geometry3D* 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); GetTimeSlicedGeometry()->ImageGeometryOn(); } void mitk::Image::PrintSelf(std::ostream& os, itk::Indent indent) const { unsigned char i; if(m_Initialized) { 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.GetTypeId().name() << std::endl; os << indent << " BitsPerElement: " << chPixelType.GetSize() << 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::Geometry3D* geo = this->GetGeometry(); bool ret = false; if(geo) { const vnl_matrix_fixed & mx = geo->GetIndexToWorldTransform()->GetMatrix().GetVnlMatrix(); float 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; } #include "mitkImageStatisticsHolder.h" //##Documentation mitk::ScalarType mitk::Image::GetScalarValueMin(int t) const { return m_ImageStatistics->GetScalarValueMin(t); } //##Documentation //## \brief Get the maximum for scalar images mitk::ScalarType mitk::Image::GetScalarValueMax(int t) const { return m_ImageStatistics->GetScalarValueMax(t); } //##Documentation //## \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); } diff --git a/Core/Code/DataManagement/mitkImage.h b/Core/Code/DataManagement/mitkImage.h index a86ca25be5..6521128fa6 100644 --- a/Core/Code/DataManagement/mitkImage.h +++ b/Core/Code/DataManagement/mitkImage.h @@ -1,673 +1,674 @@ /*=================================================================== 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 #include "mitkSlicedData.h" #include "mitkBaseData.h" #include "mitkLevelWindow.h" #include "mitkPlaneGeometry.h" #include "mitkImageDataItem.h" #include "mitkImageDescriptor.h" #include "mitkImageAccessorBase.h" #include "mitkImageVtkAccessor.h" #ifndef __itkHistogram_h #include #endif class vtkImageData; 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 Geometry3D of the Image, see //## Geometry3D::WorldToItkPhysicalPoint. //## @ingroup Data class MITK_CORE_EXPORT Image : public SlicedData { friend class SubImageSelector; friend class ImageAccessorBase; friend class ImageVtkAccessor; friend class ImageReadAccessor; friend class ImageWriteAccessor; // friend class ImageWriteAccessor; public: mitkClassMacro(Image, SlicedData); itkNewMacro(Self); mitkCloneMacro(Image); /** 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 mitk::Index3D& position, unsigned int timestep = 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)); //##Documentation //## @brief Get a volume at a specific time @a t of channel @a n as a vtkImageData. virtual ImageVtkAccessor* GetVtkImageData(int t = 0, int n = 0); //##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; //##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; //##Documentation //## @brief Check whether the channel @a n is set virtual bool IsChannelSet(int n = 0) const; //##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 Geometry3D //## //## @param tDim override time dimension (@a n[3]) if @a geometry is a TimeSlicedGeometry (if >0) virtual void Initialize(const mitk::PixelType& type, const mitk::Geometry3D& geometry, unsigned int channels = 1, int tDim=-1); //##Documentation //## initialize new (or re-initialize) image information by a Geometry2D and number of slices //## //## Initializes the bounding box according to the width/height of the //## Geometry2D and @a sDim via SlicedGeometry3D::InitializeEvenlySpaced. //## The spacing is calculated from the Geometry2D. //## @param tDim override time dimension (@a n[3]) if @a geometry is a TimeSlicedGeometry (if >0) //## \sa SlicedGeometry3D::InitializeEvenlySpaced virtual void Initialize(const mitk::PixelType& type, int sDim, const mitk::Geometry2D& 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 Geometry2D 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 <) - virtual void Initialize(vtkImageData* vtkimagedata, int channels = 1, int tDim = -1, int sDim = -1); + //## @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==NULL) 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;iGetLargestPossibleRegion().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(), 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)->GetGeometry2D(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 TimeSlicedGeometry GetTimeSlicedGeometry()->InitializeEvenlyTimed(slicedGeometry, m_Dimensions[3]); // 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(Geometry3D* aGeometry3D); /** * @warning for internal use only */ virtual ImageDataItemPointer GetSliceData(int s = 0, int t = 0, int n = 0, void *data = NULL, ImportMemoryManagementType importMemoryManagement = CopyMemory); /** * @warning for internal use only */ virtual ImageDataItemPointer GetVolumeData(int t = 0, int n = 0, void *data = NULL, ImportMemoryManagementType importMemoryManagement = CopyMemory); /** * @warning for internal use only */ virtual ImageDataItemPointer GetChannelData(int n = 0, void *data = NULL, ImportMemoryManagementType importMemoryManagement = CopyMemory); /** \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: 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 ); virtual ImageDataItemPointer AllocateSliceData(int s = 0, int t = 0, int n = 0, void *data = NULL, ImportMemoryManagementType importMemoryManagement = CopyMemory); virtual ImageDataItemPointer AllocateVolumeData(int t = 0, int n = 0, void *data = NULL, ImportMemoryManagementType importMemoryManagement = CopyMemory); virtual ImageDataItemPointer AllocateChannelData(int n = 0, void *data = NULL, ImportMemoryManagementType importMemoryManagement = CopyMemory); Image(); Image(const Image &other); virtual ~Image(); virtual void Clear(); //## @warning Has to be called by every Initialize method! virtual void Initialize(); virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; mutable ImageDataItemPointerArray m_Channels; mutable ImageDataItemPointerArray m_Volumes; mutable ImageDataItemPointerArray m_Slices; 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: /** Stores all existing ImageReadAccessors */ std::vector m_Readers; /** Stores all existing ImageWriteAccessors */ std::vector m_Writers; /** Stores all existing ImageVtkAccessors */ 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; }; //##Documentation //## @brief Cast an itk::Image (with a specific type) to an mitk::Image. //## //## CastToMitkImage does not cast pixel types etc., just image data //## Needs "mitkImage.h" header included. //## If you get a compile error, try image.GetPointer(); //## @ingroup Adaptor //## \sa mitkITKImageImport template void CastToMitkImage(const itk::SmartPointer& itkimage, itk::SmartPointer& mitkoutputimage) { if(mitkoutputimage.IsNull()) { mitkoutputimage = mitk::Image::New(); } mitkoutputimage->InitializeByItk(itkimage.GetPointer()); mitkoutputimage->SetChannel(itkimage->GetBufferPointer()); } //##Documentation //## @brief Cast an itk::Image (with a specific type) to an mitk::Image. //## //## CastToMitkImage does not cast pixel types etc., just image data //## Needs "mitkImage.h" header included. //## If you get a compile error, try image.GetPointer(); //## @ingroup Adaptor //## \sa mitkITKImageImport template void CastToMitkImage(const ItkOutputImageType* itkimage, itk::SmartPointer& mitkoutputimage) { if(mitkoutputimage.IsNull()) { mitkoutputimage = mitk::Image::New(); } mitkoutputimage->InitializeByItk(itkimage); mitkoutputimage->SetChannel(itkimage->GetBufferPointer()); } } // namespace mitk #endif /* MITKIMAGE_H_HEADER_INCLUDED_C1C2FCD2 */ diff --git a/Modules/Segmentation/Interactions/mitkPaintbrushTool.cpp b/Modules/Segmentation/Interactions/mitkPaintbrushTool.cpp index b0cad5853f..48a4306a0f 100644 --- a/Modules/Segmentation/Interactions/mitkPaintbrushTool.cpp +++ b/Modules/Segmentation/Interactions/mitkPaintbrushTool.cpp @@ -1,389 +1,487 @@ /*=================================================================== 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 "mitkPaintbrushTool.h" #include "mitkToolManager.h" #include "mitkOverwriteSliceImageFilter.h" #include "mitkBaseRenderer.h" #include "mitkImageDataItem.h" #include "ipSegmentation.h" #include "mitkLevelWindowProperty.h" #define ROUND(a) ((a)>0 ? (int)((a)+0.5) : -(int)(0.5-(a))) int mitk::PaintbrushTool::m_Size = 1; mitk::PaintbrushTool::PaintbrushTool(int paintingPixelValue) :FeedbackContourTool("PressMoveReleaseWithCTRLInversionAllMouseMoves"), m_PaintingPixelValue(paintingPixelValue), m_LastContourSize(0) // other than initial mitk::PaintbrushTool::m_Size (around l. 28) { // great magic numbers CONNECT_ACTION( 80, OnMousePressed ); CONNECT_ACTION( 90, OnMouseMoved ); CONNECT_ACTION( 42, OnMouseReleased ); CONNECT_ACTION( 49014, OnInvertLogic ); m_MasterContour = Contour::New(); m_MasterContour->Initialize(); m_CurrentPlane = NULL; m_WorkingNode = DataNode::New(); m_WorkingNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( mitk::LevelWindow(0, 1) ) ); m_WorkingNode->SetProperty( "binary", mitk::BoolProperty::New(true) ); } mitk::PaintbrushTool::~PaintbrushTool() { } void mitk::PaintbrushTool::Activated() { Superclass::Activated(); FeedbackContourTool::SetFeedbackContourVisible(true); SizeChanged.Send(m_Size); } void mitk::PaintbrushTool::Deactivated() { FeedbackContourTool::SetFeedbackContourVisible(false); if (m_ToolManager->GetDataStorage()->Exists(m_WorkingNode)) m_ToolManager->GetDataStorage()->Remove(m_WorkingNode); Superclass::Deactivated(); m_WorkingSlice = NULL; m_CurrentPlane = NULL; } void mitk::PaintbrushTool::SetSize(int value) { m_Size = value; } +mitk::Point2D mitk::PaintbrushTool::upperLeft(mitk::Point2D p) +{ + p[0] -= 0.5; + p[1] += 0.5; + return p; +} + + void mitk::PaintbrushTool::UpdateContour(const StateEvent* stateEvent) { - //MITK_INFO<<"Update..."; + //MITK_INFO<<"Update..."; // examine stateEvent and create a contour that matches the pixel mask that we are going to draw const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return; - // create a copy of this slice (at least match the pixel sizes/spacings), - // then draw the desired mask on it and create a contour from it - - // Convert to ipMITKSegmentationTYPE (because getting pixels relys on that data type) - itk::Image< ipMITKSegmentationTYPE, 2 >::Pointer correctPixelTypeImage; - CastToItkImage(m_WorkingSlice/*dynamic_cast(m_WorkingNode->GetData())*/, correctPixelTypeImage ); - assert (correctPixelTypeImage.IsNotNull() ); + // Get Spacing of current Slice + //mitk::Vector3D vSpacing = m_WorkingSlice->GetSlicedGeometry()->GetGeometry2D(0)->GetSpacing(); - itk::Image< ipMITKSegmentationTYPE, 2 >::DirectionType imageDirection; - imageDirection.SetIdentity(); - correctPixelTypeImage->SetDirection(imageDirection); + // + // Draw a contour in Square according to selected brush size + // + int radius = (m_Size)/2; + float fradius = static_cast(m_Size) / 2.0f; - Image::Pointer temporarySlice = Image::New(); - CastToMitkImage( correctPixelTypeImage, temporarySlice ); + Contour::Pointer contourInImageIndexCoordinates = Contour::New(); + contourInImageIndexCoordinates->Initialize(); - //mitkIpPicDescriptor* stupidClone = mitkIpPicClone( temporarySlice->GetSliceData()->GetPicDescriptor() ); - mitkIpPicDescriptor* stupidClone = mitkIpPicNew(); - CastToIpPicDescriptor( temporarySlice->GetSliceData(), stupidClone ); - unsigned int pixelWidth = m_Size + 1; - unsigned int pixelHeight = m_Size + 1; + // estimate center point of the brush ( relative to the pixel the mouse points on ) + // -- left upper corner for even sizes, + // -- midpoint for uneven sizes + mitk::Point2D centerCorrection; + centerCorrection.Fill(0); - if ( stupidClone->n[0] <= pixelWidth || stupidClone->n[1] <= pixelHeight ) + // even --> correction of [+0.5, +0.5] + bool evenSize = ((m_Size % 2) == 0); + if( evenSize ) { - MITK_INFO << "Brush size is bigger than your working image. Reconsider this...\n" - "(Or tell your progammer until (s)he fixes this message.)" << std::endl; - mitkIpPicFree( stupidClone ); - return; + centerCorrection[0] += 0.5; + centerCorrection[1] += 0.5; } - unsigned int lineLength( stupidClone->n[0] ); - unsigned int oneContourOffset(0); - float circleCenterX = (float)m_Size / 2.0; - float circleCenterY = (float)m_Size / 2.0; - for (unsigned int x = 0; x <= pixelWidth; ++x) + // we will compute the control points for the upper left quarter part of a circle contour + std::vector< mitk::Point2D > quarterCycleUpperRight; + std::vector< mitk::Point2D > quarterCycleLowerRight; + std::vector< mitk::Point2D > quarterCycleLowerLeft; + std::vector< mitk::Point2D > quarterCycleUpperLeft; + + + mitk::Point2D curPoint; + bool curPointIsInside = true; + curPoint[0] = 0; + curPoint[1] = radius; + quarterCycleUpperRight.push_back( upperLeft(curPoint) ); + + // to estimate if a pixel is inside the circle, we need to compare against the 'outer radius' + // i.e. the distance from the midpoint [0,0] to the border of the pixel [0,radius] + //const float outer_radius = static_cast(radius) + 0.5; + + while (curPoint[1] > 0) + { + // Move right until pixel is outside circle + float curPointX_squared = 0.0f; + float curPointY_squared = (curPoint[1] - centerCorrection[1] ) * (curPoint[1] - centerCorrection[1] ); + while( curPointIsInside ) + { + // increment posX and chec + curPoint[0]++; + curPointX_squared = (curPoint[0] - centerCorrection[0] ) * (curPoint[0] - centerCorrection[0] ); + const float len = sqrt( curPointX_squared + curPointY_squared); + if ( len > fradius ) + { + // found first Pixel in this horizontal line, that is outside the circle + curPointIsInside = false; + } + } + quarterCycleUpperRight.push_back( upperLeft(curPoint) ); + + // Move down until pixel is inside circle + while( !curPointIsInside ) + { + // increment posX and chec + curPoint[1]--; + curPointY_squared = (curPoint[1] - centerCorrection[1] ) * (curPoint[1] - centerCorrection[1] ); + const float len = sqrt( curPointX_squared + curPointY_squared); + if ( len <= fradius ) + { + // found first Pixel in this horizontal line, that is outside the circle + curPointIsInside = true; + quarterCycleUpperRight.push_back( upperLeft(curPoint) ); + } + + // Quarter cycle is full, when curPoint y position is 0 + if (curPoint[1] <= 0) + break; + } + + } + + // QuarterCycle is full! Now copy quarter cycle to other quarters. + + if( !evenSize ) { - for (unsigned int y = 0; y <= pixelHeight; ++y) + std::vector< mitk::Point2D >::const_iterator it = quarterCycleUpperRight.begin(); + while( it != quarterCycleUpperRight.end() ) { - unsigned int offset = lineLength * y + x; - ipMITKSegmentationTYPE* current = (ipMITKSegmentationTYPE*)stupidClone->data + offset; - - float pixelCenterX = x + 0.5; - float pixelCenterY = y + 0.5; - - float xoff = pixelCenterX - circleCenterX; - float yoff = pixelCenterY - circleCenterY; - - bool inside = xoff * xoff + yoff * yoff < (m_Size * m_Size) / 4.0; // no idea, if this would work for ellipses - if (inside) - { - *current = 1; - oneContourOffset = offset; - } - else - { - *current = 0; - } + mitk::Point2D p; + p = *it; + + std::cout << p << std::endl; + + // the contour points in the lower right corner have same position but with negative y values + p[1] *= -1; + quarterCycleLowerRight.push_back(p); + + // the contour points in the lower left corner have same position + // but with both x,y negative + p[0] *= -1; + quarterCycleLowerLeft.push_back(p); + + // the contour points in the upper left corner have same position + // but with x negative + p[1] *= -1; + quarterCycleUpperLeft.push_back(p); + + it++; } } - - int numberOfContourPoints( 0 ); - int newBufferSize( 0 ); - float* contourPoints = ipMITKSegmentationGetContour8N( stupidClone, oneContourOffset, numberOfContourPoints, newBufferSize ); // memory allocated with malloc - if (!contourPoints) + else { - mitkIpPicFree( stupidClone ); - return; + std::vector< mitk::Point2D >::const_iterator it = quarterCycleUpperRight.begin(); + while( it != quarterCycleUpperRight.end() ) + { + mitk::Point2D p,q; + p = *it; + + q = p; + // the contour points in the lower right corner have same position but with negative y values + q[1] *= -1; + // correct for moved offset if size even = the midpoint is not the midpoint of the current pixel + // but its upper rigt corner + q[1] += 1; + quarterCycleLowerRight.push_back(q); + + q = p; + // the contour points in the lower left corner have same position + // but with both x,y negative + q[1] = -1.0f * q[1] + 1; + q[0] = -1.0f * q[0] + 1; + quarterCycleLowerLeft.push_back(q); + + // the contour points in the upper left corner have same position + // but with x negative + q = p; + q[0] *= -1; + q[0] += 1; + quarterCycleUpperLeft.push_back(q); + + it++; + } } + std::cout << "---" << std::endl; - // copy point from float* to mitk::Contour - Contour::Pointer contourInImageIndexCoordinates = Contour::New(); - contourInImageIndexCoordinates->Initialize(); - Point3D newPoint; - //ipMITKSegmentationGetContour8N returns all points, which causes vtk warnings, since the first and the last points are coincident. - //leaving the last point out, the contour is still drawn correctly - for (int index = 0; index < numberOfContourPoints-1; ++index) + // fill contour with poins in right ordering, starting with the upperRight block + mitk::Point3D tempPoint; + for (unsigned int i=0; iAddVertex( newPoint ); + tempPoint[0] = quarterCycleUpperRight[i][0]; + tempPoint[1] = quarterCycleUpperRight[i][1]; + tempPoint[2] = 0; + contourInImageIndexCoordinates->AddVertex( tempPoint ); + } + // the lower right has to be parsed in reverse order + for (int i=quarterCycleLowerRight.size()-1; i>=0; i--) + { + tempPoint[0] = quarterCycleLowerRight[i][0]; + tempPoint[1] = quarterCycleLowerRight[i][1]; + tempPoint[2] = 0; + contourInImageIndexCoordinates->AddVertex( tempPoint); + } + for (unsigned int i=0; iAddVertex( tempPoint ); + } + // the upper left also has to be parsed in reverse order + for (int i=quarterCycleUpperLeft.size()-1; i>=0; i--) + { + tempPoint[0] = quarterCycleUpperLeft[i][0]; + tempPoint[1] = quarterCycleUpperLeft[i][1]; + tempPoint[2] = 0; + contourInImageIndexCoordinates->AddVertex( tempPoint ); } - - free(contourPoints); m_MasterContour = contourInImageIndexCoordinates; - // The PicDescriptor is only REFERENCING(!) the data, the temporarySlice takes care of deleting the data also the descriptor is pointing on - // because they got allocated by the ImageDataItem, not the descriptor. - stupidClone->data = NULL; - mitkIpPicFree( stupidClone ); } /** Just show the contour, get one point as the central point and add surrounding points to the contour. */ bool mitk::PaintbrushTool::OnMousePressed (Action* action, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; m_LastEventSender = positionEvent->GetSender(); m_LastEventSlice = m_LastEventSender->GetSlice(); return this->OnMouseMoved(action, stateEvent); } /** Insert the point to the feedback contour,finish to build the contour and at the same time the painting function */ bool mitk::PaintbrushTool::OnMouseMoved (Action* itkNotUsed(action), const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); CheckIfCurrentSliceHasChanged(positionEvent); if ( m_LastContourSize != m_Size ) { UpdateContour( stateEvent ); m_LastContourSize = m_Size; } bool leftMouseButtonPressed( stateEvent->GetId() == 530 || stateEvent->GetId() == 534 || stateEvent->GetId() == 1 || stateEvent->GetId() == 5 ); Point3D worldCoordinates = positionEvent->GetWorldPosition(); Point3D indexCoordinates; m_WorkingSlice->GetGeometry()->WorldToIndex( worldCoordinates, indexCoordinates ); MITK_DEBUG << "Mouse at W " << worldCoordinates << std::endl; MITK_DEBUG << "Mouse at I " << indexCoordinates << std::endl; // round to nearest voxel center (abort if this hasn't changed) if ( m_Size % 2 == 0 ) // even { - indexCoordinates[0] = ROUND( indexCoordinates[0] /*+ 0.5*/) + 0.5; - indexCoordinates[1] = ROUND( indexCoordinates[1] /*+ 0.5*/ ) + 0.5; + indexCoordinates[0] = ROUND( indexCoordinates[0]);// /*+ 0.5*/) + 0.5; + indexCoordinates[1] = ROUND( indexCoordinates[1]);// /*+ 0.5*/ ) + 0.5; } else // odd { indexCoordinates[0] = ROUND( indexCoordinates[0] ) ; indexCoordinates[1] = ROUND( indexCoordinates[1] ) ; } static Point3D lastPos; // uninitialized: if somebody finds out how this can be initialized in a one-liner, tell me if ( fabs(indexCoordinates[0] - lastPos[0]) > mitk::eps || fabs(indexCoordinates[1] - lastPos[1]) > mitk::eps || fabs(indexCoordinates[2] - lastPos[2]) > mitk::eps || leftMouseButtonPressed ) { lastPos = indexCoordinates; } else { MITK_DEBUG << "." << std::flush; return false; } MITK_DEBUG << "Mouse at C " << indexCoordinates; Contour::Pointer contour = Contour::New(); contour->Initialize(); for (unsigned int index = 0; index < m_MasterContour->GetNumberOfPoints(); ++index) { Point3D point = m_MasterContour->GetPoints()->ElementAt(index); point[0] += indexCoordinates[ 0 ]; point[1] += indexCoordinates[ 1 ]; contour->AddVertex( point ); } if (leftMouseButtonPressed) { FeedbackContourTool::FillContourInSlice( contour, m_WorkingSlice, m_PaintingPixelValue ); m_WorkingNode->SetData(m_WorkingSlice); m_WorkingNode->Modified(); } // visualize contour Contour::Pointer displayContour = Contour::New(); displayContour->Initialize(); //for (unsigned int index = 0; index < contour->GetNumberOfPoints(); ++index) //{ // Point3D point = contour->GetPoints()->ElementAt(index); // if ( m_Size % 2 == 0 ) // even // { // point[0] += 0.5; // point[1] += 0.5; // } // displayContour->AddVertex( point ); //} displayContour = FeedbackContourTool::BackProjectContourFrom2DSlice( m_WorkingSlice->GetGeometry(), /*displayContour*/contour ); SetFeedbackContour( *displayContour ); assert( positionEvent->GetSender()->GetRenderWindow() ); RenderingManager::GetInstance()->RequestUpdate( positionEvent->GetSender()->GetRenderWindow() ); return true; } bool mitk::PaintbrushTool::OnMouseReleased(Action* /*action*/, const StateEvent* stateEvent) { //When mouse is released write segmentationresult back into image const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; this->WriteBackSegmentationResult(positionEvent, m_WorkingSlice->Clone()); return true; } /** Called when the CTRL key is pressed. Will change the painting pixel value from 0 to 1 or from 1 to 0. */ -bool mitk::PaintbrushTool::OnInvertLogic(Action* itkNotUsed(action), const StateEvent* stateEvent) +bool mitk::PaintbrushTool::OnInvertLogic(Action* itkNotUsed(action), const StateEvent* /*stateEvent*/) { // Inversion only for 0 and 1 as painting values if (m_PaintingPixelValue == 1) { m_PaintingPixelValue = 0; FeedbackContourTool::SetFeedbackContourColor( 1.0, 0.0, 0.0 ); } else if (m_PaintingPixelValue == 0) { m_PaintingPixelValue = 1; FeedbackContourTool::SetFeedbackContourColorDefault(); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); return true; } void mitk::PaintbrushTool::CheckIfCurrentSliceHasChanged(const PositionEvent *event) { const PlaneGeometry* planeGeometry( dynamic_cast (event->GetSender()->GetCurrentWorldGeometry2D() ) ); DataNode* workingNode( m_ToolManager->GetWorkingData(0) ); if (!workingNode) return; Image::Pointer image = dynamic_cast(workingNode->GetData()); if ( !image || !planeGeometry ) return; if(m_CurrentPlane.IsNull()) { m_CurrentPlane = const_cast(planeGeometry); m_WorkingSlice = SegTool2D::GetAffectedImageSliceAs2DImage(event, image)->Clone(); m_WorkingNode->ReplaceProperty( "color", workingNode->GetProperty("color") ); m_WorkingNode->SetData(m_WorkingSlice); } else { bool isSameSlice (false); isSameSlice = mitk::MatrixEqualElementWise(planeGeometry->GetIndexToWorldTransform()->GetMatrix(),m_CurrentPlane->GetIndexToWorldTransform()->GetMatrix()); isSameSlice = mitk::Equal(planeGeometry->GetIndexToWorldTransform()->GetOffset(),m_CurrentPlane->GetIndexToWorldTransform()->GetOffset()); if (!isSameSlice) { m_ToolManager->GetDataStorage()->Remove(m_WorkingNode); m_CurrentPlane = NULL; m_WorkingSlice = NULL; m_WorkingNode = NULL; m_CurrentPlane = const_cast(planeGeometry); m_WorkingSlice = SegTool2D::GetAffectedImageSliceAs2DImage(event, image)->Clone(); m_WorkingNode = mitk::DataNode::New(); m_WorkingNode->SetProperty( "levelwindow", mitk::LevelWindowProperty::New( mitk::LevelWindow(0, 1) ) ); m_WorkingNode->SetProperty( "binary", mitk::BoolProperty::New(true) ); m_WorkingNode->SetData(m_WorkingSlice); //So that the paintbrush contour vanished in the previous render window RenderingManager::GetInstance()->RequestUpdateAll(); } } if(!m_ToolManager->GetDataStorage()->Exists(m_WorkingNode)) { m_WorkingNode->SetProperty( "outline binary", mitk::BoolProperty::New(true) ); m_WorkingNode->SetProperty( "color", workingNode->GetProperty("color") ); m_WorkingNode->SetProperty( "name", mitk::StringProperty::New("Paintbrush_Node") ); m_WorkingNode->SetProperty( "helper object", mitk::BoolProperty::New(true) ); m_WorkingNode->SetProperty( "opacity", mitk::FloatProperty::New(0.8) ); m_WorkingNode->SetProperty( "includeInBoundingBox", mitk::BoolProperty::New(false)); m_WorkingNode->SetVisibility(false, mitk::BaseRenderer::GetInstance( mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget4"))); m_ToolManager->GetDataStorage()->Add(m_WorkingNode); } } diff --git a/Modules/Segmentation/Interactions/mitkPaintbrushTool.h b/Modules/Segmentation/Interactions/mitkPaintbrushTool.h index afd28140fc..4f1c66cbcf 100644 --- a/Modules/Segmentation/Interactions/mitkPaintbrushTool.h +++ b/Modules/Segmentation/Interactions/mitkPaintbrushTool.h @@ -1,96 +1,102 @@ /*=================================================================== 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 mitkPaintbrushTool_h_Included #define mitkPaintbrushTool_h_Included #include "mitkCommon.h" #include "SegmentationExports.h" #include "mitkFeedbackContourTool.h" #include "mitkPointSet.h" #include "mitkPointOperation.h" #include "mitkLegacyAdaptors.h" namespace mitk { /** \brief Paintbrush tool for InteractiveSegmentation \sa FeedbackContourTool \sa ExtractImageFilter \sa OverwriteSliceImageFilter \ingroup Interaction \ingroup ToolManagerEtAl Simple paintbrush drawing tool. Right now there are only circular pens of varying size. \todo Bug #1727: Should be modified, so that the contour is always visible, i.e.without pressing a mouse button. \warning Only to be instantiated by mitk::ToolManager. $Author: maleike $ */ class Segmentation_EXPORT PaintbrushTool : public FeedbackContourTool { public: // sent when the pen size is changed or should be updated in a GUI. Message1 SizeChanged; mitkClassMacro(PaintbrushTool, FeedbackContourTool); void SetSize(int value); protected: PaintbrushTool(int paintingPixelValue = 1); // purposely hidden virtual ~PaintbrushTool(); virtual void Activated(); virtual void Deactivated(); virtual bool OnMousePressed (Action*, const StateEvent*); virtual bool OnMouseMoved (Action*, const StateEvent*); virtual bool OnMouseReleased(Action*, const StateEvent*); virtual bool OnInvertLogic (Action*, const StateEvent*); /** * \todo This is a possible place where to introduce * different types of pens */ void UpdateContour(const StateEvent* stateEvent); + + /** + * Little helper function. Returns the upper left corner of the given pixel. + */ + mitk::Point2D upperLeft(mitk::Point2D p); + /** * Checks if the current slice has changed */ void CheckIfCurrentSliceHasChanged(const PositionEvent* event); int m_PaintingPixelValue; static int m_Size; Contour::Pointer m_MasterContour; int m_LastContourSize; Image::Pointer m_WorkingSlice; PlaneGeometry::Pointer m_CurrentPlane; DataNode::Pointer m_WorkingNode; }; } // namespace #endif