diff --git a/Modules/AlgorithmsExt/mitkAutoCropImageFilter.cpp b/Modules/AlgorithmsExt/mitkAutoCropImageFilter.cpp index 72cba5d31e..4f9a7ca59a 100644 --- a/Modules/AlgorithmsExt/mitkAutoCropImageFilter.cpp +++ b/Modules/AlgorithmsExt/mitkAutoCropImageFilter.cpp @@ -1,363 +1,363 @@ /*=================================================================== 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 "mitkAutoCropImageFilter.h" #include "mitkImageCast.h" #include "mitkImageAccessByItk.h" #include "mitkGeometry3D.h" #include "mitkStatusBar.h" #include "mitkImageReadAccessor.h" #include "mitkPlaneGeometry.h" #include #include #include mitk::AutoCropImageFilter::AutoCropImageFilter() : m_BackgroundValue(0), m_MarginFactor(1.0), m_TimeSelector(NULL), m_OverrideCroppingRegion(false) { } mitk::AutoCropImageFilter::~AutoCropImageFilter() { } template < typename TPixel, unsigned int VImageDimension> void mitk::AutoCropImageFilter::ITKCrop3DImage( itk::Image< TPixel, VImageDimension >* inputItkImage, unsigned int timestep) { if (inputItkImage == NULL) { mitk::StatusBar::GetInstance()->DisplayErrorText ("An internal error occurred. Can't convert Image. Please report to bugs@mitk.org"); MITK_ERROR << "image is NULL...returning" << std::endl; return; } typedef itk::Image< TPixel, VImageDimension > InternalImageType; typedef typename InternalImageType::Pointer InternalImagePointer; typedef itk::RegionOfInterestImageFilter < InternalImageType, InternalImageType > ROIFilterType; typedef typename itk::RegionOfInterestImageFilter < InternalImageType, InternalImageType >::Pointer ROIFilterPointer; InternalImagePointer outputItk = InternalImageType::New(); ROIFilterPointer roiFilter = ROIFilterType::New(); roiFilter->SetInput(0,inputItkImage); roiFilter->SetRegionOfInterest(this->GetCroppingRegion()); roiFilter->Update(); outputItk = roiFilter->GetOutput(); outputItk->DisconnectPipeline(); mitk::Image::Pointer newMitkImage = mitk::Image::New(); mitk::CastToMitkImage( outputItk, newMitkImage ); MITK_INFO << "Crop-Output dimension: " << (newMitkImage->GetDimension() == 3) << " Filter-Output dimension: "<GetOutput()->GetDimension()<< " Timestep: " << timestep; mitk::ImageReadAccessor newMitkImgAcc(newMitkImage); this->GetOutput()->SetVolume( newMitkImgAcc.GetData(), timestep); } void mitk::AutoCropImageFilter::GenerateOutputInformation() { mitk::Image::Pointer input = const_cast (this->GetInput()); mitk::Image::Pointer output = this->GetOutput(); if(input->GetDimension() <= 2) { MITK_ERROR << "Only 3D any 4D images are supported." << std::endl; return; } ComputeNewImageBounds(); if ((output->IsInitialized()) && (output->GetPipelineMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; itkDebugMacro(<<"GenerateOutputInformation()"); // PART I: initialize input requested region. We do this already here (and not // later when GenerateInputRequestedRegion() is called), because we // also need the information to setup the output. // pre-initialize input-requested-region to largest-possible-region // and correct time-region; spatial part will be cropped by // bounding-box of bounding-object below m_InputRequestedRegion = input->GetLargestPossibleRegion(); // build region out of index and size calculated in ComputeNewImageBounds() mitk::SlicedData::IndexType index; index[0] = m_RegionIndex[0]; index[1] = m_RegionIndex[1]; index[2] = m_RegionIndex[2]; index[3] = m_InputRequestedRegion.GetIndex()[3]; index[4] = m_InputRequestedRegion.GetIndex()[4]; mitk::SlicedData::SizeType size; size[0] = m_RegionSize[0]; size[1] = m_RegionSize[1]; size[2] = m_RegionSize[2]; size[3] = m_InputRequestedRegion.GetSize()[3]; size[4] = m_InputRequestedRegion.GetSize()[4]; mitk::SlicedData::RegionType cropRegion(index, size); // crop input-requested-region with cropping region computed from the image data if(m_InputRequestedRegion.Crop(cropRegion)==false) { // crop not possible => do nothing: set time size to 0. size.Fill(0); m_InputRequestedRegion.SetSize(size); return; } // set input-requested-region, because we access it later in // GenerateInputRequestedRegion (there we just set the time) input->SetRequestedRegion(&m_InputRequestedRegion); // PART II: initialize output image unsigned int dimension = input->GetDimension(); unsigned int *dimensions = new unsigned int [dimension]; itk2vtk(m_InputRequestedRegion.GetSize(), dimensions); if(dimension>3) memcpy(dimensions+3, input->GetDimensions()+3, (dimension-3)*sizeof(unsigned int)); // create basic slicedGeometry that will be initialized below output->Initialize(mitk::PixelType( GetOutputPixelType() ), dimension, dimensions); delete [] dimensions; //clone the IndexToWorldTransform from the input, otherwise we will overwrite it, when adjusting the origin of the output image!! itk::ScalableAffineTransform< mitk::ScalarType,3 >::Pointer cloneTransform = itk::ScalableAffineTransform< mitk::ScalarType,3 >::New(); cloneTransform->Compose(input->GetGeometry()->GetIndexToWorldTransform()); output->GetGeometry()->SetIndexToWorldTransform( cloneTransform.GetPointer() ); // Position the output Image to match the corresponding region of the input image mitk::SlicedGeometry3D* slicedGeometry = output->GetSlicedGeometry(); mitk::SlicedGeometry3D::Pointer inputGeometry = input->GetSlicedGeometry(); const mitk::SlicedData::IndexType& start = m_InputRequestedRegion.GetIndex(); mitk::Point3D origin; vtk2itk(start, origin); input->GetSlicedGeometry()->IndexToWorld(origin, origin); slicedGeometry->SetOrigin(origin); // get the PlaneGeometry for the first slice of the original image mitk::PlaneGeometry::Pointer plane = dynamic_cast( inputGeometry->GetGeometry2D( 0 )->Clone().GetPointer() ); assert( plane ); // re-initialize the plane according to the new requirements: // dimensions of the cropped image // right- and down-vector as well as spacing do not change, so use the ones from // input image ScalarType dimX = output->GetDimensions()[0]; ScalarType dimY = output->GetDimensions()[1]; mitk::Vector3D right = plane->GetAxisVector(0); mitk::Vector3D down = plane->GetAxisVector(1); mitk::Vector3D spacing = plane->GetSpacing(); plane->InitializeStandardPlane( dimX, dimY, right, down, &spacing ); // set the new origin on the PlaneGeometry as well plane->SetOrigin(origin); // re-initialize the slicedGeometry with the correct planeGeometry // in order to get a fully initialized SlicedGeometry3D slicedGeometry->InitializeEvenlySpaced( plane, inputGeometry->GetSpacing()[2], output->GetSlicedGeometry()->GetSlices() ); mitk::TimeGeometry* timeSlicedGeometry = output->GetTimeGeometry(); mitk::ProportionalTimeGeometry* propTimeGeometry = dynamic_cast(timeSlicedGeometry); propTimeGeometry->Initialize(slicedGeometry, output->GetDimension(3)); m_TimeOfHeaderInitialization.Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); } void mitk::AutoCropImageFilter::GenerateData() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if(input.IsNull()) return; if(input->GetDimension() <= 2) { MITK_ERROR << "Only 3D and 4D images supported"; return; } if((output->IsInitialized()==false) ) return; if( m_TimeSelector.IsNull() ) m_TimeSelector = mitk::ImageTimeSelector::New(); m_TimeSelector->SetInput(input); mitk::SlicedData::RegionType outputRegion = input->GetRequestedRegion(); int tstart = outputRegion.GetIndex(3); int tmax = tstart + outputRegion.GetSize(3); for( int timestep=tstart;timestepSetTimeNr(timestep); m_TimeSelector->UpdateLargestPossibleRegion(); AccessFixedDimensionByItk_1( m_TimeSelector->GetOutput(), ITKCrop3DImage, 3, timestep ); } // this->GetOutput()->Update(); // Not sure if this is necessary... m_TimeOfHeaderInitialization.Modified(); } void mitk::AutoCropImageFilter::ComputeNewImageBounds() { mitk::Image::ConstPointer inputMitk = this->GetInput(); if (m_OverrideCroppingRegion) { for (unsigned int i=0; i<3; ++i) { m_RegionIndex[i] = m_CroppingRegion.GetIndex()[i]; m_RegionSize[i] = m_CroppingRegion.GetSize()[i]; - if (m_RegionIndex[i] >= inputMitk->GetDimension(i)) + if (m_RegionIndex[i] >= static_cast(inputMitk->GetDimension(i))) { itkExceptionMacro("Cropping index is not inside the image. " << std::endl << "Index:" << std::endl << m_CroppingRegion.GetIndex() << std::endl << "Size:" << std::endl << m_CroppingRegion.GetSize()); } if (m_RegionIndex[i] + m_RegionSize[i] >= inputMitk->GetDimension(i)) { m_RegionSize[i] = inputMitk->GetDimension(i) - m_RegionIndex[i]; } } for (unsigned int i=0; i<3; ++i) { m_RegionIndex[i] = m_CroppingRegion.GetIndex()[i]; m_RegionSize[i] = m_CroppingRegion.GetSize()[i]; } } else { // Check if a 3D or 4D image is present unsigned int timeSteps = 1; if (inputMitk->GetDimension() == 4 ) timeSteps = inputMitk->GetDimension(3); ImageType::IndexType minima,maxima; if (inputMitk->GetDimension() == 4) { // initialize with time step 0 m_TimeSelector = mitk::ImageTimeSelector::New(); m_TimeSelector->SetInput( inputMitk ); m_TimeSelector->SetTimeNr( 0 ); m_TimeSelector->UpdateLargestPossibleRegion(); inputMitk = m_TimeSelector->GetOutput(); } ImagePointer inputItk = ImageType::New(); mitk::CastToItkImage( inputMitk , inputItk ); // it is assumed that all volumes in a time series have the same 3D dimensions ImageType::RegionType origRegion = inputItk->GetLargestPossibleRegion(); // Initialize min and max on the first (or only) time step maxima = inputItk->GetLargestPossibleRegion().GetIndex(); minima[0] = inputItk->GetLargestPossibleRegion().GetSize()[0]; minima[1] = inputItk->GetLargestPossibleRegion().GetSize()[1]; minima[2] = inputItk->GetLargestPossibleRegion().GetSize()[2]; typedef itk::ImageRegionConstIterator< ImageType > ConstIteratorType; for(unsigned int idx = 0; idx < timeSteps; ++idx) { // if 4D image, update time step and itk image if( idx > 0) { m_TimeSelector->SetTimeNr( idx ); m_TimeSelector->UpdateLargestPossibleRegion(); inputMitk = m_TimeSelector->GetOutput(); mitk::CastToItkImage( inputMitk , inputItk ); } ConstIteratorType inIt( inputItk, origRegion ); for ( inIt.GoToBegin(); !inIt.IsAtEnd(); ++inIt) { float pix_val = inIt.Get(); if ( fabs(pix_val - m_BackgroundValue) > mitk::eps ) { for (int i=0; i < 3; i++) { minima[i] = vnl_math_min((int)minima[i],(int)(inIt.GetIndex()[i])); maxima[i] = vnl_math_max((int)maxima[i],(int)(inIt.GetIndex()[i])); } } } } typedef ImageType::RegionType::SizeType::SizeValueType SizeValueType; m_RegionSize[0] = (SizeValueType)(m_MarginFactor * (maxima[0] - minima[0] + 1 )); m_RegionSize[1] = (SizeValueType)(m_MarginFactor * (maxima[1] - minima[1] + 1 )); m_RegionSize[2] = (SizeValueType)(m_MarginFactor * (maxima[2] - minima[2] + 1 )); m_RegionIndex = minima; m_RegionIndex[0] -= (m_RegionSize[0] - maxima[0] + minima[0] - 1 )/2; m_RegionIndex[1] -= (m_RegionSize[1] - maxima[1] + minima[1] - 1 )/2; m_RegionIndex[2] -= (m_RegionSize[2] - maxima[2] + minima[2] - 1 )/2; ImageType::RegionType cropRegion(m_RegionIndex,m_RegionSize); origRegion.Crop(cropRegion); m_RegionSize[0] = origRegion.GetSize()[0]; m_RegionSize[1] = origRegion.GetSize()[1]; m_RegionSize[2] = origRegion.GetSize()[2]; m_RegionIndex[0] = origRegion.GetIndex()[0]; m_RegionIndex[1] = origRegion.GetIndex()[1]; m_RegionIndex[2] = origRegion.GetIndex()[2]; m_CroppingRegion = origRegion; } } void mitk::AutoCropImageFilter::GenerateInputRequestedRegion() { } const mitk::PixelType mitk::AutoCropImageFilter::GetOutputPixelType() { return this->GetInput()->GetPixelType(); } void mitk::AutoCropImageFilter::SetCroppingRegion(RegionType overrideRegion) { m_CroppingRegion = overrideRegion; m_OverrideCroppingRegion = true; } diff --git a/Modules/DataTypesExt/mitkAffineDataInteractor3D.cpp b/Modules/DataTypesExt/mitkAffineDataInteractor3D.cpp index 021fa2db0c..715854ca22 100644 --- a/Modules/DataTypesExt/mitkAffineDataInteractor3D.cpp +++ b/Modules/DataTypesExt/mitkAffineDataInteractor3D.cpp @@ -1,329 +1,329 @@ /*=================================================================== 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 "mitkAffineDataInteractor3D.h" #include "mitkInteractionConst.h" #include "mitkInteractionPositionEvent.h" #include "mitkRotationOperation.h" #include "mitkSurface.h" #include #include #include #include #include mitk::AffineDataInteractor3D::AffineDataInteractor3D() { m_OriginalGeometry = Geometry3D::New(); // Initialize vector arithmetic m_ObjectNormal[0] = 0.0; m_ObjectNormal[1] = 0.0; m_ObjectNormal[2] = 1.0; } mitk::AffineDataInteractor3D::~AffineDataInteractor3D() { } void mitk::AffineDataInteractor3D::ConnectActionsAndFunctions() { // **Conditions** that can be used in the state machine, to ensure that certain conditions are met, before actually executing an action CONNECT_CONDITION("isOverObject", CheckOverObject); // **Function** in the statmachine patterns also referred to as **Actions** CONNECT_FUNCTION("selectObject",SelectObject); CONNECT_FUNCTION("deselectObject",DeselectObject); CONNECT_FUNCTION("initTranslate",InitTranslate); CONNECT_FUNCTION("initRotate",InitRotate); CONNECT_FUNCTION("translateObject",TranslateObject); CONNECT_FUNCTION("rotateObject",RotateObject); } void mitk::AffineDataInteractor3D::DataNodeChanged() { } bool mitk::AffineDataInteractor3D::CheckOverObject(const InteractionEvent* interactionEvent) { const InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; Point3D currentPickedPoint; if(interactionEvent->GetSender()->PickObject(positionEvent->GetPointerPositionOnScreen(), currentPickedPoint) == this->GetDataNode().GetPointer()) return true; return false; } bool mitk::AffineDataInteractor3D::SelectObject(StateMachineAction*, InteractionEvent* interactionEvent) { DataNode::Pointer node = this->GetDataNode(); if (node.IsNull()) return false; node->SetColor(1.0, 0.0, 0.0); // Colorize surface / wireframe dependend on distance from picked point this->ColorizeSurface(interactionEvent->GetSender(), 0.0); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::AffineDataInteractor3D::DeselectObject(StateMachineAction*, InteractionEvent* interactionEvent) { DataNode::Pointer node = this->GetDataNode(); if (node.IsNull()) return false; node->SetColor( 1.0, 1.0, 1.0 ); // Colorize surface / wireframe as inactive this->ColorizeSurface(interactionEvent->GetSender(), -1.0); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::AffineDataInteractor3D::InitTranslate(StateMachineAction*, InteractionEvent* interactionEvent) { InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; m_InitialPickedPoint = positionEvent->GetPositionInWorld(); m_InitialPickedDisplayPoint = positionEvent->GetPointerPositionOnScreen(); // Get the timestep to also support 3D+t int timeStep = 0; if ((interactionEvent->GetSender()) != NULL) timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); // Make deep copy of current Geometry3D of the plane this->GetDataNode()->GetData()->UpdateOutputInformation(); // make sure that the Geometry is up-to-date m_OriginalGeometry = static_cast(this->GetDataNode()->GetData()->GetGeometry(timeStep)->Clone().GetPointer()); return true; } bool mitk::AffineDataInteractor3D::InitRotate(StateMachineAction*, InteractionEvent* interactionEvent) { InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; m_InitialPickedPoint = positionEvent->GetPositionInWorld(); m_InitialPickedDisplayPoint = positionEvent->GetPointerPositionOnScreen(); // Get the timestep to also support 3D+t int timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); // Make deep copy of current Geometry3D of the plane this->GetDataNode()->GetData()->UpdateOutputInformation(); // make sure that the Geometry is up-to-date m_OriginalGeometry = static_cast(this->GetDataNode()->GetData()->GetGeometry(timeStep)->Clone().GetPointer()); return true; } bool mitk::AffineDataInteractor3D::TranslateObject (StateMachineAction*, InteractionEvent* interactionEvent) { InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; Point3D currentPickedPoint = positionEvent->GetPositionInWorld(); Vector3D interactionMove; interactionMove[0] = currentPickedPoint[0] - m_InitialPickedPoint[0]; interactionMove[1] = currentPickedPoint[1] - m_InitialPickedPoint[1]; interactionMove[2] = currentPickedPoint[2] - m_InitialPickedPoint[2]; Point3D origin = m_OriginalGeometry->GetOrigin(); // Get the timestep to also support 3D+t int timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); // If data is an mitk::Surface, extract it Surface::Pointer surface = dynamic_cast< Surface* >(this->GetDataNode()->GetData()); vtkPolyData* polyData = NULL; if (surface.IsNotNull()) { polyData = surface->GetVtkPolyData( timeStep ); // Extract surface normal from surface (if existent, otherwise use default) vtkPointData* pointData = polyData->GetPointData(); if (pointData != NULL) { vtkDataArray* normal = polyData->GetPointData()->GetVectors("planeNormal"); if (normal != NULL) { m_ObjectNormal[0] = normal->GetComponent( 0, 0 ); m_ObjectNormal[1] = normal->GetComponent( 0, 1 ); m_ObjectNormal[2] = normal->GetComponent( 0, 2 ); } } } Vector3D transformedObjectNormal; this->GetDataNode()->GetData()->GetGeometry( timeStep )->IndexToWorld(m_ObjectNormal, transformedObjectNormal); this->GetDataNode()->GetData()->GetGeometry( timeStep )->SetOrigin(origin + transformedObjectNormal * (interactionMove * transformedObjectNormal)); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::AffineDataInteractor3D::RotateObject (StateMachineAction*, InteractionEvent* interactionEvent) { InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; Point2D currentPickedDisplayPoint = positionEvent->GetPointerPositionOnScreen(); Point3D currentPickedPoint = positionEvent->GetPositionInWorld(); vtkCamera* camera = NULL; vtkRenderer* currentVtkRenderer = NULL; if ((interactionEvent->GetSender()) != NULL) { vtkRenderWindow* renderWindow = interactionEvent->GetSender()->GetRenderWindow(); if (renderWindow != NULL) { vtkRenderWindowInteractor* renderWindowInteractor = renderWindow->GetInteractor(); if ( renderWindowInteractor != NULL ) { currentVtkRenderer = renderWindowInteractor->GetInteractorStyle()->GetCurrentRenderer(); if (currentVtkRenderer != NULL) camera = currentVtkRenderer->GetActiveCamera(); } } } if ( camera ) { double vpn[3]; camera->GetViewPlaneNormal( vpn ); Vector3D viewPlaneNormal; viewPlaneNormal[0] = vpn[0]; viewPlaneNormal[1] = vpn[1]; viewPlaneNormal[2] = vpn[2]; Vector3D interactionMove; interactionMove[0] = currentPickedPoint[0] - m_InitialPickedPoint[0]; interactionMove[1] = currentPickedPoint[1] - m_InitialPickedPoint[1]; interactionMove[2] = currentPickedPoint[2] - m_InitialPickedPoint[2]; if (interactionMove[0] == 0 && interactionMove[1] == 0 && interactionMove[2] == 0) return true; Vector3D rotationAxis = itk::CrossProduct(viewPlaneNormal, interactionMove); rotationAxis.Normalize(); int* size = currentVtkRenderer->GetSize(); double l2 = (currentPickedDisplayPoint[0] - m_InitialPickedDisplayPoint[0]) * (currentPickedDisplayPoint[0] - m_InitialPickedDisplayPoint[0]) + (currentPickedDisplayPoint[1] - m_InitialPickedDisplayPoint[1]) * (currentPickedDisplayPoint[1] - m_InitialPickedDisplayPoint[1]); double rotationAngle = 360.0 * sqrt(l2 / (size[0] * size[0] + size[1] * size[1])); // Use center of data bounding box as center of rotation Point3D rotationCenter = m_OriginalGeometry->GetCenter(); int timeStep = 0; if ((interactionEvent->GetSender()) != NULL) timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); // Reset current Geometry3D to original state (pre-interaction) and // apply rotation RotationOperation op( OpROTATE, rotationCenter, rotationAxis, rotationAngle ); Geometry3D::Pointer newGeometry = static_cast(m_OriginalGeometry->Clone().GetPointer()); newGeometry->ExecuteOperation( &op ); mitk::TimeGeometry::Pointer timeGeometry = this->GetDataNode()->GetData()->GetTimeGeometry(); if (timeGeometry.IsNotNull()) timeGeometry->SetTimeStepGeometry(newGeometry, timeStep); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } else return false; } bool mitk::AffineDataInteractor3D::ColorizeSurface(BaseRenderer::Pointer renderer, double scalar) { BaseData::Pointer data = this->GetDataNode()->GetData(); if(data.IsNull()) { MITK_ERROR << "AffineInteractor3D: No data object present!"; return false; } // Get the timestep to also support 3D+t int timeStep = 0; if (renderer.IsNotNull()) timeStep = renderer->GetTimeStep(data); // If data is an mitk::Surface, extract it Surface::Pointer surface = dynamic_cast< Surface* >(data.GetPointer()); vtkPolyData* polyData = NULL; if (surface.IsNotNull()) polyData = surface->GetVtkPolyData(timeStep); if (polyData == NULL) { MITK_ERROR << "AffineInteractor3D: No poly data present!"; return false; } vtkPointData* pointData = polyData->GetPointData(); if (pointData == NULL) { MITK_ERROR << "AffineInteractor3D: No point data present!"; return false; } vtkDataArray* scalars = pointData->GetScalars(); if (scalars == NULL) { MITK_ERROR << "AffineInteractor3D: No scalars for point data present!"; return false; } - for (unsigned int i = 0; i < pointData->GetNumberOfTuples(); ++i) + for (vtkIdType i = 0; i < pointData->GetNumberOfTuples(); ++i) { scalars->SetComponent(i, 0, scalar); } polyData->Modified(); pointData->Update(); return true; } diff --git a/Modules/DataTypesExt/mitkSurfaceDeformationDataInteractor3D.cpp b/Modules/DataTypesExt/mitkSurfaceDeformationDataInteractor3D.cpp index 2cb9543eb1..ad946c4722 100644 --- a/Modules/DataTypesExt/mitkSurfaceDeformationDataInteractor3D.cpp +++ b/Modules/DataTypesExt/mitkSurfaceDeformationDataInteractor3D.cpp @@ -1,302 +1,302 @@ /*=================================================================== 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 "mitkSurfaceDeformationDataInteractor3D.h" #include "mitkMouseWheelEvent.h" #include #include mitk::SurfaceDeformationDataInteractor3D::SurfaceDeformationDataInteractor3D() :m_GaussSigma(30.0) { m_OriginalPolyData = vtkPolyData::New(); // Initialize vector arithmetic m_ObjectNormal[0] = 0.0; m_ObjectNormal[1] = 0.0; m_ObjectNormal[2] = 1.0; } mitk::SurfaceDeformationDataInteractor3D::~SurfaceDeformationDataInteractor3D() { m_OriginalPolyData->Delete(); } void mitk::SurfaceDeformationDataInteractor3D::ConnectActionsAndFunctions() { // **Conditions** that can be used in the state machine, to ensure that certain conditions are met, before // actually executing an action CONNECT_CONDITION("isOverObject", CheckOverObject); // **Function** in the statmachine patterns also referred to as **Actions** CONNECT_FUNCTION("selectObject",SelectObject); CONNECT_FUNCTION("deselectObject",DeselectObject); CONNECT_FUNCTION("initDeformation",InitDeformation); CONNECT_FUNCTION("deformObject",DeformObject); CONNECT_FUNCTION("scaleRadius", ScaleRadius); } void mitk::SurfaceDeformationDataInteractor3D::DataNodeChanged() { if(this->GetDataNode().IsNotNull()) { m_Surface = dynamic_cast(this->GetDataNode()->GetData()); if (m_Surface == NULL) MITK_ERROR << "SurfaceDeformationDataInteractor3D::DataNodeChanged(): DataNode has to contain a surface."; } else m_Surface = NULL; } bool mitk::SurfaceDeformationDataInteractor3D::CheckOverObject(const InteractionEvent* interactionEvent) { const InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; Point2D currentPickedDisplayPoint = positionEvent->GetPointerPositionOnScreen(); Point3D currentPickedPoint; if(interactionEvent->GetSender()->PickObject(currentPickedDisplayPoint, currentPickedPoint) == this->GetDataNode().GetPointer()) { // Colorized surface at current picked position m_SurfaceColorizationCenter = currentPickedPoint; return true; } return false; } bool mitk::SurfaceDeformationDataInteractor3D::SelectObject(StateMachineAction*, InteractionEvent* interactionEvent) { const InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; int timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); vtkPolyData* polyData = m_Surface->GetVtkPolyData(timeStep); this->GetDataNode()->SetColor(1.0, 0.0, 0.0); // Colorize surface / wireframe dependend on distance from picked point this->ColorizeSurface(polyData, timeStep, m_SurfaceColorizationCenter, COLORIZATION_GAUSS); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::SurfaceDeformationDataInteractor3D::DeselectObject(StateMachineAction*, InteractionEvent* interactionEvent) { const InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; int timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); vtkPolyData* polyData = m_Surface->GetVtkPolyData(timeStep); this->GetDataNode()->SetColor(1.0, 1.0, 1.0); // Colorize surface / wireframe as inactive this->ColorizeSurface(polyData, timeStep, m_SurfaceColorizationCenter, COLORIZATION_CONSTANT, -1.0); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::SurfaceDeformationDataInteractor3D::InitDeformation(StateMachineAction*, InteractionEvent* interactionEvent) { const InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; int timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); vtkPolyData* polyData = m_Surface->GetVtkPolyData(timeStep); // Store current picked point interactionEvent->GetSender()->PickObject(positionEvent->GetPointerPositionOnScreen(), m_InitialPickedPoint); // Make deep copy of vtkPolyData interacted on m_OriginalPolyData->DeepCopy(polyData); return true; } bool mitk::SurfaceDeformationDataInteractor3D::DeformObject (StateMachineAction*, InteractionEvent* interactionEvent) { const InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if(positionEvent == NULL) return false; int timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); vtkPolyData* polyData = m_Surface->GetVtkPolyData(timeStep); Geometry3D::Pointer geometry = this->GetDataNode()->GetData()->GetGeometry(timeStep); Point3D currentPickedPoint = positionEvent->GetPositionInWorld(); // Calculate mouse move in 3D space Vector3D interactionMove; interactionMove[0] = currentPickedPoint[0] - m_InitialPickedPoint[0]; interactionMove[1] = currentPickedPoint[1] - m_InitialPickedPoint[1]; interactionMove[2] = currentPickedPoint[2] - m_InitialPickedPoint[2]; // Transform mouse move into geometry space this->GetDataNode()->GetData()->UpdateOutputInformation();// make sure that the Geometry is up-to-date Vector3D interactionMoveIndex; geometry->WorldToIndex(interactionMove, interactionMoveIndex); // Get picked point and transform into local coordinates Point3D pickedPoint; geometry->WorldToIndex(m_InitialPickedPoint, pickedPoint); Vector3D v1 = pickedPoint.GetVectorFromOrigin(); vtkDataArray* normal = polyData->GetPointData()->GetVectors("planeNormal"); if (normal != NULL) { m_ObjectNormal[0] = normal->GetComponent(0, 0); m_ObjectNormal[1] = normal->GetComponent(0, 1); m_ObjectNormal[2] = normal->GetComponent(0, 2); } Vector3D v2 = m_ObjectNormal * (interactionMoveIndex * m_ObjectNormal); vtkPoints* originalPoints = m_OriginalPolyData->GetPoints(); vtkPoints* deformedPoints = polyData->GetPoints(); double denom = m_GaussSigma * m_GaussSigma * 2; double point[3]; - for (unsigned int i = 0; i < deformedPoints->GetNumberOfPoints(); ++i) + for (vtkIdType i = 0; i < deformedPoints->GetNumberOfPoints(); ++i) { // Get original point double* originalPoint = originalPoints->GetPoint( i ); Vector3D v0; v0[0] = originalPoint[0]; v0[1] = originalPoint[1]; v0[2] = originalPoint[2]; // Calculate distance of this point from line through picked point double d = itk::CrossProduct(m_ObjectNormal, (v1 - v0)).GetNorm(); Vector3D t = v2 * exp(- d * d / denom); point[0] = originalPoint[0] + t[0]; point[1] = originalPoint[1] + t[1]; point[2] = originalPoint[2] + t[2]; deformedPoints->SetPoint(i, point); } // Make sure that surface is colorized at initial picked position as long as we are in deformation state m_SurfaceColorizationCenter = m_InitialPickedPoint; polyData->Modified(); m_Surface->Modified(); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::SurfaceDeformationDataInteractor3D::ScaleRadius(StateMachineAction*, InteractionEvent* interactionEvent) { const MouseWheelEvent* wheelEvent = dynamic_cast(interactionEvent); if(wheelEvent == NULL) return false; m_GaussSigma += (double) (wheelEvent->GetWheelDelta()) / 20; if ( m_GaussSigma < 10.0 ) { m_GaussSigma = 10.0; } else if ( m_GaussSigma > 128.0 ) { m_GaussSigma = 128.0; } int timeStep = interactionEvent->GetSender()->GetTimeStep(this->GetDataNode()->GetData()); vtkPolyData* polyData = m_Surface->GetVtkPolyData(timeStep); // Colorize surface / wireframe dependend on sigma and distance from picked point this->ColorizeSurface( polyData, timeStep, m_SurfaceColorizationCenter, COLORIZATION_GAUSS ); interactionEvent->GetSender()->GetRenderingManager()->RequestUpdateAll(); return true; } bool mitk::SurfaceDeformationDataInteractor3D::ColorizeSurface(vtkPolyData* polyData, int timeStep, const Point3D &pickedPoint, int mode, double scalar) { if (polyData == NULL) return false; vtkPoints* points = polyData->GetPoints(); vtkPointData* pointData = polyData->GetPointData(); if ( pointData == NULL ) return false; vtkDataArray* scalars = pointData->GetScalars(); if (scalars == NULL) return false; if (mode == COLORIZATION_GAUSS) { // Get picked point and transform into local coordinates Point3D localPickedPoint; Geometry3D::Pointer geometry = this->GetDataNode()->GetData()->GetGeometry(timeStep); geometry->WorldToIndex( pickedPoint, localPickedPoint ); Vector3D v1 = localPickedPoint.GetVectorFromOrigin(); vtkDataArray* normal = polyData->GetPointData()->GetVectors("planeNormal"); if (normal != NULL) { m_ObjectNormal[0] = normal->GetComponent(0, 0); m_ObjectNormal[1] = normal->GetComponent(0, 1); m_ObjectNormal[2] = normal->GetComponent(0, 2); } double denom = m_GaussSigma * m_GaussSigma * 2; - for (unsigned int i = 0; i < points->GetNumberOfPoints(); ++i) + for (vtkIdType i = 0; i < points->GetNumberOfPoints(); ++i) { // Get original point double* point = points->GetPoint(i); Vector3D v0; v0[0] = point[0]; v0[1] = point[1]; v0[2] = point[2]; // Calculate distance of this point from line through picked point double d = itk::CrossProduct(m_ObjectNormal, (v1 - v0)).GetNorm(); double t = exp(- d * d / denom); scalars->SetComponent(i, 0, t); } } else if (mode == COLORIZATION_CONSTANT) { - for (unsigned int i = 0; i < pointData->GetNumberOfTuples(); ++i) + for (vtkIdType i = 0; i < pointData->GetNumberOfTuples(); ++i) { scalars->SetComponent(i, 0, scalar); } } polyData->Modified(); pointData->Update(); return true; } diff --git a/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp b/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp index 9c6ab2bdbd..c0430028d5 100644 --- a/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkComputeContourSetNormalsFilter.cpp @@ -1,321 +1,321 @@ /*=================================================================== 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 "mitkComputeContourSetNormalsFilter.h" #include "mitkImagePixelReadAccessor.h" mitk::ComputeContourSetNormalsFilter::ComputeContourSetNormalsFilter() { m_MaxSpacing = 5; this->m_UseProgressBar = false; this->m_ProgressStepSize = 1; mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); } mitk::ComputeContourSetNormalsFilter::~ComputeContourSetNormalsFilter() { } void mitk::ComputeContourSetNormalsFilter::GenerateData() { unsigned int numberOfInputs = this->GetNumberOfIndexedInputs(); this->CreateOutputsForAllInputs(numberOfInputs); //Iterating over each input for(unsigned int i = 0; i < numberOfInputs; i++) { //Getting the inputs polydata and polygons Surface* currentSurface = const_cast( this->GetInput(i) ); vtkPolyData* polyData = currentSurface->GetVtkPolyData(); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType* cell (NULL); vtkIdType cellSize (0); //The array that contains all the vertex normals of the current polygon vtkSmartPointer normals = vtkSmartPointer::New(); normals->SetNumberOfComponents(3); normals->SetNumberOfTuples(polyData->GetNumberOfPoints()); //If the current contour is an inner contour then the direction is -1 //A contour lies inside another one if the pixel values in the direction of the normal is 1 m_NegativeNormalCounter = 0; m_PositiveNormalCounter = 0; vtkIdType offSet (0); //Iterating over each polygon for( existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { if(cellSize < 3)continue; //First we calculate the current polygon's normal double polygonNormal[3] = {0.0}; double p1[3]; double p2[3]; double v1[3]; double v2[3]; existingPoints->GetPoint(cell[0], p1); unsigned int index = cellSize*0.5; existingPoints->GetPoint(cell[index], p2); v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; - for (unsigned int k = 2; k < cellSize; k++) + for (vtkIdType k = 2; k < cellSize; k++) { index = cellSize*0.25; existingPoints->GetPoint(cell[index], p1); index = cellSize*0.75; existingPoints->GetPoint(cell[index], p2); v2[0] = p2[0]-p1[0]; v2[1] = p2[1]-p1[1]; v2[2] = p2[2]-p1[2]; vtkMath::Cross(v1,v2,polygonNormal); if (vtkMath::Norm(polygonNormal) != 0) break; } vtkMath::Normalize(polygonNormal); //Now we start computing the normal for each vertex double vertexNormalTemp[3]; existingPoints->GetPoint(cell[0], p1); existingPoints->GetPoint(cell[1], p2); v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; vtkMath::Cross(v1,polygonNormal,vertexNormalTemp); vtkMath::Normalize(vertexNormalTemp); double vertexNormal[3]; - for (unsigned j = 0; j < cellSize-2; j++) + for (vtkIdType j = 0; j < cellSize-2; j++) { existingPoints->GetPoint(cell[j+1], p1); existingPoints->GetPoint(cell[j+2], p2); v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; vtkMath::Cross(v1,polygonNormal,vertexNormal); vtkMath::Normalize(vertexNormal); double finalNormal[3]; finalNormal[0] = (vertexNormal[0] + vertexNormalTemp[0])*0.5; finalNormal[1] = (vertexNormal[1] + vertexNormalTemp[1])*0.5; finalNormal[2] = (vertexNormal[2] + vertexNormalTemp[2])*0.5; //Here we determine the direction of the normal if (m_SegmentationBinaryImage) { Point3D worldCoord; worldCoord[0] = p1[0]+finalNormal[0]*m_MaxSpacing; worldCoord[1] = p1[1]+finalNormal[1]*m_MaxSpacing; worldCoord[2] = p1[2]+finalNormal[2]*m_MaxSpacing; double val = 0.0; mitk::ImagePixelReadAccessor readAccess(m_SegmentationBinaryImage); mitk::Index3D idx; m_SegmentationBinaryImage->GetGeometry()->WorldToIndex(worldCoord, idx); val = readAccess.GetPixelByIndexSafe(idx); if (val == 0.0) { ++m_PositiveNormalCounter; } else { ++m_NegativeNormalCounter; } } vertexNormalTemp[0] = vertexNormal[0]; vertexNormalTemp[1] = vertexNormal[1]; vertexNormalTemp[2] = vertexNormal[2]; vtkIdType id = cell[j+1]; normals->SetTuple(id,finalNormal); } existingPoints->GetPoint(cell[0], p1); existingPoints->GetPoint(cell[1], p2); v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; vtkMath::Cross(v1,polygonNormal,vertexNormal); vtkMath::Normalize(vertexNormal); vertexNormal[0] = (vertexNormal[0] + vertexNormalTemp[0])*0.5; vertexNormal[1] = (vertexNormal[1] + vertexNormalTemp[1])*0.5; vertexNormal[2] = (vertexNormal[2] + vertexNormalTemp[2])*0.5; vtkIdType id = cell[0]; normals->SetTuple(id,vertexNormal); id = cell[cellSize-1]; normals->SetTuple(id,vertexNormal); if(m_NegativeNormalCounter > m_PositiveNormalCounter) { for(vtkIdType n = 0; n < cellSize; n++) { double normal[3]; normals->GetTuple(offSet+n, normal); normal[0] = (-1)*normal[0]; normal[1] = (-1)*normal[1]; normal[2] = (-1)*normal[2]; normals->SetTuple(offSet+n, normal); } } m_NegativeNormalCounter = 0; m_PositiveNormalCounter = 0; offSet += cellSize; }//end for all cells Surface::Pointer surface = this->GetOutput(i); surface->GetVtkPolyData()->GetCellData()->SetNormals(normals); }//end for all inputs //Setting progressbar if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize); } mitk::Surface::Pointer mitk::ComputeContourSetNormalsFilter::GetNormalsAsSurface() { //Just for debugging: vtkSmartPointer newPolyData = vtkSmartPointer::New(); vtkSmartPointer newLines = vtkSmartPointer::New(); vtkSmartPointer newPoints = vtkSmartPointer::New(); unsigned int idCounter (0); //Debug end for (unsigned int i = 0; i < this->GetNumberOfIndexedOutputs(); i++) { Surface* currentSurface = const_cast( this->GetOutput(i) ); vtkPolyData* polyData = currentSurface->GetVtkPolyData(); vtkSmartPointer currentCellNormals = vtkDoubleArray::SafeDownCast(polyData->GetCellData()->GetNormals()); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType* cell (NULL); vtkIdType cellSize (0); for( existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { - for ( unsigned int j = 0; j < cellSize; j++ ) + for ( vtkIdType j = 0; j < cellSize; j++ ) { double currentNormal[3]; currentCellNormals->GetTuple(cell[j], currentNormal); vtkSmartPointer line = vtkSmartPointer::New(); line->GetPointIds()->SetNumberOfIds(2); double newPoint[3]; double p0[3]; existingPoints->GetPoint(cell[j], p0); newPoint[0] = p0[0] + currentNormal[0]; newPoint[1] = p0[1] + currentNormal[1]; newPoint[2] = p0[2] + currentNormal[2]; line->GetPointIds()->SetId(0, idCounter); newPoints->InsertPoint(idCounter, p0); idCounter++; line->GetPointIds()->SetId(1, idCounter); newPoints->InsertPoint(idCounter, newPoint); idCounter++; newLines->InsertNextCell(line); }//end for all points }//end for all cells }//end for all outputs newPolyData->SetPoints(newPoints); newPolyData->SetLines(newLines); newPolyData->BuildCells(); mitk::Surface::Pointer surface = mitk::Surface::New(); surface->SetVtkPolyData(newPolyData); return surface; } void mitk::ComputeContourSetNormalsFilter::SetMaxSpacing(double maxSpacing) { m_MaxSpacing = maxSpacing; } void mitk::ComputeContourSetNormalsFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } void mitk::ComputeContourSetNormalsFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { this->PopBackInput(); } this->SetNumberOfIndexedInputs(0); this->SetNumberOfIndexedOutputs(0); mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); } void mitk::ComputeContourSetNormalsFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::ComputeContourSetNormalsFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; } diff --git a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp index 522698602d..6d4a64d6c3 100644 --- a/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkCreateDistanceImageFromSurfaceFilter.cpp @@ -1,629 +1,631 @@ /*=================================================================== 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 "mitkCreateDistanceImageFromSurfaceFilter.h" mitk::CreateDistanceImageFromSurfaceFilter::CreateDistanceImageFromSurfaceFilter() { m_DistanceImageVolume = 50000; this->m_UseProgressBar = false; this->m_ProgressStepSize = 5; mitk::Image::Pointer output = mitk::Image::New(); this->SetNthOutput(0, output.GetPointer()); } mitk::CreateDistanceImageFromSurfaceFilter::~CreateDistanceImageFromSurfaceFilter() { } void mitk::CreateDistanceImageFromSurfaceFilter::GenerateData() { //First of all we have to build the equation-system from the existing contour-edge-points this->CreateSolutionMatrixAndFunctionValues(); //Then we solve the equation-system via QR - decomposition. The interpolation weights are obtained in that way vnl_qr solver (m_SolutionMatrix); m_Weights = solver.solve(m_FunctionValues); //Setting progressbar if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(2); //The last step is to create the distance map with the interpolated distance function this->CreateDistanceImage(); m_Centers.clear(); m_FunctionValues.clear(); m_Normals.clear(); m_Weights.clear(); m_SolutionMatrix.clear(); //Setting progressbar if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(3); } void mitk::CreateDistanceImageFromSurfaceFilter::CreateSolutionMatrixAndFunctionValues() { unsigned int numberOfInputs = this->GetNumberOfIndexedInputs(); if (numberOfInputs == 0) { MITK_ERROR << "mitk::CreateDistanceImageFromSurfaceFilter: No input available. Please set an input!" << std::endl; itkExceptionMacro("mitk::CreateDistanceImageFromSurfaceFilter: No input available. Please set an input!"); return; } //First of all we have to extract the nomals and the surface points. //Duplicated points can be eliminated Surface* currentSurface; vtkSmartPointer polyData; vtkSmartPointer currentCellNormals; vtkSmartPointer existingPolys; vtkSmartPointer existingPoints; double p[3]; PointType currentPoint; PointType normal; for (unsigned int i = 0; i < numberOfInputs; i++) { currentSurface = const_cast( this->GetInput(i) ); polyData = currentSurface->GetVtkPolyData(); if (polyData->GetNumberOfPolys() == 0) { MITK_INFO << "mitk::CreateDistanceImageFromSurfaceFilter: No input-polygons available. Please be sure the input surface consists of polygons!" << std::endl; } currentCellNormals = vtkDoubleArray::SafeDownCast(polyData->GetCellData()->GetNormals()); existingPolys = polyData->GetPolys(); existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType* cell (NULL); vtkIdType cellSize (0); for( existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { - for ( unsigned int j = 0; j < cellSize; j++ ) + for ( vtkIdType j = 0; j < cellSize; j++ ) { existingPoints->GetPoint(cell[j], p); currentPoint.copy_in(p); int count = std::count(m_Centers.begin() ,m_Centers.end(),currentPoint); if (count == 0) { double currentNormal[3]; currentCellNormals->GetTuple(cell[j], currentNormal); normal.copy_in(currentNormal); m_Normals.push_back(normal); m_Centers.push_back(currentPoint); } }//end for all points }//end for all cells }//end for all outputs //For we can now calculate the exact size of the centers we initialize the data structures unsigned int numberOfCenters = m_Centers.size(); m_Centers.reserve(numberOfCenters*3); m_FunctionValues.set_size(numberOfCenters*3); m_FunctionValues.fill(0); //Create inner points for (unsigned int i = 0; i < numberOfCenters; i++) { currentPoint = m_Centers.at(i); normal = m_Normals.at(i); currentPoint[0] = currentPoint[0] - normal[0]; currentPoint[1] = currentPoint[1] - normal[1]; currentPoint[2] = currentPoint[2] - normal[2]; m_Centers.push_back(currentPoint); m_FunctionValues.put(numberOfCenters+i, -1); } //Create outer points for (unsigned int i = 0; i < numberOfCenters; i++) { currentPoint = m_Centers.at(i); normal = m_Normals.at(i); currentPoint[0] = currentPoint[0] + normal[0]; currentPoint[1] = currentPoint[1] + normal[1]; currentPoint[2] = currentPoint[2] + normal[2]; m_Centers.push_back(currentPoint); m_FunctionValues.put(numberOfCenters*2+i, 1); } //Now we have created all centers and all function values. Next step is to create the solution matrix numberOfCenters = m_Centers.size(); m_SolutionMatrix.set_size(numberOfCenters, numberOfCenters); m_Weights.set_size(numberOfCenters); PointType p1; PointType p2; double norm; for (unsigned int i = 0; i < numberOfCenters; i++) { for (unsigned int j = 0; j < numberOfCenters; j++) { //Calculate the RBF value. Currently using Phi(r) = r with r is the euclidian distance between two points p1 = m_Centers.at(i); p2 = m_Centers.at(j); p1 = p1 - p2; norm = p1.two_norm(); m_SolutionMatrix(i,j) = norm; } } } void mitk::CreateDistanceImageFromSurfaceFilter::CreateDistanceImage() { DistanceImageType::Pointer distanceImg = DistanceImageType::New(); // Determine the bounds of the input points in index- and world-coordinates DistanceImageType::PointType minPointInWorldCoordinates, maxPointInWorldCoordinates; DistanceImageType::IndexType minPointInIndexCoordinates, maxPointInIndexCoordinates; DetermineBounds( minPointInWorldCoordinates, maxPointInWorldCoordinates, minPointInIndexCoordinates, maxPointInIndexCoordinates ); // Calculate the extent of the region that contains all given points in MM. // To do this, we take the difference between the maximal and minimal // index-coordinates (must not be less than 1) and multiply it with the // spacing of the reference-image. Vector3D extentMM; for (unsigned int dim = 0; dim < 3; ++dim) { extentMM[dim] = (int) ( (std::max( std::abs(maxPointInIndexCoordinates[dim] - minPointInIndexCoordinates[dim]), (DistanceImageType::IndexType::IndexValueType) 1 ) + 1.0) // (max-index - min-index)+1 because the pixels between index 3 and 5 cover 2+1=3 pixels (pixel 3,4, and 5) * m_ReferenceImage->GetSpacing()[dim] ) + 1; // (int) ((...) + 1) -> we round up to the next BIGGER int value } /* * Now create an empty distance image. The create image will always have the same sizeOfRegion, independent from * the original image (e.g. always consists of 500000 pixels) and will have an isotropic spacing. * The spacing is calculated like the following: * The image's volume = 500000 Pixels = extentX*spacing*extentY*spacing*extentZ*spacing * So the spacing is: spacing = ( 500000 / extentX*extentY*extentZ )^(1/3) */ double basis = (extentMM[0]*extentMM[1]*extentMM[2]) / m_DistanceImageVolume; double exponent = 1.0/3.0; double distImgSpacing = pow(basis, exponent); int tempSpacing = (distImgSpacing+0.05)*10; m_DistanceImageSpacing = (double)tempSpacing/10.0; // calculate the number of pixels of the distance image for each direction unsigned int numberOfXPixel = extentMM[0] / m_DistanceImageSpacing; unsigned int numberOfYPixel = extentMM[1] / m_DistanceImageSpacing; unsigned int numberOfZPixel = extentMM[2] / m_DistanceImageSpacing; // We increase the sizeOfRegion by 4 as we decrease the origin by 2 later. // This expansion of the region is necessary to achieve a complete // interpolation. DistanceImageType::SizeType sizeOfRegion; sizeOfRegion[0] = numberOfXPixel + 4; sizeOfRegion[1] = numberOfYPixel + 4; sizeOfRegion[2] = numberOfZPixel + 4; // The region starts at index 0,0,0 DistanceImageType::IndexType initialOriginAsIndex; initialOriginAsIndex.Fill(0); DistanceImageType::PointType originAsWorld = minPointInWorldCoordinates; DistanceImageType::RegionType lpRegion; lpRegion.SetSize(sizeOfRegion); lpRegion.SetIndex(initialOriginAsIndex); // We initialize the itk::Image with // * origin and direction to have it correctly placed and rotated in the world // * the largest possible region to set the extent to be calculated // * the isotropic spacing that we have calculated above distanceImg->SetOrigin( originAsWorld ); distanceImg->SetDirection( m_ReferenceImage->GetDirection() ); distanceImg->SetRegions( lpRegion ); distanceImg->SetSpacing( m_DistanceImageSpacing ); distanceImg->Allocate(); //First of all the image is initialized with the value 10 for each pixel distanceImg->FillBuffer(10); // Now we move the origin of the distanceImage 2 index-Coordinates // in all directions DistanceImageType::IndexType originAsIndex; distanceImg->TransformPhysicalPointToIndex( originAsWorld, originAsIndex ); originAsIndex[0] -= 2; originAsIndex[1] -= 2; originAsIndex[2] -= 2; distanceImg->TransformIndexToPhysicalPoint( originAsIndex, originAsWorld ); distanceImg->SetOrigin( originAsWorld ); /* * Now we must calculate the distance for each pixel. But instead of calculating the distance value * for all of the image's pixels we proceed similar to the region growing algorithm: * * 1. Take the first pixel from the narrowband_point_list and calculate the distance for each neighbor (6er) * 2. If the current index's distance value is below a certain threshold push it into the list * 3. Next iteration take the next index from the list and originAsIndex with 1. again * * This is done until the narrowband_point_list is empty. */ std::queue narrowbandPoints; PointType currentPoint = m_Centers.at(0); double distance = this->CalculateDistanceValue(currentPoint); // create itk::Point from vnl_vector DistanceImageType::PointType currentPointAsPoint; currentPointAsPoint[0] = currentPoint[0]; currentPointAsPoint[1] = currentPoint[1]; currentPointAsPoint[2] = currentPoint[2]; // Transform the input point in world-coordinates to index-coordinates DistanceImageType::IndexType currentIndex; distanceImg->TransformPhysicalPointToIndex( currentPointAsPoint, currentIndex ); assert( lpRegion.IsInside(currentIndex) ); // we are quite certain this should hold narrowbandPoints.push(currentIndex); distanceImg->SetPixel(currentIndex, distance); NeighborhoodImageIterator::RadiusType radius; radius.Fill(1); NeighborhoodImageIterator nIt(radius, distanceImg, distanceImg->GetLargestPossibleRegion()); unsigned int relativeNbIdx[] = {4, 10, 12, 14, 16, 22}; bool isInBounds = false; while ( !narrowbandPoints.empty() ) { nIt.SetLocation(narrowbandPoints.front()); narrowbandPoints.pop(); unsigned int* relativeNb = &relativeNbIdx[0]; for (int i = 0; i < 6; i++) { nIt.GetPixel(*relativeNb, isInBounds); if( isInBounds && nIt.GetPixel(*relativeNb) == 10) { currentIndex = nIt.GetIndex(*relativeNb); // Transform the currently checked point from index-coordinates to // world-coordinates distanceImg->TransformIndexToPhysicalPoint( currentIndex, currentPointAsPoint ); // create a vnl_vector currentPoint[0] = currentPointAsPoint[0]; currentPoint[1] = currentPointAsPoint[1]; currentPoint[2] = currentPointAsPoint[2]; // and check the distance distance = this->CalculateDistanceValue(currentPoint); if ( abs(distance) <= m_DistanceImageSpacing ) { nIt.SetPixel(*relativeNb, distance); narrowbandPoints.push(currentIndex); } } relativeNb++; } } ImageIterator imgRegionIterator (distanceImg, distanceImg->GetLargestPossibleRegion()); imgRegionIterator.GoToBegin(); double prevPixelVal = 1; - unsigned int _size[3] = { (unsigned int)(sizeOfRegion[0] - 1), (unsigned int)(sizeOfRegion[1] - 1), (unsigned int)(sizeOfRegion[2] - 1) }; + DistanceImageType::IndexType _size; + _size.Fill(-1); + _size += sizeOfRegion; double center [3] = {_size[0]/2.0, _size[1]/2.0, _size[2]/2.0}; MITK_INFO<<"Size: ["<<_size[0]<<","<<_size[1]<<","<<_size[2]<<"] Center: ["<GetOutput(); // Cast the created distance-Image from itk::Image to the mitk::Image // that is our output. CastToMitkImage(distanceImg, resultImage); } void mitk::CreateDistanceImageFromSurfaceFilter::FillImageRegion(DistanceImageType::RegionType reqRegion, DistanceImageType::PixelType pixelValue, DistanceImageType::Pointer image) { image->SetRequestedRegion(reqRegion); ImageIterator it (image, image->GetRequestedRegion()); while (!it.IsAtEnd()) { it.Set(pixelValue); ++it; } } double mitk::CreateDistanceImageFromSurfaceFilter::CalculateDistanceValue(PointType p) { double distanceValue (0); PointType p1; PointType p2; double norm; CenterList::iterator centerIter; InterpolationWeights::iterator weightsIter; for ( centerIter=m_Centers.begin(), weightsIter=m_Weights.begin(); centerIter!=m_Centers.end() && weightsIter!=m_Weights.end(); centerIter++, weightsIter++ ) { p1 = *centerIter; p2 = p-p1; norm = p2.two_norm(); distanceValue = distanceValue + norm* (*weightsIter); } return distanceValue; } void mitk::CreateDistanceImageFromSurfaceFilter::GenerateOutputInformation() { } void mitk::CreateDistanceImageFromSurfaceFilter::PrintEquationSystem() { std::ofstream esfile; esfile.open("C:/Users/fetzer/Desktop/equationSystem/es.txt"); esfile<<"Nummber of rows: "<SetInput( 0, const_cast( surface ) ); } void mitk::CreateDistanceImageFromSurfaceFilter::SetInput( unsigned int idx, const mitk::Surface* surface ) { if ( this->GetInput(idx) != surface ) { this->SetNthInput( idx, const_cast( surface ) ); this->Modified(); } } const mitk::Surface* mitk::CreateDistanceImageFromSurfaceFilter::GetInput() { if (this->GetNumberOfIndexedInputs() < 1) return NULL; return static_cast(this->ProcessObject::GetInput(0)); } const mitk::Surface* mitk::CreateDistanceImageFromSurfaceFilter::GetInput( unsigned int idx) { if (this->GetNumberOfIndexedInputs() < 1) return NULL; return static_cast(this->ProcessObject::GetInput(idx)); } void mitk::CreateDistanceImageFromSurfaceFilter::RemoveInputs(mitk::Surface* input) { DataObjectPointerArraySizeType nb = this->GetNumberOfIndexedInputs(); for(DataObjectPointerArraySizeType i = 0; i < nb; i++) { if( this->GetInput(i) == input ) { this->RemoveInput(i); return; } } } void mitk::CreateDistanceImageFromSurfaceFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { this->PopBackInput(); } this->SetNumberOfIndexedInputs(0); this->SetNumberOfIndexedOutputs(1); mitk::Image::Pointer output = mitk::Image::New(); this->SetNthOutput(0, output.GetPointer()); } void mitk::CreateDistanceImageFromSurfaceFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::CreateDistanceImageFromSurfaceFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; } void mitk::CreateDistanceImageFromSurfaceFilter::SetReferenceImage( itk::ImageBase<3>::Pointer referenceImage ) { m_ReferenceImage = referenceImage; } void mitk::CreateDistanceImageFromSurfaceFilter::DetermineBounds( DistanceImageType::PointType &minPointInWorldCoordinates, DistanceImageType::PointType &maxPointInWorldCoordinates, DistanceImageType::IndexType &minPointInIndexCoordinates, DistanceImageType::IndexType &maxPointInIndexCoordinates ) { PointType firstCenter = m_Centers.at(0); DistanceImageType::PointType tmpPoint; tmpPoint[0] = firstCenter[0]; tmpPoint[1] = firstCenter[1]; tmpPoint[2] = firstCenter[2]; // transform the first point from world-coordinates to index-coordinates DistanceImageType::IndexType tmpIndex; m_ReferenceImage->TransformPhysicalPointToIndex( tmpPoint, tmpIndex ); // initialize the variables with this first point int xmin = tmpIndex[0]; int ymin = tmpIndex[1]; int zmin = tmpIndex[2]; int xmax = tmpIndex[0]; int ymax = tmpIndex[1]; int zmax = tmpIndex[2]; // iterate over the rest of the points CenterList::iterator centerIter = m_Centers.begin(); for ( ++centerIter; centerIter!=m_Centers.end(); centerIter++) { tmpPoint[0] = (*centerIter)[0]; tmpPoint[1] = (*centerIter)[1]; tmpPoint[2] = (*centerIter)[2]; // transform each point from world-coordinates to index-coordinates m_ReferenceImage->TransformPhysicalPointToIndex( tmpPoint, tmpIndex ); // and set the variables accordingly to find the minimum // and maximum in all directions in index-coordinates if (xmin > tmpIndex[0]) { xmin = tmpIndex[0]; } if (ymin > tmpIndex[1]) { ymin = tmpIndex[1]; } if (zmin > tmpIndex[2]) { zmin = tmpIndex[2]; } if (xmax < tmpIndex[0]) { xmax = tmpIndex[0]; } if (ymax < tmpIndex[1]) { ymax = tmpIndex[1]; } if (zmax < tmpIndex[2]) { zmax = tmpIndex[2]; } } // put the found coordinates into Index-Points minPointInIndexCoordinates[0] = xmin; minPointInIndexCoordinates[1] = ymin; minPointInIndexCoordinates[2] = zmin; maxPointInIndexCoordinates[0] = xmax; maxPointInIndexCoordinates[1] = ymax; maxPointInIndexCoordinates[2] = zmax; // and transform them into world-coordinates m_ReferenceImage->TransformIndexToPhysicalPoint( minPointInIndexCoordinates, minPointInWorldCoordinates ); m_ReferenceImage->TransformIndexToPhysicalPoint( maxPointInIndexCoordinates, maxPointInWorldCoordinates ); } diff --git a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp index 6fc087e02c..487bbf7af5 100644 --- a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp @@ -1,502 +1,502 @@ /*=================================================================== 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 "mitkReduceContourSetFilter.h" mitk::ReduceContourSetFilter::ReduceContourSetFilter() { m_MaxSegmentLenght = 0; m_StepSize = 10; m_Tolerance = -1; m_ReductionType = DOUGLAS_PEUCKER; m_MaxSpacing = -1; m_MinSpacing = -1; this->m_UseProgressBar = false; this->m_ProgressStepSize = 1; m_NumberOfPointsAfterReduction = 0; mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); } mitk::ReduceContourSetFilter::~ReduceContourSetFilter() { } void mitk::ReduceContourSetFilter::GenerateData() { unsigned int numberOfInputs = this->GetNumberOfIndexedInputs(); unsigned int numberOfOutputs (0); vtkSmartPointer newPolyData; vtkSmartPointer newPolygons; vtkSmartPointer newPoints; //For the purpose of evaluation // unsigned int numberOfPointsBefore (0); m_NumberOfPointsAfterReduction=0; for(unsigned int i = 0; i < numberOfInputs; i++) { mitk::Surface* currentSurface = const_cast( this->GetInput(i) ); vtkSmartPointer polyData = currentSurface->GetVtkPolyData(); newPolyData = vtkPolyData::New(); newPolygons = vtkCellArray::New(); newPoints = vtkPoints::New(); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType* cell (NULL); vtkIdType cellSize (0); for( existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { bool incorporatePolygon = this->CheckForIntersection(cell,cellSize,existingPoints, /*numberOfIntersections, intersectionPoints, */i); if ( !incorporatePolygon ) continue; vtkSmartPointer newPolygon = vtkSmartPointer::New(); if(m_ReductionType == NTH_POINT) { this->ReduceNumberOfPointsByNthPoint(cellSize, cell, existingPoints, newPolygon, newPoints); if (newPolygon->GetPointIds()->GetNumberOfIds() != 0) { newPolygons->InsertNextCell(newPolygon); } } else if (m_ReductionType == DOUGLAS_PEUCKER) { this->ReduceNumberOfPointsByDouglasPeucker(cellSize, cell, existingPoints, newPolygon, newPoints); if (newPolygon->GetPointIds()->GetNumberOfIds() > 3) { newPolygons->InsertNextCell(newPolygon); } } //Again for evaluation // numberOfPointsBefore += cellSize; m_NumberOfPointsAfterReduction += newPolygon->GetPointIds()->GetNumberOfIds(); } if (newPolygons->GetNumberOfCells() != 0) { newPolyData->SetPolys(newPolygons); newPolyData->SetPoints(newPoints); newPolyData->BuildLinks(); this->SetNumberOfIndexedOutputs(numberOfOutputs + 1); mitk::Surface::Pointer surface = mitk::Surface::New(); this->SetNthOutput(numberOfOutputs, surface.GetPointer()); surface->SetVtkPolyData(newPolyData); numberOfOutputs++; } } // MITK_INFO<<"Points before: "<SetNumberOfIndexedOutputs(numberOfOutputs); //Setting progressbar if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize); } void mitk::ReduceContourSetFilter::ReduceNumberOfPointsByNthPoint (vtkIdType cellSize, vtkIdType* cell, vtkPoints* points, vtkPolygon* reducedPolygon, vtkPoints* reducedPoints) { unsigned int newNumberOfPoints (0); unsigned int mod = cellSize%m_StepSize; if(mod == 0) { newNumberOfPoints = cellSize/m_StepSize; } else { newNumberOfPoints = ( (cellSize-mod)/m_StepSize )+1; } if (newNumberOfPoints <= 3) { return; } reducedPolygon->GetPointIds()->SetNumberOfIds(newNumberOfPoints); reducedPolygon->GetPoints()->SetNumberOfPoints(newNumberOfPoints); - for (unsigned int i = 0; i < cellSize; i++) + for (vtkIdType i = 0; i < cellSize; i++) { if (i%m_StepSize == 0) { double point[3]; points->GetPoint(cell[i], point); vtkIdType id = reducedPoints->InsertNextPoint(point); reducedPolygon->GetPointIds()->SetId(i/m_StepSize, id); } } vtkIdType id = cell[0]; double point[3]; points->GetPoint(id, point); id = reducedPoints->InsertNextPoint(point); reducedPolygon->GetPointIds()->SetId(newNumberOfPoints-1, id); } void mitk::ReduceContourSetFilter::ReduceNumberOfPointsByDouglasPeucker(vtkIdType cellSize, vtkIdType* cell, vtkPoints* points, vtkPolygon* reducedPolygon, vtkPoints* reducedPoints) { //If the cell is too small to obtain a reduced polygon with the given stepsize return - if (cellSize <= m_StepSize*3)return; + if (cellSize <= static_cast(m_StepSize*3))return; /* What we do now is (see the Douglas Peucker Algorithm): 1. Divide the current contour in two line segments (start - middle; middle - end), put them into the stack 2. Fetch first line segment and create the following vectors: - v1 = (start;end) - v2 = (start;currentPoint) -> for each point of the current line segment! 3. Calculate the distance from the currentPoint to v1: a. Determine the length of the orthogonal projection of v2 to v1 by: l = v2 * (normalized v1) b. There a three possibilities for the distance then: d = sqrt(lenght(v2)^2 - l^2) if l > 0 and l < length(v1) d = lenght(v2-v1) if l > 0 and l > lenght(v1) d = length(v2) if l < 0 because v2 is then pointing in a different direction than v1 4. Memorize the point with the biggest distance and create two new line segments with it at the end of the iteration and put it into the stack 5. If the distance value D <= m_Tolerance, then add the start and end index and the corresponding points to the reduced ones */ //First of all set tolerance if none is specified if(m_Tolerance < 0) { if(m_MaxSpacing > 0) { m_Tolerance = m_MinSpacing; } else { m_Tolerance = 1.5; } } std::stack lineSegments; //1. Divide in line segments LineSegment ls2; ls2.StartIndex = cell[cellSize/2]; ls2.EndIndex = cell[cellSize-1]; lineSegments.push(ls2); LineSegment ls1; ls1.StartIndex = cell[0]; ls1.EndIndex = cell[cellSize/2]; lineSegments.push(ls1); LineSegment currentSegment; double v1[3]; double v2[3]; double tempV[3]; double lenghtV1; double currentMaxDistance (0); vtkIdType currentMaxDistanceIndex (0); double l; double d; vtkIdType pointId (0); //Add the start index to the reduced points. From now on just the end indices will be added pointId = reducedPoints->InsertNextPoint(points->GetPoint(cell[0])); reducedPolygon->GetPointIds()->InsertNextId(pointId); while (!lineSegments.empty()) { currentSegment = lineSegments.top(); lineSegments.pop(); //2. Create vectors points->GetPoint(currentSegment.EndIndex, tempV); points->GetPoint(currentSegment.StartIndex, v1); v1[0] = tempV[0]-v1[0]; v1[1] = tempV[1]-v1[1]; v1[2] = tempV[2]-v1[2]; lenghtV1 = vtkMath::Norm(v1); vtkMath::Normalize(v1); int range = currentSegment.EndIndex - currentSegment.StartIndex; for (int i = 1; i < abs(range); ++i) { points->GetPoint(currentSegment.StartIndex+i, tempV); points->GetPoint(currentSegment.StartIndex, v2); v2[0] = tempV[0]-v2[0]; v2[1] = tempV[1]-v2[1]; v2[2] = tempV[2]-v2[2]; //3. Calculate the distance l = vtkMath::Dot(v2, v1); d = vtkMath::Norm(v2); if (l > 0 && l < lenghtV1) { d = sqrt((d*d-l*l)); } else if (l > 0 && l > lenghtV1) { tempV[0] = lenghtV1*v1[0] - v2[0]; tempV[1] = lenghtV1*v1[1] - v2[1]; tempV[2] = lenghtV1*v1[2] - v2[2]; d = vtkMath::Norm(tempV); } //4. Memorize maximum distance if (d > currentMaxDistance) { currentMaxDistance = d; currentMaxDistanceIndex = currentSegment.StartIndex+i; } } //4. & 5. if (currentMaxDistance <= m_Tolerance) { //double temp[3]; int segmentLenght = currentSegment.EndIndex - currentSegment.StartIndex; if (segmentLenght > (int)m_MaxSegmentLenght) { m_MaxSegmentLenght = (unsigned int)segmentLenght; } // MITK_INFO<<"Lenght: "< 25) { unsigned int newLenght(segmentLenght); while (newLenght > 25) { newLenght = newLenght*0.5; } unsigned int divisions = abs(segmentLenght)/newLenght; // MITK_INFO<<"Divisions: "<InsertNextPoint(points->GetPoint(currentSegment.StartIndex + newLenght*i)); reducedPolygon->GetPointIds()->InsertNextId(pointId); } } // MITK_INFO<<"Inserting END: "<InsertNextPoint(points->GetPoint(currentSegment.EndIndex)); reducedPolygon->GetPointIds()->InsertNextId(pointId); } else { ls2.StartIndex = currentMaxDistanceIndex; ls2.EndIndex = currentSegment.EndIndex; lineSegments.push(ls2); ls1.StartIndex = currentSegment.StartIndex; ls1.EndIndex = currentMaxDistanceIndex; lineSegments.push(ls1); } currentMaxDistance = 0; } } bool mitk::ReduceContourSetFilter::CheckForIntersection (vtkIdType* currentCell, vtkIdType currentCellSize, vtkPoints* currentPoints,/* vtkIdType numberOfIntersections, vtkIdType* intersectionPoints,*/ unsigned int currentInputIndex) { /* If we check the current cell for intersections then we have to consider three possibilies: 1. There is another cell among all the other input surfaces which intersects the current polygon: - That means we have to save the intersection points because these points should not be eliminated 2. There current polygon exists just because of an intersection of another polygon with the current plane defined by the current polygon - That means the current polygon should not be incorporated and all of its points should be eliminated 3. There is no intersection - That mean we can just reduce the current polygons points without considering any intersections */ for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { //Don't check for intersection with the polygon itself if (i == currentInputIndex) continue; //Get the next polydata to check for intersection vtkSmartPointer poly = const_cast( this->GetInput(i) )->GetVtkPolyData(); vtkSmartPointer polygonArray = poly->GetPolys(); polygonArray->InitTraversal(); vtkIdType anotherInputPolygonSize (0); vtkIdType* anotherInputPolygonIDs(NULL); /* The procedure is: - Create the equation of the plane, defined by the points of next input - Calculate the distance of each point of the current polygon to the plane - If the maximum distance is not bigger than 1.5 of the maximum spacing AND the minimal distance is not bigger than 0.5 of the minimum spacing then the current contour is an intersection contour */ for( polygonArray->InitTraversal(); polygonArray->GetNextCell(anotherInputPolygonSize, anotherInputPolygonIDs);) { //Choosing three plane points to calculate the plane vectors double p1[3]; double p2[3]; double p3[3]; //The plane vectors double v1[3]; double v2[3] = { 0 }; //The plane normal double normal[3]; //Create first Vector poly->GetPoint(anotherInputPolygonIDs[0], p1); poly->GetPoint(anotherInputPolygonIDs[1], p2); v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; //Find 3rd point for 2nd vector (The angle between the two plane vectors should be bigger than 30 degrees) double maxDistance (0); double minDistance (10000); - for (unsigned int j = 2; j < anotherInputPolygonSize; j++) + for (vtkIdType j = 2; j < anotherInputPolygonSize; j++) { poly->GetPoint(anotherInputPolygonIDs[j], p3); v2[0] = p3[0]-p1[0]; v2[1] = p3[1]-p1[1]; v2[2] = p3[2]-p1[2]; //Calculate the angle between the two vector for the current point double dotV1V2 = vtkMath::Dot(v1,v2); double absV1 = sqrt(vtkMath::Dot(v1,v1)); double absV2 = sqrt(vtkMath::Dot(v2,v2)); double cosV1V2 = dotV1V2/(absV1*absV2); double arccos = acos(cosV1V2); double degree = vtkMath::DegreesFromRadians(arccos); //If angle is bigger than 30 degrees break if (degree > 30) break; }//for (to find 3rd point) //Calculate normal of the plane by taking the cross product of the two vectors vtkMath::Cross(v1,v2,normal); vtkMath::Normalize(normal); //Determine position of the plane double lambda = vtkMath::Dot(normal, p1); /* Calculate the distance to the plane for each point of the current polygon If the distance is zero then save the currentPoint as intersection point */ - for (unsigned int k = 0; k < currentCellSize; k++) + for (vtkIdType k = 0; k < currentCellSize; k++) { double currentPoint[3]; currentPoints->GetPoint(currentCell[k], currentPoint); double tempPoint[3]; tempPoint[0] = normal[0]*currentPoint[0]; tempPoint[1] = normal[1]*currentPoint[1]; tempPoint[2] = normal[2]*currentPoint[2]; double temp = tempPoint[0]+tempPoint[1]+tempPoint[2]-lambda; double distance = fabs(temp); if (distance > maxDistance) { maxDistance = distance; } if (distance < minDistance) { minDistance = distance; } }//for (to calculate distance and intersections with currentPolygon) if (maxDistance < 1.5*m_MaxSpacing && minDistance < 0.5*m_MinSpacing) { return false; } //Because we are considering the plane defined by the acual input polygon only one iteration is sufficient //We do not need to consider each cell of the plane break; }//for (to traverse through all cells of actualInputPolyData) }//for (to iterate through all inputs) return true; } void mitk::ReduceContourSetFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } void mitk::ReduceContourSetFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { this->PopBackInput(); } this->SetNumberOfIndexedInputs(0); this->SetNumberOfIndexedOutputs(0); mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); m_NumberOfPointsAfterReduction = 0; } void mitk::ReduceContourSetFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::ReduceContourSetFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; }