diff --git a/Modules/Classification/CLActiveLearning/include/mitkActiveLearningInteractor.h b/Modules/Classification/CLActiveLearning/include/mitkActiveLearningInteractor.h index 2c84ced023..5f0a33534c 100644 --- a/Modules/Classification/CLActiveLearning/include/mitkActiveLearningInteractor.h +++ b/Modules/Classification/CLActiveLearning/include/mitkActiveLearningInteractor.h @@ -1,64 +1,73 @@ /*=================================================================== 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 mitkActiveLearningInteractor_h #define mitkActiveLearningInteractor_h #include #include #include #include #include namespace mitk { class MITKCLACTIVELEARNING_EXPORT ActiveLearningInteractor : public DataInteractor { public: typedef unsigned short AnnotationPixelType; mitkClassMacro(ActiveLearningInteractor, DataInteractor) itkFactorylessNewMacro(Self) void SetPaintingPixelValue(AnnotationPixelType value){m_PaintingPixelValue = static_cast(value);} + void SetSize(unsigned int size){m_Size = size;} + bool IsUsed(){return m_Used;} private: ActiveLearningInteractor(); ~ActiveLearningInteractor(); void ConnectActionsAndFunctions() override; void DataNodeChanged() override; void Paint(mitk::StateMachineAction* action, mitk::InteractionEvent* event); void PaintInterpolate(mitk::StateMachineAction* action, mitk::InteractionEvent* event); - int m_Size; - itk::Index<3> m_LastPixelIndex; + void SetWorkingSlice(mitk::StateMachineAction* action, mitk::InteractionEvent* event); + + void WriteBackWorkingSlice(mitk::StateMachineAction* action, mitk::InteractionEvent* event); + + mitk::Image::Pointer m_WorkingSlice; + mitk::PlaneGeometry::Pointer m_WorkingPlane; + + unsigned int m_Size; + mitk::Point3D m_LastPosition; AnnotationPixelType m_PaintingPixelValue; bool m_Used; }; } #endif diff --git a/Modules/Classification/CLActiveLearning/resource/Interactions/Paint.xml b/Modules/Classification/CLActiveLearning/resource/Interactions/Paint.xml index 871d3ff218..ef1d27da29 100644 --- a/Modules/Classification/CLActiveLearning/resource/Interactions/Paint.xml +++ b/Modules/Classification/CLActiveLearning/resource/Interactions/Paint.xml @@ -1,25 +1,29 @@ + + - + + + diff --git a/Modules/Classification/CLActiveLearning/src/mitkActiveLearningInteractor.cpp b/Modules/Classification/CLActiveLearning/src/mitkActiveLearningInteractor.cpp index db205205f2..8e4fe02c60 100644 --- a/Modules/Classification/CLActiveLearning/src/mitkActiveLearningInteractor.cpp +++ b/Modules/Classification/CLActiveLearning/src/mitkActiveLearningInteractor.cpp @@ -1,637 +1,551 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include +#include #include #include #include #include #include #include - -#define ROUND(a) ((a)>0 ? (int)((a)+0.5) : -(int)(0.5-(a))) +#include +#include +#include +#include +#include +#include // Helper function to get an image from a data node. static mitk::Image::Pointer GetImage(mitk::DataNode::Pointer dataNode) { if (dataNode.IsNull()) mitkThrow(); mitk::Image::Pointer image = dynamic_cast(dataNode->GetData()); if (image.IsNull()) mitkThrow(); return image; } // Helper function to get a geometry of an image for a specific time step. static mitk::BaseGeometry::Pointer GetGeometry(mitk::Image* image, unsigned int timeStep) { if (image == nullptr) mitkThrow(); mitk::TimeGeometry::Pointer timeGeometry = image->GetTimeGeometry(); if (timeGeometry.IsNull()) mitkThrow(); auto geometry = timeGeometry->GetGeometryForTimeStep(timeStep); if (geometry.IsNull()) mitkThrow(); return geometry; } -static mitk::Point2D upperLeft(mitk::Point2D p) -{ - p[0] -= 0.5; - p[1] += 0.5; - return p; -} - static double EuclideanDistance(const mitk::Point3D p1, const mitk::Point3D p2) { return std::sqrt( (p1[0] - p2[0])*(p1[0] - p2[0]) + (p1[1] - p2[1])*(p1[1] - p2[1]) + (p1[2] - p2[2])*(p1[2] - p2[2]) ); } +static double EuclideanDistance2D(const mitk::Point2D p1, const mitk::Point2D p2) +{ + return std::sqrt( (p1[0] - p2[0])*(p1[0] - p2[0]) + + (p1[1] - p2[1])*(p1[1] - p2[1]) ); +} + static double ScalarProduct(const mitk::Point3D p1, const mitk::Point3D p2) { return p1[0]*p2[0] + p1[1]*p2[1] + p1[2]*p2[2]; } -static std::vector> InterpolateIndices(itk::Index<3> startIndex, itk::Index<3> endIndex, mitk::BaseGeometry* geometry, unsigned int size) +static double ScalarProduct2D(const mitk::Point2D p1, const mitk::Point2D p2) +{ + return p1[0]*p2[0] + p1[1]*p2[1]; +} + +static std::vector> InterpolateIndices2D(itk::Index<2> startIndex, itk::Index<2> endIndex, const mitk::PlaneGeometry* geometry, unsigned int size) { if (geometry == nullptr) mitkThrow(); - MITK_INFO << "Interpolating between " << startIndex << " and " << endIndex; + std::vector> resultIndices; + mitk::Point2D startPoint; + mitk::Point2D endPoint; + startPoint[0] = startIndex[0]; + startPoint[1] = startIndex[1]; + endPoint[0] = endIndex[0]; + endPoint[1] = endIndex[1]; + geometry->IndexToWorld(startPoint, startPoint); + geometry->IndexToWorld(endPoint, endPoint); + + // Distance between end points + double dist = EuclideanDistance2D(startPoint, endPoint); + + // Define region between startIndex and endIndex, padded by (size - 1) + int regionBounds[4]; + for (int i=0; i<2; ++i) + { + regionBounds[2*i] = std::min(startIndex[i], endIndex[i]) - (size - 1); + regionBounds[2*i+1] = std::max(startIndex[i], endIndex[i]) + (size - 1); + } + + // We only have a threshold given in pixels (size), but image can be spaced in different units. + // To get the corresponding distances, transform unit vectors and get their lengths. + // The minimum spacing will be what we compare to. + double spacingInIndexSystem[2]; + for (int i=0; i<2; ++i) + { + mitk::Point2D origin; + origin.Fill(0); + mitk::Point2D index; + index.Fill(0); + index[i] = 1; + mitk::Point2D p_origin; + mitk::Point2D p_index; + geometry->IndexToWorld(origin, p_origin); + geometry->IndexToWorld(index, p_index); + spacingInIndexSystem[i] = EuclideanDistance2D(p_origin, p_index); + } + double minSpacing = (spacingInIndexSystem[0] > spacingInIndexSystem[1]) ? spacingInIndexSystem[0] : spacingInIndexSystem[1]; + + double t = 0; + double d = 0; + for (int x = regionBounds[0]; x<=regionBounds[1]; ++x) + { + for (int y = regionBounds[2]; y<=regionBounds[3]; ++y) + { + mitk::Point2D p; + mitk::Point2D index; + index[0] = x; + index[1] = y; + geometry->IndexToWorld(index, p); + + // if there is not distance between start and end, just get distance to start + if (dist < mitk::eps) + { + d = EuclideanDistance2D(startPoint, p); + } + else + { + t = -1./(dist*dist) * ScalarProduct2D(startPoint - p, endPoint - startPoint); + if (t > 0. && t < 1.) + { + d = std::sqrt( ScalarProduct2D(startPoint - p, startPoint - p) + + 2. * t * ScalarProduct2D(endPoint - startPoint, startPoint - p) + + t * t * dist * dist ); + } + else if (t <= 0.) + { + d = EuclideanDistance2D(startPoint, p); + } + else + { + d = EuclideanDistance2D(endPoint, p); + } + } + + if (d <= minSpacing) + { + resultIndices.push_back(itk::Index<2>({index[0], index[1]})); + } + } + } + + return resultIndices; +} + +static std::vector> InterpolateIndices(itk::Index<3> startIndex, itk::Index<3> endIndex, const mitk::BaseGeometry* geometry, unsigned int size) +{ + if (geometry == nullptr) + mitkThrow(); std::vector> resultIndices; mitk::Point3D startPoint; mitk::Point3D endPoint; geometry->IndexToWorld(startIndex, startPoint); geometry->IndexToWorld(endIndex, endPoint); // Distance between end points double dist = EuclideanDistance(startPoint, endPoint); - MITK_INFO << "Distance between points is " << dist; - // Define region between startIndex and endIndex, padded by (size - 1) int regionBounds[6]; for (int i=0; i<3; ++i) { regionBounds[2*i] = std::min(startIndex[i], endIndex[i]) - (size - 1); regionBounds[2*i+1] = std::max(startIndex[i], endIndex[i]) + (size - 1); } - MITK_INFO << "Region size is " << regionBounds[0] << ", " << regionBounds[1] << ", " << regionBounds[2] << ", " << regionBounds[3] << ", " << regionBounds[4] << ", " << regionBounds[5]; - // We only have a threshold given in pixels (size), but image can be spaced in different units. // To get the corresponding distances, transform unit vectors and get their lengths. // The minimum spacing will be what we compare to. + // Could potentially use the spacing to normalize. double spacingInIndexSystem[3]; double minSpacing = -1.; for (int i=0; i<3; ++i) { itk::Index<3> origin; origin.Fill(0); itk::Index<3> index; index.Fill(0); index[i] = 1; mitk::Point3D p_origin; mitk::Point3D p_index; geometry->IndexToWorld(origin, p_origin); geometry->IndexToWorld(index, p_index); double spacing = EuclideanDistance(p_origin, p_index); if ( (minSpacing > 0. && spacing < minSpacing) || minSpacing < 0. ) { minSpacing = spacing; } spacingInIndexSystem[i] = spacing; } - MITK_INFO << "Spacings are " << spacingInIndexSystem[0] << ", " << spacingInIndexSystem[1] << ", " << spacingInIndexSystem[2]; - // Iterate over all indices in the given region and get distance to the line between startPoint and endPoint. // If distance is smaller than size, add to resultIndices. // // Let (x1,y1,z1) = startPoint, (x2,y2,z2) = endPoint, (x0,y0,z0) = p a point. // // Line is defined by: // [x1 + (x2-x1) * t] // v = [y1 + (y2-y1) * t] // [z1 + (z2-z1) * t] // // Then (with * dot product): // t(p) = - (startPoint - p) * (endPoint - startPoint) / |endPoint - startPoint|^2 // // And (with x cross product): // d(p) = |(p - startPoint) x (p - endPoint)| / |endPoint - startPoint| double t = 0; double d = 0; for (int x = regionBounds[0]; x<=regionBounds[1]; ++x) { for (int y = regionBounds[2]; y<=regionBounds[3]; ++y) { for (int z = regionBounds[4]; z<=regionBounds[5]; ++z) { mitk::Point3D p; itk::Index<3> index = {x,y,z}; geometry->IndexToWorld(index, p); // if there is not distance between start and end, just get distance to start if (dist < mitk::eps) { d = EuclideanDistance(startPoint, p); } else { t = -1./(dist*dist) * ScalarProduct(startPoint - p, endPoint - startPoint); if (t > 0. && t < 1.) { d = std::sqrt( ScalarProduct(startPoint - p, startPoint - p) + 2. * t * ScalarProduct(endPoint - startPoint, startPoint - p) + t * t * dist * dist ); } else if (t <= 0.) { d = EuclideanDistance(startPoint, p); } else { d = EuclideanDistance(endPoint, p); } } if (d <= minSpacing) { resultIndices.push_back(index); } } } } - MITK_INFO << "Found " << resultIndices.size() << " indices"; - return resultIndices; +// ======= OLD INTERPOLATION ======== + // std::vector> resultIndices; // mitk::Point3D startPoint; // mitk::Point3D endPoint; // geometry->IndexToWorld(startIndex, startPoint); // geometry->IndexToWorld(endIndex, endPoint); // itk::Index<3> indexDelta; // int indexDeltaInc[3]; // for (int i=0; i<3; i++) // { // indexDelta[i] = endIndex[i] - startIndex[i]; // indexDeltaInc[i] = (indexDelta[i] > 0) ? 1 : (indexDelta[i] < 0) ? -1 : 0; // } // int argm[3] = {0, 1, 2}; // if (abs(indexDelta[1]) > abs(indexDelta[0])) // { // argm[0] = 1; // argm[1] = 0; // } // if (abs(indexDelta[2]) > abs(indexDelta[argm[1]])) // { // argm[2] = argm[1]; // argm[1] = 2; // } // if (abs(indexDelta[2]) > abs(indexDelta[argm[0]])) // { // argm[1] = argm[0]; // argm[0] = 2; // } // double slopes[2]; // slopes[0] = (endPoint[argm[1]] - startPoint[argm[1]]) / (endPoint[argm[0]] - startPoint[argm[0]]); // slopes[1] = (endPoint[argm[2]] - startPoint[argm[2]]) / sqrt((endPoint[argm[1]] - startPoint[argm[1]]) * (endPoint[argm[1]] - startPoint[argm[1]]) + (endPoint[argm[0]] - startPoint[argm[0]]) * (endPoint[argm[0]] - startPoint[argm[0]])); // itk::Index<3> currentIndex = startIndex; // mitk::Point3D currentPoint = startPoint; // while (currentIndex != endIndex) // { // currentIndex[argm[0]] += indexDeltaInc[argm[0]]; // geometry->IndexToWorld(currentIndex, currentPoint); // currentPoint[argm[1]] = startPoint[argm[1]] + slopes[0] * (currentPoint[argm[0]] - startPoint[argm[0]]); // currentPoint[argm[2]] = startPoint[argm[2]] + slopes[1] * sqrt((currentPoint[argm[1]] - startPoint[argm[1]]) * (currentPoint[argm[1]] - startPoint[argm[1]]) + (currentPoint[argm[0]] - startPoint[argm[0]]) * (currentPoint[argm[0]] - startPoint[argm[0]])); // geometry->WorldToIndex(currentPoint, currentIndex); // resultIndices.push_back(currentIndex); // } // return resultIndices; } mitk::ActiveLearningInteractor::ActiveLearningInteractor() : - m_Size(1), + m_WorkingSlice(nullptr), + m_WorkingPlane(nullptr), + m_Size(5), m_PaintingPixelValue(0) { } mitk::ActiveLearningInteractor::~ActiveLearningInteractor() { } void mitk::ActiveLearningInteractor::ConnectActionsAndFunctions() { CONNECT_FUNCTION("paint", Paint) - CONNECT_FUNCTION("paint_interpolate", PaintInterpolate) + CONNECT_FUNCTION("paint_interpolate", PaintInterpolate) + CONNECT_FUNCTION("set_workingslice", SetWorkingSlice) + CONNECT_FUNCTION("writeback_workingslice", WriteBackWorkingSlice) } void mitk::ActiveLearningInteractor::DataNodeChanged() { this->ResetToStartState(); } -void mitk::ActiveLearningInteractor::Paint(mitk::StateMachineAction* /*action*/, mitk::InteractionEvent* event) +void mitk::ActiveLearningInteractor::SetWorkingSlice(mitk::StateMachineAction* /*action*/, + mitk::InteractionEvent *event) +{ + try + { + auto renderer = event->GetSender(); + auto image = GetImage(this->GetDataNode()); + auto timeStep = renderer->GetTimeStep(); + auto geometry = GetGeometry(image, timeStep); + + // Set current plane geometry + m_WorkingPlane = renderer->GetCurrentWorldPlaneGeometry()->Clone(); + + // Extract corresponding slice + vtkSmartPointer reslice = vtkSmartPointer::New(); + reslice->SetOverwriteMode(false); + reslice->Modified(); + auto extract = mitk::ExtractSliceFilter::New(reslice); + extract->SetInput(image); + extract->SetTimeStep(timeStep); + extract->SetWorldGeometry(m_WorkingPlane); + extract->SetResliceTransformByGeometry(geometry); + extract->Update(); + m_WorkingSlice = extract->GetOutput()->Clone(); + MITK_INFO << ""; + } + catch (itk::ExceptionObject& e) + { + mitkThrow() << "Could not set working slice, because:"; + mitkThrow() << e.GetDescription(); + } + +} + +void mitk::ActiveLearningInteractor::WriteBackWorkingSlice(mitk::StateMachineAction* /*action*/, + mitk::InteractionEvent *event) +{ + try + { + auto renderer = event->GetSender(); + auto image = GetImage(this->GetDataNode()); + auto timeStep = renderer->GetTimeStep(); + auto geometry = GetGeometry(image, timeStep); + + // Write back + vtkSmartPointer reslice = vtkSmartPointer::New(); + reslice->SetInputSlice(m_WorkingSlice->GetVtkImageData()); + reslice->SetOverwriteMode(true); + reslice->Modified(); + mitk::ExtractSliceFilter::Pointer insert = mitk::ExtractSliceFilter::New(reslice); + insert->SetInput(image); + insert->SetTimeStep(timeStep); + insert->SetWorldGeometry(m_WorkingPlane); + insert->SetResliceTransformByGeometry(geometry); + insert->Modified(); + insert->Update(); + + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } + catch (itk::ExceptionObject& e) + { + mitkThrow() << "Could not set working slice, because:"; + mitkThrow() << e.GetDescription(); + } + +} + +void mitk::ActiveLearningInteractor::Paint(mitk::StateMachineAction* /*action*/, + mitk::InteractionEvent* event) { - if (m_PaintingPixelValue == -1) return; + if (m_PaintingPixelValue == -1) mitkThrow() << "Cannot paint negative values"; try { auto renderer = event->GetSender(); auto image = GetImage(this->GetDataNode()); auto timeStep = renderer->GetTimeStep(); auto geometry = GetGeometry(image, timeStep); auto positionEvent = dynamic_cast(event); auto position = positionEvent->GetPositionInWorld(); if (!geometry->IsInside(position)) return; // Okay, we're safe. Convert the mouse position to the index of the pixel // we're pointing at. itk::Index<3> index; + itk::Index<3> oldIndex; geometry->WorldToIndex(position, index); + geometry->WorldToIndex(m_LastPosition, oldIndex); // We don't need to paint over and over again while moving the mouse // pointer inside the same pixel. That's especially relevant when operating // on zoomed images. - if (index != m_LastPixelIndex) + if (index != oldIndex) { - mitk::ImagePixelWriteAccessor writeAccessor(image, image->GetVolumeData(timeStep)); - writeAccessor.SetPixelByIndexSafe(index, m_PaintingPixelValue); + // Convert index to slice geometry + m_WorkingSlice->GetGeometry()->WorldToIndex(position, index); + itk::Index<2> indexInPlane2D; + indexInPlane2D[0] = index[0]; + indexInPlane2D[1] = index[1]; + + // Get indices + auto indices = InterpolateIndices2D(indexInPlane2D, indexInPlane2D, m_WorkingPlane, m_Size); + + // Fill indices + mitk::ImagePixelWriteAccessor writeAccessor(m_WorkingSlice, m_WorkingSlice->GetVolumeData(timeStep)); + for (auto i : indices) + { + writeAccessor.SetPixelByIndexSafe(i, m_PaintingPixelValue); + } + + // Update image->Modified(); this->GetDataNode()->Modified(); - mitk::RenderingManager::GetInstance()->RequestUpdateAll(); - m_LastPixelIndex = index; + m_LastPosition = position; m_Used = true; } } - catch (...) + catch (itk::ExceptionObject& e) { - return; + mitkThrow() << "Could not paint, because:"; + mitkThrow() << e.GetDescription(); } } void mitk::ActiveLearningInteractor::PaintInterpolate(mitk::StateMachineAction* /*action*/, mitk::InteractionEvent* event) { - if (m_PaintingPixelValue == -1) return; + if (m_PaintingPixelValue == -1) mitkThrow() << "Cannot paint negative values"; try { auto renderer = event->GetSender(); auto image = GetImage(this->GetDataNode()); auto timeStep = renderer->GetTimeStep(); auto geometry = GetGeometry(image, timeStep); - auto planeGeometry = renderer->GetCurrentWorldPlaneGeometry(); auto positionEvent = dynamic_cast(event); auto position = positionEvent->GetPositionInWorld(); -// mitk::Point3D index; -// geometry->WorldToIndex(position, index); if (!geometry->IsInside(position)) return; -// // Extract current slice -// vtkSmartPointer reslice = vtkSmartPointer::New(); -// //set to false to extract a slice -// reslice->SetOverwriteMode(false); -// reslice->Modified(); -// mitk::ExtractSliceFilter::Pointer extractor = mitk::ExtractSliceFilter::New(reslice); -// extractor->SetInput(image); -// extractor->SetTimeStep(timeStep); -// extractor->SetWorldGeometry(planeGeometry); -// extractor->SetVtkOutputRequest(false); -// extractor->SetResliceTransformByGeometry(image->GetTimeGeometry()->GetGeometryForTimeStep(timeStep)); -// extractor->Modified(); -// extractor->Update(); -// auto workingSlice = extractor->GetOutput(); - -// // Draw a contour in Square according to selected brush size -// int radius = m_Size / 2; -// double fradius = static_cast(m_Size) / 2.0; - -// ContourModel::Pointer contourInImageIndexCoordinates = ContourModel::New(); - -// // estimate center point of the brush ( relative to the pixel the mouse points on ) -// // -- left upper corner for even sizes, -// // -- midpoint for uneven sizes -// mitk::Point2D centerCorrection; -// centerCorrection.Fill(0); - -// // even --> correction of [+0.5, +0.5] -// bool evenSize = ((m_Size % 2) == 0); -// if( evenSize ) -// { -// centerCorrection[0] += 0.5; -// centerCorrection[1] += 0.5; -// } - -// // we will compute the control points for the upper left quarter part of a circle contour -// std::vector< mitk::Point2D > quarterCycleUpperRight; -// std::vector< mitk::Point2D > quarterCycleLowerRight; -// std::vector< mitk::Point2D > quarterCycleLowerLeft; -// std::vector< mitk::Point2D > quarterCycleUpperLeft; - - -// mitk::Point2D curPoint; -// bool curPointIsInside = true; -// curPoint[0] = 0; -// curPoint[1] = radius; -// quarterCycleUpperRight.push_back(upperLeft(curPoint)); - -// // to estimate if a pixel is inside the circle, we need to compare against the 'outer radius' -// // i.e. the distance from the midpoint [0,0] to the border of the pixel [0,radius] -// //const float outer_radius = static_cast(radius) + 0.5; - -// while (curPoint[1] > 0) -// { -// // Move right until pixel is outside circle -// float curPointX_squared = 0.0f; -// float curPointY_squared = (curPoint[1] - centerCorrection[1] ) * (curPoint[1] - centerCorrection[1] ); -// while( curPointIsInside ) -// { -// // increment posX and chec -// curPoint[0]++; -// curPointX_squared = (curPoint[0] - centerCorrection[0] ) * (curPoint[0] - centerCorrection[0] ); -// const float len = sqrt( curPointX_squared + curPointY_squared); -// if ( len > fradius ) -// { -// // found first Pixel in this horizontal line, that is outside the circle -// curPointIsInside = false; -// } -// } -// quarterCycleUpperRight.push_back( upperLeft(curPoint) ); - -// // Move down until pixel is inside circle -// while( !curPointIsInside ) -// { -// // increment posX and chec -// curPoint[1]--; -// curPointY_squared = (curPoint[1] - centerCorrection[1] ) * (curPoint[1] - centerCorrection[1] ); -// const float len = sqrt( curPointX_squared + curPointY_squared); -// if ( len <= fradius ) -// { -// // found first Pixel in this horizontal line, that is outside the circle -// curPointIsInside = true; -// quarterCycleUpperRight.push_back( upperLeft(curPoint) ); -// } - -// // Quarter cycle is full, when curPoint y position is 0 -// if (curPoint[1] <= 0) -// break; -// } - -// } - -// // QuarterCycle is full! Now copy quarter cycle to other quarters. - -// if( !evenSize ) -// { -// std::vector< mitk::Point2D >::const_iterator it = quarterCycleUpperRight.begin(); -// while( it != quarterCycleUpperRight.end() ) -// { -// mitk::Point2D p; -// p = *it; - -// // the contour points in the lower right corner have same position but with negative y values -// p[1] *= -1; -// quarterCycleLowerRight.push_back(p); - -// // the contour points in the lower left corner have same position -// // but with both x,y negative -// p[0] *= -1; -// quarterCycleLowerLeft.push_back(p); - -// // the contour points in the upper left corner have same position -// // but with x negative -// p[1] *= -1; -// quarterCycleUpperLeft.push_back(p); - -// it++; -// } -// } -// else -// { -// std::vector< mitk::Point2D >::const_iterator it = quarterCycleUpperRight.begin(); -// while( it != quarterCycleUpperRight.end() ) -// { -// mitk::Point2D p,q; -// p = *it; - -// q = p; -// // the contour points in the lower right corner have same position but with negative y values -// q[1] *= -1; -// // correct for moved offset if size even = the midpoint is not the midpoint of the current pixel -// // but its upper rigt corner -// q[1] += 1; -// quarterCycleLowerRight.push_back(q); - -// q = p; -// // the contour points in the lower left corner have same position -// // but with both x,y negative -// q[1] = -1.0f * q[1] + 1; -// q[0] = -1.0f * q[0] + 1; -// quarterCycleLowerLeft.push_back(q); - -// // the contour points in the upper left corner have same position -// // but with x negative -// q = p; -// q[0] *= -1; -// q[0] += 1; -// quarterCycleUpperLeft.push_back(q); - -// it++; -// } -// } - -// // fill contour with points in right ordering, starting with the upperRight block -// mitk::Point3D tempPoint; -// for (unsigned int i=0; iAddVertex( tempPoint ); -// } -// // the lower right has to be parsed in reverse order -// for (int i=quarterCycleLowerRight.size()-1; i>=0; i--) -// { -// tempPoint[0] = quarterCycleLowerRight[i][0]; -// tempPoint[1] = quarterCycleLowerRight[i][1]; -// tempPoint[2] = 0; -// contourInImageIndexCoordinates->AddVertex( tempPoint ); -// } -// for (unsigned int i=0; iAddVertex( tempPoint ); -// } -// // the upper left also has to be parsed in reverse order -// for (int i=quarterCycleUpperLeft.size()-1; i>=0; i--) -// { -// tempPoint[0] = quarterCycleUpperLeft[i][0]; -// tempPoint[1] = quarterCycleUpperLeft[i][1]; -// tempPoint[2] = 0; -// contourInImageIndexCoordinates->AddVertex( tempPoint ); -// } - -// // round to nearest voxel center (abort if this hasn't changed) -// if ( m_Size % 2 == 0 ) // even -// { -// index[0] = ROUND(index[0]); -// index[1] = ROUND(index[1]); -// } -// else // odd -// { -// index[0] = ROUND(index[0]); -// index[1] = ROUND(index[1]); -// } - -// ContourModel::Pointer contour = ContourModel::New(); -// contour->Expand(timeStep + 1); -// contour->SetClosed(true, timeStep); - -// ContourModel::VertexIterator it = contourInImageIndexCoordinates->Begin(); -// ContourModel::VertexIterator end = contourInImageIndexCoordinates->End(); - -// while(it != end) -// { -// mitk::Point3D point = (*it)->Coordinates; -// point[0] += index[0]; -// point[1] += index[1]; - -// contour->AddVertex(point, timeStep); -// it++; -// } - -// double dist = index.EuclideanDistanceTo(m_LastPixelIndex); - -// // draw it the old way -// FeedbackContourTool::FillContourInSlice(contour, timeStep, workingSlice, m_PaintingPixelValue); - -// // if points are >= radius away draw rectangle to fill empty holes -// // in between the 2 points -// if (dist > radius) -// { -// const mitk::Point3D& currentPos = index; -// mitk::Point3D direction; -// mitk::Point3D vertex; -// mitk::Point3D normal; - -// direction[0] = index[0] - m_LastPixelIndex[0]; -// direction[1] = index[1] - m_LastPixelIndex[1]; -// direction[2] = index[2] - m_LastPixelIndex[2]; - -// direction[0] = direction.GetVnlVector().normalize()[0]; -// direction[1] = direction.GetVnlVector().normalize()[1]; -// direction[2] = direction.GetVnlVector().normalize()[2]; - -// // 90 degrees rotation of direction -// normal[0] = -1.0 * direction[1]; -// normal[1] = direction[0]; - -// contour->Clear(); - -// // upper left corner -// vertex[0] = m_LastPixelIndex[0] + (normal[0] * radius); -// vertex[1] = m_LastPixelIndex[1] + (normal[1] * radius); - -// contour->AddVertex(vertex); - -// // upper right corner -// vertex[0] = currentPos[0] + (normal[0] * radius); -// vertex[1] = currentPos[1] + (normal[1] * radius); - -// contour->AddVertex(vertex); - -// // lower right corner -// vertex[0] = currentPos[0] - (normal[0] * radius); -// vertex[1] = currentPos[1] - (normal[1] * radius); - -// contour->AddVertex(vertex); - -// // lower left corner -// vertex[0] = m_LastPixelIndex[0] - (normal[0] * radius); -// vertex[1] = m_LastPixelIndex[1] - (normal[1] * radius); - -// contour->AddVertex(vertex); - -// FeedbackContourTool::FillContourInSlice(contour, timeStep, workingSlice, m_PaintingPixelValue); -// } - -// m_LastPixelIndex = index; -// m_Used = true; - -// RenderingManager::GetInstance()->RequestUpdateAll(); - - // Okay, we're safe. Convert the mouse position to the index of the pixel - // we're pointing at. - itk::Index<3> index; - geometry->WorldToIndex<3>(position, index); - - // We don't need to paint over and over again while moving the mouse - // pointer inside the same pixel. That's especially relevant when operating - // on zoomed images. - if (index != m_LastPixelIndex) - { - // And finally... - mitk::ImagePixelWriteAccessor writeAccessor(image, image->GetVolumeData(timeStep)); - - // Paint all points between current and last pixel - auto indices = InterpolateIndices(m_LastPixelIndex, index, geometry, m_Size); - for (auto i : indices) - { - if (geometry->IsIndexInside(i)) - { - writeAccessor.SetPixelByIndexSafe(i, m_PaintingPixelValue); - } - } + // Okay, we're safe. Convert the mouse position to the index of the pixel + // we're pointing at. + itk::Index<3> index; + itk::Index<3> oldIndex; + geometry->WorldToIndex(position, index); + geometry->WorldToIndex(m_LastPosition, oldIndex); - image->Modified(); - this->GetDataNode()->Modified(); + // We don't need to paint over and over again while moving the mouse + // pointer inside the same pixel. That's especially relevant when operating + // on zoomed images. + if (index != oldIndex) + { + // Convert index to slice geometry + m_WorkingSlice->GetGeometry()->WorldToIndex(position, index); + m_WorkingSlice->GetGeometry()->WorldToIndex(m_LastPosition, oldIndex); + itk::Index<2> indexInPlane2D; + itk::Index<2> oldIndexInPlane2D; + indexInPlane2D[0] = index[0]; + indexInPlane2D[1] = index[1]; + oldIndexInPlane2D[0] = oldIndex[0]; + oldIndexInPlane2D[1] = oldIndex[1]; + + // Get indices + auto indices = InterpolateIndices2D(oldIndexInPlane2D, indexInPlane2D, m_WorkingPlane, m_Size); + + // Fill indices + mitk::ImagePixelWriteAccessor writeAccessor(m_WorkingSlice, m_WorkingSlice->GetVolumeData(timeStep)); + for (auto i : indices) + { + writeAccessor.SetPixelByIndexSafe(i, m_PaintingPixelValue); + } - mitk::RenderingManager::GetInstance()->RequestUpdateAll(); - m_LastPixelIndex = index; - m_Used = true; - } + // Update + image->Modified(); + this->GetDataNode()->Modified(); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + m_LastPosition = position; + m_Used = true; + } } - catch (...) + catch (itk::ExceptionObject& e) { - return; + mitkThrow() << "Could not paint with interpolation, because:"; + mitkThrow() << e.GetDescription(); } }