diff --git a/Modules/Core/include/mitkDisplayInteractor.h b/Modules/Core/include/mitkDisplayInteractor.h index 78eb244f1e..b9c170d102 100644 --- a/Modules/Core/include/mitkDisplayInteractor.h +++ b/Modules/Core/include/mitkDisplayInteractor.h @@ -1,267 +1,281 @@ /*=================================================================== 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 mitkDisplayInteractor_h #define mitkDisplayInteractor_h #include #include "mitkInteractionEventObserver.h" namespace mitk { /** *\class DisplayInteractor *@brief Observer that manages the interaction with the display. * * This includes the interaction of Zooming, Panning, Scrolling and adjusting the LevelWindow. * * @ingroup Interaction **/ /** * Inherits from mitk::InteractionEventObserver since it doesn't alter any data (only their representation), * and its actions cannot be associated with a DataNode. Also inherits from EventStateMachine */ class MITKCORE_EXPORT DisplayInteractor: public EventStateMachine, public InteractionEventObserver { public: mitkClassMacro(DisplayInteractor, EventStateMachine) itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** * By this function the Observer gets notified about new events. * Here it is adapted to pass the events to the state machine in order to use * its infrastructure. * It also checks if event is to be accepted when it already has been processed by a DataInteractor. */ virtual void Notify(InteractionEvent* interactionEvent, bool isHandled) override; protected: DisplayInteractor(); virtual ~DisplayInteractor(); /** * Derived function. * Connects the action names used in the state machine pattern with functions implemented within * this InteractionEventObserver. This is only necessary here because the events are processed by the state machine. */ void ConnectActionsAndFunctions() override; /** * Derived function. * Is executed when config object is set / changed. * Here it is used to read out the parameters set in the configuration file, * and set the member variables accordingly. */ virtual void ConfigurationChanged() override; /** * Derived function. * Is executed when config object is set / changed. * Here it is used to read out the parameters set in the configuration file, * and set the member variables accordingly. */ virtual bool FilterEvents(InteractionEvent* interactionEvent, DataNode* dataNode) override; virtual bool CheckPositionEvent( const InteractionEvent* interactionEvent ); virtual bool CheckRotationPossible( const InteractionEvent* interactionEvent ); virtual bool CheckSwivelPossible( const InteractionEvent* interactionEvent ); /** * \brief Initializes an interaction, saves the pointers start position for further reference. */ virtual void Init(StateMachineAction*, InteractionEvent*); /** * \brief Performs panning of the data set in the render window. */ virtual void Move(StateMachineAction*, InteractionEvent*); /** * \brief Sets crosshair at clicked position* */ virtual void SetCrosshair(StateMachineAction*, InteractionEvent*); /** * \brief Performs zooming relative to mouse/pointer movement. * * Behavior is determined by \see m_ZoomDirection and \see m_ZoomFactor. * */ virtual void Zoom(StateMachineAction*, InteractionEvent*); /** * \brief Performs scrolling relative to mouse/pointer movement. * * Behavior is determined by \see m_ScrollDirection and \see m_AutoRepeat. * */ virtual void Scroll(StateMachineAction*, InteractionEvent*); /** * \brief Scrolls one layer up */ virtual void ScrollOneDown(StateMachineAction*, InteractionEvent*); /** * \brief Scrolls one layer down */ virtual void ScrollOneUp(StateMachineAction*, InteractionEvent*); /** * \brief Adjusts the level windows relative to mouse/pointer movement. */ virtual void AdjustLevelWindow(StateMachineAction*, InteractionEvent*); /** * \brief Starts crosshair rotation */ virtual void StartRotation(StateMachineAction*, InteractionEvent*); /** * \brief Ends crosshair rotation */ virtual void EndRotation(StateMachineAction*, InteractionEvent*); /** * \brief */ virtual void Rotate(StateMachineAction*, InteractionEvent* event); virtual void Swivel(StateMachineAction*, InteractionEvent* event); /** * \brief Updates the Statusbar information with the information about the clicked position */ - virtual void UpdateStatusbar(StateMachineAction*, InteractionEvent* event); + virtual void SetStatusBarInfoFromClickedPosition(StateMachineAction*, InteractionEvent* event); /** * \brief Method to retrieve bool-value for given property from string-property * in given propertylist. */ bool GetBoolProperty( mitk::PropertyList::Pointer propertyList, const char* propertyName, bool defaultValue ); // Typedefs typedef std::vector SNCVector; private: mitk::DataNode::Pointer GetTopLayerNode(mitk::DataStorage::SetOfObjects::ConstPointer nodes,mitk::Point3D worldposition, BaseRenderer* ren); + /** + * @brief UpdateStatusBar + * @param image3D + * @param idx + * @param time + * @param component If the PixelType of image3D is a vector (for example a 2D velocity vector), then only one of the vector components can be + * displayed at once. Setting this parameter will determine which of the vector's components will be used to determine the displayed PixelValue. + * Set this to 0 for scalar images + */ + void UpdateStatusBar(itk::SmartPointer image3D, itk::Index<3> idx, TimeStepType time=0, int component=0); + + /** * \brief Coordinate of the pointer at begin of an interaction translated to mm unit */ mitk::Point2D m_StartCoordinateInMM; /** * \brief Coordinate of the pointer in the last step within an interaction. */ mitk::Point2D m_LastDisplayCoordinate; /** * \brief Coordinate of the pointer in the last step within an interaction translated to mm unit */ mitk::Point2D m_LastCoordinateInMM; /** * \brief Current coordinates of the pointer. */ mitk::Point2D m_CurrentDisplayCoordinate; + + /** * \brief Modifier that defines how many slices are scrolled per pixel that the mouse has moved * * This modifier defines how many slices the scene is scrolled per pixel that the mouse cursor has moved. * By default the modifier is 4. This means that when the user moves the cursor by 4 pixels in Y-direction * the scene is scrolled by one slice. If the user has moved the the cursor by 20 pixels, the scene is * scrolled by 5 slices. * * If the cursor has moved less than m_IndexToSliceModifier pixels the scene is scrolled by one slice. */ int m_IndexToSliceModifier; /** Defines behavior at end of data set. * If set to true it will restart at end of data set from the beginning. */ bool m_AutoRepeat; /** * Defines scroll behavior. * Default is up/down movement of pointer performs scrolling */ std::string m_ScrollDirection; /** * Defines how the axis of interaction influences scroll behavior. */ bool m_InvertScrollDirection; /** * Defines scroll behavior. * Default is up/down movement of pointer performs zooming */ std::string m_ZoomDirection; /** * Defines how the axis of interaction influences zoom behavior. */ bool m_InvertZoomDirection; /** * Defines how the axis of interaction influences move behavior. */ bool m_InvertMoveDirection; /** * Defines level/window behavior. * Default is left/right movement of pointer modifies the level. */ std::string m_LevelDirection; /** * Defines how the axis of interaction influences level/window behavior. */ bool m_InvertLevelWindowDirection; /** * Determines if the Observer reacts to events that already have been processed by a DataInteractor. * The default value is false. */ bool m_AlwaysReact; /** * Factor to adjust zooming speed. */ float m_ZoomFactor; ///// Members to deal with rotating slices /** * @brief m_LinkPlanes Determines if angle between crosshair remains fixed when rotating */ bool m_LinkPlanes; SNCVector m_RotatableSNCs; /// all SNCs that currently have CreatedWorldGeometries, that can be rotated. SNCVector m_SNCsToBeRotated; /// all SNCs that will be rotated (exceptions are the ones parallel to the one being clicked) Point3D m_LastCursorPosition; /// used for calculation of the rotation angle Point3D m_CenterOfRotation; /// used for calculation of the rotation angle Point2D m_ReferenceCursor; Vector3D m_RotationPlaneNormal; Vector3D m_RotationPlaneXVector; Vector3D m_RotationPlaneYVector; Vector3D m_PreviousRotationAxis; ScalarType m_PreviousRotationAngle; }; } #endif diff --git a/Modules/Core/src/Interactions/mitkDisplayInteractor.cpp b/Modules/Core/src/Interactions/mitkDisplayInteractor.cpp index f994e33ccf..2257483832 100644 --- a/Modules/Core/src/Interactions/mitkDisplayInteractor.cpp +++ b/Modules/Core/src/Interactions/mitkDisplayInteractor.cpp @@ -1,944 +1,996 @@ /*=================================================================== 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 "mitkDisplayInteractor.h" #include "mitkBaseRenderer.h" #include "mitkInteractionPositionEvent.h" #include "mitkCameraController.h" #include "mitkPropertyList.h" #include #include #include // level window #include "mitkStandaloneDataStorage.h" #include "mitkNodePredicateDataType.h" #include "mitkLevelWindowProperty.h" #include "mitkLevelWindow.h" #include "vtkRenderWindowInteractor.h" #include "mitkLine.h" // Rotation #include #include #include "rotate_cursor.xpm" #include "mitkInteractionConst.h" // #include "mitkStatusBar.h" #include "mitkImage.h" #include "mitkPixelTypeMultiplex.h" #include "mitkImagePixelReadAccessor.h" void mitk::DisplayInteractor::Notify(InteractionEvent* interactionEvent, bool isHandled) { // to use the state machine pattern, // the event is passed to the state machine interface to be handled if (!isHandled || m_AlwaysReact) { this->HandleEvent(interactionEvent, NULL); } } void mitk::DisplayInteractor::ConnectActionsAndFunctions() { CONNECT_CONDITION( "check_position_event", CheckPositionEvent ); CONNECT_CONDITION( "check_can_rotate", CheckRotationPossible ); CONNECT_CONDITION( "check_can_swivel", CheckSwivelPossible ); CONNECT_FUNCTION("init", Init); CONNECT_FUNCTION("move", Move); CONNECT_FUNCTION("zoom", Zoom); CONNECT_FUNCTION("scroll", Scroll); CONNECT_FUNCTION("ScrollOneDown", ScrollOneDown); CONNECT_FUNCTION("ScrollOneUp", ScrollOneUp); CONNECT_FUNCTION("levelWindow", AdjustLevelWindow); CONNECT_FUNCTION("setCrosshair", SetCrosshair); - CONNECT_FUNCTION ("updateStatusbar", UpdateStatusbar) + CONNECT_FUNCTION ("updateStatusbar", SetStatusBarInfoFromClickedPosition) CONNECT_FUNCTION("startRotation", StartRotation); CONNECT_FUNCTION("endRotation", EndRotation); CONNECT_FUNCTION("rotate", Rotate); CONNECT_FUNCTION("swivel", Swivel); } mitk::DisplayInteractor::DisplayInteractor() : m_IndexToSliceModifier(4) , m_AutoRepeat(false) , m_InvertScrollDirection( false ) , m_InvertZoomDirection( false ) , m_InvertMoveDirection( false ) , m_InvertLevelWindowDirection( false ) , m_AlwaysReact(false) , m_ZoomFactor(2) , m_LinkPlanes(true) { m_StartCoordinateInMM.Fill(0); m_LastDisplayCoordinate.Fill(0); m_LastCoordinateInMM.Fill(0); m_CurrentDisplayCoordinate.Fill(0); } mitk::DisplayInteractor::~DisplayInteractor() { } bool mitk::DisplayInteractor::CheckPositionEvent( const InteractionEvent* interactionEvent ) { const InteractionPositionEvent* positionEvent = dynamic_cast(interactionEvent); if (positionEvent == NULL) { return false; } return true; } bool mitk::DisplayInteractor::CheckRotationPossible(const mitk::InteractionEvent *interactionEvent) { // Decide between moving and rotation slices. /* Detailed logic: 1. Find the SliceNavigationController that has sent the event: this one defines our rendering plane and will NOT be rotated. Needs not even be counted or checked. 2. Inspect every other SliceNavigationController - calculate the line intersection of this SliceNavigationController's plane with our rendering plane - if there is NO interesection, ignore and continue - IF there is an intersection - check the mouse cursor's distance from that line. 0. if the line is NOT near the cursor, remember the plane as "one of the other planes" (which can be rotated in "locked" mode) 1. on first line near the cursor, just remember this intersection line as THE other plane that we want to rotate 2. on every consecutive line near the cursor, check if the line is geometrically identical to the line that we want to rotate - if yes, we just push this line to the "other" lines and rotate it along - if no, then we have a situation where the mouse is near two other lines (e.g. crossing point) and don't want to rotate */ const InteractionPositionEvent* posEvent = dynamic_cast(interactionEvent); if (posEvent == nullptr) return false; BaseRenderer* clickedRenderer = posEvent->GetSender(); const PlaneGeometry* ourViewportGeometry = (clickedRenderer->GetCurrentWorldPlaneGeometry()); if (!ourViewportGeometry) return false; Point3D cursorPosition = posEvent->GetPositionInWorld(); const PlaneGeometry* geometryToBeRotated = NULL; // this one is under the mouse cursor const PlaneGeometry* anyOtherGeometry = NULL; // this is also visible (for calculation of intersection ONLY) Line3D intersectionLineWithGeometryToBeRotated; bool hitMultipleLines(false); m_SNCsToBeRotated.clear(); const double threshholdDistancePixels = 12.0; auto renWindows = interactionEvent->GetSender()->GetRenderingManager()->GetAllRegisteredRenderWindows(); for(auto renWin : renWindows) { SliceNavigationController* snc = BaseRenderer::GetInstance(renWin)->GetSliceNavigationController(); // If the mouse cursor is in 3D Renderwindow, do not check for intersecting planes. if ( BaseRenderer::GetInstance(renWin)->GetMapperID() == BaseRenderer::Standard3D) continue; const PlaneGeometry* otherRenderersRenderPlane = snc->GetCurrentPlaneGeometry(); if (otherRenderersRenderPlane == NULL) continue; // ignore, we don't see a plane // check if there is an intersection Line3D intersectionLine; // between rendered/clicked geometry and the one being analyzed if (!ourViewportGeometry->IntersectionLine(otherRenderersRenderPlane, intersectionLine)) { continue; // we ignore this plane, it's parallel to our plane } // check distance from intersection line double distanceFromIntersectionLine = intersectionLine.Distance(cursorPosition); // far away line, only remember for linked rotation if necessary if (distanceFromIntersectionLine > threshholdDistancePixels) { anyOtherGeometry = otherRenderersRenderPlane; // we just take the last one, so overwrite each iteration (we just need some crossing point) // TODO what about multiple crossings? NOW we have undefined behavior / random crossing point is used if (m_LinkPlanes) { m_SNCsToBeRotated.push_back(snc); } } else // close to cursor { if (geometryToBeRotated == NULL) // first one close to the cursor { geometryToBeRotated = otherRenderersRenderPlane; intersectionLineWithGeometryToBeRotated = intersectionLine; m_SNCsToBeRotated.push_back(snc); } else { // compare to the line defined by geometryToBeRotated: if identical, just rotate this otherRenderersRenderPlane together with the primary one // if different, DON'T rotate if (intersectionLine.IsParallel(intersectionLineWithGeometryToBeRotated) && intersectionLine.Distance(intersectionLineWithGeometryToBeRotated.GetPoint1()) < mitk::eps) { m_SNCsToBeRotated.push_back(snc); } else { hitMultipleLines = true; } } } } bool moveSlices(true); if (geometryToBeRotated && anyOtherGeometry && ourViewportGeometry && !hitMultipleLines) { // assure all three are valid, so calculation of center of rotation can be done moveSlices = false; } // question in state machine is: "rotate?" if (moveSlices) // i.e. NOT rotate { return false; } else { // we DO have enough information for rotation m_LastCursorPosition = intersectionLineWithGeometryToBeRotated.Project(cursorPosition); // remember where the last cursor position ON THE LINE has been observed if (anyOtherGeometry->IntersectionPoint(intersectionLineWithGeometryToBeRotated, m_CenterOfRotation)) // find center of rotation by intersection with any of the OTHER lines { return true; } else { return false; } } return false; } bool mitk::DisplayInteractor::CheckSwivelPossible(const mitk::InteractionEvent *interactionEvent) { const ScalarType ThresholdDistancePixels = 6.0; // Decide between moving and rotation: if we're close to the crossing // point of the planes, moving mode is entered, otherwise // rotation/swivel mode const InteractionPositionEvent *posEvent = dynamic_cast(interactionEvent); BaseRenderer *renderer = interactionEvent->GetSender(); if ( !posEvent || !renderer ) return false; const Point3D &cursor = posEvent->GetPositionInWorld(); m_SNCsToBeRotated.clear(); const PlaneGeometry *clickedGeometry( NULL ); const PlaneGeometry *otherGeometry1( NULL ); const PlaneGeometry *otherGeometry2( NULL ); auto renWindows = interactionEvent->GetSender()->GetRenderingManager()->GetAllRegisteredRenderWindows(); for(auto renWin : renWindows) { SliceNavigationController* snc = BaseRenderer::GetInstance(renWin)->GetSliceNavigationController(); // If the mouse cursor is in 3D Renderwindow, do not check for intersecting planes. if ( BaseRenderer::GetInstance(renWin)->GetMapperID() == BaseRenderer::Standard3D) continue; //unsigned int slice = (*iter)->GetSlice()->GetPos(); //unsigned int time = (*iter)->GetTime()->GetPos(); const PlaneGeometry *planeGeometry = snc->GetCurrentPlaneGeometry(); if ( !planeGeometry ) continue; if ( snc == renderer->GetSliceNavigationController() ) { clickedGeometry = planeGeometry; m_SNCsToBeRotated.push_back(snc); } else { if ( otherGeometry1 == NULL ) { otherGeometry1 = planeGeometry; } else { otherGeometry2 = planeGeometry; } if ( m_LinkPlanes ) { // If planes are linked, apply rotation to all planes m_SNCsToBeRotated.push_back(snc); } } } mitk::Line3D line; mitk::Point3D point; if ( (clickedGeometry != NULL) && (otherGeometry1 != NULL) && (otherGeometry2 != NULL) && clickedGeometry->IntersectionLine( otherGeometry1, line ) && otherGeometry2->IntersectionPoint( line, point )) { m_CenterOfRotation = point; if ( m_CenterOfRotation.EuclideanDistanceTo( cursor ) < ThresholdDistancePixels ) { return false; } else { m_ReferenceCursor = posEvent->GetPointerPositionOnScreen(); // Get main axes of rotation plane and store it for rotation step m_RotationPlaneNormal = clickedGeometry->GetNormal(); ScalarType xVector[] = { 1.0, 0.0, 0.0 }; ScalarType yVector[] = { 0.0, 1.0, 0.0 }; clickedGeometry->BaseGeometry::IndexToWorld( Vector3D( xVector), m_RotationPlaneXVector ); clickedGeometry->BaseGeometry::IndexToWorld( Vector3D( yVector), m_RotationPlaneYVector ); m_RotationPlaneNormal.Normalize(); m_RotationPlaneXVector.Normalize(); m_RotationPlaneYVector.Normalize(); m_PreviousRotationAxis.Fill( 0.0 ); m_PreviousRotationAxis[2] = 1.0; m_PreviousRotationAngle = 0.0; return true; } } else { return false; } return false; } void mitk::DisplayInteractor::Init(StateMachineAction*, InteractionEvent* interactionEvent) { InteractionPositionEvent* positionEvent = static_cast(interactionEvent); m_LastDisplayCoordinate = positionEvent->GetPointerPositionOnScreen(); m_CurrentDisplayCoordinate = m_LastDisplayCoordinate; positionEvent->GetSender()->DisplayToPlane(m_LastDisplayCoordinate,m_StartCoordinateInMM); m_LastCoordinateInMM = m_StartCoordinateInMM; } void mitk::DisplayInteractor::Move(StateMachineAction*, InteractionEvent* interactionEvent) { BaseRenderer* sender = interactionEvent->GetSender(); InteractionPositionEvent* positionEvent = static_cast(interactionEvent); float invertModifier = -1.0; if ( m_InvertMoveDirection ) { invertModifier = 1.0; } // perform translation Vector2D moveVector = (positionEvent->GetPointerPositionOnScreen() - m_LastDisplayCoordinate) * invertModifier; moveVector*=sender->GetScaleFactorMMPerDisplayUnit(); sender->GetCameraController()->MoveBy(moveVector); sender->GetRenderingManager()->RequestUpdate(sender->GetRenderWindow()); m_LastDisplayCoordinate = positionEvent->GetPointerPositionOnScreen(); } void mitk::DisplayInteractor::SetCrosshair(mitk::StateMachineAction *, mitk::InteractionEvent *interactionEvent) { const BaseRenderer::Pointer sender = interactionEvent->GetSender(); auto renWindows = sender->GetRenderingManager()->GetAllRegisteredRenderWindows(); InteractionPositionEvent* positionEvent = static_cast(interactionEvent); Point3D pos = positionEvent->GetPositionInWorld(); for(auto renWin : renWindows) { if (BaseRenderer::GetInstance(renWin)->GetMapperID() == BaseRenderer::Standard2D && renWin != sender->GetRenderWindow()) BaseRenderer::GetInstance(renWin)->GetSliceNavigationController()->SelectSliceByPoint(pos); } } void mitk::DisplayInteractor::Zoom(StateMachineAction*, InteractionEvent* interactionEvent) { const BaseRenderer::Pointer sender = interactionEvent->GetSender(); InteractionPositionEvent* positionEvent = static_cast(interactionEvent); float factor = 1.0; float distance = 0; if (m_ZoomDirection == "updown") { distance = m_CurrentDisplayCoordinate[1] - m_LastDisplayCoordinate[1]; } else { distance = m_CurrentDisplayCoordinate[0] - m_LastDisplayCoordinate[0]; } if ( m_InvertZoomDirection ) { distance *= -1.0; } // set zooming speed if (distance < 0.0) { factor = 1.0 / m_ZoomFactor; } else if (distance > 0.0) { factor = 1.0 * m_ZoomFactor; } if (factor != 1.0) { sender->GetCameraController()->Zoom(factor, m_StartCoordinateInMM); sender->GetRenderingManager()->RequestUpdate(sender->GetRenderWindow()); } m_LastDisplayCoordinate = m_CurrentDisplayCoordinate; m_CurrentDisplayCoordinate = positionEvent->GetPointerPositionOnScreen(); } void mitk::DisplayInteractor::Scroll(StateMachineAction*, InteractionEvent* interactionEvent) { InteractionPositionEvent* positionEvent = static_cast(interactionEvent); mitk::SliceNavigationController::Pointer sliceNaviController = interactionEvent->GetSender()->GetSliceNavigationController(); if (sliceNaviController) { int delta = 0; // Scrolling direction if (m_ScrollDirection == "updown") { delta = static_cast(m_LastDisplayCoordinate[1] - positionEvent->GetPointerPositionOnScreen()[1]); } else { delta = static_cast(m_LastDisplayCoordinate[0] - positionEvent->GetPointerPositionOnScreen()[0]); } if ( m_InvertScrollDirection ) { delta *= -1; } // Set how many pixels the mouse has to be moved to scroll one slice // if we moved less than 'm_IndexToSliceModifier' pixels slice ONE slice only if (delta > 0 && delta < m_IndexToSliceModifier) { delta = m_IndexToSliceModifier; } else if (delta < 0 && delta > -m_IndexToSliceModifier) { delta = -m_IndexToSliceModifier; } delta /= m_IndexToSliceModifier; int newPos = sliceNaviController->GetSlice()->GetPos() + delta; // if auto repeat is on, start at first slice if you reach the last slice and vice versa int maxSlices = sliceNaviController->GetSlice()->GetSteps(); if (m_AutoRepeat) { while (newPos < 0) { newPos += maxSlices; } while (newPos >= maxSlices) { newPos -= maxSlices; } } else { // if the new slice is below 0 we still show slice 0 // due to the stepper using unsigned int we have to do this ourselves if (newPos < 1) { newPos = 0; } } // set the new position sliceNaviController->GetSlice()->SetPos(newPos); m_LastDisplayCoordinate = m_CurrentDisplayCoordinate; m_CurrentDisplayCoordinate = positionEvent->GetPointerPositionOnScreen(); + MITK_INFO << "using mitk::DisplayInteractor::Scroll"; } } void mitk::DisplayInteractor::ScrollOneDown(StateMachineAction*, InteractionEvent* interactionEvent) { mitk::SliceNavigationController::Pointer sliceNaviController = interactionEvent->GetSender()->GetSliceNavigationController(); if (!sliceNaviController->GetSliceLocked()) { mitk::Stepper* stepper = sliceNaviController->GetSlice(); if (stepper->GetSteps() <= 1) { stepper = sliceNaviController->GetTime(); } stepper->Next(); + MITK_INFO << "using mitk::DisplayInteractor::ScrollOneDown"; } } void mitk::DisplayInteractor::ScrollOneUp(StateMachineAction*, InteractionEvent* interactionEvent) { mitk::SliceNavigationController::Pointer sliceNaviController = interactionEvent->GetSender()->GetSliceNavigationController(); if (!sliceNaviController->GetSliceLocked()) { mitk::Stepper* stepper = sliceNaviController->GetSlice(); if (stepper->GetSteps() <= 1) { stepper = sliceNaviController->GetTime(); } stepper->Previous(); + MITK_INFO << "using mitk::DisplayInteractor::ScrollOneUp"; } } void mitk::DisplayInteractor::AdjustLevelWindow(StateMachineAction*, InteractionEvent* interactionEvent) { BaseRenderer::Pointer sender = interactionEvent->GetSender(); InteractionPositionEvent* positionEvent = static_cast(interactionEvent); m_LastDisplayCoordinate = m_CurrentDisplayCoordinate; m_CurrentDisplayCoordinate = positionEvent->GetPointerPositionOnScreen(); // search for active image mitk::DataStorage::Pointer storage = sender->GetDataStorage(); mitk::DataNode::Pointer node = NULL; mitk::DataStorage::SetOfObjects::ConstPointer allImageNodes = storage->GetSubset(mitk::NodePredicateDataType::New("Image")); for (unsigned int i = 0; i < allImageNodes->size(); i++) { bool isActiveImage = false; bool propFound = allImageNodes->at(i)->GetBoolProperty("imageForLevelWindow", isActiveImage); if (propFound && isActiveImage) { node = allImageNodes->at(i); continue; } } if (node.IsNull()) { node = storage->GetNode(mitk::NodePredicateDataType::New("Image")); } if (node.IsNull()) { return; } mitk::LevelWindow lv = mitk::LevelWindow(); node->GetLevelWindow(lv); ScalarType level = lv.GetLevel(); ScalarType window = lv.GetWindow(); int levelIndex = 0; int windowIndex = 1; if ( m_LevelDirection != "leftright" ) { levelIndex = 1; windowIndex = 0; } int directionModifier = 1; if ( m_InvertLevelWindowDirection ) { directionModifier = -1; } // calculate adjustments from mouse movements level += (m_CurrentDisplayCoordinate[levelIndex] - m_LastDisplayCoordinate[levelIndex]) * static_cast(2) * directionModifier; window += (m_CurrentDisplayCoordinate[windowIndex] - m_LastDisplayCoordinate[windowIndex]) * static_cast(2) * directionModifier; lv.SetLevelWindow(level, window); dynamic_cast(node->GetProperty("levelwindow"))->SetLevelWindow(lv); sender->GetRenderingManager()->RequestUpdateAll(); } void mitk::DisplayInteractor::StartRotation(mitk::StateMachineAction *, mitk::InteractionEvent *) { this->SetMouseCursor(rotate_cursor_xpm, 0, 0); } void mitk::DisplayInteractor::EndRotation(mitk::StateMachineAction *, mitk::InteractionEvent *) { this->ResetMouseCursor(); } void mitk::DisplayInteractor::Rotate(mitk::StateMachineAction *, mitk::InteractionEvent * event) { const InteractionPositionEvent* posEvent = dynamic_cast(event); if (posEvent == nullptr) return; Point3D cursor = posEvent->GetPositionInWorld(); Vector3D toProjected = m_LastCursorPosition - m_CenterOfRotation; Vector3D toCursor = cursor - m_CenterOfRotation; // cross product: | A x B | = |A| * |B| * sin(angle) Vector3D axisOfRotation; vnl_vector_fixed< ScalarType, 3 > vnlDirection = vnl_cross_3d(toCursor.GetVnlVector(), toProjected.GetVnlVector()); axisOfRotation.SetVnlVector(vnlDirection); // scalar product: A * B = |A| * |B| * cos(angle) // tan = sin / cos ScalarType angle = -atan2((double)(axisOfRotation.GetNorm()), (double)(toCursor * toProjected)); angle *= 180.0 / vnl_math::pi; m_LastCursorPosition = cursor; // create RotationOperation and apply to all SNCs that should be rotated RotationOperation rotationOperation(OpROTATE, m_CenterOfRotation, axisOfRotation, angle); // iterate the OTHER slice navigation controllers: these are filled in DoDecideBetweenRotationAndSliceSelection for (SNCVector::iterator iter = m_SNCsToBeRotated.begin(); iter != m_SNCsToBeRotated.end(); ++iter) { TimeGeometry* timeGeometry = (*iter)->GetCreatedWorldGeometry(); if (!timeGeometry) continue; timeGeometry->ExecuteOperation(&rotationOperation); (*iter)->SendCreatedWorldGeometryUpdate(); } RenderingManager::GetInstance()->RequestUpdateAll(); } void mitk::DisplayInteractor::Swivel(mitk::StateMachineAction *, mitk::InteractionEvent *event) { const InteractionPositionEvent *posEvent = dynamic_cast(event); if (!posEvent) return; // Determine relative mouse movement projected onto world space Point2D cursor = posEvent->GetPointerPositionOnScreen(); Vector2D relativeCursor = cursor - m_ReferenceCursor; Vector3D relativeCursorAxis = m_RotationPlaneXVector * relativeCursor[0] + m_RotationPlaneYVector * relativeCursor[1]; // Determine rotation axis (perpendicular to rotation plane and cursor // movement) Vector3D rotationAxis = itk::CrossProduct( m_RotationPlaneNormal, relativeCursorAxis ); ScalarType rotationAngle = relativeCursor.GetNorm() / 2.0; // Restore the initial plane pose by undoing the previous rotation // operation RotationOperation op( OpROTATE, m_CenterOfRotation, m_PreviousRotationAxis, -m_PreviousRotationAngle ); SNCVector::iterator iter; for ( iter = m_SNCsToBeRotated.begin(); iter != m_SNCsToBeRotated.end(); ++iter ) { if ( !(*iter)->GetSliceRotationLocked() ) { TimeGeometry* timeGeometry = (*iter)->GetCreatedWorldGeometry(); if (!timeGeometry) continue; timeGeometry->ExecuteOperation(&op); (*iter)->SendCreatedWorldGeometryUpdate(); } } // Apply new rotation operation to all relevant SNCs RotationOperation op2( OpROTATE, m_CenterOfRotation, rotationAxis, rotationAngle ); for ( iter = m_SNCsToBeRotated.begin(); iter != m_SNCsToBeRotated.end(); ++iter) { if ( !(*iter)->GetSliceRotationLocked() ) { // Retrieve the TimeGeometry of this SliceNavigationController TimeGeometry* timeGeometry = (*iter)->GetCreatedWorldGeometry(); if (!timeGeometry) continue; // Execute the new rotation timeGeometry->ExecuteOperation(&op2); // Notify listeners (*iter)->SendCreatedWorldGeometryUpdate(); } } m_PreviousRotationAxis = rotationAxis; m_PreviousRotationAngle = rotationAngle; RenderingManager::GetInstance()->RequestUpdateAll(); return; } -void mitk::DisplayInteractor::UpdateStatusbar(mitk::StateMachineAction *, mitk::InteractionEvent *event) +void mitk::DisplayInteractor::SetStatusBarInfoFromClickedPosition(mitk::StateMachineAction *, mitk::InteractionEvent *event) { const InteractionPositionEvent *posEvent = dynamic_cast(event); if (!posEvent) return; std::string statusText; TNodePredicateDataType::Pointer isImageData = TNodePredicateDataType::New(); mitk::DataStorage::SetOfObjects::ConstPointer nodes = posEvent->GetSender()->GetDataStorage()->GetSubset(isImageData).GetPointer(); mitk::Point3D worldposition = posEvent->GetPositionInWorld(); mitk::Image::Pointer image3D; mitk::DataNode::Pointer node; mitk::DataNode::Pointer topSourceNode; bool isBinary (false); int component = 0; node = this->GetTopLayerNode(nodes,worldposition,posEvent->GetSender()); if(node.IsNotNull()) { node->GetBoolProperty("binary", isBinary); if(isBinary) { mitk::DataStorage::SetOfObjects::ConstPointer sourcenodes = posEvent->GetSender()->GetDataStorage()->GetSources(node, NULL, true); if(!sourcenodes->empty()) { topSourceNode = this->GetTopLayerNode(sourcenodes,worldposition,posEvent->GetSender()); } if(topSourceNode.IsNotNull()) { image3D = dynamic_cast(topSourceNode->GetData()); topSourceNode->GetIntProperty("Image.Displayed Component", component); } else { image3D = dynamic_cast(node->GetData()); node->GetIntProperty("Image.Displayed Component", component); } } else { image3D = dynamic_cast(node->GetData()); node->GetIntProperty("Image.Displayed Component", component); } } - std::stringstream stream; - stream.imbue(std::locale::classic()); - // get the position and gray value from the image and build up status bar text +// std::stringstream stream; +// stream.imbue(std::locale::classic()); + +// // get the position and gray value from the image and build up status bar text if(image3D.IsNotNull()) { itk::Index<3> p; image3D->GetGeometry()->WorldToIndex(worldposition, p); - stream.precision(2); - stream<<"Position: <" << std::fixed < mm"; - stream<<"; Index: <"< "; +// stream.precision(2); +// stream<<"Position: <" << std::fixed < mm"; +// stream<<"; Index: <"< "; + +// mitk::ScalarType pixelValue; + +// mitkPixelTypeMultiplex5( +// mitk::FastSinglePixelAccess, +// image3D->GetChannelDescriptor().GetPixelType(), +// image3D, +// image3D->GetVolumeData(posEvent->GetSender()->GetTimeStep()), +// p, +// pixelValue, +// component); + + + +// if (fabs(pixelValue)>1000000 || fabs(pixelValue) < 0.01) +// { +// stream<<"; Time: " << posEvent->GetSender()->GetTime() << " ms; Pixelvalue: " << std::scientific<< pixelValue <<" "; +// } +// else +// { +// stream<<"; Time: " << posEvent->GetSender()->GetTime() << " ms; Pixelvalue: "<< pixelValue <<" "; +// } +// } +// else +// { +// stream << "No image information at this position!"; +// } + +// statusText = stream.str(); +// mitk::StatusBar::GetInstance()->DisplayGreyValueText(statusText.c_str()); + TimeStepType time = posEvent->GetSender()->GetTimeStep(); + UpdateStatusBar(image3D, p, time, component); + } +} - mitk::ScalarType pixelValue; +void mitk::DisplayInteractor::UpdateStatusBar(itk::SmartPointer image3D, itk::Index<3> idx, TimeStepType time, int component) +{ + std::stringstream stream; + stream.imbue(std::locale::classic()); - mitkPixelTypeMultiplex5( - mitk::FastSinglePixelAccess, - image3D->GetChannelDescriptor().GetPixelType(), - image3D, - image3D->GetVolumeData(posEvent->GetSender()->GetTimeStep()), - p, - pixelValue, - component); + // get the position and gray value from the image and build up status bar text + if(image3D.IsNotNull()) + { + mitk::Point3D worldposition; + image3D->GetGeometry()->IndexToWorld(idx, worldposition); + stream.precision(2); + stream<<"Position: <" << std::fixed < mm"; + stream<<"; Index: <"< "; + mitk::ScalarType pixelValue; + mitkPixelTypeMultiplex5( + mitk::FastSinglePixelAccess, + image3D->GetChannelDescriptor().GetPixelType(), + image3D, + image3D->GetVolumeData(time), + idx, + pixelValue, + component); - if (fabs(pixelValue)>1000000 || fabs(pixelValue) < 0.01) - { - stream<<"; Time: " << posEvent->GetSender()->GetTime() << " ms; Pixelvalue: " << std::scientific<< pixelValue <<" "; + + + if (fabs(pixelValue)>1000000 || fabs(pixelValue) < 0.01) + { + stream<<"; Pixelvalue: " << std::scientific<< pixelValue <<" "; + } + else + { + stream<<"; Pixelvalue: "<< pixelValue <<" "; + } } else { - stream<<"; Time: " << posEvent->GetSender()->GetTime() << " ms; Pixelvalue: "<< pixelValue <<" "; + stream << "No image information at this position!"; } - } - else - { - stream << "No image information at this position!"; - } - statusText = stream.str(); - mitk::StatusBar::GetInstance()->DisplayGreyValueText(statusText.c_str()); + std::string statusText = stream.str(); + mitk::StatusBar::GetInstance()->DisplayGreyValueText(statusText.c_str()); } void mitk::DisplayInteractor::ConfigurationChanged() { mitk::PropertyList::Pointer properties = GetAttributes(); // auto repeat std::string strAutoRepeat = ""; if (properties->GetStringProperty("autoRepeat", strAutoRepeat)) { if (strAutoRepeat == "true") { m_AutoRepeat = true; } else { m_AutoRepeat = false; } } // pixel movement for scrolling one slice std::string strPixelPerSlice = ""; if (properties->GetStringProperty("pixelPerSlice", strPixelPerSlice)) { m_IndexToSliceModifier = atoi(strPixelPerSlice.c_str()); } else { m_IndexToSliceModifier = 4; } // scroll direction if (!properties->GetStringProperty("scrollDirection", m_ScrollDirection)) { m_ScrollDirection = "updown"; } m_InvertScrollDirection = GetBoolProperty( properties, "invertScrollDirection", false ); // zoom direction if (!properties->GetStringProperty("zoomDirection", m_ZoomDirection)) { m_ZoomDirection = "updown"; } m_InvertZoomDirection = GetBoolProperty( properties, "invertZoomDirection", false ); m_InvertMoveDirection = GetBoolProperty( properties, "invertMoveDirection", false ); if (!properties->GetStringProperty("levelWindowDirection", m_LevelDirection)) { m_LevelDirection = "leftright"; } m_InvertLevelWindowDirection = GetBoolProperty( properties, "invertLevelWindowDirection", false ); // coupled rotation std::string strCoupled= ""; if (properties->GetStringProperty("coupled", strCoupled)) { if (strCoupled == "true") m_LinkPlanes = true; else m_LinkPlanes = false; } // zoom factor std::string strZoomFactor = ""; properties->GetStringProperty("zoomFactor", strZoomFactor); m_ZoomFactor = .05; if (atoi(strZoomFactor.c_str()) > 0) { m_ZoomFactor = 1.0 + (atoi(strZoomFactor.c_str()) / 100.0); } // allwaysReact std::string strAlwaysReact = ""; if (properties->GetStringProperty("alwaysReact", strAlwaysReact)) { if (strAlwaysReact == "true") { m_AlwaysReact = true; } else { m_AlwaysReact = false; } } else { m_AlwaysReact = false; } } bool mitk::DisplayInteractor::FilterEvents(InteractionEvent* interactionEvent, DataNode* /*dataNode*/) { if (interactionEvent->GetSender() == nullptr) return false; if (interactionEvent->GetSender()->GetMapperID() == BaseRenderer::Standard3D) return false; return true; } bool mitk::DisplayInteractor::GetBoolProperty( mitk::PropertyList::Pointer propertyList, const char* propertyName, bool defaultValue ) { std::string valueAsString; if ( !propertyList->GetStringProperty( propertyName, valueAsString ) ) { return defaultValue; } else { if ( valueAsString == "true" ) { return true; } else { return false; } } } mitk::DataNode::Pointer mitk::DisplayInteractor::GetTopLayerNode(mitk::DataStorage::SetOfObjects::ConstPointer nodes, mitk::Point3D worldposition, BaseRenderer *ren) { mitk::DataNode::Pointer node; int maxlayer = -32768; bool isHelper (false); if(nodes.IsNotNull()) { for (unsigned int x = 0; x < nodes->size(); x++) { nodes->at(x)->GetBoolProperty("helper object", isHelper); if(nodes->at(x)->GetData()->GetGeometry()->IsInside(worldposition) && isHelper == false) { int layer = 0; if(!(nodes->at(x)->GetIntProperty("layer", layer))) continue; if(layer > maxlayer) { if(static_cast(nodes->at(x))->IsVisible(ren)) { node = nodes->at(x); maxlayer = layer; } } } } } return node; } diff --git a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt index 5736f18a80..e4564aff4f 100755 --- a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt +++ b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt @@ -1,117 +1,120 @@ option(BUILD_DiffusionMiniApps "Build commandline tools for diffusion" OFF) if(BUILD_DiffusionMiniApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion miniapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionminiapps DwiDenoising^^ ImageResampler^^ NetworkCreation^^MitkFiberTracking_MitkConnectomics NetworkStatistics^^MitkConnectomics ExportShImage^^ + ImageStatisticsMiniApp^^MitkImageStatistics + PlanarFigureMaskCreatorMiniapp^^MitkImageStatistics + HotspotMiniApp^^MitkImageStatistics Fiberfox^^MitkFiberTracking MultishellMethods^^MitkFiberTracking PeaksAngularError^^MitkFiberTracking PeakExtraction^^MitkFiberTracking FiberExtraction^^MitkFiberTracking FiberProcessing^^MitkFiberTracking FiberDirectionExtraction^^MitkFiberTracking LocalDirectionalFiberPlausibility^^MitkFiberTracking StreamlineTracking^^MitkFiberTracking GibbsTracking^^MitkFiberTracking CopyGeometry^^ DiffusionIndices^^ TractometerMetrics^^MitkFiberTracking QballReconstruction^^ Registration^^ FileFormatConverter^^MitkFiberTracking TensorReconstruction^^ TensorDerivedMapsExtraction^^ DiffusionDICOMLoader^^ DFTraining^^MitkFiberTracking DFTracking^^MitkFiberTracking ) foreach(diffusionminiapp ${diffusionminiapps}) # extract mini app name and dependencies string(REPLACE "^^" "\\;" miniapp_info ${diffusionminiapp}) set(miniapp_info_list ${miniapp_info}) list(GET miniapp_info_list 0 appname) list(GET miniapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitk_create_executable(${appname} DEPENDS MitkCore MitkDiffusionCore MitkCommandLine ${dependencies_list} PACKAGE_DEPENDS ITK CPP_FILES ${appname}.cpp ) if(EXECUTABLE_IS_ENABLED) # On Linux, create a shell script to start a relocatable application if(UNIX AND NOT APPLE) install(PROGRAMS "${MITK_SOURCE_DIR}/CMake/RunInstalledApp.sh" DESTINATION "." RENAME ${EXECUTABLE_TARGET}.sh) endif() get_target_property(_is_bundle ${EXECUTABLE_TARGET} MACOSX_BUNDLE) if(APPLE) if(_is_bundle) set(_target_locations ${EXECUTABLE_TARGET}.app) set(${_target_locations}_qt_plugins_install_dir ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_bundle_dest_dir ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_qt_plugins_for_current_bundle ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_qt_conf_install_dirs ${EXECUTABLE_TARGET}.app/Contents/Resources) install(TARGETS ${EXECUTABLE_TARGET} BUNDLE DESTINATION . ) else() if(NOT MACOSX_BUNDLE_NAMES) set(_qt_conf_install_dirs bin) set(_target_locations bin/${EXECUTABLE_TARGET}) set(${_target_locations}_qt_plugins_install_dir bin) install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) else() foreach(bundle_name ${MACOSX_BUNDLE_NAMES}) list(APPEND _qt_conf_install_dirs ${bundle_name}.app/Contents/Resources) set(_current_target_location ${bundle_name}.app/Contents/MacOS/${EXECUTABLE_TARGET}) list(APPEND _target_locations ${_current_target_location}) set(${_current_target_location}_qt_plugins_install_dir ${bundle_name}.app/Contents/MacOS) message( " set(${_current_target_location}_qt_plugins_install_dir ${bundle_name}.app/Contents/MacOS) ") install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION ${bundle_name}.app/Contents/MacOS/) endforeach() endif() endif() else() set(_target_locations bin/${EXECUTABLE_TARGET}${CMAKE_EXECUTABLE_SUFFIX}) set(${_target_locations}_qt_plugins_install_dir bin) set(_qt_conf_install_dirs bin) install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) endif() endif() endforeach() # This mini app does not depend on mitkDiffusionImaging at all mitk_create_executable(Dicom2Nrrd DEPENDS MitkCore MitkCommandLine CPP_FILES Dicom2Nrrd.cpp ) # On Linux, create a shell script to start a relocatable application if(UNIX AND NOT APPLE) install(PROGRAMS "${MITK_SOURCE_DIR}/CMake/RunInstalledApp.sh" DESTINATION "." RENAME ${EXECUTABLE_TARGET}.sh) endif() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/MiniApps/ExtractImageStatistics.cpp b/Modules/DiffusionImaging/MiniApps/ExtractImageStatistics.cpp index e30bbfa158..b89d0995b8 100644 --- a/Modules/DiffusionImaging/MiniApps/ExtractImageStatistics.cpp +++ b/Modules/DiffusionImaging/MiniApps/ExtractImageStatistics.cpp @@ -1,127 +1,129 @@ /*=================================================================== 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 "MiniAppManager.h" +//#include "MiniAppManager.h" -#include "ctkCommandLineParser.h" +//#include "ctkCommandLineParser.h" +#include "mitkCommandLineParser.cpp" #include "mitkImage.h" #include "mitkImageStatisticsCalculator.h" #include "mitkIOUtil.h" #include #include #include int ExtractImageStatistics(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Extract Image Statistics"); parser.setCategory("Preprocessing Tools"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),false); parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: filenameOfRoi.nrrd_statistics.txt)", us::Any()); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << "\n\n MiniApp Description: \nCalculates statistics on the supplied image using given mask." << endl; std::cout << "Output is written to the designated output file in this order:" << endl; std::cout << "Mean, Standard Deviation, RMS, Max, Min, Number of Voxels, Volume [mm3]" << endl; std::cout << "\n\n Parameters:"<< endl; std::cout << parser.helpText(); return EXIT_SUCCESS; } + // Parameters: bool ignoreZeroValues = false; unsigned int timeStep = 0; std::string inputImageFile = us::any_cast(parsedArgs["input"]); std::string maskImageFile = us::any_cast(parsedArgs["mask"]); std::string outFile; if (parsedArgs.count("out") || parsedArgs.count("o") ) outFile = us::any_cast(parsedArgs["out"]); else outFile = inputImageFile + "_statistics.txt"; // Load image and mask mitk::Image::Pointer maskImage = mitk::IOUtil::LoadImage(maskImageFile); mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); // Calculate statistics mitk::ImageStatisticsCalculator::Statistics statisticsStruct; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); try { calculator->SetImage(inputImage); calculator->SetImageMask(maskImage); calculator->SetMaskingModeToImage(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Statistic Calculation Failed - ITK Exception:" << e.what(); return -1; } calculator->SetDoIgnorePixelValue(ignoreZeroValues); calculator->SetIgnorePixelValue(0); try { calculator->ComputeStatistics(timeStep); } catch ( mitk::Exception& e) { MITK_ERROR<< "MITK Exception: " << e.what(); return -1; } statisticsStruct = calculator->GetStatistics(timeStep); // Calculate Volume double volume = 0; const mitk::BaseGeometry *geometry = inputImage->GetGeometry(); if ( geometry != NULL ) { const mitk::Vector3D &spacing = inputImage->GetGeometry()->GetSpacing(); volume = spacing[0] * spacing[1] * spacing[2] * (double) statisticsStruct.GetN(); } // Write Results to file std::ofstream output; output.open(outFile.c_str()); output << statisticsStruct.GetMean() << " , "; output << statisticsStruct.GetSigma() << " , "; output << statisticsStruct.GetRMS() << " , "; output << statisticsStruct.GetMax() << " , "; output << statisticsStruct.GetMin() << " , "; output << statisticsStruct.GetN() << " , "; output << volume << "\n"; output.flush(); output.close(); return 0; } RegisterDiffusionMiniApp(ExtractImageStatistics); diff --git a/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp new file mode 100644 index 0000000000..2f76e385e9 --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp @@ -0,0 +1,104 @@ +/*=================================================================== + +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 "mitkCommandLineParser.h" +#include "mitkImage.h" +#include "mitkIOUtil.h" +#include +#include +#include +#include + + +int main( int argc, char* argv[] ) +{ + mitkCommandLineParser parser; + + parser.setTitle("Test basic hotspot functionality"); + parser.setCategory("Preprocessing Tools"); + parser.setDescription(""); + parser.setContributor("MBI"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); + parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); + parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),true); + parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: hotspotMiniApp_output.nrrd)", us::Any()); + + std::cout << "test...." << std::endl; + + std::map parsedArgs = parser.parseArguments(argc, argv); + std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; + if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) + { +// std::cout << "\n\n MiniApp Description: \nCalculates and saves the hotspot of an image." << endl; +// std::cout << "Output is written to the designated output file" << endl; +// std::cout << parser.helpText(); +// return EXIT_SUCCESS; + } + + // Parameters: + std::string inputImageFile; + if (!parsedArgs.count("i") && !parsedArgs.count("input")) + { + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + } + else + { + inputImageFile = us::any_cast(parsedArgs["input"]); + } + + mitk::Image::Pointer maskImage; + if (parsedArgs.count("mask") || parsedArgs.count("m")) + { + std::string maskImageFile = us::any_cast(parsedArgs["mask"]); + maskImage = mitk::IOUtil::LoadImage(maskImageFile); + } + + std::string outFile; + if (parsedArgs.count("out") || parsedArgs.count("o") ) + outFile = us::any_cast(parsedArgs["out"]); + else + outFile = "hotspotMiniApp_output.nrrd"; + + // Load image and mask + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); + + // Calculate statistics + mitk::HotspotCalculator hotspotCalc; + hotspotCalc.setInputImage(inputImage); + hotspotCalc.setMaskImage(maskImage); + hotspotCalc.setHotspotRadiusInMM(10); + + mitk::Image::Pointer outImage; + try + { + outImage = hotspotCalc.getHotspotMask(0, 0); + } + catch( const itk::ExceptionObject& e) + { + MITK_ERROR << "Failed - ITK Exception:" << e.what(); + return -1; + } + + if (outImage != nullptr) + { + mitk::IOUtil::SaveImage(outImage, outFile); + } + + + return EXIT_SUCCESS; +} diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniApp.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniApp.cpp new file mode 100644 index 0000000000..fe950e639f --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniApp.cpp @@ -0,0 +1,210 @@ +/*=================================================================== + +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 "mitkCommandLineParser.h" +#include "mitkImage.h" +#include "mitkImageStatisticsCalculator.h" +#include "mitkIOUtil.h" +#include +#include +#include +#include +#include "mitkImageAccessByItk.h" +#include +#include +#include +#include +#include +#include +#include +#include + + +struct statistics_res{ + double mean, variance, min, max, count, moment; +}; + +void printstats(statistics_res s) +{ + std::cout << "mean: " << s.mean << std::endl + << "variance: " << s.variance << std::endl + << "min: " << s.min << std::endl + << "max: " << s.max << std::endl + << "count: " << s.count << std::endl + << "moment: " << s.moment << std::endl; +} + +template < typename TPixel, unsigned int VImageDimension > +void get_statistics_boost(itk::Image* itkImage, statistics_res& res){ + typedef itk::Image ImageType; + + itk::ImageRegionConstIterator it(itkImage, itkImage->GetLargestPossibleRegion()); + + TPixel currentPixel; + int ctr=0; + double sum=0; + + boost::accumulators::accumulator_set> > acc; + + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + acc(it.Get()); +// currentPixel = it.Get(); +// sum+=currentPixel; +// ctr+=1; + } + +// res.mean=(double)sum/(double)ctr; + res.mean = boost::accumulators::mean(acc); + res.variance = boost::accumulators::variance(acc); + res.min = boost::accumulators::min(acc); + res.max = boost::accumulators::max(acc); + res.count = boost::accumulators::count(acc); + res.moment = boost::accumulators::moment<2>(acc); + + std::cout << "sum: " << sum << " N: " << ctr << " mean: " << res.mean << std::endl; +} + +int main( int argc, char* argv[] ) +{ + mitkCommandLineParser parser; + + parser.setTitle("Extract Image Statistics"); + parser.setCategory("Preprocessing Tools"); + parser.setDescription(""); + parser.setContributor("MBI"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); + parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); + parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),true); + parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: filenameOfRoi.nrrd_statistics.txt)", us::Any()); + + std::cout << "test...." << std::endl; + + std::map parsedArgs = parser.parseArguments(argc, argv); + std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; + if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) + { + std::cout << "\n\n MiniApp Description: \nCalculates statistics on the supplied image using given mask." << endl; + std::cout << "Output is written to the designated output file in this order:" << endl; + std::cout << "Mean, Standard Deviation, RMS, Max, Min, Number of Voxels, Volume [mm3]" << endl; + std::cout << "\n\n Parameters:"<< endl; + std::cout << parser.helpText(); + return EXIT_SUCCESS; + } + + + // Parameters: + bool ignoreZeroValues = false; + unsigned int timeStep = 0; + + std::string inputImageFile = us::any_cast(parsedArgs["input"]); + mitk::Image::Pointer maskImage; + if (parsedArgs.count("mask") || parsedArgs.count("m")) + { + std::string maskImageFile = us::any_cast(parsedArgs["mask"]); + maskImage = mitk::IOUtil::LoadImage(maskImageFile); + } + + std::string outFile; + if (parsedArgs.count("out") || parsedArgs.count("o") ) + outFile = us::any_cast(parsedArgs["out"]); + else + outFile = inputImageFile + "_statistics.txt"; + + // Load image and mask + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); + + // Calculate statistics + mitk::ImageStatisticsCalculator::Statistics statisticsStruct; + mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); + try + { + calculator->SetImage(inputImage); + if (parsedArgs.count("mask") || parsedArgs.count("m")) + { + calculator->SetImageMask(maskImage); + calculator->SetMaskingModeToImage(); + } + else + { + calculator->SetMaskingModeToNone(); + } + + } + catch( const itk::ExceptionObject& e) + { + MITK_ERROR << "Statistic Calculation Failed - ITK Exception:" << e.what(); + return -1; + } + + + calculator->SetDoIgnorePixelValue(ignoreZeroValues); + calculator->SetIgnorePixelValue(0); + std::cout << "calculating statistics itk: " << std::endl; + try + { + calculator->ComputeStatistics(timeStep); + } + catch ( mitk::Exception& e) + { + MITK_ERROR<< "MITK Exception: " << e.what(); + return -1; + } + + + statisticsStruct = calculator->GetStatistics(timeStep); + + + // Calculate Volume + double volume = 0; + const mitk::BaseGeometry *geometry = inputImage->GetGeometry(); + if ( geometry != NULL ) + { + const mitk::Vector3D &spacing = inputImage->GetGeometry()->GetSpacing(); + volume = spacing[0] * spacing[1] * spacing[2] * (double) statisticsStruct.GetN(); + } + + // Write Results to file + std::ofstream output; + output.open(outFile.c_str()); + output << statisticsStruct.GetMean() << " , "; + output << statisticsStruct.GetSigma() << " , "; + output << statisticsStruct.GetRMS() << " , "; + output << statisticsStruct.GetMax() << " , "; + output << statisticsStruct.GetMin() << " , "; + output << statisticsStruct.GetN() << " , "; + output << volume << "\n"; + + output.flush(); + output.close(); + + std::cout << "calculating statistics boost: " << std::endl; + + statistics_res res; + AccessByItk_n(inputImage, get_statistics_boost, (res)); + + printstats(res); + + return EXIT_SUCCESS; +} diff --git a/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp b/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp new file mode 100644 index 0000000000..c099bdf1e9 --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp @@ -0,0 +1,131 @@ +/*=================================================================== + +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 "mitkCommandLineParser.h" +#include "mitkImage.h" +#include "mitkIOUtil.h" +#include +#include +#include +#include +#include +#include + + +int main( int argc, char* argv[] ) +{ + mitkCommandLineParser parser; + + parser.setTitle("Test basic hotspot functionality"); + parser.setCategory("Preprocessing Tools"); + parser.setDescription(""); + parser.setContributor("MBI"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); + parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); + parser.addArgument("planarfigure", "p", mitkCommandLineParser::InputFile, "PlanarFigure:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),false); + parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: hotspotMiniApp_output.nrrd)", us::Any()); + + std::cout << "test...." << std::endl; + + std::map parsedArgs = parser.parseArguments(argc, argv); + std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; + if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) + { +// std::cout << "\n\n MiniApp Description: \nCalculates and saves the hotspot of an image." << endl; +// std::cout << "Output is written to the designated output file" << endl; +// std::cout << parser.helpText(); +// return EXIT_SUCCESS; + } + + // Parameters: + std::string inputImageFile; + if (!parsedArgs.count("i") && !parsedArgs.count("input")) + { + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + } + else + { + inputImageFile = us::any_cast(parsedArgs["input"]); + } + + std::string inputPFFile; + if (parsedArgs.count("planarfigure") || parsedArgs.count("p")) + { + inputPFFile = us::any_cast(parsedArgs["planarfigure"]); + } + else + { + inputPFFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; + } + + + mitk::PlanarFigure::Pointer planarFigure; + try + { + std::vector loadedObjects; + // try except + loadedObjects = mitk::IOUtil::Load(inputPFFile); + mitk::BaseData::Pointer pf = loadedObjects[0]; + planarFigure = dynamic_cast(pf.GetPointer()); + if (planarFigure.IsNull()) + { + MITK_ERROR << "something went wrong"; + return -1; + } + } + catch (const itk::ExceptionObject& e) + { + MITK_ERROR << "Failed: " << e.what(); + return -1; + } + + std::string outFile; + if (parsedArgs.count("out") || parsedArgs.count("o") ) + outFile = us::any_cast(parsedArgs["out"]); + else + outFile = "planarFigureExtraction_output.nrrd"; + + // Load image and mask + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); + + // Calculate statistics + mitk::PlanarFigureMaskGenerator planarFigMaskExtr; + planarFigMaskExtr.SetPlanarFigure(planarFigure); + planarFigMaskExtr.SetImage(inputImage); + + mitk::Image::Pointer outImage; + try + { + outImage = planarFigMaskExtr.GetMask(); + } + catch( const itk::ExceptionObject& e) + { + MITK_ERROR << "Failed - ITK Exception:" << e.what(); + return -1; + } + + mitk::BaseGeometry *imageGeo = outImage->GetGeometry(); + std::cout << "origin: " << imageGeo->GetOrigin()[0] << " " << imageGeo->GetOrigin()[1] << " " << imageGeo->GetOrigin()[2] << std::endl; + if (outImage != nullptr) + { + mitk::IOUtil::SaveImage(outImage, outFile); + } + + + return EXIT_SUCCESS; +} diff --git a/Modules/ImageStatistics/CMakeLists.txt b/Modules/ImageStatistics/CMakeLists.txt index 5264c1c320..fa0a7cf973 100644 --- a/Modules/ImageStatistics/CMakeLists.txt +++ b/Modules/ImageStatistics/CMakeLists.txt @@ -1,13 +1,13 @@ MITK_CREATE_MODULE( - DEPENDS MitkImageExtraction MitkPlanarFigure + DEPENDS MitkImageExtraction MitkPlanarFigure MitkMultilabel PACKAGE_DEPENDS PUBLIC ITK|ITKIOXML PRIVATE ITK|ITKVTK+ITKConvolution # WARNINGS_AS_ERRORS ) if(BUILD_TESTING) add_subdirectory(Testing) endif(BUILD_TESTING) diff --git a/Modules/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index 29d3c3547a..be094a5f7a 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,14 +1,24 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp + mitkImageStatisticsCalculator2.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp + mitkHotspotCalculator.cpp + mitkMaskGenerator.cpp + mitkPlanarFigureMaskGenerator.cpp + mitkMultiLabelMaskGenerator.cpp ) set(H_FILES mitkImageStatisticsCalculator.h + mitkImageStatisticsCalculator2.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkExtendedStatisticsImageFilter.h mitkExtendedLabelStatisticsImageFilter.h + mitkHotspotCalculator.h + mitkMaskGenerator.h + mitkPlanarFigureMaskGenerator.h + mitkMultiLabelMaskGenerator.h ) diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx.autosave b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx.autosave new file mode 100644 index 0000000000..36fd1884f4 --- /dev/null +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx.autosave @@ -0,0 +1,324 @@ +/*=================================================================== + +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 _mitkExtendedLabelStatisticsImageFilter_hxx +#define _mitkExtendedLabelStatisticsImageFilter_hxx + +#include "mitkExtendedLabelStatisticsImageFilter.h" + +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageRegionConstIterator.h" +#include +#include +#include "mitkNumericConstants.h" +#include "mitkLogMacros.h" + +namespace itk +{ + template< class TInputImage , class TLabelImage> + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::ExtendedLabelStatisticsImageFilter() + : LabelStatisticsImageFilter< TInputImage, TLabelImage >() + { + CalculateSettingsForLabels(); + } + + + template< class TInputImage, class TLabelImage > + std::list + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::GetRelevantLabels() const + { + return m_RelevantLabels; + } + + + template< class TInputImage, class TLabelImage > + bool + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::GetMaskingNonEmpty() const + { + return m_MaskNonEmpty; + } + + + template< class TInputImage, class TLabelImage > + typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::GetUniformity(LabelPixelType label) const + { + CoefficientsMapConstIterator mapIt; + + mapIt = m_LabelStatisticsCoefficients.find(label); + if ( mapIt == m_LabelStatisticsCoefficients.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::Zero; + } + else + { + return ( *mapIt ).second.m_Uniformity; + } + } + + template< class TInputImage, class TLabelImage > + typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::GetEntropy(LabelPixelType label) const + { + CoefficientsMapConstIterator mapIt; + + mapIt = m_LabelStatisticsCoefficients.find(label); + if ( mapIt == m_LabelStatisticsCoefficients.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::Zero; + } + else + { + return ( *mapIt ).second.m_Entropy; + } + } + + template< class TInputImage, class TLabelImage > + typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::GetUPP(LabelPixelType label) const + { + CoefficientsMapConstIterator mapIt; + + mapIt = m_LabelStatisticsCoefficients.find(label); + if ( mapIt == m_LabelStatisticsCoauch efficients.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::Zero; + } + else + { + return ( *mapIt ).second.m_UPP; + } + } + + template< class TInputImage, class TLabelImage > + typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::GetMPP(LabelPixelType label) const + { + CoefficientsMapConstIterator mapIt; + + mapIt = m_LabelStatisticsCoefficients.find(label); + if ( mapIt == m_LabelStatisticsCoefficients.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::Zero; + } + else + { + return ( *mapIt ).second.m_MPP; + } + } + + template< class TInputImage, class TLabelImage > + typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::GetKurtosis(LabelPixelType label) const + { + CoefficientsMapConstIterator mapIt; + + mapIt = m_LabelStatisticsCoefficients.find(label); + if ( mapIt == m_LabelStatisticsCoefficients.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::Zero; + } + else + { + return ( *mapIt ).second.m_Kurtosis; + } + } + + + template< class TInputImage, class TLabelImage > + typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > + ::GetSkewness(LabelPixelType label) const + { + CoefficientsMapConstIterator mapIt; + + mapIt = m_LabelStatisticsCoefficients.find(label); + if ( mapIt == m_LabelStatisticsCoefficients.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::Zero; + } + else + { + return ( *mapIt ).second.m_Skewness; + } + } + + + template< class TInputImage, class TLabelImage > + void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: + CalculateSettingsForLabels() + { + LabelPixelType i; + m_MaskNonEmpty = false; + for ( i = 1; i < 4096; ++i ) + { + if ( this->HasLabel( i ) ) + { + m_RelevantLabels.push_back( i ); + m_MaskNonEmpty = true; + m_LabelStatisticsCoefficients.insert( std::make_pair(i, CoefficientsClass()) ); + } + } + } + + + template< class TInputImage, class TLabelImage > + void + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: + ComputeEntropyUniformityAndUPP() + { + double baseChange = std::log10(2); + RealType partialProbability( 0.0 ); + RealType uniformity( 0.0 ); + RealType entropy( 0.0 ); + RealType upp( 0.0 ); + + LabelPixelType i; + if ( m_MaskNonEmpty ) + { + typename std::list< int >::const_iterator it; + for ( it = m_RelevantLabels.cbegin(), i = 0; + it != m_RelevantLabels.cend(); + ++it, ++i ) + { + HistogramType::Pointer histogramForEntropy = this->GetHistogram(*it); + for (int i = 0; i < histogramForEntropy->Size(); i++) + { + partialProbability = histogramForEntropy->GetFrequency(i,0) / double ( histogramForEntropy->GetTotalFrequency() ) ; + + if( partialProbability != 0) + { + entropy -= partialProbability *( std::log10(partialProbability) / std::log10(2) ) ; + uniformity += std::pow(partialProbability,2); + + + if(histogramForEntropy->GetMeasurement(i,0) > 0) + { + upp += std::pow(partialProbability,2); + } + } + } + m_LabelStatisticsCoefficients[*it].m_Entropy = entropy; + m_LabelStatisticsCoefficients[*it].m_Uniformity = uniformity; + m_LabelStatisticsCoefficients[*it].m_UPP = upp; + } + } + } + + template< class TInputImage, class TLabelImage > + void + ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: + ComputeSkewnessKurtosisAndMPP() + { + typename TLabelImage::RegionType Subregion; + + RealType baseOfSkewnessAndCurtosis( 0.0 ); + RealType kurtosis( 0.0 ); + RealType skewness( 0.0 ); + RealType mpp( 0.0 ); + RealType currentPixel( 0.0 ); + + std::list< LabelPixelType> relevantLabels; + LabelPixelType i; + if ( m_MaskNonEmpty ) + { + typename std::list< int >::const_iterator it; + for ( it = m_RelevantLabels.cbegin(), i = 0; + it != m_RelevantLabels.cend(); + ++it ) + { + RealType sigma = this->GetSigma( *it ); + RealType mean = this->GetMean( *it ); + Subregion = Superclass::GetRegion(*it); + + int count( this->GetCount(*it) ); + if ( count == 0 || sigma < mitk::eps) + { + throw std::logic_error( "Empty segmentation" ); + } + + if ( fabs( sigma ) < mitk::sqrteps ) + { + throw std::logic_error( "Sigma == 0" ); + } + + + ImageRegionConstIteratorWithIndex< TInputImage > it1 (this->GetInput(), + Subregion); + ImageRegionConstIterator< TLabelImage > labelIt (this->GetLabelInput(), + Subregion); + + for (it1.GoToBegin(); !it1.IsAtEnd(); ++it1, ++labelIt) + { + if (labelIt.Get() == *it) + { + currentPixel = it1.Get(); + baseOfSkewnessAndCurtosis = (currentPixel -mean) / sigma; + kurtosis += std::pow( baseOfSkewnessAndCurtosis, 4.0 ); + skewness += std::pow( baseOfSkewnessAndCurtosis, 3.0 ); + + if(currentPixel > 0) + { + mpp+= currentPixel; + } + + } + } + m_LabelStatisticsCoefficients[*it].m_Skewness = RealType(skewness/count); + m_LabelStatisticsCoefficients[*it].m_Kurtosis = RealType(kurtosis/count); + m_LabelStatisticsCoefficients[*it].m_MPP = RealType(mpp/count); + } + } + } + + + + template< class TInputImage, class TLabelImage > + void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: + AfterThreadedGenerateData() + { + Superclass::AfterThreadedGenerateData(); + + CalculateSettingsForLabels(); + + ComputeSkewnessKurtosisAndMPP(); + + if(this->GetUseHistograms()) + { + ComputeEntropyUniformityAndUPP(); + } + else + { + MITK_WARN << "Cannot compute coefficients UPP,Entropy,Uniformity because of missing histogram"; + } + } + +} // end namespace itk + +#endif diff --git a/Modules/ImageStatistics/mitkHotspotCalculator.cpp b/Modules/ImageStatistics/mitkHotspotCalculator.cpp new file mode 100644 index 0000000000..e385b36903 --- /dev/null +++ b/Modules/ImageStatistics/mitkHotspotCalculator.cpp @@ -0,0 +1,524 @@ +#include +#include +#include +#include +#include +#include "mitkImageAccessByItk.h" +#include +#include +#include + +namespace mitk +{ + HotspotCalculator::HotspotCalculator(): + m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm + m_HotspotMustBeCompletelyInsideImage(true), + m_HotspotParamsChanged(true) + { + } + + void HotspotCalculator::setInputImage(mitk::Image::Pointer inputImage) + { + m_Image = inputImage; + } + + void HotspotCalculator::setMaskImage(mitk::Image::Pointer maskImage) + { + m_Mask = maskImage; + } + + HotspotCalculator::~HotspotCalculator() + { + } + + void HotspotCalculator::setHotspotRadiusInMM(double radiusInMillimeter) + { + if(radiusInMillimeter != m_HotspotRadiusinMM) + { + m_HotspotRadiusinMM = radiusInMillimeter; + m_HotspotParamsChanged = true; + } + } + + const double& HotspotCalculator::getHotspotRadiusinMM() const + { + return m_HotspotRadiusinMM; + } + + bool HotspotCalculator::getHotspotMustBeCompletelyInsideImage() const + { + return m_HotspotMustBeCompletelyInsideImage; + } + + mitk::Image::Pointer HotspotCalculator::getHotspotMask(const unsigned int timeStep, + unsigned int label) + { + if ( m_Image.IsNull() ) + { + throw std::runtime_error( "Error: image empty!" ); + } + + if ( timeStep >= m_Image->GetTimeSteps() ) + { + throw std::runtime_error( "Error: invalid time step!" ); + } + + mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); + imageTimeSelector->SetInput( m_Image ); + imageTimeSelector->SetTimeNr( timeStep ); + imageTimeSelector->UpdateLargestPossibleRegion(); + mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); + + m_internalImage = timeSliceImage; + m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? + m_internalMask3D = nullptr; + + if ( m_Mask != nullptr ) + { + unsigned int timeStep_mask = timeStep; + if ( timeStep >= m_Mask->GetTimeSteps() ) + { + timeStep_mask = m_Mask->GetTimeSteps(); + } + + mitk::ImageTimeSelector::Pointer maskTimeSelector = mitk::ImageTimeSelector::New(); + maskTimeSelector->SetInput( m_Mask ); + maskTimeSelector->SetTimeNr( timeStep_mask ); + maskTimeSelector->UpdateLargestPossibleRegion(); + mitk::Image::Pointer timeSliceMask = maskTimeSelector->GetOutput(); + + if ( m_internalImage->GetDimension() == 3 ) + { + CastToItkImage(timeSliceMask, m_internalMask3D); + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, label); + } + else if ( m_internalImage->GetDimension() == 2 ) + { + CastToItkImage(timeSliceMask, m_internalMask2D); + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, label); + } + else + { + throw std::runtime_error( "Error: invalid image dimension" ); + } + } + else + { + + if ( m_internalImage->GetDimension() == 3 ) + { + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, label); + } + else if ( m_internalImage->GetDimension() == 2 ) + { + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, label); + } + else + { + throw std::runtime_error( "Error: invalid image dimension" ); + } + } + + return m_hotspotMask; + } + + template + HotspotCalculator::ImageExtrema + HotspotCalculator::CalculateExtremaWorld( const itk::Image* inputImage, + typename itk::Image::Pointer maskImage, + double neccessaryDistanceToImageBorderInMM, + unsigned int label ) + { + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskImageType; + + typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; + typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; + + typename ImageType::SpacingType spacing = inputImage->GetSpacing(); + + ImageExtrema minMax; + minMax.Defined = false; + minMax.MaxIndex.set_size(VImageDimension); + minMax.MaxIndex.set_size(VImageDimension); + + typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); + + bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); + if (keepDistanceToImageBorders) + { + long distanceInPixels[VImageDimension]; + for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) + { + // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as + // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: + // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. + // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 + distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); + } + + allowedExtremaRegion.ShrinkByRadius(distanceInPixels); + } + + InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); + + float maxValue = itk::NumericTraits::min(); + float minValue = itk::NumericTraits::max(); + + typename ImageType::IndexType maxIndex; + typename ImageType::IndexType minIndex; + + for(unsigned short i = 0; i < VImageDimension; ++i) + { + maxIndex[i] = 0; + minIndex[i] = 0; + } + + if (maskImage != nullptr) + { + MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); + typename ImageType::IndexType imageIndex; + typename ImageType::PointType worldPosition; + typename ImageType::IndexType maskIndex; + + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + imageIndex = maskIndex = maskIt.GetIndex(); + + if(maskIt.Get() == label) + { + if( allowedExtremaRegion.IsInside(imageIndex) ) + { + imageIndexIt.SetIndex( imageIndex ); + double value = imageIndexIt.Get(); + minMax.Defined = true; + + //Calculate minimum, maximum and corresponding index-values + if( value > maxValue ) + { + maxIndex = imageIndexIt.GetIndex(); + maxValue = value; + } + + if(value < minValue ) + { + minIndex = imageIndexIt.GetIndex(); + minValue = value; + } + } + } + } + } + else + { + for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) + { + double value = imageIndexIt.Get(); + minMax.Defined = true; + + //Calculate minimum, maximum and corresponding index-values + if( value > maxValue ) + { + maxIndex = imageIndexIt.GetIndex(); + maxValue = value; + } + + if(value < minValue ) + { + minIndex = imageIndexIt.GetIndex(); + minValue = value; + } + } + } + + minMax.MaxIndex.set_size(VImageDimension); + minMax.MinIndex.set_size(VImageDimension); + + for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) + { + minMax.MaxIndex[i] = maxIndex[i]; + } + + for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) + { + minMax.MinIndex[i] = minIndex[i]; + } + + minMax.Max = maxValue; + minMax.Min = minValue; + + return minMax; + } + + template + itk::Size + HotspotCalculator::CalculateConvolutionKernelSize( double spacing[VImageDimension], + double radiusInMM ) + { + typedef itk::Image< float, VImageDimension > KernelImageType; + typedef typename KernelImageType::SizeType SizeType; + SizeType maskSize; + + for(unsigned int i = 0; i < VImageDimension; ++i) + { + maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); + + // We always want an uneven size to have a clear center point in the convolution mask + if(maskSize[i] % 2 == 0 ) + { + ++maskSize[i]; + } + } + return maskSize; + } + + template + itk::SmartPointer< itk::Image > + HotspotCalculator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], + double radiusInMM ) + { + std::stringstream ss; + for (unsigned int i = 0; i < VImageDimension; ++i) + { + ss << mmPerPixel[i]; + if (i < VImageDimension -1) + ss << ","; + } + MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; + + + double radiusInMMSquared = radiusInMM * radiusInMM; + typedef itk::Image< float, VImageDimension > KernelImageType; + typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); + + // Calculate size and allocate mask image + typedef typename KernelImageType::SizeType SizeType; + SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); + + mitk::Point3D convolutionMaskCenterIndex; + convolutionMaskCenterIndex.Fill(0.0); + for(unsigned int i = 0; i < VImageDimension; ++i) + { + convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); + } + + typedef typename KernelImageType::IndexType IndexType; + IndexType maskIndex; + maskIndex.Fill(0); + + typedef typename KernelImageType::RegionType RegionType; + RegionType maskRegion; + maskRegion.SetSize(maskSize); + maskRegion.SetIndex(maskIndex); + + convolutionKernel->SetRegions(maskRegion); + convolutionKernel->SetSpacing(mmPerPixel); + convolutionKernel->Allocate(); + + // Fill mask image values by subsampling the image grid + typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; + MaskIteratorType maskIt(convolutionKernel,maskRegion); + + int numberOfSubVoxelsPerDimension = 2; // per dimension! + int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); + double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; + double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; + double maskValue = 0.0; + mitk::Point3D subVoxelIndexPosition; + double distanceSquared = 0.0; + + typedef itk::ContinuousIndex ContinuousIndexType; + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + ContinuousIndexType indexPoint(maskIt.GetIndex()); + mitk::Point3D voxelPosition; + for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) + { + voxelPosition[dimension] = indexPoint[dimension]; + } + + maskValue = 0.0; + mitk::Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); + // iterate sub-voxels by iterating all possible offsets + for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; + subVoxelOffset[0] < +0.5; + subVoxelOffset[0] += subVoxelSizeInPixels) + { + for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; + subVoxelOffset[1] < +0.5; + subVoxelOffset[1] += subVoxelSizeInPixels) + { + for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; + subVoxelOffset[2] < +0.5; + subVoxelOffset[2] += subVoxelSizeInPixels) + { + subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) + distanceSquared = + (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; + + if (distanceSquared <= radiusInMMSquared) + { + maskValue += valueOfOneSubVoxel; + } + } + } + } + maskIt.Set( maskValue ); + } + + return convolutionKernel; + } + + template + itk::SmartPointer > + HotspotCalculator::GenerateConvolutionImage( const itk::Image* inputImage ) + { + double mmPerPixel[VImageDimension]; + for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) + { + mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; + } + + // update convolution kernel + typedef itk::Image< float, VImageDimension > KernelImageType; + typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusinMM); + + // update convolution image + typedef itk::Image< TPixel, VImageDimension > InputImageType; + typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; + typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; + + typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); + typedef itk::ConstantBoundaryCondition BoundaryConditionType; + BoundaryConditionType boundaryCondition; + boundaryCondition.SetConstant(0.0); + + if (m_HotspotMustBeCompletelyInsideImage) + { + // overwrite default boundary condition + convolutionFilter->SetBoundaryCondition(&boundaryCondition); + } + + convolutionFilter->SetInput(inputImage); + convolutionFilter->SetKernelImage(convolutionKernel); + convolutionFilter->SetNormalize(true); + MITK_DEBUG << "Update Convolution image for hotspot search"; + convolutionFilter->UpdateLargestPossibleRegion(); + + typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); + convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image + + m_modified = false; + return convolutionImage; + } + + template < typename TPixel, unsigned int VImageDimension> + void + HotspotCalculator::FillHotspotMaskPixels( itk::Image* maskImage, + itk::Point sphereCenter, + double sphereRadiusInMM ) + { + typedef itk::Image< TPixel, VImageDimension > MaskImageType; + typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; + + MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); + + typename MaskImageType::IndexType maskIndex; + typename MaskImageType::PointType worldPosition; + + // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + maskIndex = maskIt.GetIndex(); + maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); + maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); + } + } + + template + void + HotspotCalculator::CalculateHotspotMask(itk::Image* inputImage, + typename itk::Image::Pointer maskImage, + unsigned int label) + { + typedef itk::Image< TPixel, VImageDimension > InputImageType; + typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; + typedef itk::Image< float, VImageDimension > KernelImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskImageType; + typedef itk::ImageDuplicator< InputImageType > DuplicatorType; + typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType; + + typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); + + if (convolutionImage.IsNull()) + { + MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; + throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); + } + + // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's + // there is maybe a better way to do this!? + if (maskImage == nullptr) + { + maskImage = MaskImageType::New(); + typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); + typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); + typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); + typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); + maskImage->SetRegions(maskRegion); + maskImage->Allocate(); + maskImage->SetOrigin(maskOrigin); + maskImage->SetSpacing(maskSpacing); + maskImage->SetDirection(maskDirection); + + maskImage->FillBuffer(1); + + label = 1; + } + + // find maximum in convolution image, given the current mask + double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusinMM : -1.0; + ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); + + bool isHotspotDefined = convolutionImageInformation.Defined; + + if (!isHotspotDefined) + { + MITK_ERROR << "No origin of hotspot-sphere was calculated!"; + m_hotspotMask = nullptr; + } + else + { + // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center + typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); + copyMachine->SetInputImage(inputImage); + copyMachine->Update(); + + typename CastFilterType::Pointer caster = CastFilterType::New(); + caster->SetInput( copyMachine->GetOutput() ); + caster->Update(); + typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput(); + + typedef typename InputImageType::IndexType IndexType; + IndexType maskCenterIndex; + for (unsigned int d =0; d< VImageDimension;++d) + { + maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; + } + + typename ConvolutionImageType::PointType maskCenter; + inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); + + FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM); + + //obtain mitk::Image::Pointer from itk::Image + mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); + + m_hotspotMask = hotspotMaskAsMITKImage; + } + } +} diff --git a/Modules/ImageStatistics/mitkHotspotCalculator.h b/Modules/ImageStatistics/mitkHotspotCalculator.h new file mode 100644 index 0000000000..1dfd7514b6 --- /dev/null +++ b/Modules/ImageStatistics/mitkHotspotCalculator.h @@ -0,0 +1,110 @@ +#ifndef MITKHOTSPOTCALCULATOR_H +#define MITKHOTSPOTCALCULATOR_H + +#include +#include +#include +#include +#include +#include +#include + + +namespace mitk +{ + class MITKIMAGESTATISTICS_EXPORT HotspotCalculator: itk::Object + { + public: + HotspotCalculator(); + + ~HotspotCalculator(); + + void setInputImage(mitk::Image::Pointer inputImage); + + void setMaskImage(mitk::Image::Pointer maskImage); + + void setHotspotRadiusInMM(double radiusInMillimeter); + + const double& getHotspotRadiusinMM() const; + + void setHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage); + + bool getHotspotMustBeCompletelyInsideImage() const; + + mitk::Image::Pointer getHotspotMask(unsigned int timeStep, unsigned int label=0); + + protected: + class ImageExtrema + { + public: + bool Defined; + double Max; + double Min; + vnl_vector MaxIndex; + vnl_vector MinIndex; + + ImageExtrema() + :Defined(false) + ,Max(itk::NumericTraits::min()) + ,Min(itk::NumericTraits::max()) + { + } + }; + + private: + /** \brief Returns size of convolution kernel depending on spacing and radius. */ + template + itk::Size + CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); + + /** \brief Generates image of kernel which is needed for convolution. */ + template + itk::SmartPointer< itk::Image > + GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); + + /** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */ + template + itk::SmartPointer< itk::Image > + GenerateConvolutionImage( const itk::Image* inputImage ); + + + /** \brief Fills pixels of the spherical hotspot mask. */ + template < typename TPixel, unsigned int VImageDimension> + void + FillHotspotMaskPixels( itk::Image* maskImage, + itk::Point sphereCenter, + double sphereRadiusInMM); + + + /** \brief */ + template + void + CalculateHotspotMask(itk::Image* inputImage, + typename itk::Image::Pointer maskImage, + unsigned int label); + + + template + ImageExtrema CalculateExtremaWorld( const itk::Image* inputImage, + typename itk::Image::Pointer maskImage, + double neccessaryDistanceToImageBorderInMM, + unsigned int label); + + HotspotCalculator(const HotspotCalculator &); + HotspotCalculator & operator=(const HotspotCalculator &); + + mitk::Image::Pointer m_Image; + mitk::Image::Pointer m_Mask = nullptr; + mitk::Image::Pointer m_internalImage; + itk::Image::Pointer m_internalMask2D; + itk::Image::Pointer m_internalMask3D; + mitk::Image::Pointer m_hotspotMask; + double m_HotspotRadiusinMM; + bool m_HotspotMustBeCompletelyInsideImage; + bool m_HotspotParamsChanged; + bool m_modified = false; + }; +} +#endif // MITKHOTSPOTCALCULATOR + + diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index f55900fafc..03128441fc 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,2240 +1,2246 @@ /*=================================================================== 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 "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkExtractImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkITKImageImport.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkImage.h" //#define DEBUG_HOTSPOTSEARCH #define _USE_MATH_DEFINES #include #include "vtkLassoStencilSource.h" namespace mitk { ImageStatisticsCalculator::ImageStatisticsCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_IgnorePixelValue(0.0), m_DoIgnorePixelValue(false), m_IgnorePixelValueChanged(false), m_PlanarFigureAxis (0), m_PlanarFigureSlice (0), m_PlanarFigureCoordinate0 (0), m_PlanarFigureCoordinate1 (0), m_HistogramBinSize(1.0), m_UseDefaultBinSize(true), m_UseBinSizeBasedOnVOIRegion(false), m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_CalculateHotspot(false), m_HotspotRadiusInMMChanged(false), m_HotspotMustBeCompletelyInsideImage(true) { m_EmptyHistogram = HistogramType::New(); m_EmptyHistogram->SetMeasurementVectorSize(1); HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } ImageStatisticsCalculator::~ImageStatisticsCalculator() { } void ImageStatisticsCalculator::SetUseDefaultBinSize(bool useDefault) { m_UseDefaultBinSize = useDefault; } ImageStatisticsCalculator::Statistics::Statistics(bool withHotspotStatistics) :m_HotspotStatistics(withHotspotStatistics ? new Statistics(false) : nullptr) { Reset(); } ImageStatisticsCalculator::Statistics::Statistics(const Statistics& other) :m_HotspotStatistics( nullptr) { this->SetLabel( other.GetLabel() ); this->SetN( other.GetN() ); this->SetMin( other.GetMin() ); this->SetMax( other.GetMax() ); this->SetMedian( other.GetMedian() ); this->SetMean( other.GetMean() ); this->SetVariance( other.GetVariance() ); this->SetKurtosis( other.GetKurtosis() ); this->SetSkewness( other.GetSkewness() ); this->SetUniformity( other.GetUniformity() ); this->SetEntropy( other.GetEntropy() ); this->SetUPP( other.GetUPP() ); this->SetMPP( other.GetMPP() ); this->SetSigma( other.GetSigma() ); this->SetRMS( other.GetRMS() ); this->SetMaxIndex( other.GetMaxIndex() ); this->SetMinIndex( other.GetMinIndex() ); this->SetHotspotIndex( other.GetHotspotIndex() ); if (other.m_HotspotStatistics) { this->m_HotspotStatistics = new Statistics(false); *this->m_HotspotStatistics = *other.m_HotspotStatistics; } } bool ImageStatisticsCalculator::Statistics::HasHotspotStatistics() const { return m_HotspotStatistics != nullptr; } void ImageStatisticsCalculator::Statistics::SetHasHotspotStatistics(bool hasHotspotStatistics) { m_HasHotspotStatistics = hasHotspotStatistics; } ImageStatisticsCalculator::Statistics::~Statistics() { delete m_HotspotStatistics; } double ImageStatisticsCalculator::Statistics::GetVariance() const { return this->Variance; } void ImageStatisticsCalculator::Statistics::SetVariance( const double value ) { if( this->Variance != value ) { if( value < 0.0 ) { this->Variance = 0.0; // if given value is negative set variance to 0.0 } else { this->Variance = value; } } } double ImageStatisticsCalculator::Statistics::GetSigma() const { return this->Sigma; } void ImageStatisticsCalculator::Statistics::SetSigma( const double value ) { if( this->Sigma != value ) { // for some compiler the value != value works to check for NaN but not for all // but we can always be sure that the standard deviation is a positive value if( value != value || value < 0.0 ) { // if standard deviation is NaN we just assume 0.0 this->Sigma = 0.0; } else { this->Sigma = value; } } } void ImageStatisticsCalculator::Statistics::Reset(unsigned int dimension) { SetLabel(0); SetN( 0 ); SetMin( 0.0 ); SetMax( 0.0 ); SetMedian( 0.0 ); SetVariance( 0.0 ); SetMean( 0.0 ); SetSigma( 0.0 ); SetRMS( 0.0 ); SetSkewness( 0.0 ); SetKurtosis( 0.0 ); SetUniformity( 0.0 ); SetEntropy( 0.0 ); SetMPP( 0.0 ); SetUPP( 0.0 ); vnl_vector zero; zero.set_size(dimension); for(unsigned int i = 0; i < dimension; ++i) { zero[i] = 0; } SetMaxIndex(zero); SetMinIndex(zero); SetHotspotIndex(zero); if (m_HotspotStatistics != nullptr) { m_HotspotStatistics->Reset(dimension); } } const ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::GetHotspotStatistics() const { if (m_HotspotStatistics) { return *m_HotspotStatistics; } else { throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()"); } } ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::GetHotspotStatistics() { if (m_HotspotStatistics) { return *m_HotspotStatistics; } else { throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()"); } } ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::operator=(ImageStatisticsCalculator::Statistics const& other) { if (this == &other) return *this; this->SetLabel( other.GetLabel() ); this->SetN( other.GetN() ); this->SetMin( other.GetMin() ); this->SetMax( other.GetMax() ); this->SetMean( other.GetMean() ); this->SetMedian( other.GetMedian() ); this->SetVariance( other.GetVariance() ); this->SetSigma( other.GetSigma() ); this->SetRMS( other.GetRMS() ); this->SetMinIndex( other.GetMinIndex() ); this->SetMaxIndex( other.GetMaxIndex() ); this->SetHotspotIndex( other.GetHotspotIndex() ); this->SetSkewness( other.GetSkewness() ); this->SetKurtosis( other.GetKurtosis() ); this->SetUniformity( other.GetUniformity() ); this->SetEntropy( other.GetEntropy() ); this->SetUPP( other.GetUPP() ); this->SetMPP( other.GetMPP() ); delete this->m_HotspotStatistics; this->m_HotspotStatistics = nullptr; if (other.m_HotspotStatistics) { this->m_HotspotStatistics = new Statistics(false); *this->m_HotspotStatistics = *other.m_HotspotStatistics; } return *this; } void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); unsigned int numberOfTimeSteps = image->GetTimeSteps(); // Initialize vectors to time-size of this image m_ImageHistogramVector.resize( numberOfTimeSteps ); m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); m_ImageStatisticsVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) { m_ImageStatisticsTimeStampVector[t].Modified(); m_ImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_MaskedImageStatisticsTimeStampVector[t].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_PlanarFigureStatisticsTimeStampVector[t].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetIgnorePixelValue(double value) { if ( m_IgnorePixelValue != value ) { m_IgnorePixelValue = value; if(m_DoIgnorePixelValue) { m_IgnorePixelValueChanged = true; } this->Modified(); } } double ImageStatisticsCalculator::GetIgnorePixelValue() { return m_IgnorePixelValue; } void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) { if ( m_DoIgnorePixelValue != value ) { m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; this->Modified(); } } bool ImageStatisticsCalculator::GetDoIgnorePixelValue() { return m_DoIgnorePixelValue; } void ImageStatisticsCalculator::SetHistogramBinSize(double size) { this->m_HistogramBinSize = size; } double ImageStatisticsCalculator::GetHistogramBinSize() { return this->m_HistogramBinSize; } void ImageStatisticsCalculator::SetHotspotRadiusInMM(double value) { if ( m_HotspotRadiusInMM != value ) { m_HotspotRadiusInMM = value; if(m_CalculateHotspot) { m_HotspotRadiusInMMChanged = true; //MITK_INFO <<"Hotspot radius changed, new convolution required"; } this->Modified(); } } double ImageStatisticsCalculator::GetHotspotRadiusInMM() { return m_HotspotRadiusInMM; } void ImageStatisticsCalculator::SetCalculateHotspot(bool on) { if ( m_CalculateHotspot != on ) { m_CalculateHotspot = on; m_HotspotRadiusInMMChanged = true; //MITK_INFO <<"Hotspot calculation changed, new convolution required"; this->Modified(); } } bool ImageStatisticsCalculator::IsHotspotCalculated() { return m_CalculateHotspot; } void ImageStatisticsCalculator::SetHotspotMustBeCompletlyInsideImage(bool hotspotMustBeCompletelyInsideImage, bool warn) { m_HotspotMustBeCompletelyInsideImage = hotspotMustBeCompletelyInsideImage; if (!m_HotspotMustBeCompletelyInsideImage && warn) { MITK_WARN << "Hotspot calculation will extrapolate pixels at image borders. Be aware of the consequences for the hotspot location."; } } bool ImageStatisticsCalculator::GetHotspotMustBeCompletlyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } /* Implementation of the min max values for setting the range of the histogram */ template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::GetMinAndMaxValue( double &min, double &max, int &counter, double &sigma, const itk::Image< TPixel, VImageDimension > *InputImage, itk::Image< unsigned short, VImageDimension > *MaskedImage ) { typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::ImageRegionConstIteratorWithIndex Imageie; typedef itk::ImageRegionConstIteratorWithIndex Imageie2; Imageie2 labelIterator2( MaskedImage, MaskedImage->GetRequestedRegion() ); Imageie labelIterator3( InputImage, InputImage->GetRequestedRegion() ); max = 0; min = 0; counter = 0; sigma = 0; double SumOfSquares = 0; double sumSquared = 0; double actualPielValue = 0; int counterOfPixelsInROI = 0; for( labelIterator2.GoToBegin(); !labelIterator2.IsAtEnd(); ++labelIterator2, ++labelIterator3) { if( labelIterator2.Value()== 1.0) { counter++; counterOfPixelsInROI++; actualPielValue = labelIterator3.Value(); sumSquared = sumSquared + actualPielValue; SumOfSquares = SumOfSquares + std::pow(actualPielValue,2); if(counterOfPixelsInROI == 1) { max = actualPielValue; min = actualPielValue; } if(actualPielValue >= max) { max = actualPielValue; } else if(actualPielValue <= min) { min = actualPielValue; } } } if (counter > 1) { sigma = ( SumOfSquares - std::pow( sumSquared, 2) / counter ) / ( counter-1 ); } else { sigma = 0; } } bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) { if (m_Image.IsNull() ) { mitkThrow() << "Image not set!"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if ( m_Image->GetReferenceCount() == 1 ) { // Image no longer valid; we are the only ones to still hold a reference on it return false; } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = nullptr; } // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; if ( !m_IgnorePixelValueChanged && !m_HotspotRadiusInMMChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; m_IgnorePixelValueChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( timeStep ); StatisticsContainer *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } ImageStatisticsCalculator::BinFrequencyType ImageStatisticsCalculator::GetBinsAndFreuqencyForHistograms( unsigned int timeStep , unsigned int label ) const { const HistogramType *binsAndFrequencyToCalculate = this->GetHistogram(0); // ToDo: map should be created on stack not on heap std::map returnedHistogramMap; unsigned int size = binsAndFrequencyToCalculate->Size(); for( unsigned int bin=0; bin < size; ++bin ) { double frequency = binsAndFrequencyToCalculate->GetFrequency( bin, 0 ); //if( frequency > mitk::eps ) { returnedHistogramMap.insert( std::pair(binsAndFrequencyToCalculate->GetMeasurement( bin, 0 ), binsAndFrequencyToCalculate->GetFrequency( bin, 0 ) ) ); } } return returnedHistogramMap; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return nullptr; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep][label]; return m_ImageHistogramVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep][label]; } } const ImageStatisticsCalculator::HistogramContainer & ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyHistogramContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep]; return m_ImageHistogramVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep][label]; return m_ImageStatisticsVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep][label]; } } const ImageStatisticsCalculator::StatisticsContainer & ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatisticsContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep]; return m_ImageStatisticsVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep]; } } void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( timeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; m_InternalImageMask2D = nullptr; m_InternalImageMask3D = nullptr; if(m_DoIgnorePixelValue) { if( m_InternalImage->GetDimension() == 3 ) { if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType()) CastToItkImage( timeSliceImage, m_InternalImageMask3D ); else CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } if( m_InternalImage->GetDimension() == 2 ) { if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType()) CastToItkImage( timeSliceImage, m_InternalImageMask2D ); else CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask2D ); m_InternalImageMask2D->FillBuffer(1); } } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { if ( timeStep >= m_ImageMask->GetTimeSteps() ) { // Use the last mask time step in case the current time step is bigger than the total // number of mask time steps. // It makes more sense setting this to the last mask time step than to 0. // For instance if you have a mask with 2 time steps and an image with 5: // If time step 0 is selected, the mask will use time step 0. // If time step 1 is selected, the mask will use time step 1. // If time step 2+ is selected, the mask will use time step 1. // If you have a mask with only one time step instead, this will always default to 0. timeStep = m_ImageMask->GetTimeSteps() - 1; } ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( timeStep ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); m_InternalImage = timeSliceImage; CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D ); } else { throw std::runtime_error( "Error: image mask empty!" ); } break; } case MASKING_MODE_PLANARFIGURE: { m_InternalImageMask2D = nullptr; if ( m_PlanarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !m_PlanarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const BaseGeometry *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == nullptr ) { throw std::runtime_error( "Image geometry invalid!" ); } const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == nullptr ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); if ( planarFigureGeometry == nullptr ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image MaskImage3DType::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // Extract slice with given position and direction from image unsigned int dimension = timeSliceImage->GetDimension(); if (dimension != 2) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( timeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); m_InternalImage = imageExtractor->GetOutput(); } else { m_InternalImage = timeSliceImage; } // Compute mask from PlanarFigure AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 2, axis ); } } if(m_DoIgnorePixelValue) { if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); } } } bool ImageStatisticsCalculator::GetPrincipalAxis( const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } unsigned int ImageStatisticsCalculator::calcNumberOfBins(mitk::ScalarType min, mitk::ScalarType max) { return std::ceil( ( (max - min ) / m_HistogramBinSize) ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer *statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef typename ImageType::IndexType IndexType; typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; statisticsContainer->clear(); histogramContainer->clear(); // Progress listening... typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator // does not (yet?) support progress reporting this->InvokeEvent( itk::StartEvent() ); for ( unsigned int i = 0; i < 100; ++i ) { this->UnmaskedStatisticsProgressUpdate(); } // Calculate statistics (separate filter) typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); statisticsFilter->SetBinSize( 100 ); statisticsFilter->SetCoordinateTolerance( 0.001 ); statisticsFilter->SetDirectionTolerance( 0.001 ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } statisticsFilter->RemoveObserver( observerTag ); this->InvokeEvent( itk::EndEvent() ); // Calculate minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( image ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); Statistics statistics; statistics.Reset(); statistics.SetLabel(1); statistics.SetN(image->GetBufferedRegion().GetNumberOfPixels()); statistics.SetMin(statisticsFilter->GetMinimum()); statistics.SetMax(statisticsFilter->GetMaximum()); statistics.SetMean(statisticsFilter->GetMean()); statistics.SetMedian(statisticsFilter->GetMedian()); statistics.SetVariance(statisticsFilter->GetVariance()); statistics.SetSkewness(statisticsFilter->GetSkewness()); statistics.SetKurtosis(statisticsFilter->GetKurtosis()); statistics.SetUniformity( statisticsFilter->GetUniformity()); statistics.SetEntropy( statisticsFilter->GetEntropy()); statistics.SetUPP( statisticsFilter->GetUPP()); statistics.SetMPP( statisticsFilter->GetMPP()); statistics.SetSigma(statisticsFilter->GetSigma()); statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() )); statistics.GetMinIndex().set_size(image->GetImageDimension()); statistics.GetMaxIndex().set_size(image->GetImageDimension()); vnl_vector tmpMaxIndex; vnl_vector tmpMinIndex; tmpMaxIndex.set_size(image->GetImageDimension() ); tmpMinIndex.set_size(image->GetImageDimension() ); for (unsigned int i=0; iGetIndexOfMaximum()[i]; tmpMinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } statistics.SetMinIndex(tmpMaxIndex); statistics.SetMinIndex(tmpMinIndex); if( IsHotspotCalculated() && VImageDimension == 3 ) { typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typename MaskImageType::Pointer nullMask; bool isHotspotDefined(false); Statistics hotspotStatistics = this->CalculateHotspotStatistics(image, nullMask.GetPointer(), m_HotspotRadiusInMM, isHotspotDefined, 0 ); if (isHotspotDefined) { statistics.SetHasHotspotStatistics(true); statistics.GetHotspotStatistics() = hotspotStatistics; } else { statistics.SetHasHotspotStatistics(false); } if(statistics.GetHotspotStatistics().HasHotspotStatistics() ) { MITK_DEBUG << "Hotspot statistics available"; statistics.SetHotspotIndex(hotspotStatistics.GetHotspotIndex()); } else { MITK_ERROR << "No hotspot statistics available!"; } } statisticsContainer->push_back( statistics ); // Calculate histogram // calculate bin size or number of bins unsigned int numberOfBins = 200; // default number of bins if (m_UseDefaultBinSize) { m_HistogramBinSize = std::ceil( (statistics.GetMax() - statistics.GetMin() + 1)/numberOfBins ); } else { numberOfBins = calcNumberOfBins(statistics.GetMin(), statistics.GetMax()); } typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); histogramGenerator->SetNumberOfBins( numberOfBins ); histogramGenerator->SetHistogramMin( statistics.GetMin() ); histogramGenerator->SetHistogramMax( statistics.GetMax() ); histogramGenerator->Compute(); histogramContainer->push_back( histogramGenerator->GetOutput() ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(m_IgnorePixelValue == itimage.Get()) { itmask.Set(0); } ++itmask; ++itimage; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::SpacingType SpacingType; typedef typename ImageType::Pointer ImagePointer; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; statisticsContainer->clear(); histogramContainer->clear(); // Make sure that mask is set if ( maskImage == nullptr ) { itkExceptionMacro( << "Mask image needs to be set!" ); } // Make sure that spacing of mask and image are the same //SpacingType imageSpacing = image->GetSpacing(); //SpacingType maskSpacing = maskImage->GetSpacing(); //PointType zeroPoint; zeroPoint.Fill( 0.0 ); //if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps ) //{ // itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" ); //} // Make sure that orientation of mask and image are the same typedef typename ImageType::DirectionType DirectionType; DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > 0.001 /*mitk::eps*/ ) // TODO: temp fix (bug 17121) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } } // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); long offset[ImageType::ImageDimension]; typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); if ( fabs( misalignment ) > mitk::eps ) { itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; offset[i] = int( indexCoordDistance + image->GetBufferedRegion().GetIndex()[i] + 0.5 ); } // Adapt the origin and region (index/size) of the mask so that the origin of both are the same typename ChangeInformationFilterType::Pointer adaptMaskFilter; adaptMaskFilter = ChangeInformationFilterType::New(); adaptMaskFilter->ChangeOriginOn(); adaptMaskFilter->ChangeRegionOn(); adaptMaskFilter->SetInput( maskImage ); adaptMaskFilter->SetOutputOrigin( image->GetOrigin() ); adaptMaskFilter->SetOutputOffset( offset ); adaptMaskFilter->SetCoordinateTolerance( 0.001 ); adaptMaskFilter->SetDirectionTolerance( 0.001 ); typename MaskImageType::Pointer adaptedMaskImage; try { adaptMaskFilter->Update(); adaptedMaskImage = adaptMaskFilter->GetOutput(); } catch( const itk::ExceptionObject &e) { mitkThrow() << "Attempt to adapt shifted origin of the mask image failed due to ITK Exception: \n" << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // Make sure that mask region is contained within image region if ( adaptedMaskImage.IsNotNull() && !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) { itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); } // If mask region is smaller than image region, extract the sub-sampled region from the original image typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } typename ImageType::ConstPointer adaptedImage; if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); extractImageFilter->SetCoordinateTolerance( 0.001 ); extractImageFilter->SetDirectionTolerance( 0.001 ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); try { statisticsFilter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Image statistics initialization computation failed with ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // Calculate bin size or number of bins unsigned int numberOfBins = 200; // default number of bins double maximum = 0.0; double minimum = 0.0; if (m_UseBinSizeBasedOnVOIRegion) { maximum = statisticsFilter->GetMaximum(); minimum = statisticsFilter->GetMinimum(); if (m_UseDefaultBinSize) { m_HistogramBinSize = std::ceil( static_cast((statisticsFilter->GetMaximum() - statisticsFilter->GetMinimum() + 1)/numberOfBins) ); } else { numberOfBins = calcNumberOfBins(statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum()); } } else { double sig = 0.0; int counter = 0; //Find the min and max values for the Roi to set the range for the histogram GetMinAndMaxValue( minimum, maximum, counter, sig, image, maskImage); numberOfBins = maximum - minimum; if(maximum - minimum <= 10) { numberOfBins = 100; } } typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->SetCoordinateTolerance( 0.001 ); labelStatisticsFilter->SetDirectionTolerance( 0.001 ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( numberOfBins, floor(minimum), ceil(maximum) ); //statisticsFilter->GetMinimum() statisticsFilter->GetMaximum() // Add progress listening typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate ); unsigned long observerTag = labelStatisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); // Execute filter this->InvokeEvent( itk::StartEvent() ); // Make sure that only the mask region is considered (otherwise, if the mask region is smaller // than the image region, the Update() would result in an exception). labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); // Execute the filter try { labelStatisticsFilter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } this->InvokeEvent( itk::EndEvent() ); if( observerTag ) labelStatisticsFilter->RemoveObserver( observerTag ); // Find all relevant labels of mask (other than 0) std::list< int > relevantLabels = labelStatisticsFilter->GetRelevantLabels(); unsigned int i; if ( labelStatisticsFilter->GetMaskingNonEmpty() ) { std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { Statistics statistics; // restore previous code labelStatisticsFilter->GetHistogram(*it) ; histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); statistics.SetLabel (*it); statistics.SetN(labelStatisticsFilter->GetCount( *it )); statistics.SetMin(labelStatisticsFilter->GetMinimum( *it )); statistics.SetMax(labelStatisticsFilter->GetMaximum( *it )); statistics.SetMean(labelStatisticsFilter->GetMean( *it )); statistics.SetMedian(labelStatisticsFilter->GetMedian( *it)); statistics.SetMedian(labelStatisticsFilter->GetMedian( *it )); statistics.SetVariance(labelStatisticsFilter->GetVariance( *it )); statistics.SetSigma(labelStatisticsFilter->GetSigma( *it )); statistics.SetSkewness(labelStatisticsFilter->GetSkewness( *it )); statistics.SetKurtosis(labelStatisticsFilter->GetKurtosis( *it )); statistics.SetUniformity( labelStatisticsFilter->GetUniformity( *it )); statistics.SetEntropy( labelStatisticsFilter->GetEntropy( *it )); statistics.SetUPP( labelStatisticsFilter->GetUPP( *it)); statistics.SetMPP( labelStatisticsFilter->GetMPP( *it)); statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() )); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); bool isMinAndMaxSameValue = (statistics.GetMin() == statistics.GetMax()); // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here double outsideValue = (isMinAndMaxSameValue ? (statistics.GetMax()/2) : (statistics.GetMin()+statistics.GetMax())/2); masker->SetOutsideValue( outsideValue ); masker->SetInput1(adaptedImage); masker->SetInput2(adaptedMaskImage); masker->SetCoordinateTolerance( 0.001 ); masker->SetDirectionTolerance( 0.001 ); masker->Update(); // get index of minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( masker->GetOutput() ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here typename MinMaxFilterType::IndexType tempMinIndex = (isMinAndMaxSameValue ? minMaxFilter->GetIndexOfMaximum() : minMaxFilter->GetIndexOfMinimum()); // FIX BUG 14644 //If a PlanarFigure is used for segmentation the //adaptedImage is a single slice (2D). Adding the // 3. dimension. vnl_vector maxIndex; vnl_vector minIndex; maxIndex.set_size(m_Image->GetDimension()); minIndex.set_size(m_Image->GetDimension()); if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { maxIndex[m_PlanarFigureCoordinate0] = tempMaxIndex[0]; maxIndex[m_PlanarFigureCoordinate1] = tempMaxIndex[1]; maxIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; minIndex[m_PlanarFigureCoordinate0] = tempMinIndex[0] ; minIndex[m_PlanarFigureCoordinate1] = tempMinIndex[1]; minIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; } else { for (unsigned int i = 0; ipush_back( statistics ); } } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); statisticsContainer->push_back( Statistics() ); } } template ImageStatisticsCalculator::ImageExtrema ImageStatisticsCalculator::CalculateExtremaWorld( const itk::Image *inputImage, itk::Image *maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { long distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::PointType worldPosition; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size ImageStatisticsCalculator ::CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > ImageStatisticsCalculator ::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; double maskValue = 0.0; Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } maskValue = 0.0; Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > ImageStatisticsCalculator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusInMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (GetHotspotMustBeCompletlyInsideImage()) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image m_HotspotRadiusInMMChanged = false; return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void ImageStatisticsCalculator ::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template < typename TPixel, unsigned int VImageDimension> ImageStatisticsCalculator::Statistics ImageStatisticsCalculator::CalculateHotspotStatistics( const itk::Image* inputImage, itk::Image* maskImage, double radiusInMM, bool& isHotspotDefined, unsigned int label) { // get convolution image (updated in GenerateConvolutionImage()) typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< float, VImageDimension > KernelImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; //typename ConvolutionImageType::Pointer convolutionImage = dynamic_cast(this->GenerateConvolutionImage(inputImage)); typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } // find maximum in convolution image, given the current mask double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusInMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { m_EmptyStatistics.Reset(VImageDimension); MITK_ERROR << "No origin of hotspot-sphere was calculated! Returning empty statistics"; return m_EmptyStatistics; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center typedef itk::ImageDuplicator< InputImageType > DuplicatorType; typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); copyMachine->SetInputImage(inputImage); copyMachine->Update(); typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType; typename CastFilterType::Pointer caster = CastFilterType::New(); caster->SetInput( copyMachine->GetOutput() ); caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput(); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); this->FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, radiusInMM); // calculate statistics within the binary mask typedef itk::ExtendedLabelStatisticsImageFilter< InputImageType, MaskImageType> LabelStatisticsFilterType; typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( inputImage ); labelStatisticsFilter->SetLabelInput( hotspotMaskITK ); labelStatisticsFilter->SetCoordinateTolerance( 0.001 ); labelStatisticsFilter->SetDirectionTolerance( 0.001 ); labelStatisticsFilter->Update(); Statistics hotspotStatistics; hotspotStatistics.SetHotspotIndex(convolutionImageInformation.MaxIndex); hotspotStatistics.SetMean(convolutionImageInformation.Max); if ( labelStatisticsFilter->HasLabel( 1 ) ) { hotspotStatistics.SetLabel (1); hotspotStatistics.SetN(labelStatisticsFilter->GetCount(1)); hotspotStatistics.SetMin(labelStatisticsFilter->GetMinimum(1)); hotspotStatistics.SetMax(labelStatisticsFilter->GetMaximum(1)); hotspotStatistics.SetMedian(labelStatisticsFilter->GetMedian(1)); hotspotStatistics.SetVariance(labelStatisticsFilter->GetVariance(1)); hotspotStatistics.SetSigma(labelStatisticsFilter->GetSigma(1)); hotspotStatistics.SetRMS(sqrt( hotspotStatistics.GetMean() * hotspotStatistics.GetMean() + hotspotStatistics.GetSigma() * hotspotStatistics.GetSigma() )); MITK_DEBUG << "Statistics for inside hotspot: Mean " << hotspotStatistics.GetMean() << ", SD " << hotspotStatistics.GetSigma() << ", Max " << hotspotStatistics.GetMax() << ", Min " << hotspotStatistics.GetMin(); } else { MITK_ERROR << "Uh oh! Unable to calculate statistics for hotspot region..."; return m_EmptyStatistics; } return hotspotStatistics; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with 1. typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput( image ); castFilter->Update(); castFilter->GetOutput()->FillBuffer( 1 ); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::BaseGeometry *imageGeometry3D = m_Image->GetGeometry( 0 ); // If there is a second poly line in a closed planar figure, treat it as a hole. PlanarFigure::PolyLineType planarFigureHolePolyline; if (m_PlanarFigure->GetPolyLinesSize() == 2) planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1); // Determine x- and y-dimensions depending on principal axis + // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } m_PlanarFigureCoordinate0= i0; m_PlanarFigureCoordinate1= i1; // store the polyline contour as vtkPoints object bool outOfBounds = false; vtkSmartPointer points = vtkSmartPointer::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image + // Fabian: From PlaneGeometry documentation: + // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm). + // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld. planarFigurePlaneGeometry->Map( *it, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } + // Fabian: Why convert to index coordinates? imageGeometry3D->WorldToIndex( point3D, point3D ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } vtkSmartPointer holePoints = nullptr; if (!planarFigureHolePolyline.empty()) { holePoints = vtkSmartPointer::New(); Point3D point3D; PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end(); for (it = planarFigureHolePolyline.begin(); it != end; ++it) { + // Fabian: same as above planarFigurePlaneGeometry->Map(*it, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0); } } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); vtkSmartPointer holeLassoStencil = nullptr; if (holePoints.GetPointer() != nullptr) { holeLassoStencil = vtkSmartPointer::New(); holeLassoStencil->SetShapeToPolygon(); holeLassoStencil->SetPoints(holePoints); } // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); vtkSmartPointer holeStencilFilter = nullptr; if (holeLassoStencil.GetPointer() != nullptr) { holeStencilFilter = vtkSmartPointer::New(); holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); holeStencilFilter->ReverseStencilOn(); holeStencilFilter->SetBackgroundValue(0); holeStencilFilter->Update(); } // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkSmartPointer::New(); vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr ? imageStencilFilter->GetOutputPort() : holeStencilFilter->GetOutputPort()); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage( itkImporter->GetOutput() ); duplicator->Update(); // Store mask m_InternalImageMask2D = duplicator->GetOutput(); } void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp new file mode 100644 index 0000000000..10e5f1e4cb --- /dev/null +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -0,0 +1,5 @@ +#ifndef MITKIMAGESTATISTICSCALCULATOR2 +#define MITKIMAGESTATISTICSCALCULATOR2 + +#endif // MITKIMAGESTATISTICSCALCULATOR2 + diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h new file mode 100644 index 0000000000..10e5f1e4cb --- /dev/null +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -0,0 +1,5 @@ +#ifndef MITKIMAGESTATISTICSCALCULATOR2 +#define MITKIMAGESTATISTICSCALCULATOR2 + +#endif // MITKIMAGESTATISTICSCALCULATOR2 + diff --git a/Modules/ImageStatistics/mitkMaskGenerator.cpp b/Modules/ImageStatistics/mitkMaskGenerator.cpp new file mode 100644 index 0000000000..6f34324163 --- /dev/null +++ b/Modules/ImageStatistics/mitkMaskGenerator.cpp @@ -0,0 +1,33 @@ +#include + +namespace mitk +{ +void MaskGenerator::SetTimeStep(unsigned int timeStep) +{ + if (m_TimeStep != timeStep) + { + m_TimeStep = timeStep; + m_Modified = true; + } +} + +MaskGenerator::MaskGenerator(): + m_TimeStep(0), + m_Modified(false) +{ + +} + +//typename itk::Region<3>::Pointer MaskGenerator::GetImageRegionOfMask(Image::Pointer image) +//{ +// if (m_InternalMask.IsNull() || m_Modified) +// { +// MITK_ERROR << "Update MaskGenerator first!"; +// } + +// mitk::BaseGeometry::Pointer imageGeometry = image->GetGeometry(); +// mitk::BaseGeometry::Pointer maskGeometry = m_InternalMask->GetGeometry(); + + +//} +} diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h new file mode 100644 index 0000000000..7c19bf6b6b --- /dev/null +++ b/Modules/ImageStatistics/mitkMaskGenerator.h @@ -0,0 +1,33 @@ +#ifndef MITKMASKGENERATOR +#define MITKMASKGENERATOR + +#include +#include +#include + +namespace mitk +{ +class MITKIMAGESTATISTICS_EXPORT MaskGenerator +{ +public: + MaskGenerator(); + + //~MaskGenerator(); + + virtual mitk::Image::Pointer GetMask() = 0; + + void SetTimeStep(unsigned int timeStep); + +protected: + unsigned int m_TimeStep; + mitk::Image::Pointer m_InternalMask; + bool m_Modified; + +private: + void CheckMaskSanity(mitk::Image::Pointer); + +}; +} + +#endif // MITKMASKGENERATOR + diff --git a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp new file mode 100644 index 0000000000..e3d467c857 --- /dev/null +++ b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp @@ -0,0 +1 @@ +#include diff --git a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h new file mode 100644 index 0000000000..4a0a0be26f --- /dev/null +++ b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h @@ -0,0 +1,25 @@ +#ifndef MITKMULTILABELMASKGENERATOR +#define MITKMULTILABELMASKGENERATOR + +#include +#include +#include +#include + +namespace mitk +{ + +class MITKIMAGESTATISTICS_EXPORT MultiLabelMaskGenerator: public MaskGenerator +{ +public: + void SetLabelSetImage(mitk::LabelSetImage::Pointer labelSetImage); + +protected: + +private: + +}; + +} + +#endif diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp new file mode 100644 index 0000000000..2df2feb9e8 --- /dev/null +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp @@ -0,0 +1,388 @@ +#include +#include +#include +#include "mitkImageAccessByItk.h" +#include +#include + +#include +#include +#include +#include + + +#include +#include +#include +#include +#include +#include + + + + +namespace mitk +{ + +void PlanarFigureMaskGenerator::SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure) +{ + if ( planarFigure.IsNull() ) + { + throw std::runtime_error( "Error: planar figure empty!" ); + } + if ( !planarFigure->IsClosed() ) + { + throw std::runtime_error( "Masking not possible for non-closed figures" ); + } + + const PlaneGeometry *planarFigurePlaneGeometry = planarFigure->GetPlaneGeometry(); + if ( planarFigurePlaneGeometry == nullptr ) + { + throw std::runtime_error( "Planar-Figure not yet initialized!" ); + } + + const PlaneGeometry *planarFigureGeometry = + dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); + if ( planarFigureGeometry == nullptr ) + { + throw std::runtime_error( "Non-planar planar figures not supported!" ); + } + + if (planarFigure != m_PlanarFigure) + { + m_Modified = true; + m_PlanarFigure = planarFigure; + } + +} + +void PlanarFigureMaskGenerator::SetImage(mitk::Image::Pointer image) +{ + // check dimension + unsigned int dimension = image->GetDimension(); + if (dimension > 3) + { + MITK_ERROR << "unsupported image dimension"; + } + + const BaseGeometry *imageGeometry = image->GetGeometry(); + if ( imageGeometry == nullptr ) + { + throw std::runtime_error( "Image geometry invalid!" ); + } + + if (image != m_InternalInputImage) + { + m_Modified = true; + m_InternalInputImage = image; + } +} + + +template < typename TPixel, unsigned int VImageDimension > +void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure( + const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) +{ + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, 2 > MaskImage2DType; + typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; + + // Generate mask image as new image with same header as input image and + // initialize with 1. + typename CastFilterType::Pointer castFilter = CastFilterType::New(); + castFilter->SetInput( image ); + castFilter->Update(); + castFilter->GetOutput()->FillBuffer( 1 ); +// typename itk::Image::Pointer maskImage = itk::Image::New(); +// maskImage->SetDirection(image->GetDirection()); +// maskImage->SetOrigin(image->GetOrigin()); +// maskImage->SetSpacing(image->GetSpacing()); +// maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); +// maskImage->Allocate(); +// maskImage->FillBuffer(1); + + // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. + // These points are used by the vtkLassoStencilSource to create + // a vtkImageStencil. + const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); + const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); + const mitk::BaseGeometry *imageGeometry3D = m_InternalInputImage->GetGeometry( 0 ); + // If there is a second poly line in a closed planar figure, treat it as a hole. + PlanarFigure::PolyLineType planarFigureHolePolyline; + + if (m_PlanarFigure->GetPolyLinesSize() == 2) + planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1); + + + // Determine x- and y-dimensions depending on principal axis + // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis + int i0, i1; + switch ( axis ) + { + case 0: + i0 = 1; + i1 = 2; + break; + + case 1: + i0 = 0; + i1 = 2; + break; + + case 2: + default: + i0 = 0; + i1 = 1; + break; + } + + // store the polyline contour as vtkPoints object + bool outOfBounds = false; + vtkSmartPointer points = vtkSmartPointer::New(); + typename PlanarFigure::PolyLineType::const_iterator it; + for ( it = planarFigurePolyline.begin(); + it != planarFigurePolyline.end(); + ++it ) + { + Point3D point3D; + + // Convert 2D point back to the local index coordinates of the selected + // image + // Fabian: From PlaneGeometry documentation: + // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm). + // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld. + planarFigurePlaneGeometry->Map( *it, point3D ); + + // Polygons (partially) outside of the image bounds can not be processed + // further due to a bug in vtkPolyDataToImageStencil + if ( !imageGeometry3D->IsInside( point3D ) ) + { + outOfBounds = true; + } + + imageGeometry3D->WorldToIndex( point3D, point3D ); + + points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); + } + + vtkSmartPointer holePoints = nullptr; + + if (!planarFigureHolePolyline.empty()) + { + holePoints = vtkSmartPointer::New(); + + Point3D point3D; + PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end(); + + for (it = planarFigureHolePolyline.begin(); it != end; ++it) + { + // Fabian: same as above + planarFigurePlaneGeometry->Map(*it, point3D); + imageGeometry3D->WorldToIndex(point3D, point3D); + holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0); + } + } + + // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds + // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero + double bounds[6] = {0, 0, 0, 0, 0, 0}; + points->GetBounds( bounds ); + bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; + bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; + bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; + + // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent + if ( m_PlanarFigure->IsClosed() && + ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) + { + mitkThrow() << "Figure has a zero area and cannot be used for masking."; + } + + if ( outOfBounds ) + { + throw std::runtime_error( "Figure at least partially outside of image bounds!" ); + } + + // create a vtkLassoStencilSource and set the points of the Polygon + vtkSmartPointer lassoStencil = vtkSmartPointer::New(); + lassoStencil->SetShapeToPolygon(); + lassoStencil->SetPoints( points ); + + vtkSmartPointer holeLassoStencil = nullptr; + + if (holePoints.GetPointer() != nullptr) + { + holeLassoStencil = vtkSmartPointer::New(); + holeLassoStencil->SetShapeToPolygon(); + holeLassoStencil->SetPoints(holePoints); + } + + // Export from ITK to VTK (to use a VTK filter) + typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; + typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; + + typename ImageExportType::Pointer itkExporter = ImageExportType::New(); + //itkExporter->SetInput( maskImage ); + itkExporter->SetInput( castFilter->GetOutput() ); + + vtkSmartPointer vtkImporter = vtkSmartPointer::New(); + this->ConnectPipelines( itkExporter, vtkImporter ); + + // Apply the generated image stencil to the input image + vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); + imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); + imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); + imageStencilFilter->ReverseStencilOff(); + imageStencilFilter->SetBackgroundValue( 0 ); + imageStencilFilter->Update(); + + vtkSmartPointer holeStencilFilter = nullptr; + + if (holeLassoStencil.GetPointer() != nullptr) + { + holeStencilFilter = vtkSmartPointer::New(); + holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); + holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); + holeStencilFilter->ReverseStencilOn(); + holeStencilFilter->SetBackgroundValue(0); + holeStencilFilter->Update(); + } + + // Export from VTK back to ITK + vtkSmartPointer vtkExporter = vtkSmartPointer::New(); + vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr + ? imageStencilFilter->GetOutputPort() + : holeStencilFilter->GetOutputPort()); + vtkExporter->Update(); + + typename ImageImportType::Pointer itkImporter = ImageImportType::New(); + this->ConnectPipelines( vtkExporter, itkImporter ); + itkImporter->Update(); + + typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType; + DuplicatorType::Pointer duplicator = DuplicatorType::New(); + duplicator->SetInputImage( itkImporter->GetOutput() ); + duplicator->Update(); + + // Store mask + m_InternalITKImageMask2D = duplicator->GetOutput(); +} + +bool PlanarFigureMaskGenerator::GetPrincipalAxis( + const BaseGeometry *geometry, Vector3D vector, + unsigned int &axis ) +{ + vector.Normalize(); + for ( unsigned int i = 0; i < 3; ++i ) + { + Vector3D axisVector = geometry->GetAxisVector( i ); + axisVector.Normalize(); + + if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) + { + axis = i; + return true; + } + } + + return false; +} + +void PlanarFigureMaskGenerator::CalculateMask() +{ + m_InternalITKImageMask2D = nullptr; + const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); + const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); + const BaseGeometry *imageGeometry = m_InternalInputImage->GetGeometry(); + + // Find principal direction of PlanarFigure in input image + unsigned int axis; + if ( !this->GetPrincipalAxis( imageGeometry, + planarFigureGeometry->GetNormal(), axis ) ) + { + throw std::runtime_error( "Non-aligned planar figures not supported!" ); + } + m_PlanarFigureAxis = axis; + + // Find slice number corresponding to PlanarFigure in input image + typename itk::Image< unsigned short, 3 >::IndexType index; + imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); + + unsigned int slice = index[axis]; + + // extract image slice which corresponds to the planarFigure and sotre it in m_InternalImageSlice + mitk::Image::Pointer inputImageSlice = extract2DImageSlice(axis, slice); + + // Compute mask from PlanarFigure + AccessFixedDimensionByItk_1(inputImageSlice, + InternalCalculateMaskFromPlanarFigure, + 2, axis) + + //convert itk mask to mitk::Image::Pointer and return it + mitk::Image::Pointer planarFigureMaskImage; + planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D); + planarFigureMaskImage->SetGeometry(inputImageSlice->GetGeometry()); + + Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New(); + sliceTo3DImageConverter->SetInput(planarFigureMaskImage); + sliceTo3DImageConverter->Update(); + + m_InternalMask = sliceTo3DImageConverter->GetOutput(); +} + +mitk::Image::Pointer PlanarFigureMaskGenerator::GetMask() +{ + if (m_Modified) + { + if (m_InternalInputImage.IsNull()) + { + MITK_ERROR << "Image is not set."; + } + + if (m_PlanarFigure.IsNull()) + { + MITK_ERROR << "PlanarFigure is not set."; + } + + if (m_TimeStep != 0) + { + MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet)."; + } + + this->CalculateMask(); + } + + m_Modified = false; + + return m_InternalMask; +} + +mitk::Image::Pointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice) +{ + // Extract slice with given position and direction from image + unsigned int dimension = m_InternalInputImage->GetDimension(); + mitk::Image::Pointer imageSlice = mitk::Image::New(); + + if (dimension == 3) + { + ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); + imageExtractor->SetInput( m_InternalInputImage ); + imageExtractor->SetSliceDimension( axis ); + imageExtractor->SetSliceIndex( slice ); + imageExtractor->Update(); + imageSlice = imageExtractor->GetOutput(); + } + else if(dimension == 2) + { + imageSlice = m_InternalInputImage; + } + else + { + MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported."; + } + + return imageSlice; +} + +} + diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h new file mode 100644 index 0000000000..38c42f8cac --- /dev/null +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -0,0 +1,95 @@ +#ifndef MITKPLANARFIGUREMASKGENERATOR +#define MITKPLANARFIGUREMASKGENERATOR + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mitk +{ +class MITKIMAGESTATISTICS_EXPORT PlanarFigureMaskGenerator: public MaskGenerator + { + public: + + mitk::Image::Pointer GetMask(); + + void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure); + + /*We need an image because we need its geometry in order to give the points of the planar figure the correct world coordinates*/ + void SetImage(const mitk::Image::Pointer inputImage); + + protected: + + + private: + void CalculateMask(); + + template < typename TPixel, unsigned int VImageDimension > + void InternalCalculateMaskFromPlanarFigure( + const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); + + mitk::Image::Pointer extract2DImageSlice(unsigned int axis, unsigned int slice); + + bool GetPrincipalAxis(const BaseGeometry *geometry, Vector3D vector, + unsigned int &axis ); + + /** Connection from ITK to VTK */ + template + void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) + { + importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); + + importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); + importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); + importer->SetSpacingCallback(exporter->GetSpacingCallback()); + importer->SetOriginCallback(exporter->GetOriginCallback()); + importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); + + importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); + + importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); + importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); + importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); + importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); + importer->SetCallbackUserData(exporter->GetCallbackUserData()); + } + + /** Connection from VTK to ITK */ + template + void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) + { + importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); + + importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); + importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); + importer->SetSpacingCallback(exporter->GetSpacingCallback()); + importer->SetOriginCallback(exporter->GetOriginCallback()); + importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); + + importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); + + importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); + importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); + importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); + importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); + importer->SetCallbackUserData(exporter->GetCallbackUserData()); + } + + + mitk::PlanarFigure::Pointer m_PlanarFigure; + typename itk::Image::Pointer m_InternalITKImageMask2D; + mitk::Image::Pointer m_InternalInputImage; + + unsigned int m_PlanarFigureAxis; + }; +} + +#endif // MITKPLANARFIGUREMASKGENERATOR + diff --git a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp index 46e692d7c3..cc725b7871 100644 --- a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp @@ -1,519 +1,520 @@ /*=================================================================== 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::SetInput( unsigned int idx, const mitk::Surface* surface ) { this->SetNthInput( idx, const_cast( surface ) ); this->Modified(); } void mitk::ReduceContourSetFilter::SetInput( const mitk::Surface* surface ) { this->SetInput( 0, const_cast( surface ) ); } 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 = vtkSmartPointer::New(); newPolygons = vtkSmartPointer::New(); newPoints = vtkSmartPointer::New(); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType* cell (nullptr); 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); if (numberOfOutputs == 0) { mitk::Surface::Pointer tmp_output = mitk::Surface::New(); tmp_output->SetVtkPolyData(vtkPolyData::New()); this->SetNthOutput(0, tmp_output.GetPointer()); } //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 (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 <= 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(nullptr); /* 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 (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 (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); + // BUG XXXXX Fix 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; }