diff --git a/Modules/Classification/CLActiveLearning/src/mitkActiveLearningInteractor.cpp b/Modules/Classification/CLActiveLearning/src/mitkActiveLearningInteractor.cpp index 8e4fe02c60..3abcaba547 100644 --- a/Modules/Classification/CLActiveLearning/src/mitkActiveLearningInteractor.cpp +++ b/Modules/Classification/CLActiveLearning/src/mitkActiveLearningInteractor.cpp @@ -1,551 +1,599 @@ /*=================================================================== 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 #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 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 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(); 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); // 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); } // 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; } // 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); } } } } 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_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("set_workingslice", SetWorkingSlice) CONNECT_FUNCTION("writeback_workingslice", WriteBackWorkingSlice) } void mitk::ActiveLearningInteractor::DataNodeChanged() { this->ResetToStartState(); } 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(); - } - +// 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(); +// } +// 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(); - } +// 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) 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(); if (!geometry->IsInside(position)) return; + // 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(planeGeometry); +// extract->SetResliceTransformByGeometry(geometry); +// extract->Update(); +// auto workingSlice = extract->GetOutput()->Clone(); + + vtkSmartPointer reslice = vtkSmartPointer::New(); + //set to false to extract a slice + reslice->SetOverwriteMode(false); + reslice->Modified(); + + //use ExtractSliceFilter with our specific vtkImageReslice for overwriting and extracting + 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()->Clone(); + + mitk::IOUtil::Save(workingSlice, "/home/jens/Desktop/emptyslice.nrrd"); + // 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 != oldIndex) { // Convert index to slice geometry - m_WorkingSlice->GetGeometry()->WorldToIndex(position, index); + 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); + auto indices = InterpolateIndices2D(indexInPlane2D, indexInPlane2D, planeGeometry, m_Size); + + MITK_INFO << indices.size(); // Fill indices - mitk::ImagePixelWriteAccessor writeAccessor(m_WorkingSlice, m_WorkingSlice->GetVolumeData(timeStep)); + mitk::ImagePixelWriteAccessor writeAccessor(workingSlice, workingSlice->GetSliceData(0)); + + MITK_INFO << writeAccessor.GetPixelByIndexSafe({0,0}); + for (auto i : indices) { writeAccessor.SetPixelByIndexSafe(i, m_PaintingPixelValue); } + vtkSmartPointer reslice = vtkSmartPointer::New(); + reslice->SetInputSlice(workingSlice->GetVtkImageData()); + reslice->SetOverwriteMode(true); + reslice->Modified(); + mitk::ExtractSliceFilter::Pointer insert = mitk::ExtractSliceFilter::New(reslice); + insert->SetInput(image); + insert->SetTimeStep(timeStep); + insert->SetWorldGeometry(planeGeometry); + insert->SetResliceTransformByGeometry(geometry); + insert->Modified(); + insert->Update(); + // Update image->Modified(); this->GetDataNode()->Modified(); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); m_LastPosition = position; m_Used = true; } } catch (itk::ExceptionObject& e) { mitkThrow() << "Could not paint, because:"; mitkThrow() << e.GetDescription(); } } void mitk::ActiveLearningInteractor::PaintInterpolate(mitk::StateMachineAction* /*action*/, mitk::InteractionEvent* event) { 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 != 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); } // Update image->Modified(); this->GetDataNode()->Modified(); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); m_LastPosition = position; m_Used = true; } } catch (itk::ExceptionObject& e) { mitkThrow() << "Could not paint with interpolation, because:"; mitkThrow() << e.GetDescription(); } }