diff --git a/Modules/Core/include/mitkPointSetVtkMapper2D.h b/Modules/Core/include/mitkPointSetVtkMapper2D.h index 56abd70e00..adcb44f8c8 100644 --- a/Modules/Core/include/mitkPointSetVtkMapper2D.h +++ b/Modules/Core/include/mitkPointSetVtkMapper2D.h @@ -1,224 +1,224 @@ /*=================================================================== 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 mitkPointSetVtkMapper2D_h #define mitkPointSetVtkMapper2D_h #include #include #include "mitkVtkMapper.h" #include "mitkBaseRenderer.h" #include "mitkLocalStorageHandler.h" //VTK #include class vtkActor; class vtkPropAssembly; class vtkPolyData; class vtkPolyDataMapper; class vtkGlyphSource2D; class vtkGlyph3D; class vtkFloatArray; class vtkCellArray; namespace mitk { class PointSet; /** * @brief Vtk-based 2D mapper for PointSet * * Due to the need of different colors for selected * and unselected points and the facts, that we also have a contour and * labels for the points, the vtk structure is build up the following way: * * We have three PolyData, one selected, and one unselected and one * for a contour between the points. Each one is connected to an own * PolyDataMapper and an Actor. The different color for the unselected and * selected state and for the contour is read from properties. * * This mapper has several additional functionalities, such as rendering * a contour between points, calculating and displaying distances or angles * between points. * * Then the three Actors are combined inside a vtkPropAssembly and this * object is returned in GetProp() and so hooked up into the rendering * pipeline. * Properties that can be set for point sets and influence the PointSetVTKMapper2D are: * * - \b "line width": (IntProperty 2) // line width of the line from one point to another * - \b "point line width": (IntProperty 1) // line width of the cross marking a point - * - \b "point 2D size": (IntProperty 6) // size of the glyph marking a point + * - \b "point 2D size": (FloatProperty 6) // size of the glyph marking a point * - \b "show contour": (BoolProperty false) // enable contour rendering between points (lines) * - \b "close contour": (BoolProperty false) // if enabled, the open strip is closed (first point connected with last point) * - \b "show points": (BoolProperty true) // show or hide points * - \b "show distances": (BoolProperty false) // show or hide distance measure * - \b "distance decimal digits": (IntProperty 2) // set the number of decimal digits to be shown when rendering the distance information * - \b "show angles": (BoolProperty false) // show or hide angle measurement * - \b "show distant lines": (BoolProperty false) // show the line between to points from a distant view (equals "always on top" option) * - \b "layer": (IntProperty 1) // default is drawing pointset above images (they have a default layer of 0) * - \b "PointSet.2D.shape" (EnumerationProperty Cross) // provides different shapes marking a point * 0 = "None", 1 = "Vertex", 2 = "Dash", 3 = "Cross", 4 = "ThickCross", 5 = "Triangle", 6 = "Square", 7 = "Circle", * 8 = "Diamond", 9 = "Arrow", 10 = "ThickArrow", 11 = "HookedArrow", 12 = "Cross" * - \b "PointSet.2D.fill shape": (BoolProperty false) // fill or do not fill the glyph shape * - \b "Pointset.2D.distance to plane": (FloatProperty 4.0) //In the 2D render window, points are rendered which lie within a certain distance * to the current plane. They are projected on the current plane and scalled according to their distance. * Point markers appear smaller as the plane moves away from their true location. * The distance threshold can be adjusted by this float property, which ables the user to delineate the points * that lie exactly on the plane. (+/- rounding error) * * Other Properties used here but not defined in this class: * * - \b "selectedcolor": (ColorProperty (1.0f, 0.0f, 0.0f)) // default color of the selected pointset e.g. the current point is red * - \b "contourcolor" : (ColorProperty (1.0f, 0.0f, 0.0f)) // default color for the contour is red * - \b "color": (ColorProperty (1.0f, 1.0f, 0.0f)) // default color of the (unselected) pointset is yellow * - \b "opacity": (FloatProperty 1.0) // opacity of point set, contours * - \b "label": (StringProperty NULL) // a label can be defined for each point, which is rendered in proximity to the point * * @ingroup Mapper */ class MITKCORE_EXPORT PointSetVtkMapper2D : public VtkMapper { public: mitkClassMacro(PointSetVtkMapper2D, VtkMapper); itkFactorylessNewMacro(Self) itkCloneMacro(Self) virtual const mitk::PointSet* GetInput() const; /** \brief returns the a prop assembly */ virtual vtkProp* GetVtkProp(mitk::BaseRenderer* renderer) override; /** \brief set the default properties for this mapper */ static void SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer = NULL, bool overwrite = false); /** \brief Internal class holding the mapper, actor, etc. for each of the 3 2D render windows */ class LocalStorage : public mitk::Mapper::BaseLocalStorage { public: /* constructor */ LocalStorage(); /* destructor */ ~LocalStorage(); // points vtkSmartPointer m_UnselectedPoints; vtkSmartPointer m_SelectedPoints; vtkSmartPointer m_ContourPoints; // scales vtkSmartPointer m_UnselectedScales; vtkSmartPointer m_SelectedScales; // distances vtkSmartPointer m_DistancesBetweenPoints; // lines vtkSmartPointer m_ContourLines; // glyph source (provides different shapes for the points) vtkSmartPointer m_UnselectedGlyphSource2D; vtkSmartPointer m_SelectedGlyphSource2D; // glyph vtkSmartPointer m_UnselectedGlyph3D; vtkSmartPointer m_SelectedGlyph3D; // polydata vtkSmartPointer m_VtkUnselectedPointListPolyData; vtkSmartPointer m_VtkSelectedPointListPolyData; vtkSmartPointer m_VtkContourPolyData; // actor vtkSmartPointer m_UnselectedActor; vtkSmartPointer m_SelectedActor; vtkSmartPointer m_ContourActor; vtkSmartPointer m_VtkTextActor; std::vector < vtkSmartPointer > m_VtkTextLabelActors; std::vector < vtkSmartPointer > m_VtkTextDistanceActors; std::vector < vtkSmartPointer > m_VtkTextAngleActors; // mappers vtkSmartPointer m_VtkUnselectedPolyDataMapper; vtkSmartPointer m_VtkSelectedPolyDataMapper; vtkSmartPointer m_VtkContourPolyDataMapper; // propassembly vtkSmartPointer m_PropAssembly; }; /** \brief The LocalStorageHandler holds all (three) LocalStorages for the three 2D render windows. */ mitk::LocalStorageHandler m_LSH; protected: /* constructor */ PointSetVtkMapper2D(); /* destructor */ virtual ~PointSetVtkMapper2D(); /* \brief Applies the color and opacity properties and calls CreateVTKRenderObjects */ virtual void GenerateDataForRenderer(mitk::BaseRenderer* renderer) override; /* \brief Called in mitk::Mapper::Update * If TimeSlicedGeometry or time step is not valid of point set: reset mapper so that nothing is * displayed e.g. toggle visiblity of the propassembly */ virtual void ResetMapper( BaseRenderer* renderer ) override; /* \brief Fills the vtk objects, thus it is only called when the point set has been changed. * This function iterates over the input point set and determines the glyphs which lie in a specific * range around the current slice. Those glyphs are rendered using a specific shape defined in vtk glyph source * to mark each point. The shape can be changed in MITK using the property "PointSet.2D.shape". * * There were issues when rendering vtk glyphs in the 2D-render windows. By default, the glyphs are * rendered within the x-y plane in each 2D-render window, so you would only see them from the * side in the saggital and coronal 2D-render window. The solution to this is to rotate the glyphs in order * to be ortogonal to the current view vector. To achieve this, the rotation (vtktransform) of the current * PlaneGeometry is applied to the orienation of the glyphs. */ virtual void CreateVTKRenderObjects(mitk::BaseRenderer* renderer); // member variables holding the current value of the properties used in this mapper bool m_ShowContour; // "show contour" property bool m_CloseContour; // "close contour" property bool m_ShowPoints; // "show points" property bool m_ShowDistances; // "show distances" property int m_DistancesDecimalDigits; // "distance decimal digits" property bool m_ShowAngles; // "show angles" property bool m_ShowDistantLines; // "show distant lines" property int m_LineWidth; // "line width" property int m_PointLineWidth; // "point line width" property - int m_Point2DSize; // "point 2D size" property + float m_Point2DSize; // "point 2D size" property int m_IDShapeProperty; // ID for mitkPointSetShape Enumeration Property "Pointset.2D.shape" bool m_FillShape; // "Pointset.2D.fill shape" property float m_DistanceToPlane; // "Pointset.2D.distance to plane" property }; } // namespace mitk #endif /* mitkPointSetVtkMapper2D_h */ diff --git a/Modules/Core/src/Rendering/mitkPointSetVtkMapper2D.cpp b/Modules/Core/src/Rendering/mitkPointSetVtkMapper2D.cpp index 833c828976..7afb1a805b 100644 --- a/Modules/Core/src/Rendering/mitkPointSetVtkMapper2D.cpp +++ b/Modules/Core/src/Rendering/mitkPointSetVtkMapper2D.cpp @@ -1,713 +1,718 @@ /*=================================================================== 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 "mitkPointSetVtkMapper2D.h" //mitk includes #include "mitkVtkPropRenderer.h" #include #include #include #include //vtk includes #include #include #include #include #include #include #include #include #include #include #include #include #include #include // constructor LocalStorage mitk::PointSetVtkMapper2D::LocalStorage::LocalStorage() { // points m_UnselectedPoints = vtkSmartPointer::New(); m_SelectedPoints = vtkSmartPointer::New(); m_ContourPoints = vtkSmartPointer::New(); // scales m_UnselectedScales = vtkSmartPointer::New(); m_SelectedScales = vtkSmartPointer::New(); // distances m_DistancesBetweenPoints = vtkSmartPointer::New(); // lines m_ContourLines = vtkSmartPointer::New(); // glyph source (provides the different shapes) m_UnselectedGlyphSource2D = vtkSmartPointer::New(); m_SelectedGlyphSource2D = vtkSmartPointer::New(); // glyphs m_UnselectedGlyph3D = vtkSmartPointer::New(); m_SelectedGlyph3D = vtkSmartPointer::New(); // polydata m_VtkUnselectedPointListPolyData = vtkSmartPointer::New(); m_VtkSelectedPointListPolyData = vtkSmartPointer ::New(); m_VtkContourPolyData = vtkSmartPointer::New(); // actors m_UnselectedActor = vtkSmartPointer ::New(); m_SelectedActor = vtkSmartPointer ::New(); m_ContourActor = vtkSmartPointer ::New(); // mappers m_VtkUnselectedPolyDataMapper = vtkSmartPointer::New(); m_VtkSelectedPolyDataMapper = vtkSmartPointer::New(); m_VtkContourPolyDataMapper = vtkSmartPointer::New(); // propassembly m_PropAssembly = vtkSmartPointer ::New(); } // destructor LocalStorage mitk::PointSetVtkMapper2D::LocalStorage::~LocalStorage() { } // input for this mapper ( = point set) const mitk::PointSet* mitk::PointSetVtkMapper2D::GetInput() const { return static_cast ( GetDataNode()->GetData() ); } // constructor PointSetVtkMapper2D mitk::PointSetVtkMapper2D::PointSetVtkMapper2D() : m_ShowContour(false), m_CloseContour(false), m_ShowPoints(true), m_ShowDistances(false), m_DistancesDecimalDigits(1), m_ShowAngles(false), m_ShowDistantLines(false), m_LineWidth(1), m_PointLineWidth(1), m_Point2DSize(6), m_IDShapeProperty(mitk::PointSetShapeProperty::CROSS), m_FillShape(false), m_DistanceToPlane(4.0f) { } // destructor mitk::PointSetVtkMapper2D::~PointSetVtkMapper2D() { } // reset mapper so that nothing is displayed e.g. toggle visiblity of the propassembly void mitk::PointSetVtkMapper2D::ResetMapper( BaseRenderer* renderer ) { LocalStorage *ls = m_LSH.GetLocalStorage(renderer); ls->m_PropAssembly->VisibilityOff(); } // returns propassembly vtkProp* mitk::PointSetVtkMapper2D::GetVtkProp(mitk::BaseRenderer * renderer) { LocalStorage *ls = m_LSH.GetLocalStorage(renderer); return ls->m_PropAssembly; } static bool makePerpendicularVector2D(const mitk::Vector2D& in, mitk::Vector2D& out) { // The dot product of orthogonal vectors is zero. // In two dimensions the slopes of perpendicular lines are negative reciprocals. if((fabs(in[0])>0) && ( (fabs(in[0])>fabs(in[1])) || (in[1] == 0) ) ) { // negative reciprocal out[0]=-in[1]/in[0]; out[1]=1; out.Normalize(); return true; } else if(fabs(in[1])>0) { out[0]=1; // negative reciprocal out[1]=-in[0]/in[1]; out.Normalize(); return true; } else return false; } void mitk::PointSetVtkMapper2D::CreateVTKRenderObjects(mitk::BaseRenderer* renderer) { LocalStorage *ls = m_LSH.GetLocalStorage(renderer); unsigned i = 0; // The vtk text actors need to be removed manually from the propassembly // since the same vtk text actors are not overwriten within this function, // but new actors are added to the propassembly each time this function is executed. // Thus, the actors from the last call must be removed in the beginning. for(i=0; i< ls->m_VtkTextLabelActors.size(); i++) { if(ls->m_PropAssembly->GetParts()->IsItemPresent(ls->m_VtkTextLabelActors.at(i))) ls->m_PropAssembly->RemovePart(ls->m_VtkTextLabelActors.at(i)); } for(i=0; i< ls->m_VtkTextDistanceActors.size(); i++) { if(ls->m_PropAssembly->GetParts()->IsItemPresent(ls->m_VtkTextDistanceActors.at(i))) ls->m_PropAssembly->RemovePart(ls->m_VtkTextDistanceActors.at(i)); } for(i=0; i< ls->m_VtkTextAngleActors.size(); i++) { if(ls->m_PropAssembly->GetParts()->IsItemPresent(ls->m_VtkTextAngleActors.at(i))) ls->m_PropAssembly->RemovePart(ls->m_VtkTextAngleActors.at(i)); } // initialize polydata here, otherwise we have update problems when // executing this function again ls->m_VtkUnselectedPointListPolyData = vtkSmartPointer::New(); ls->m_VtkSelectedPointListPolyData = vtkSmartPointer ::New(); ls->m_VtkContourPolyData = vtkSmartPointer::New(); // get input point set and update the PointSet mitk::PointSet::Pointer input = const_cast(this->GetInput()); // only update the input data, if the property tells us to bool update = true; this->GetDataNode()->GetBoolProperty("updateDataOnRender", update); if (update == true) input->Update(); int timestep = this->GetTimestep(); mitk::PointSet::DataType::Pointer itkPointSet = input->GetPointSet( timestep ); if ( itkPointSet.GetPointer() == NULL) { ls->m_PropAssembly->VisibilityOff(); return; } //iterator for point set mitk::PointSet::PointsContainer::Iterator pointsIter = itkPointSet->GetPoints()->Begin(); // PointDataContainer has additional information to each point, e.g. whether // it is selected or not mitk::PointSet::PointDataContainer::Iterator pointDataIter; pointDataIter = itkPointSet->GetPointData()->Begin(); //check if the list for the PointDataContainer is the same size as the PointsContainer. //If not, then the points were inserted manually and can not be visualized according to the PointData (selected/unselected) bool pointDataBroken = (itkPointSet->GetPointData()->Size() != itkPointSet->GetPoints()->Size()); if( itkPointSet->GetPointData()->size() == 0 || pointDataBroken) { ls->m_PropAssembly->VisibilityOff(); return; } ls->m_PropAssembly->VisibilityOn(); // empty point sets, cellarrays, scalars ls->m_UnselectedPoints->Reset(); ls->m_SelectedPoints->Reset(); ls->m_ContourPoints->Reset(); ls->m_ContourLines->Reset(); ls->m_UnselectedScales->Reset(); ls->m_SelectedScales->Reset(); ls->m_DistancesBetweenPoints->Reset(); ls->m_VtkTextLabelActors.clear(); ls->m_VtkTextDistanceActors.clear(); ls->m_VtkTextAngleActors.clear(); ls->m_UnselectedScales->SetNumberOfComponents(3); ls->m_SelectedScales->SetNumberOfComponents(3); int NumberContourPoints = 0; bool pointsOnSameSideOfPlane = false; const int text2dDistance = 10; // initialize points with a random start value // current point in point set itk::Point point = pointsIter->Value(); mitk::Point3D p = point; // currently visited point mitk::Point3D lastP = point; // last visited point (predecessor in point set of "point") mitk::Vector3D vec; // p - lastP mitk::Vector3D lastVec; // lastP - point before lastP vec.Fill(0.0); lastVec.Fill(0.0); mitk::Point3D projected_p = point; // p projected on viewplane mitk::Point2D pt2d; pt2d[0] = point[0]; // projected_p in display coordinates pt2d[1] = point[1]; mitk::Point2D lastPt2d = pt2d; // last projected_p in display coordinates (predecessor in point set of "pt2d") mitk::Point2D preLastPt2d = pt2d ; // projected_p in display coordinates before lastPt2 mitk::DisplayGeometry::Pointer displayGeometry = renderer->GetDisplayGeometry(); const mitk::PlaneGeometry* geo2D = renderer->GetCurrentWorldPlaneGeometry(); vtkLinearTransform* dataNodeTransform = input->GetGeometry()->GetVtkTransform(); int count = 0; for (pointsIter=itkPointSet->GetPoints()->Begin(); pointsIter!=itkPointSet->GetPoints()->End(); pointsIter++) { lastP = p; // valid for number of points count > 0 preLastPt2d = lastPt2d; // valid only for count > 1 lastPt2d = pt2d; // valid for number of points count > 0 lastVec = vec; // valid only for counter > 1 // get current point in point set point = pointsIter->Value(); // transform point { float vtkp[3]; itk2vtk(point, vtkp); dataNodeTransform->TransformPoint(vtkp, vtkp); vtk2itk(vtkp,point); } p[0] = point[0]; p[1] = point[1]; p[2] = point[2]; displayGeometry->Project(p, projected_p); displayGeometry->Map(projected_p, pt2d); displayGeometry->WorldToDisplay(pt2d, pt2d); vec = p-lastP; // valid only for counter > 0 // compute distance to current plane - float diff = geo2D->Distance(point); - diff = diff * diff; + float diff = fabs(geo2D->Distance(point)); +// diff = diff * diff; //draw markers on slices a certain distance away from the points //location according to the tolerance threshold (m_DistanceToPlane) if(diff < m_DistanceToPlane) { // is point selected or not? if (pointDataIter->Value().selected) { ls->m_SelectedPoints->InsertNextPoint(point[0],point[1],point[2]); // point is scaled according to its distance to the plane - ls->m_SelectedScales->InsertNextTuple3(m_Point2DSize - (2*diff),0,0); + ls->m_SelectedScales->InsertNextTuple3((double)m_Point2DSize*(1-diff/m_DistanceToPlane),0,0); } else { ls->m_UnselectedPoints->InsertNextPoint(point[0],point[1],point[2]); // point is scaled according to its distance to the plane - ls->m_UnselectedScales->InsertNextTuple3(m_Point2DSize - (2*diff),0,0); + ls->m_UnselectedScales->InsertNextTuple3((double)m_Point2DSize*(1-diff/m_DistanceToPlane),0,0); } //---- LABEL -----// //paint label for each point if available if (dynamic_cast(this->GetDataNode()->GetProperty("label")) != NULL) { const char * pointLabel = dynamic_cast( this->GetDataNode()->GetProperty("label"))->GetValue(); std::string l = pointLabel; if (input->GetSize()>1) { std::stringstream ss; ss << pointsIter->Index(); l.append(ss.str()); } ls->m_VtkTextActor = vtkSmartPointer::New(); ls->m_VtkTextActor->SetPosition(pt2d[0] + text2dDistance, pt2d[1] + text2dDistance); ls->m_VtkTextActor->SetInput(l.c_str()); ls->m_VtkTextActor->GetTextProperty()->SetOpacity( 100 ); float unselectedColor[4] = {1.0, 1.0, 0.0, 1.0}; //check if there is a color property GetDataNode()->GetColor(unselectedColor); ls->m_VtkTextActor->GetTextProperty()->SetColor(unselectedColor[0], unselectedColor[1], unselectedColor[2]); ls->m_VtkTextLabelActors.push_back(ls->m_VtkTextActor); } } // draw contour, distance text and angle text in render window // lines between points, which intersect the current plane, are drawn if( m_ShowContour && count > 0 ) { ScalarType distance = displayGeometry->GetWorldGeometry()->SignedDistance(point); ScalarType lastDistance = displayGeometry->GetWorldGeometry()->SignedDistance(lastP); pointsOnSameSideOfPlane = (distance * lastDistance) > 0.5; // Points must be on different side of plane in order to draw a contour. // If "show distant lines" is enabled this condition is disregarded. if ( !pointsOnSameSideOfPlane || m_ShowDistantLines) { vtkSmartPointer line = vtkSmartPointer::New(); ls->m_ContourPoints->InsertNextPoint(lastP[0],lastP[1],lastP[2]); line->GetPointIds()->SetId(0, NumberContourPoints); NumberContourPoints++; ls->m_ContourPoints->InsertNextPoint(point[0], point[1], point[2]); line->GetPointIds()->SetId(1, NumberContourPoints); NumberContourPoints++; ls->m_ContourLines->InsertNextCell(line); if(m_ShowDistances) // calculate and print distance between adjacent points { float distancePoints = point.EuclideanDistanceTo(lastP); std::stringstream buffer; buffer<m_VtkTextActor = vtkSmartPointer::New(); ls->m_VtkTextActor->SetPosition(pos2d[0],pos2d[1]); ls->m_VtkTextActor->SetInput(buffer.str().c_str()); ls->m_VtkTextActor->GetTextProperty()->SetColor(0.0, 1.0, 0.0); ls->m_VtkTextDistanceActors.push_back(ls->m_VtkTextActor); } if(m_ShowAngles && count > 1) // calculate and print angle between connected lines { std::stringstream buffer; buffer << angle(vec.GetVnlVector(), -lastVec.GetVnlVector())*180/vnl_math::pi << "°"; //compute desired display position of text Vector2D vec2d = pt2d-lastPt2d; // first arm enclosing the angle vec2d.Normalize(); Vector2D lastVec2d = lastPt2d-preLastPt2d; // second arm enclosing the angle lastVec2d.Normalize(); vec2d=vec2d-lastVec2d; // vector connecting both arms vec2d.Normalize(); // middle between two vectors that enclose the angle Vector2D pos2d = lastPt2d.GetVectorFromOrigin() + vec2d * text2dDistance * text2dDistance; ls->m_VtkTextActor = vtkSmartPointer::New(); ls->m_VtkTextActor->SetPosition(pos2d[0],pos2d[1]); ls->m_VtkTextActor->SetInput(buffer.str().c_str()); ls->m_VtkTextActor->GetTextProperty()->SetColor(0.0, 1.0, 0.0); ls->m_VtkTextAngleActors.push_back(ls->m_VtkTextActor); } } } if(pointDataIter != itkPointSet->GetPointData()->End()) { pointDataIter++; count++; } } // add each single text actor to the assembly for(i=0; i< ls->m_VtkTextLabelActors.size(); i++) { ls->m_PropAssembly->AddPart(ls->m_VtkTextLabelActors.at(i)); } for(i=0; i< ls->m_VtkTextDistanceActors.size(); i++) { ls->m_PropAssembly->AddPart(ls->m_VtkTextDistanceActors.at(i)); } for(i=0; i< ls->m_VtkTextAngleActors.size(); i++) { ls->m_PropAssembly->AddPart(ls->m_VtkTextAngleActors.at(i)); } //---- CONTOUR -----// //create lines between the points which intersect the plane if (m_ShowContour) { // draw line between first and last point which is rendered if(m_CloseContour && NumberContourPoints > 1){ vtkSmartPointer closingLine = vtkSmartPointer::New(); closingLine->GetPointIds()->SetId(0, 0); // index of first point closingLine->GetPointIds()->SetId(1, NumberContourPoints-1); // index of last point ls->m_ContourLines->InsertNextCell(closingLine); } ls->m_VtkContourPolyData->SetPoints(ls->m_ContourPoints); ls->m_VtkContourPolyData->SetLines(ls->m_ContourLines); ls->m_VtkContourPolyDataMapper->SetInputData(ls->m_VtkContourPolyData); ls->m_ContourActor->SetMapper(ls->m_VtkContourPolyDataMapper); ls->m_ContourActor->GetProperty()->SetLineWidth(m_LineWidth); ls->m_PropAssembly->AddPart(ls->m_ContourActor); } // the point set must be transformed in order to obtain the appropriate glyph orientation // according to the current view vtkSmartPointer transform = vtkSmartPointer::New(); vtkSmartPointer a,b = vtkSmartPointer::New(); a = geo2D->GetVtkTransform()->GetMatrix(); b->DeepCopy( a ); // delete transformation from matrix, only take orientation b->SetElement(3,3,1); b->SetElement(2,3,0); b->SetElement(1,3,0); b->SetElement(0,3,0); b->SetElement(3,2,0); b->SetElement(3,1,0); b->SetElement(3,0,0); transform->SetMatrix( b ); //---- UNSELECTED POINTS -----// // apply properties to glyph ls->m_UnselectedGlyphSource2D->SetGlyphType(m_IDShapeProperty); if(m_FillShape) ls->m_UnselectedGlyphSource2D->FilledOn(); else ls->m_UnselectedGlyphSource2D->FilledOff(); // apply transform vtkSmartPointer transformFilterU = vtkSmartPointer::New(); transformFilterU->SetInputConnection(ls->m_UnselectedGlyphSource2D->GetOutputPort()); transformFilterU->SetTransform(transform); ls->m_VtkUnselectedPointListPolyData->SetPoints(ls->m_UnselectedPoints); ls->m_VtkUnselectedPointListPolyData->GetPointData()->SetVectors(ls->m_UnselectedScales); // apply transform of current plane to glyphs ls->m_UnselectedGlyph3D->SetSourceConnection(transformFilterU->GetOutputPort()); ls->m_UnselectedGlyph3D->SetInputData(ls->m_VtkUnselectedPointListPolyData); ls->m_UnselectedGlyph3D->SetScaleModeToScaleByVector(); ls->m_UnselectedGlyph3D->SetVectorModeToUseVector(); ls->m_VtkUnselectedPolyDataMapper->SetInputConnection(ls->m_UnselectedGlyph3D->GetOutputPort()); ls->m_UnselectedActor->SetMapper(ls->m_VtkUnselectedPolyDataMapper); ls->m_UnselectedActor->GetProperty()->SetLineWidth(m_PointLineWidth); ls->m_PropAssembly->AddPart(ls->m_UnselectedActor); //---- SELECTED POINTS -----// ls->m_SelectedGlyphSource2D->SetGlyphTypeToDiamond(); ls->m_SelectedGlyphSource2D->CrossOn(); ls->m_SelectedGlyphSource2D->FilledOff(); // apply transform vtkSmartPointer transformFilterS = vtkSmartPointer::New(); transformFilterS->SetInputConnection(ls->m_SelectedGlyphSource2D->GetOutputPort()); transformFilterS->SetTransform(transform); ls->m_VtkSelectedPointListPolyData->SetPoints(ls->m_SelectedPoints); ls->m_VtkSelectedPointListPolyData->GetPointData()->SetVectors(ls->m_SelectedScales); // apply transform of current plane to glyphs ls->m_SelectedGlyph3D->SetSourceConnection(transformFilterS->GetOutputPort()); ls->m_SelectedGlyph3D->SetInputData(ls->m_VtkSelectedPointListPolyData); ls->m_SelectedGlyph3D->SetScaleModeToScaleByVector(); ls->m_SelectedGlyph3D->SetVectorModeToUseVector(); ls->m_VtkSelectedPolyDataMapper->SetInputConnection(ls->m_SelectedGlyph3D->GetOutputPort()); ls->m_SelectedActor->SetMapper(ls->m_VtkSelectedPolyDataMapper); ls->m_SelectedActor->GetProperty()->SetLineWidth(m_PointLineWidth); ls->m_PropAssembly->AddPart(ls->m_SelectedActor); } void mitk::PointSetVtkMapper2D::GenerateDataForRenderer( mitk::BaseRenderer *renderer ) { const mitk::DataNode* node = GetDataNode(); if( node == NULL ) return; LocalStorage *ls = m_LSH.GetLocalStorage(renderer); // check whether the input data has been changed bool needGenerateData = ls->IsGenerateDataRequired( renderer, this, GetDataNode() ); // toggle visibility bool visible = true; node->GetVisibility(visible, renderer, "visible"); if(!visible) { ls->m_UnselectedActor->VisibilityOff(); ls->m_SelectedActor->VisibilityOff(); ls->m_ContourActor->VisibilityOff(); ls->m_PropAssembly->VisibilityOff(); return; }else{ ls->m_PropAssembly->VisibilityOn(); } node->GetBoolProperty("show contour", m_ShowContour, renderer); node->GetBoolProperty("close contour", m_CloseContour, renderer); node->GetBoolProperty("show points", m_ShowPoints, renderer); node->GetBoolProperty("show distances", m_ShowDistances, renderer); node->GetIntProperty("distance decimal digits", m_DistancesDecimalDigits, renderer); node->GetBoolProperty("show angles", m_ShowAngles, renderer); node->GetBoolProperty("show distant lines", m_ShowDistantLines, renderer); node->GetIntProperty("line width", m_LineWidth, renderer); node->GetIntProperty("point line width", m_PointLineWidth, renderer); - node->GetIntProperty("point 2D size", m_Point2DSize, renderer); + if(!node->GetFloatProperty("point 2D size", m_Point2DSize, renderer)) + { + int size = 0; + if(node->GetIntProperty("point 2D size", size, renderer)) + m_Point2DSize = size; + } node->GetBoolProperty("Pointset.2D.fill shape", m_FillShape, renderer); node->GetFloatProperty("Pointset.2D.distance to plane", m_DistanceToPlane, renderer ); mitk::PointSetShapeProperty::Pointer shape = dynamic_cast(this->GetDataNode()->GetProperty( "Pointset.2D.shape", renderer )); if(shape.IsNotNull()) { m_IDShapeProperty = shape->GetPointSetShape(); } //check for color props and use it for rendering of selected/unselected points and contour //due to different params in VTK (double/float) we have to convert float unselectedColor[4]; double selectedColor[4]={1.0f,0.0f,0.0f,1.0f}; //red double contourColor[4]={1.0f,0.0f,0.0f,1.0f}; //red float opacity = 1.0; GetDataNode()->GetOpacity(opacity, renderer); // apply color and opacity if(m_ShowPoints) { ls->m_UnselectedActor->VisibilityOn(); ls->m_SelectedActor->VisibilityOn(); //check if there is a color property GetDataNode()->GetColor(unselectedColor); //get selected color property if (dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("selectedcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("selectedcolor"))->GetValue(); selectedColor[0] = tmpColor[0]; selectedColor[1] = tmpColor[1]; selectedColor[2] = tmpColor[2]; selectedColor[3] = 1.0f; // alpha value } else if (dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("selectedcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("selectedcolor"))->GetValue(); selectedColor[0] = tmpColor[0]; selectedColor[1] = tmpColor[1]; selectedColor[2] = tmpColor[2]; selectedColor[3] = 1.0f; // alpha value } ls->m_SelectedActor->GetProperty()->SetColor(selectedColor); ls->m_SelectedActor->GetProperty()->SetOpacity(opacity); ls->m_UnselectedActor->GetProperty()->SetColor(unselectedColor[0],unselectedColor[1],unselectedColor[2]); ls->m_UnselectedActor->GetProperty()->SetOpacity(opacity); } else { ls->m_UnselectedActor->VisibilityOff(); ls-> m_SelectedActor->VisibilityOff(); } if (m_ShowContour) { ls->m_ContourActor->VisibilityOn(); //get contour color property if (dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("contourcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("contourcolor"))->GetValue(); contourColor[0] = tmpColor[0]; contourColor[1] = tmpColor[1]; contourColor[2] = tmpColor[2]; contourColor[3] = 1.0f; } else if (dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("contourcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("contourcolor"))->GetValue(); contourColor[0] = tmpColor[0]; contourColor[1] = tmpColor[1]; contourColor[2] = tmpColor[2]; contourColor[3] = 1.0f; } ls->m_ContourActor->GetProperty()->SetColor(contourColor); ls->m_ContourActor->GetProperty()->SetOpacity(opacity); } else { ls->m_ContourActor->VisibilityOff(); } if(needGenerateData) { // create new vtk render objects (e.g. a circle for a point) this->CreateVTKRenderObjects(renderer); } } void mitk::PointSetVtkMapper2D::SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer, bool overwrite) { node->AddProperty( "line width", mitk::IntProperty::New(2), renderer, overwrite ); node->AddProperty( "point line width", mitk::IntProperty::New(1), renderer, overwrite ); - node->AddProperty( "point 2D size", mitk::IntProperty::New(6), renderer, overwrite ); + node->AddProperty( "point 2D size", mitk::FloatProperty::New(6), renderer, overwrite ); node->AddProperty( "show contour", mitk::BoolProperty::New(false), renderer, overwrite ); node->AddProperty( "close contour", mitk::BoolProperty::New(false), renderer, overwrite ); node->AddProperty( "show points", mitk::BoolProperty::New(true), renderer, overwrite ); node->AddProperty( "show distances", mitk::BoolProperty::New(false), renderer, overwrite ); node->AddProperty( "distance decimal digits", mitk::IntProperty::New(2), renderer, overwrite ); node->AddProperty( "show angles", mitk::BoolProperty::New(false), renderer, overwrite ); node->AddProperty( "show distant lines", mitk::BoolProperty::New(false), renderer, overwrite ); node->AddProperty( "layer", mitk::IntProperty::New(1), renderer, overwrite ); node->AddProperty( "Pointset.2D.fill shape", mitk::BoolProperty::New(false), renderer, overwrite); // fill or do not fill the glyph shape mitk::PointSetShapeProperty::Pointer pointsetShapeProperty = mitk::PointSetShapeProperty::New(); node->AddProperty( "Pointset.2D.shape", pointsetShapeProperty, renderer, overwrite); node->AddProperty( "Pointset.2D.distance to plane", mitk::FloatProperty::New(4.0f), renderer, overwrite ); //show the point at a certain distance above/below the 2D imaging plane. Superclass::SetDefaultProperties(node, renderer, overwrite); } diff --git a/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp b/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp index e492d58937..afee3da252 100644 --- a/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp +++ b/Modules/DiffusionImaging/DiffusionIO/mitkFiberBundleVtkReader.cpp @@ -1,277 +1,277 @@ /*=================================================================== 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 "mitkFiberBundleVtkReader.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkDiffusionIOMimeTypes.h" mitk::FiberBundleVtkReader::FiberBundleVtkReader() : mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::FIBERBUNDLE_VTK_MIMETYPE_NAME(), "VTK Fiber Bundle Reader" ) { m_ServiceReg = this->RegisterService(); } mitk::FiberBundleVtkReader::FiberBundleVtkReader(const FiberBundleVtkReader &other) :mitk::AbstractFileReader(other) { } mitk::FiberBundleVtkReader * mitk::FiberBundleVtkReader::Clone() const { return new FiberBundleVtkReader(*this); } std::vector > mitk::FiberBundleVtkReader::Read() { std::vector > result; const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); setlocale(LC_ALL, locale.c_str()); std::string filename = this->GetInputLocation(); std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); try { MITK_INFO << "Trying to load fiber file as VTK format."; vtkSmartPointer reader = vtkSmartPointer::New(); reader->SetFileName( this->GetInputLocation().c_str() ); if (reader->IsFilePolyData()) { reader->Update(); if ( reader->GetOutput() != NULL ) { vtkSmartPointer fiberPolyData = reader->GetOutput(); FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData); vtkSmartPointer weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS")); if (weights!=NULL) { -// float weight=0; -// for (int i=0; iGetSize(); i++) -// if (!mitk::Equal(weights->GetValue(i),weight,0.00001)) -// { -// MITK_INFO << "Weight: " << weights->GetValue(i); -// weight = weights->GetValue(i); -// } + float weight=0; + for (int i=0; iGetSize(); i++) + if (!mitk::Equal(weights->GetValue(i),weight,0.00001)) + { + MITK_INFO << "Weight: " << weights->GetValue(i); + weight = weights->GetValue(i); + } fiberBundle->SetFiberWeights(weights); } vtkSmartPointer fiberColors = vtkUnsignedCharArray::SafeDownCast(fiberPolyData->GetPointData()->GetArray("FIBER_COLORS")); if (fiberColors!=NULL) fiberBundle->SetFiberColors(fiberColors); result.push_back(fiberBundle.GetPointer()); return result; } } else MITK_INFO << "File is not VTK format."; } catch(...) { throw; } try { MITK_INFO << "Trying to load fiber file as VTP format."; vtkSmartPointer reader = vtkSmartPointer::New(); reader->SetFileName( this->GetInputLocation().c_str() ); if ( reader->CanReadFile(this->GetInputLocation().c_str()) ) { reader->Update(); if ( reader->GetOutput() != NULL ) { vtkSmartPointer fiberPolyData = reader->GetOutput(); FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData); vtkSmartPointer weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS")); if (weights!=NULL) { // float weight=0; // for (int i=0; iGetSize(); i++) // if (!mitk::Equal(weights->GetValue(i),weight,0.00001)) // { // MITK_INFO << "Weight: " << weights->GetValue(i); // weight = weights->GetValue(i); // } fiberBundle->SetFiberWeights(weights); } vtkSmartPointer fiberColors = vtkUnsignedCharArray::SafeDownCast(fiberPolyData->GetPointData()->GetArray("FIBER_COLORS")); if (fiberColors!=NULL) fiberBundle->SetFiberColors(fiberColors); result.push_back(fiberBundle.GetPointer()); return result; } } else MITK_INFO << "File is not VTP format."; } catch(...) { throw; } try { MITK_INFO << "Trying to load fiber file as deprecated XML format."; vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); vtkSmartPointer cellArray = vtkSmartPointer::New(); vtkSmartPointer points = vtkSmartPointer::New(); TiXmlDocument doc( this->GetInputLocation().c_str() ); if(doc.LoadFile()) { TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); // save this for later hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement("geometry").Element(); // read geometry mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); // read origin mitk::Point3D origin; double temp = 0; pElem->Attribute("origin_x", &temp); origin[0] = temp; pElem->Attribute("origin_y", &temp); origin[1] = temp; pElem->Attribute("origin_z", &temp); origin[2] = temp; geometry->SetOrigin(origin); // read spacing ScalarType spacing[3]; pElem->Attribute("spacing_x", &temp); spacing[0] = temp; pElem->Attribute("spacing_y", &temp); spacing[1] = temp; pElem->Attribute("spacing_z", &temp); spacing[2] = temp; geometry->SetSpacing(spacing); // read transform vtkMatrix4x4* m = vtkMatrix4x4::New(); pElem->Attribute("xx", &temp); m->SetElement(0,0,temp); pElem->Attribute("xy", &temp); m->SetElement(1,0,temp); pElem->Attribute("xz", &temp); m->SetElement(2,0,temp); pElem->Attribute("yx", &temp); m->SetElement(0,1,temp); pElem->Attribute("yy", &temp); m->SetElement(1,1,temp); pElem->Attribute("yz", &temp); m->SetElement(2,1,temp); pElem->Attribute("zx", &temp); m->SetElement(0,2,temp); pElem->Attribute("zy", &temp); m->SetElement(1,2,temp); pElem->Attribute("zz", &temp); m->SetElement(2,2,temp); m->SetElement(0,3,origin[0]); m->SetElement(1,3,origin[1]); m->SetElement(2,3,origin[2]); m->SetElement(3,3,1); geometry->SetIndexToWorldTransformByVtkMatrix(m); // read bounds float bounds[] = {0, 0, 0, 0, 0, 0}; pElem->Attribute("size_x", &temp); bounds[1] = temp; pElem->Attribute("size_y", &temp); bounds[3] = temp; pElem->Attribute("size_z", &temp); bounds[5] = temp; geometry->SetFloatBounds(bounds); geometry->SetImageGeometry(true); pElem = hRoot.FirstChildElement("fiber_bundle").FirstChild().Element(); for( ; pElem ; pElem=pElem->NextSiblingElement()) { TiXmlElement* pElem2 = pElem->FirstChildElement(); vtkSmartPointer container = vtkSmartPointer::New(); for( ; pElem2; pElem2=pElem2->NextSiblingElement()) { Point3D point; pElem2->Attribute("pos_x", &temp); point[0] = temp; pElem2->Attribute("pos_y", &temp); point[1] = temp; pElem2->Attribute("pos_z", &temp); point[2] = temp; geometry->IndexToWorld(point, point); vtkIdType id = points->InsertNextPoint(point.GetDataPointer()); container->GetPointIds()->InsertNextId(id); } cellArray->InsertNextCell(container); } fiberPolyData->SetPoints(points); fiberPolyData->SetLines(cellArray); vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInputData(fiberPolyData); cleaner->Update(); fiberPolyData = cleaner->GetOutput(); FiberBundle::Pointer image = FiberBundle::New(fiberPolyData); result.push_back(image.GetPointer()); return result; } else { MITK_INFO << "File is not deprectaed XML format."; } setlocale(LC_ALL, currLocale.c_str()); MITK_INFO << "Fiber bundle read"; } catch(...) { throw; } throw "Selected file is no vtk readable fiber format (binary or ascii vtk or vtp file)."; return result; } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.cpp index da4a53b7ed..6e66836263 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.cpp @@ -1,793 +1,816 @@ /*=================================================================== 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 __itkMLBSTrackingFilter_txx #define __itkMLBSTrackingFilter_txx #include #include #include #include "itkMLBSTrackingFilter.h" #include #include #include #include //#include #define _USE_MATH_DEFINES #include namespace itk { template< int NumImageFeatures > MLBSTrackingFilter< NumImageFeatures > ::MLBSTrackingFilter() : m_FiberPolyData(NULL) , m_Points(NULL) , m_Cells(NULL) , m_AngularThreshold(0.7) , m_StepSize(0) , m_MaxLength(10000) , m_MinTractLength(20.0) , m_MaxTractLength(400.0) , m_SeedsPerVoxel(1) - , m_UseDirection(true) , m_NumberOfSamples(50) , m_SamplingDistance(-1) , m_SeedImage(NULL) , m_MaskImage(NULL) , m_DecisionForest(NULL) , m_StoppingRegions(NULL) , m_DemoMode(false) , m_PauseTracking(false) , m_AbortTracking(false) , m_RemoveWmEndFibers(false) , m_AposterioriCurvCheck(false) , m_AvoidStop(true) + , m_RandomSampling(false) + , m_Verbose(false) { this->SetNumberOfRequiredInputs(1); } template< int NumImageFeatures > double MLBSTrackingFilter< NumImageFeatures > ::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< int NumImageFeatures > void MLBSTrackingFilter< NumImageFeatures >::BeforeThreadedGenerateData() { m_InputImage = const_cast(this->GetInput(0)); PreprocessRawData(); m_FiberPolyData = PolyDataType::New(); m_Points = vtkSmartPointer< vtkPoints >::New(); m_Cells = vtkSmartPointer< vtkCellArray >::New(); m_ImageSize.resize(3); m_ImageSize[0] = m_FeatureImage->GetLargestPossibleRegion().GetSize()[0]; m_ImageSize[1] = m_FeatureImage->GetLargestPossibleRegion().GetSize()[1]; m_ImageSize[2] = m_FeatureImage->GetLargestPossibleRegion().GetSize()[2]; m_ImageSpacing.resize(3); m_ImageSpacing[0] = m_FeatureImage->GetSpacing()[0]; m_ImageSpacing[1] = m_FeatureImage->GetSpacing()[1]; m_ImageSpacing[2] = m_FeatureImage->GetSpacing()[2]; double minSpacing; if(m_ImageSpacing[0]GetNumberOfThreads(); i++) { PolyDataType poly = PolyDataType::New(); m_PolyDataContainer.push_back(poly); } m_NotWmImage = ItkDoubleImgType::New(); m_NotWmImage->SetSpacing( m_FeatureImage->GetSpacing() ); m_NotWmImage->SetOrigin( m_FeatureImage->GetOrigin() ); m_NotWmImage->SetDirection( m_FeatureImage->GetDirection() ); m_NotWmImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() ); m_NotWmImage->Allocate(); m_NotWmImage->FillBuffer(0); m_WmImage = ItkDoubleImgType::New(); m_WmImage->SetSpacing( m_FeatureImage->GetSpacing() ); m_WmImage->SetOrigin( m_FeatureImage->GetOrigin() ); m_WmImage->SetDirection( m_FeatureImage->GetDirection() ); m_WmImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() ); m_WmImage->Allocate(); m_WmImage->FillBuffer(0); m_AvoidStopImage = ItkDoubleImgType::New(); m_AvoidStopImage->SetSpacing( m_FeatureImage->GetSpacing() ); m_AvoidStopImage->SetOrigin( m_FeatureImage->GetOrigin() ); m_AvoidStopImage->SetDirection( m_FeatureImage->GetDirection() ); m_AvoidStopImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() ); m_AvoidStopImage->Allocate(); m_AvoidStopImage->FillBuffer(0); if (m_StoppingRegions.IsNull()) { m_StoppingRegions = ItkUcharImgType::New(); m_StoppingRegions->SetSpacing( m_FeatureImage->GetSpacing() ); m_StoppingRegions->SetOrigin( m_FeatureImage->GetOrigin() ); m_StoppingRegions->SetDirection( m_FeatureImage->GetDirection() ); m_StoppingRegions->SetRegions( m_FeatureImage->GetLargestPossibleRegion() ); m_StoppingRegions->Allocate(); m_StoppingRegions->FillBuffer(0); } if (m_SeedImage.IsNull()) { m_SeedImage = ItkUcharImgType::New(); m_SeedImage->SetSpacing( m_FeatureImage->GetSpacing() ); m_SeedImage->SetOrigin( m_FeatureImage->GetOrigin() ); m_SeedImage->SetDirection( m_FeatureImage->GetDirection() ); m_SeedImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() ); m_SeedImage->Allocate(); m_SeedImage->FillBuffer(1); } if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( m_FeatureImage->GetSpacing() ); m_MaskImage->SetOrigin( m_FeatureImage->GetOrigin() ); m_MaskImage->SetDirection( m_FeatureImage->GetDirection() ); m_MaskImage->SetRegions( m_FeatureImage->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } else std::cout << "MLBSTrackingFilter: using mask image" << std::endl; if (m_AngularThreshold<0.0) m_AngularThreshold = 0.5*minSpacing; m_BuildFibersReady = 0; m_BuildFibersFinished = false; m_Threads = 0; m_Tractogram.clear(); + m_SamplingPointset = mitk::PointSet::New(); + m_AlternativePointset = mitk::PointSet::New(); std::cout << "MLBSTrackingFilter: Angular threshold: " << m_AngularThreshold << std::endl; std::cout << "MLBSTrackingFilter: Stepsize: " << m_StepSize << " mm" << std::endl; std::cout << "MLBSTrackingFilter: Seeds per voxel: " << m_SeedsPerVoxel << std::endl; std::cout << "MLBSTrackingFilter: Max. sampling distance: " << m_SamplingDistance << " mm" << std::endl; std::cout << "MLBSTrackingFilter: Number of samples: " << m_NumberOfSamples << std::endl; std::cout << "MLBSTrackingFilter: Max. tract length: " << m_MaxTractLength << " mm" << std::endl; std::cout << "MLBSTrackingFilter: Min. tract length: " << m_MinTractLength << " mm" << std::endl; std::cout << "MLBSTrackingFilter: Starting streamline tracking using " << this->GetNumberOfThreads() << " threads." << std::endl; } template< int NumImageFeatures > void MLBSTrackingFilter< NumImageFeatures >::PreprocessRawData() { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter InterpolationFilterType; std::cout << "MLBSTrackingFilter: Spherical signal interpolation and sampling ..." << std::endl; typename InterpolationFilterType::Pointer filter = InterpolationFilterType::New(); filter->SetGradientImage( m_GradientDirections, m_InputImage ); filter->SetBValue( m_B_Value ); filter->SetLambda(0.006); filter->SetNormalizationMethod(InterpolationFilterType::QBAR_RAW_SIGNAL); filter->Update(); // FeatureImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage(); // featureImageVector.push_back(itkFeatureImage); std::cout << "MLBSTrackingFilter: Creating feature image ..." << std::endl; vnl_vector_fixed ref; ref.fill(0); ref[0]=1; itk::OrientationDistributionFunction< double, NumImageFeatures*2 > odf; m_DirectionIndices.clear(); for (unsigned int f=0; f0) // only used directions on one hemisphere m_DirectionIndices.push_back(f); } m_FeatureImage = FeatureImageType::New(); m_FeatureImage->SetSpacing(filter->GetOutput()->GetSpacing()); m_FeatureImage->SetOrigin(filter->GetOutput()->GetOrigin()); m_FeatureImage->SetDirection(filter->GetOutput()->GetDirection()); m_FeatureImage->SetLargestPossibleRegion(filter->GetOutput()->GetLargestPossibleRegion()); m_FeatureImage->SetBufferedRegion(filter->GetOutput()->GetLargestPossibleRegion()); m_FeatureImage->SetRequestedRegion(filter->GetOutput()->GetLargestPossibleRegion()); m_FeatureImage->Allocate(); itk::ImageRegionIterator< typename InterpolationFilterType::OutputImageType > it(filter->GetOutput(), filter->GetOutput()->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { typename FeatureImageType::PixelType pix; for (unsigned int f=0; fSetPixel(it.GetIndex(), pix); ++it; } } template< int NumImageFeatures > void MLBSTrackingFilter< NumImageFeatures >::CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir) { // vnl_matrix_fixed< double, 3, 3 > rot = m_FeatureImage->GetDirection().GetTranspose(); // dir = rot*dir; dir *= m_StepSize; pos[0] += dir[0]; pos[1] += dir[1]; pos[2] += dir[2]; } template< int NumImageFeatures > bool MLBSTrackingFilter< NumImageFeatures > ::IsValidPosition(itk::Point &pos) { typename FeatureImageType::IndexType idx; m_FeatureImage->TransformPhysicalPointToIndex(pos, idx); if (!m_FeatureImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->GetPixel(idx)==0) return false; return true; } template< int NumImageFeatures > typename MLBSTrackingFilter< NumImageFeatures >::FeatureImageType::PixelType MLBSTrackingFilter< NumImageFeatures >::GetImageValues(itk::Point itkP) { itk::Index<3> idx; itk::ContinuousIndex< double, 3> cIdx; m_FeatureImage->TransformPhysicalPointToIndex(itkP, idx); m_FeatureImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx); typename FeatureImageType::PixelType pix; pix.Fill(0.0); if ( m_FeatureImage->GetLargestPossibleRegion().IsInside(idx) ) pix = m_FeatureImage->GetPixel(idx); else return pix; double frac_x = cIdx[0] - idx[0]; double frac_y = cIdx[1] - idx[1]; double frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < m_FeatureImage->GetLargestPossibleRegion().GetSize(0)-1 && idx[1] >= 0 && idx[1] < m_FeatureImage->GetLargestPossibleRegion().GetSize(1)-1 && idx[2] >= 0 && idx[2] < m_FeatureImage->GetLargestPossibleRegion().GetSize(2)-1) { vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = m_FeatureImage->GetPixel(idx) * interpWeights[0]; typename FeatureImageType::IndexType tmpIdx = idx; tmpIdx[0]++; pix += m_FeatureImage->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += m_FeatureImage->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += m_FeatureImage->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += m_FeatureImage->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += m_FeatureImage->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += m_FeatureImage->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += m_FeatureImage->GetPixel(tmpIdx) * interpWeights[7]; } return pix; } template< int NumImageFeatures > vnl_vector_fixed MLBSTrackingFilter< NumImageFeatures >::Classify(itk::Point& pos, int& candidates, vnl_vector_fixed& olddir, double angularThreshold, double& prob, bool avoidStop) { vnl_vector_fixed direction; direction.fill(0); - vigra::MultiArray<2, double> featureData; - if(m_UseDirection) - featureData = vigra::MultiArray<2, double>( vigra::Shape2(1,NumImageFeatures+3) ); - else - featureData = vigra::MultiArray<2, double>( vigra::Shape2(1,NumImageFeatures) ); + vigra::MultiArray<2, double> featureData = vigra::MultiArray<2, double>( vigra::Shape2(1,NumImageFeatures+3) ); typename FeatureImageType::PixelType featurePixel = GetImageValues(pos); // pixel values for (unsigned int f=0; f ref; ref.fill(0); ref[0]=1; for (unsigned int f=NumImageFeatures; f probs(vigra::Shape2(1, m_DecisionForest->class_count())); m_DecisionForest->predictProbabilities(featureData, probs); double outProb = 0; prob = 0; candidates = 0; // directions with probability > 0 for (int i=0; iclass_count(); i++) { if (probs(0,i)>0) { int classLabel = 0; m_DecisionForest->ext_param_.to_classlabel(i, classLabel); if (classLabel d = m_ODF.GetDirection(m_DirectionIndices.at(classLabel)); double dot = dot_product(d, olddir); if (olddir.magnitude()>0) { if (fabs(dot)>angularThreshold) { if (dot<0) d *= -1; dot = fabs(dot); direction += probs(0,i)*dot*d; prob += probs(0,i)*dot; } } else { direction += probs(0,i)*d; prob += probs(0,i); } } else outProb += probs(0,i); } } ItkDoubleImgType::IndexType idx; - m_NotWmImage->TransformPhysicalPointToIndex(pos, idx); - if (m_NotWmImage->GetLargestPossibleRegion().IsInside(idx)) + if (m_Verbose) { - m_NotWmImage->SetPixel(idx, m_NotWmImage->GetPixel(idx)+outProb); - m_WmImage->SetPixel(idx, m_WmImage->GetPixel(idx)+prob); + m_NotWmImage->TransformPhysicalPointToIndex(pos, idx); + if (m_NotWmImage->GetLargestPossibleRegion().IsInside(idx)) + { + m_NotWmImage->SetPixel(idx, m_NotWmImage->GetPixel(idx)+outProb); + m_WmImage->SetPixel(idx, m_WmImage->GetPixel(idx)+prob); + } } if (outProb>prob && prob>0) { candidates = 0; prob = 0; direction.fill(0.0); } - if (avoidStop && m_AvoidStopImage->GetLargestPossibleRegion().IsInside(idx) && candidates>0 && direction.magnitude()>0.001) + if (m_Verbose && avoidStop && m_AvoidStopImage->GetLargestPossibleRegion().IsInside(idx) && candidates>0 && direction.magnitude()>0.001) m_AvoidStopImage->SetPixel(idx, m_AvoidStopImage->GetPixel(idx)+0.1); return direction; } template< int NumImageFeatures > double MLBSTrackingFilter< NumImageFeatures >::GetRandDouble(double min, double max) { return (double)(rand()%((int)(10000*(max-min))) + 10000*min)/10000; } template< int NumImageFeatures > vnl_vector_fixed MLBSTrackingFilter< NumImageFeatures >::GetNewDirection(itk::Point &pos, vnl_vector_fixed& olddir) { + if (m_DemoMode) + { + m_SamplingPointset->Clear(); + m_AlternativePointset->Clear(); + } vnl_vector_fixed direction; direction.fill(0); ItkUcharImgType::IndexType idx; m_StoppingRegions->TransformPhysicalPointToIndex(pos, idx); if (m_StoppingRegions->GetPixel(idx)>0) return direction; if (olddir.magnitude()>0) olddir.normalize(); int candidates = 0; // number of directions with probability > 0 double prob = 0; direction = Classify(pos, candidates, olddir, m_AngularThreshold, prob); // sample neighborhood direction *= prob; itk::OrientationDistributionFunction< double, 50 > probeVecs; + itk::Point sample_pos; + int alternatives = 1; for (int i=0; i d; -// probe[0] = GetRandDouble()*m_SamplingDistance; -// probe[1] = GetRandDouble()*m_SamplingDistance; -// probe[2] = GetRandDouble()*m_SamplingDistance; - d = probeVecs.GetDirection(i)*m_SamplingDistance; + if (m_RandomSampling) + { + d[0] = GetRandDouble(); + d[1] = GetRandDouble(); + d[2] = GetRandDouble(); + d.normalize(); + d *= GetRandDouble(0,m_SamplingDistance); + } + else + { + d = probeVecs.GetDirection(i)*m_SamplingDistance; + } - itk::Point sample_pos; sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; + if(m_DemoMode) + m_SamplingPointset->InsertPoint(i, sample_pos); candidates = 0; vnl_vector_fixed tempDir = Classify(sample_pos, candidates, olddir, m_AngularThreshold, prob); // sample neighborhood if (candidates>0 && tempDir.magnitude()>0.001) { direction += tempDir*prob; } else if (m_AvoidStop && candidates==0 && olddir.magnitude()>0) // out of white matter { double dot = dot_product(d, olddir); if (dot >= 0.0) // in front of plane defined by pos and olddir d = -d + 2*dot*olddir; // reflect else d = -d; // invert // look a bit further into the other direction - sample_pos[0] = pos[0] + d[0]*2; - sample_pos[1] = pos[1] + d[1]*2; - sample_pos[2] = pos[2] + d[2]*2; + sample_pos[0] = pos[0] + d[0]; + sample_pos[1] = pos[1] + d[1]; + sample_pos[2] = pos[2] + d[2]; + if(m_DemoMode) + m_AlternativePointset->InsertPoint(alternatives, sample_pos); + alternatives++; candidates = 0; vnl_vector_fixed tempDir = Classify(sample_pos, candidates, olddir, m_AngularThreshold, prob, true); // sample neighborhood if (candidates>0 && tempDir.magnitude()>0.001) // are we back in the white matter? { direction += d; // go into the direction of the white matter direction += tempDir*prob; // go into the direction of the white matter direction at this location } } } if (direction.magnitude()>0.001) { direction.normalize(); olddir[0] = direction[0]; olddir[1] = direction[1]; olddir[2] = direction[2]; } else direction.fill(0); return direction; } template< int NumImageFeatures > double MLBSTrackingFilter< NumImageFeatures >::FollowStreamline(ThreadIdType threadId, itk::Point pos, vnl_vector_fixed dir, FiberType* fib, double tractLength, bool front) { vnl_vector_fixed dirOld = dir; dirOld = dir; for (int step=0; step< m_MaxLength/2; step++) { - while (m_PauseTracking){} - if (m_DemoMode) - { - m_Mutex.Lock(); - m_BuildFibersReady++; - m_Tractogram.push_back(*fib); - BuildFibers(true); - m_Stop = true; - m_Mutex.Unlock(); - while (m_Stop){} - } // get new position CalculateNewPosition(pos, dir); // is new position inside of image and mask if (!IsValidPosition(pos) || m_AbortTracking) // if not end streamline { return tractLength; } else // if yes, add new point to streamline { tractLength += m_StepSize; if (front) fib->push_front(pos); else fib->push_back(pos); if (m_AposterioriCurvCheck) { int curv = CheckCurvature(fib, front); // TODO: Move into classification ??? if (curv>0) { tractLength -= m_StepSize*curv; while (curv>0) { if (front) fib->pop_front(); else fib->pop_back(); curv--; } return tractLength; } } if (tractLength>m_MaxTractLength) return tractLength; } + if (m_DemoMode) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? + { + m_Mutex.Lock(); + m_BuildFibersReady++; + m_Tractogram.push_back(*fib); + BuildFibers(true); + m_Stop = true; + m_Mutex.Unlock(); + while (m_Stop){ + } + } dir = GetNewDirection(pos, dirOld); + while (m_PauseTracking){} + if (dir.magnitude()<0.0001) return tractLength; } return tractLength; } template< int NumImageFeatures > int MLBSTrackingFilter::CheckCurvature(FiberType* fib, bool front) { double m_Distance = 5; if (fib->size()<3) return 0; double dist = 0; std::vector< vnl_vector_fixed< float, 3 > > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0); double dev = 0; if (front) { int c=0; while(distsize()-1) { itk::Point p1 = fib->at(c); itk::Point p2 = fib->at(c+1); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==0) meanV += v; c++; } } else { int c=fib->size()-1; while(dist0) { itk::Point p1 = fib->at(c); itk::Point p2 = fib->at(c-1); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==fib->size()-1) meanV += v; c--; } } meanV.normalize(); for (int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/M_PI; } if (vectors.size()>0) dev /= vectors.size(); if (dev<30) return 0; else return vectors.size(); } template< int NumImageFeatures > void MLBSTrackingFilter< NumImageFeatures >::ThreadedGenerateData(const InputImageRegionType ®ionForThread, ThreadIdType threadId) { m_Mutex.Lock(); m_Threads++; m_Mutex.Unlock(); typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, regionForThread ); MaskIteratorType mit(m_MaskImage, regionForThread ); sit.GoToBegin(); mit.GoToBegin(); itk::Point worldPos; while( !sit.IsAtEnd() ) { if (sit.Value()==0 || mit.Value()==0) { ++sit; ++mit; continue; } for (int s=0; s start; unsigned int counter = 0; if (m_SeedsPerVoxel>1) { start[0] = index[0]+GetRandDouble(-0.5, 0.5); start[1] = index[1]+GetRandDouble(-0.5, 0.5); start[2] = index[2]+GetRandDouble(-0.5, 0.5); } else { start[0] = index[0]; start[1] = index[1]; start[2] = index[2]; } // get staring position m_SeedImage->TransformContinuousIndexToPhysicalPoint( start, worldPos ); // get starting direction int candidates = 0; double prob = 0; vnl_vector_fixed dirOld; dirOld.fill(0.0); vnl_vector_fixed dir = Classify(worldPos, candidates, dirOld, 0, prob); if (dir.magnitude()<0.0001) continue; // forward tracking tractLength = FollowStreamline(threadId, worldPos, dir, &fib, 0, false); fib.push_front(worldPos); if (m_RemoveWmEndFibers) { itk::Point check = fib.back(); dirOld.fill(0.0); vnl_vector_fixed check2 = GetNewDirection(check, dirOld); if (check2.magnitude()>0.001) { MITK_INFO << "Detected WM ending. Discarding fiber."; continue; } } // backward tracking tractLength = FollowStreamline(threadId, worldPos, -dir, &fib, tractLength, true); counter = fib.size(); if (m_RemoveWmEndFibers) { itk::Point check = fib.front(); dirOld.fill(0.0); vnl_vector_fixed check2 = GetNewDirection(check, dirOld); if (check2.magnitude()>0.001) { MITK_INFO << "Detected WM ending. Discarding fiber."; continue; } } if (tractLength void MLBSTrackingFilter< NumImageFeatures >::BuildFibers(bool check) { if (m_BuildFibersReady::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); for (int i=0; i container = vtkSmartPointer::New(); FiberType fib = m_Tractogram.at(i); - for (FiberType::iterator it = fib.begin(); it!=fib.end(); it++) + for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it) { vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } + if (check) for (int i=0; iSetPoints(vNewPoints); m_FiberPolyData->SetLines(vNewLines); m_BuildFibersFinished = true; } template< int NumImageFeatures > void MLBSTrackingFilter< NumImageFeatures >::AfterThreadedGenerateData() { MITK_INFO << "Generating polydata "; BuildFibers(false); MITK_INFO << "done"; } } #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.h index b3d6561252..ad08ddeb0a 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/itkMLBSTrackingFilter.h @@ -1,191 +1,195 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkMLBSTrackingFilter_h_ #define __itkMLBSTrackingFilter_h_ #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include // classification includes #include #include #include namespace itk{ /** * \brief Performes deterministic streamline tracking on the input tensor image. */ template< int NumImageFeatures=100 > class MLBSTrackingFilter : public ImageToImageFilter< VectorImage< short, 3 >, Image< double, 3 > > { public: typedef MLBSTrackingFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage< short, 3 >, Image< double, 3 > > Superclass; typedef vigra::RandomForest DecisionForestType; typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::InputImageRegionType InputImageRegionType; typedef Image< Vector< float, NumImageFeatures > , 3 > FeatureImageType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(MLBSTrackingFilter, ImageToImageFilter) typedef itk::Image ItkUcharImgType; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkFloatImgType; typedef vtkSmartPointer< vtkPolyData > PolyDataType; typedef std::deque< itk::Point > FiberType; typedef std::vector< FiberType > BundleType; - bool m_PauseTracking; + volatile bool m_PauseTracking; bool m_AbortTracking; bool m_BuildFibersFinished; int m_BuildFibersReady; - bool m_Stop; + volatile bool m_Stop; + mitk::PointSet::Pointer m_SamplingPointset; + mitk::PointSet::Pointer m_AlternativePointset; // void RequestFibers(){ m_Stop=true; m_BuildFibersReady=0; m_BuildFibersFinished=false; } itkGetMacro( FiberPolyData, PolyDataType ) ///< Output fibers itkSetMacro( SeedImage, ItkUcharImgType::Pointer) ///< Seeds are only placed inside of this mask. itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< Tracking is only performed inside of this mask image. itkSetMacro( SeedsPerVoxel, int) ///< One seed placed in the center of each voxel or multiple seeds randomly placed inside each voxel. itkSetMacro( StepSize, double) ///< Integration step size in mm itkSetMacro( MinTractLength, double ) ///< Shorter tracts are discarded. itkSetMacro( MaxTractLength, double ) itkSetMacro( AngularThreshold, double ) - itkSetMacro( UseDirection, bool ) itkSetMacro( SamplingDistance, double ) itkSetMacro( NumberOfSamples, int ) itkSetMacro( StoppingRegions, ItkUcharImgType::Pointer) itkSetMacro( B_Value, float ) itkSetMacro( GradientDirections, mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer ) itkSetMacro( DemoMode, bool ) itkSetMacro( RemoveWmEndFibers, bool ) itkSetMacro( AposterioriCurvCheck, bool ) itkSetMacro( AvoidStop, bool ) + itkSetMacro( RandomSampling, bool ) + itkSetMacro( Verbose, bool ) void SetDecisionForest( DecisionForestType* forest ) { m_DecisionForest = forest; } itkGetMacro( WmImage, ItkDoubleImgType::Pointer ) itkGetMacro( NotWmImage, ItkDoubleImgType::Pointer ) itkGetMacro( AvoidStopImage, ItkDoubleImgType::Pointer ) protected: MLBSTrackingFilter(); ~MLBSTrackingFilter() {} void CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir); ///< Calculate next integration step. double FollowStreamline(ThreadIdType threadId, itk::Point pos, vnl_vector_fixed dir, FiberType* fib, double tractLength, bool front); ///< Start streamline in one direction. bool IsValidPosition(itk::Point& pos); ///< Are we outside of the mask image? vnl_vector_fixed GetNewDirection(itk::Point& pos, vnl_vector_fixed& olddir); vnl_vector_fixed Classify(itk::Point& pos, int& candidates, vnl_vector_fixed& olddir, double angularThreshold, double& prob, bool avoidStop=false); typename FeatureImageType::PixelType GetImageValues(itk::Point itkP); double GetRandDouble(double min=-1, double max=1); double RoundToNearest(double num); void BeforeThreadedGenerateData() override; void PreprocessRawData(); void ThreadedGenerateData( const InputImageRegionType &outputRegionForThread, ThreadIdType threadId) override; void AfterThreadedGenerateData() override; PolyDataType m_FiberPolyData; vtkSmartPointer m_Points; vtkSmartPointer m_Cells; BundleType m_Tractogram; double m_AngularThreshold; double m_StepSize; int m_MaxLength; double m_MinTractLength; double m_MaxTractLength; int m_SeedsPerVoxel; - bool m_UseDirection; + bool m_RandomSampling; double m_SamplingDistance; int m_NumberOfSamples; std::vector< int > m_ImageSize; std::vector< double > m_ImageSpacing; SimpleFastMutexLock m_Mutex; ItkUcharImgType::Pointer m_StoppingRegions; ItkDoubleImgType::Pointer m_WmImage; ItkDoubleImgType::Pointer m_NotWmImage; ItkDoubleImgType::Pointer m_AvoidStopImage; ItkUcharImgType::Pointer m_SeedImage; ItkUcharImgType::Pointer m_MaskImage; typename FeatureImageType::Pointer m_FeatureImage; typename InputImageType::Pointer m_InputImage; mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer m_GradientDirections; float m_B_Value; bool m_AposterioriCurvCheck; bool m_RemoveWmEndFibers; bool m_AvoidStop; - + bool m_Verbose; int m_Threads; bool m_DemoMode; void BuildFibers(bool check); int CheckCurvature(FiberType* fib, bool front); // decision forest DecisionForestType* m_DecisionForest; itk::OrientationDistributionFunction< double, NumImageFeatures*2 > m_ODF; std::vector< int > m_DirectionIndices; std::vector< PolyDataType > m_PolyDataContainer; private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkMLBSTrackingFilter.cpp" #endif #endif //__itkMLBSTrackingFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.cpp index aa50256bdb..566b167309 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.cpp @@ -1,540 +1,477 @@ /*=================================================================== 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 _TrackingForestHandler_cpp #define _TrackingForestHandler_cpp #include "mitkTrackingForestHandler.h" #include #include namespace mitk { template< int NumberOfSignalFeatures > TrackingForestHandler< NumberOfSignalFeatures >::TrackingForestHandler() : m_GrayMatterSamplesPerVoxel(50) , m_StepSize(-1) - , m_UsePreviousDirection(true) , m_NumTrees(30) , m_MaxTreeDepth(50) , m_SampleFraction(1.0) { } template< int NumberOfSignalFeatures > TrackingForestHandler< NumberOfSignalFeatures >::~TrackingForestHandler() { } template< int NumberOfSignalFeatures > typename TrackingForestHandler< NumberOfSignalFeatures >::InterpolatedRawImageType::PixelType TrackingForestHandler< NumberOfSignalFeatures >::GetImageValues(itk::Point itkP, typename InterpolatedRawImageType::Pointer image) { itk::Index<3> idx; itk::ContinuousIndex< double, 3> cIdx; image->TransformPhysicalPointToIndex(itkP, idx); image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); typename InterpolatedRawImageType::PixelType pix; pix.Fill(0.0); if ( image->GetLargestPossibleRegion().IsInside(idx) ) pix = image->GetPixel(idx); else return pix; double frac_x = cIdx[0] - idx[0]; double frac_y = cIdx[1] - idx[1]; double frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < image->GetLargestPossibleRegion().GetSize(0)-1 && idx[1] >= 0 && idx[1] < image->GetLargestPossibleRegion().GetSize(1)-1 && idx[2] >= 0 && idx[2] < image->GetLargestPossibleRegion().GetSize(2)-1) { vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); pix = image->GetPixel(idx) * interpWeights[0]; typename InterpolatedRawImageType::IndexType tmpIdx = idx; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; pix += image->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; pix += image->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; pix += image->GetPixel(tmpIdx) * interpWeights[7]; } return pix; } template< int NumberOfSignalFeatures > void TrackingForestHandler< NumberOfSignalFeatures >::InputDataValidForTracking() { if (m_RawData.empty()) mitkThrow() << "No diffusion-weighted images set!"; } template< int NumberOfSignalFeatures > void TrackingForestHandler< NumberOfSignalFeatures >::StartTraining() { InputDataValidForTraining(); PreprocessInputData(); CalculateFeatures(); TrainForest(); } template< int NumberOfSignalFeatures > void TrackingForestHandler< NumberOfSignalFeatures >::InputDataValidForTraining() { if (m_RawData.empty()) mitkThrow() << "No diffusion-weighted images set!"; if (m_Tractograms.empty()) mitkThrow() << "No tractograms set!"; if (m_RawData.size()!=m_Tractograms.size()) mitkThrow() << "Unequal number of diffusion-weighted images and tractograms detected!"; } template< int NumberOfSignalFeatures > void TrackingForestHandler< NumberOfSignalFeatures >::PreprocessInputData() { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter InterpolationFilterType; MITK_INFO << "Spherical signal interpolation and sampling ..."; for (unsigned int i=0; iSetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(m_RawData.at(i)), mitk::DiffusionPropertyHelper::GetItkVectorImage(m_RawData.at(i)) ); qballfilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(m_RawData.at(i))); qballfilter->SetLambda(0.006); qballfilter->SetNormalizationMethod(InterpolationFilterType::QBAR_RAW_SIGNAL); qballfilter->Update(); // FeatureImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage(); // featureImageVector.push_back(itkFeatureImage); m_InterpolatedRawImages.push_back(qballfilter->GetOutput()); if (i>=m_MaskImages.size()) { ItkUcharImgType::Pointer newMask = ItkUcharImgType::New(); newMask->SetSpacing( m_InterpolatedRawImages.at(i)->GetSpacing() ); newMask->SetOrigin( m_InterpolatedRawImages.at(i)->GetOrigin() ); newMask->SetDirection( m_InterpolatedRawImages.at(i)->GetDirection() ); newMask->SetLargestPossibleRegion( m_InterpolatedRawImages.at(i)->GetLargestPossibleRegion() ); newMask->SetBufferedRegion( m_InterpolatedRawImages.at(i)->GetLargestPossibleRegion() ); newMask->SetRequestedRegion( m_InterpolatedRawImages.at(i)->GetLargestPossibleRegion() ); newMask->Allocate(); newMask->FillBuffer(1); m_MaskImages.push_back(newMask); } } MITK_INFO << "Resampling fibers and calculating number of samples ..."; m_NumberOfSamples = 0; for (unsigned int t=0; t::Pointer env = itk::TractDensityImageFilter< ItkUcharImgType >::New(); env->SetFiberBundle(m_Tractograms.at(t)); env->SetInputImage(mask); env->SetBinaryOutput(true); env->SetUseImageGeometry(true); env->Update(); wmmask = env->GetOutput(); m_WhiteMatterImages.push_back(wmmask); } itk::ImageRegionConstIterator it(wmmask, wmmask->GetLargestPossibleRegion()); int OUTOFWM = 0; while(!it.IsAtEnd()) { if (it.Get()==0 && mask->GetPixel(it.GetIndex())>0) OUTOFWM++; ++it; } m_NumberOfSamples += m_GrayMatterSamplesPerVoxel*OUTOFWM; MITK_INFO << "Samples outside of WM: " << m_NumberOfSamples; if (m_StepSize<0) { typename InterpolatedRawImageType::Pointer image = m_InterpolatedRawImages.at(t); float minSpacing = 1; if(image->GetSpacing()[0]GetSpacing()[1] && image->GetSpacing()[0]GetSpacing()[2]) minSpacing = image->GetSpacing()[0]; else if (image->GetSpacing()[1] < image->GetSpacing()[2]) minSpacing = image->GetSpacing()[1]; else minSpacing = image->GetSpacing()[2]; m_StepSize = minSpacing*0.5; } m_Tractograms.at(t)->ResampleSpline(m_StepSize); m_NumberOfSamples += m_Tractograms.at(t)->GetNumberOfPoints(); m_NumberOfSamples -= 2*m_Tractograms.at(t)->GetNumFibers(); } MITK_INFO << "Number of samples: " << m_NumberOfSamples; } template< int NumberOfSignalFeatures > void TrackingForestHandler< NumberOfSignalFeatures >::CalculateFeatures() { vnl_vector_fixed ref; ref.fill(0); ref[0]=1; itk::OrientationDistributionFunction< double, 2*NumberOfSignalFeatures > directions; std::vector< int > directionIndices; for (unsigned int f=0; f<2*NumberOfSignalFeatures; f++) { if (dot_product(ref, directions.GetDirection(f))>0) directionIndices.push_back(f); } - int numDirectionFeatures = 0; - if (m_UsePreviousDirection) - numDirectionFeatures = 3; + int numDirectionFeatures = 3; m_FeatureData.reshape( vigra::Shape2(m_NumberOfSamples, NumberOfSignalFeatures+numDirectionFeatures) ); m_LabelData.reshape( vigra::Shape2(m_NumberOfSamples,1) ); MITK_INFO << "Number of features: " << m_FeatureData.shape(1); itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); MITK_INFO << "Creating training data ..."; int sampleCounter = 0; for (unsigned int t=0; t it(wmMask, wmMask->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (it.Get()==0 && (mask.IsNull() || (mask.IsNotNull() && mask->GetPixel(it.GetIndex())>0))) { typename InterpolatedRawImageType::PixelType pix = image->GetPixel(it.GetIndex()); - if (m_UsePreviousDirection) - { + // null direction for (unsigned int f=0; f probe; probe[0] = m_RandGen->GetVariate()*2-1; probe[1] = m_RandGen->GetVariate()*2-1; probe[2] = m_RandGen->GetVariate()*2-1; probe.normalize(); if (dot_product(ref, probe)<0) probe *= -1; for (unsigned int f=NumberOfSignalFeatures; f idx; - idx[0] = it.GetIndex()[0]; - idx[1] = it.GetIndex()[1]; - idx[2] = it.GetIndex()[2]; - itk::Point itkP1; - image->TransformContinuousIndexToPhysicalPoint(idx, itkP1); - typename InterpolatedRawImageType::PixelType pix = GetImageValues(itkP1, image);; - for (unsigned int f=0; f idx; - idx[0] = it.GetIndex()[0] + m_RandGen->GetVariate()-0.5; - idx[1] = it.GetIndex()[1] + m_RandGen->GetVariate()-0.5; - idx[2] = it.GetIndex()[2] + m_RandGen->GetVariate()-0.5; - itk::Point itkP1; - image->TransformContinuousIndexToPhysicalPoint(idx, itkP1); - typename InterpolatedRawImageType::PixelType pix = GetImageValues(itkP1, image);; - for (unsigned int f=0; f polyData = fib->GetFiberPolyData(); for (int i=0; iGetNumFibers(); i++) { vtkCell* cell = polyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vnl_vector_fixed dirOld; dirOld.fill(0.0); for (int j=0; jGetPoint(j); itk::Point itkP1; itkP1[0] = p1[0]; itkP1[1] = p1[1]; itkP1[2] = p1[2]; vnl_vector_fixed dir; dir.fill(0.0); itk::Point itkP2; double* p2 = points->GetPoint(j+1); itkP2[0] = p2[0]; itkP2[1] = p2[1]; itkP2[2] = p2[2]; dir[0]=itkP2[0]-itkP1[0]; dir[1]=itkP2[1]-itkP1[1]; dir[2]=itkP2[2]-itkP1[2]; if (dir.magnitude()<0.0001) { MITK_INFO << "streamline error!"; continue; } dir.normalize(); if (dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2]) { MITK_INFO << "ERROR: NaN direction!"; continue; } if (j==0) { dirOld = dir; continue; } // get voxel values typename InterpolatedRawImageType::PixelType pix = GetImageValues(itkP1, image); for (unsigned int f=0; f0.0001) { int label = 0; for (unsigned int f=0; fangle) { m_LabelData(sampleCounter,0) = f; angle = a; label = f; } } } dirOld = dir; sampleCounter++; } } } } template< int NumberOfSignalFeatures > void TrackingForestHandler< NumberOfSignalFeatures >::TrainForest() { + MITK_INFO << "Maximum tree depths: " << m_MaxTreeDepth; MITK_INFO << "Sample fraction per tree: " << m_SampleFraction; MITK_INFO << "Number of trees: " << m_NumTrees; - bool random_split = false; - vigra::rf::visitors::OOB_Error oob_v; - MITK_INFO << "Create Split Function"; - // typedef ThresholdSplit, vigra::ClassificationTag> DefaultSplitType; - - m_Forest.set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal - m_Forest.set_options().sample_with_replacement(true); // if sampled with replacement or not - m_Forest.set_options().samples_per_tree(m_SampleFraction); // Fraction of samples that are used to train a tree - m_Forest.set_options().tree_count(1); // Number of trees that are calculated; - m_Forest.set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node - // rf.set_options().features_per_node(10); - - m_Forest.learn(m_FeatureData, m_LabelData, vigra::rf::visitors::create_visitor(oob_v)); - - // Prepare parallel VariableImportance Calculation - int numMod = m_FeatureData.shape(1); - const int numClass = 2 + 2; - - float** varImp = new float*[numMod]; - - for(int i = 0; i < numMod; i++) - varImp[i] = new float[numClass]; - - for (int i = 0; i < numMod; ++i) - for (int j = 0; j < numClass; ++j) - varImp[i][j] = 0.0; + std::vector< vigra::RandomForest > trees; + int count = 0; #pragma omp parallel for - for (int i = 0; i < m_NumTrees - 1; ++i) + for (int i = 0; i < m_NumTrees; ++i) { vigra::RandomForest lrf; vigra::rf::visitors::OOB_Error loob_v; lrf.set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal lrf.set_options().sample_with_replacement(true); // if sampled with replacement or not lrf.set_options().samples_per_tree(m_SampleFraction); // Fraction of samples that are used to train a tree lrf.set_options().tree_count(1); // Number of trees that are calculated; lrf.set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node - // lrf.set_options().features_per_node(10); - - vigra::rf::visitors::VariableImportanceVisitor lvariableImportance; - lrf.learn(m_FeatureData, m_LabelData, vigra::rf::visitors::create_visitor(loob_v)); + lrf.ext_param_.max_tree_depth = m_MaxTreeDepth; + lrf.learn(m_FeatureData, m_LabelData); #pragma omp critical { - m_Forest.trees_.push_back(lrf.trees_[0]); + count++; + MITK_INFO << "Tree " << count << " finished training."; + trees.push_back(lrf); } } + for (int i = 1; i < m_NumTrees; ++i) + trees.at(0).trees_.push_back(trees.at(i).trees_[0]); + + m_Forest = trees.at(0); m_Forest.options_.tree_count_ = m_NumTrees; MITK_INFO << "Training finsihed"; - MITK_INFO << "The out-of-bag error is: " << oob_v.oob_breiman << std::endl; } template< int NumberOfSignalFeatures > void TrackingForestHandler< NumberOfSignalFeatures >::SaveForest(std::string forestFile) { MITK_INFO << "Saving forest to " << forestFile; vigra::rf_export_HDF5( m_Forest, forestFile, "" ); } template< int NumberOfSignalFeatures > void TrackingForestHandler< NumberOfSignalFeatures >::LoadForest(std::string forestFile) { MITK_INFO << "Loading forest from " << forestFile; vigra::rf_import_HDF5(m_Forest, forestFile); } //// superclass implementations //template< int NumberOfSignalFeatures > //void TrackingForestHandler< NumberOfSignalFeatures >::UpdateOutputInformation() //{ //} //template< int NumberOfSignalFeatures > //void TrackingForestHandler< NumberOfSignalFeatures >::SetRequestedRegionToLargestPossibleRegion() //{ //} //template< int NumberOfSignalFeatures > //bool TrackingForestHandler< NumberOfSignalFeatures >::RequestedRegionIsOutsideOfTheBufferedRegion() //{ // return false; //} //template< int NumberOfSignalFeatures > //bool TrackingForestHandler< NumberOfSignalFeatures >::VerifyRequestedRegion() //{ // return true; //} //template< int NumberOfSignalFeatures > //void TrackingForestHandler< NumberOfSignalFeatures >::SetRequestedRegion(const itk::DataObject* ) //{ //} } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.h index 4aecb1ac41..5c34f6a0a8 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/MLTracking/mitkTrackingForestHandler.h @@ -1,115 +1,113 @@ /*=================================================================== 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 _TrackingForestHandler #define _TrackingForestHandler #include "mitkBaseData.h" #include #include #include #include // classification includes //#include "RegressionForestClasses.hxx" #undef DIFFERENCE #define VIGRA_STATIC_LIB #include #include #include #define _USE_MATH_DEFINES #include namespace mitk { /** * \brief */ template< int NumberOfSignalFeatures=100 > class TrackingForestHandler { public: TrackingForestHandler(); ~TrackingForestHandler(); typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::Vector< float, NumberOfSignalFeatures*2 > , 3 > InterpolatedRawImageType; void SetRawData( std::vector< Image::Pointer > images ){ m_RawData = images; } void SetTractograms( std::vector< FiberBundle::Pointer > tractograms ) { m_Tractograms.clear(); for (int i=0; iGetDeepCopy()); } } void SetMaskImages( std::vector< ItkUcharImgType::Pointer > images ){ m_MaskImages = images; } void SetWhiteMatterImages( std::vector< ItkUcharImgType::Pointer > images ){ m_WhiteMatterImages = images; } void StartTraining(); void SaveForest(std::string forestFile); void LoadForest(std::string forestFile); void SetNumTrees(int num){ m_NumTrees = num; } void SetMaxTreeDepth(int depth){ m_MaxTreeDepth = depth; } - void SetUsePreviousDirection(bool use){ m_UsePreviousDirection = use; } void SetStepSize(double step){ m_StepSize = step; } void SetGrayMatterSamplesPerVoxel(int samples){ m_GrayMatterSamplesPerVoxel = samples; } void SetSampleFraction(double fraction){ m_SampleFraction = fraction; } vigra::RandomForest GetForest(){ return m_Forest; } protected: void InputDataValidForTracking(); void InputDataValidForTraining(); void PreprocessInputData(); void CalculateFeatures(); void TrainForest(); int m_GrayMatterSamplesPerVoxel; double m_StepSize; - bool m_UsePreviousDirection; int m_NumTrees; int m_MaxTreeDepth; double m_SampleFraction; std::vector< Image::Pointer > m_RawData; std::vector< FiberBundle::Pointer > m_Tractograms; std::vector< ItkUcharImgType::Pointer > m_MaskImages; std::vector< ItkUcharImgType::Pointer > m_WhiteMatterImages; std::vector< ItkUcharImgType::Pointer > m_SeedImages; std::vector< ItkUcharImgType::Pointer > m_StopImages; int m_NumberOfSamples; vigra::RandomForest m_Forest; vigra::MultiArray<2, double> m_FeatureData; vigra::MultiArray<2, double> m_LabelData; std::vector< typename InterpolatedRawImageType::Pointer > m_InterpolatedRawImages; typename InterpolatedRawImageType::PixelType GetImageValues(itk::Point itkP, typename InterpolatedRawImageType::Pointer image); }; } #include "mitkTrackingForestHandler.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp index 56d48bbb46..a9d5174fe4 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp @@ -1,363 +1,372 @@ /*=================================================================== 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 __itkEvaluateDirectionImagesFilter_cpp #define __itkEvaluateDirectionImagesFilter_cpp #include "itkEvaluateDirectionImagesFilter.h" #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { template< class PixelType > EvaluateDirectionImagesFilter< PixelType > ::EvaluateDirectionImagesFilter(): m_ImageSet(NULL), m_ReferenceImageSet(NULL), m_IgnoreMissingDirections(false), + m_IgnoreEmptyVoxels(false), m_Eps(0.0001) { this->SetNumberOfIndexedOutputs(2); } template< class PixelType > void EvaluateDirectionImagesFilter< PixelType >::GenerateData() { if (m_ImageSet.IsNull() || m_ReferenceImageSet.IsNull()) return; DirectionImageContainerType::Pointer set1 = DirectionImageContainerType::New(); DirectionImageContainerType::Pointer set2 = DirectionImageContainerType::New(); for (unsigned int i=0; iSize(); i++) { typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New(); duplicator->SetInputImage( m_ImageSet->GetElement(i) ); duplicator->Update(); set1->InsertElement(i, dynamic_cast(duplicator->GetOutput())); } for (unsigned int i=0; iSize(); i++) { typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New(); duplicator->SetInputImage( m_ReferenceImageSet->GetElement(i) ); duplicator->Update(); set2->InsertElement(i, dynamic_cast(duplicator->GetOutput())); } m_ImageSet = set1; m_ReferenceImageSet = set2; // angular error image typename OutputImageType::Pointer outputImage = OutputImageType::New(); outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); outputImage->Allocate(); outputImage->FillBuffer(0.0); this->SetNthOutput(0, outputImage); // length error image outputImage = OutputImageType::New(); outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); outputImage->Allocate(); outputImage->FillBuffer(0.0); this->SetNthOutput(1, outputImage); if (m_MaskImage.IsNull()) { m_MaskImage = UCharImageType::New(); m_MaskImage->SetOrigin( outputImage->GetOrigin() ); m_MaskImage->SetRegions( outputImage->GetLargestPossibleRegion() ); m_MaskImage->SetSpacing( outputImage->GetSpacing() ); m_MaskImage->SetDirection( outputImage->GetDirection() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } m_MeanAngularError = 0.0; m_MedianAngularError = 0; m_MaxAngularError = 0.0; m_MinAngularError = itk::NumericTraits::max(); m_VarAngularError = 0.0; m_AngularErrorVector.clear(); m_MeanLengthError = 0.0; m_MedianLengthError = 0; m_MaxLengthError = 0.0; m_MinLengthError = itk::NumericTraits::max(); m_VarLengthError = 0.0; m_LengthErrorVector.clear(); if (m_ImageSet.IsNull() || m_ReferenceImageSet.IsNull()) return; outputImage = static_cast< OutputImageType* >(this->ProcessObject::GetOutput(0)); typename OutputImageType::Pointer outputImage2 = static_cast< OutputImageType* >(this->ProcessObject::GetOutput(1)); ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion()); ImageRegionIterator< OutputImageType > oit2(outputImage2, outputImage2->GetLargestPossibleRegion()); ImageRegionIterator< UCharImageType > mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion()); int numTestImages = m_ImageSet->Size(); int numReferenceImages = m_ReferenceImageSet->Size(); int maxNumDirections = std::max(numReferenceImages, numTestImages); // matrix containing the angular error between the directions vnl_matrix< float > diffM; diffM.set_size(maxNumDirections, maxNumDirections); boost::progress_display disp(outputImage->GetLargestPossibleRegion().GetSize()[0]*outputImage->GetLargestPossibleRegion().GetSize()[1]*outputImage->GetLargestPossibleRegion().GetSize()[2]); while( !oit.IsAtEnd() ) { ++disp; if( mit.Get()!=1 ) { ++oit; ++oit2; ++mit; continue; } typename OutputImageType::IndexType index = oit.GetIndex(); float maxAngularError = 1.0; diffM.fill(10); // initialize with invalid error value // get number of valid directions (length > 0) int numRefDirs = 0; int numTestDirs = 0; std::vector< vnl_vector_fixed< PixelType, 3 > > testDirs; std::vector< vnl_vector_fixed< PixelType, 3 > > refDirs; for (int i=0; i refDir = m_ReferenceImageSet->GetElement(i)->GetPixel(index).GetVnlVector(); if (refDir.magnitude() > m_Eps ) { refDir.normalize(); refDirs.push_back(refDir); numRefDirs++; } } for (int i=0; i testDir = m_ImageSet->GetElement(i)->GetPixel(index).GetVnlVector(); if (testDir.magnitude() > m_Eps ) { testDir.normalize(); testDirs.push_back(testDir); numTestDirs++; } } + if (m_IgnoreEmptyVoxels && (numRefDirs==0 || numTestDirs==0) ) + { + ++oit; + ++oit2; + ++mit; + continue; + } + // i: index of reference direction // j: index of test direction for (int i=0; i refDir; if (i testDir; if (j1.0) diffM[i][j] = 1.0; } } float angularError = 0.0; float lengthError = 0.0; int counter = 0; vnl_matrix< float > diffM_copy = diffM; for (int k=0; k small error) for (int i=0; ierror && diffM[i][j]<2) // found valid error entry { error = diffM[i][j]; a = i; b = j; missingDir = false; } else if (diffM[i][j]<0 && error<0) // found missing direction { a = i; b = j; missingDir = true; } } if (a<0 || b<0 || (m_IgnoreMissingDirections && missingDir)) continue; // no more directions found if (a>=numRefDirs && b>=numTestDirs) { MITK_INFO << "ERROR: missing test and reference direction. should not be possible. check code."; continue; } // remove processed directions from error matrix diffM.set_row(a, 10.0); diffM.set_column(b, 10.0); if (a>=numRefDirs) // missing reference direction (find next closest) { for (int i=0; ierror) { error = diffM_copy[i][b]; a = i; } } else if (b>=numTestDirs) // missing test direction (find next closest) { for (int i=0; ierror) { error = diffM_copy[a][i]; b = i; } } float refLength = 0; float testLength = 1; if (a>=numRefDirs || b>=numTestDirs || error<0) error = 0; else { refLength = m_ReferenceImageSet->GetElement(a)->GetPixel(index).GetVnlVector().magnitude(); testLength = m_ImageSet->GetElement(b)->GetPixel(index).GetVnlVector().magnitude(); } m_LengthErrorVector.push_back( fabs(refLength-testLength) ); m_AngularErrorVector.push_back( acos(error)*180.0/M_PI ); m_MeanAngularError += m_AngularErrorVector.back(); m_MeanLengthError += m_LengthErrorVector.back(); angularError += m_AngularErrorVector.back(); lengthError += m_LengthErrorVector.back(); counter++; } if (counter>0) { lengthError /= counter; angularError /= counter; } oit2.Set(lengthError); oit.Set(angularError); ++oit; ++oit2; ++mit; } std::sort( m_AngularErrorVector.begin(), m_AngularErrorVector.end() ); m_MeanAngularError /= m_AngularErrorVector.size(); // mean for (unsigned int i=0; im_MaxAngularError ) m_MaxAngularError = m_AngularErrorVector.at(i); if ( m_AngularErrorVector.at(i)1) { m_VarAngularError /= (m_AngularErrorVector.size()-1); // variance // median if (m_AngularErrorVector.size()%2 == 0) m_MedianAngularError = 0.5*( m_AngularErrorVector.at( m_AngularErrorVector.size()/2 ) + m_AngularErrorVector.at( m_AngularErrorVector.size()/2+1 ) ); else m_MedianAngularError = m_AngularErrorVector.at( (m_AngularErrorVector.size()+1)/2 ) ; } std::sort( m_LengthErrorVector.begin(), m_LengthErrorVector.end() ); m_MeanLengthError /= m_LengthErrorVector.size(); // mean for (unsigned int i=0; im_MaxLengthError ) m_MaxLengthError = m_LengthErrorVector.at(i); if ( m_LengthErrorVector.at(i)1) { m_VarLengthError /= (m_LengthErrorVector.size()-1); // variance // median if (m_LengthErrorVector.size()%2 == 0) m_MedianLengthError = 0.5*( m_LengthErrorVector.at( m_LengthErrorVector.size()/2 ) + m_LengthErrorVector.at( m_LengthErrorVector.size()/2+1 ) ); else m_MedianLengthError = m_LengthErrorVector.at( (m_LengthErrorVector.size()+1)/2 ) ; } } } #endif // __itkEvaluateDirectionImagesFilter_cpp diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h index 320d98de29..17a08617ca 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h @@ -1,112 +1,114 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkEvaluateDirectionImagesFilter_h_ #define __itkEvaluateDirectionImagesFilter_h_ #include #include #include namespace itk{ /** \brief Evaluates the voxel-wise angular error between two sets of directions. */ template< class PixelType > class EvaluateDirectionImagesFilter : public ImageSource< Image< PixelType, 3 > > { public: typedef EvaluateDirectionImagesFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageSource< Image< PixelType, 3 > > Superclass; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef typename Superclass::OutputImageType OutputImageType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(EvaluateDirectionImagesFilter, ImageToImageFilter) typedef Vector< float, 3 > DirectionType; typedef Image< DirectionType, 3 > DirectionImageType; typedef VectorContainer< unsigned int, DirectionImageType::Pointer > DirectionImageContainerType; typedef Image< float, 3 > FloatImageType; typedef Image< bool, 3 > BoolImageType; typedef Image< unsigned char, 3 > UCharImageType; itkSetMacro( ImageSet , DirectionImageContainerType::Pointer) ///< test image containers itkSetMacro( ReferenceImageSet , DirectionImageContainerType::Pointer) ///< reference image containers itkSetMacro( MaskImage , UCharImageType::Pointer) ///< Calculation is only performed inside of the mask image. itkSetMacro( IgnoreMissingDirections , bool) ///< If in one voxel, the number of directions differs between the test container and the reference, the excess directions are ignored. Otherwise, the error to the next closest direction is calculated. + itkSetMacro( IgnoreEmptyVoxels , bool) ///< Don't increase error if either reference or test voxel is empty. /** Output statistics of the measured angular errors. */ itkGetMacro( MeanAngularError, float) itkGetMacro( MinAngularError, float) itkGetMacro( MaxAngularError, float) itkGetMacro( VarAngularError, float) itkGetMacro( MedianAngularError, float) /** Output statistics of the measured peak length errors. */ itkGetMacro( MeanLengthError, float) itkGetMacro( MinLengthError, float) itkGetMacro( MaxLengthError, float) itkGetMacro( VarLengthError, float) itkGetMacro( MedianLengthError, float) protected: EvaluateDirectionImagesFilter(); ~EvaluateDirectionImagesFilter() {} void GenerateData(); UCharImageType::Pointer m_MaskImage; DirectionImageContainerType::Pointer m_ImageSet; DirectionImageContainerType::Pointer m_ReferenceImageSet; bool m_IgnoreMissingDirections; + bool m_IgnoreEmptyVoxels; double m_MeanAngularError; double m_MedianAngularError; double m_MaxAngularError; double m_MinAngularError; double m_VarAngularError; std::vector< double > m_AngularErrorVector; double m_MeanLengthError; double m_MedianLengthError; double m_MaxLengthError; double m_MinLengthError; double m_VarLengthError; std::vector< double > m_LengthErrorVector; double m_Eps; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkEvaluateDirectionImagesFilter.cpp" #endif #endif //__itkEvaluateDirectionImagesFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp index 4d434011cb..d844388d70 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp @@ -1,565 +1,575 @@ #include "itkTractsToVectorImageFilter.h" // VTK #include #include #include // ITK #include #include // misc #define _USE_MATH_DEFINES #include #include namespace itk{ static bool CompareVectorLengths(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2) { return (v1.magnitude()>v2.magnitude()); } template< class PixelType > TractsToVectorImageFilter< PixelType >::TractsToVectorImageFilter(): m_AngularThreshold(0.7), m_Epsilon(0.999), m_MaskImage(NULL), m_NormalizeVectors(false), m_UseWorkingCopy(true), m_MaxNumDirections(3), m_SizeThreshold(0.3), m_NumDirectionsImage(NULL), m_CreateDirectionImages(true) { this->SetNumberOfRequiredOutputs(1); } template< class PixelType > TractsToVectorImageFilter< PixelType >::~TractsToVectorImageFilter() { } template< class PixelType > vnl_vector_fixed TractsToVectorImageFilter< PixelType >::GetVnlVector(double point[]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } template< class PixelType > itk::Point TractsToVectorImageFilter< PixelType >::GetItkPoint(double point[]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class PixelType > void TractsToVectorImageFilter< PixelType >::GenerateData() { mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry(); // calculate new image parameters itk::Vector spacing; itk::Point origin; itk::Matrix direction; ImageRegion<3> imageRegion; if (!m_MaskImage.IsNull()) { spacing = m_MaskImage->GetSpacing(); imageRegion = m_MaskImage->GetLargestPossibleRegion(); origin = m_MaskImage->GetOrigin(); direction = m_MaskImage->GetDirection(); } else { spacing = geometry->GetSpacing(); origin = geometry->GetOrigin(); mitk::BaseGeometry::BoundsArrayType bounds = geometry->GetBounds(); origin[0] += bounds.GetElement(0); origin[1] += bounds.GetElement(2); origin[2] += bounds.GetElement(4); for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]; imageRegion.SetSize(0, geometry->GetExtent(0)); imageRegion.SetSize(1, geometry->GetExtent(1)); imageRegion.SetSize(2, geometry->GetExtent(2)); m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( spacing ); m_MaskImage->SetOrigin( origin ); m_MaskImage->SetDirection( direction ); m_MaskImage->SetRegions( imageRegion ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } OutputImageType::RegionType::SizeType outImageSize = imageRegion.GetSize(); m_OutImageSpacing = m_MaskImage->GetSpacing(); m_ClusteredDirectionsContainer = ContainerType::New(); // initialize num directions image m_NumDirectionsImage = ItkUcharImgType::New(); m_NumDirectionsImage->SetSpacing( spacing ); m_NumDirectionsImage->SetOrigin( origin ); m_NumDirectionsImage->SetDirection( direction ); m_NumDirectionsImage->SetRegions( imageRegion ); m_NumDirectionsImage->Allocate(); m_NumDirectionsImage->FillBuffer(0); // initialize direction images m_DirectionImageContainer = DirectionImageContainerType::New(); // resample fiber bundle double minSpacing = 1; if(m_OutImageSpacing[0]GetDeepCopy(); // resample fiber bundle for sufficient voxel coverage m_FiberBundle->ResampleSpline(minSpacing/10); // iterate over all fibers vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); int numFibers = m_FiberBundle->GetNumFibers(); m_DirectionsContainer = ContainerType::New(); + VectorContainer< unsigned int, std::vector< double > >::Pointer peakLengths = VectorContainer< unsigned int, std::vector< double > >::New(); + MITK_INFO << "Generating directions from tractogram"; boost::progress_display disp(numFibers); for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) continue; vnl_vector_fixed dir; itk::Point worldPos; vnl_vector v; + + float fiberWeight = m_FiberBundle->GetFiberWeight(i); + for( int j=0; jGetPoint(j); worldPos = GetItkPoint(temp); itk::Index<3> index; m_MaskImage->TransformPhysicalPointToIndex(worldPos, index); if (!m_MaskImage->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0) continue; // get fiber tangent direction at this position v = GetVnlVector(temp); dir = GetVnlVector(points->GetPoint(j+1))-v; if (dir.is_zero()) continue; dir.normalize(); // add direction to container unsigned int idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]); DirectionContainerType::Pointer dirCont; if (m_DirectionsContainer->IndexExists(idx)) { + peakLengths->ElementAt(idx).push_back(fiberWeight); + dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->push_back(dir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->push_back(dir); } else { dirCont = DirectionContainerType::New(); dirCont->push_back(dir); m_DirectionsContainer->InsertElement(idx, dirCont); + + std::vector< double > lengths; lengths.push_back(fiberWeight); + peakLengths->InsertElement(idx, lengths); } } } vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); itk::ImageRegionIterator dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion()); MITK_INFO << "Clustering directions"; boost::progress_display disp2(outImageSize[0]*outImageSize[1]*outImageSize[2]); while(!dirIt.IsAtEnd()) { ++disp2; OutputImageType::IndexType index = dirIt.GetIndex(); int idx = index[0]+(index[1]+index[2]*outImageSize[1])*outImageSize[0]; if (!m_DirectionsContainer->IndexExists(idx)) { ++dirIt; continue; } DirectionContainerType::Pointer dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull() || dirCont->empty()) { ++dirIt; continue; } - std::vector< double > lengths; lengths.resize(dirCont->size(), 1); // all peaks have size 1 +// std::vector< double > lengths; lengths.resize(dirCont->size(), 1); // all peaks have size 1 DirectionContainerType::Pointer directions; if (m_MaxNumDirections>0) { - directions = FastClustering(dirCont, lengths); + directions = FastClustering(dirCont, peakLengths->GetElement(idx)); std::sort( directions->begin(), directions->end(), CompareVectorLengths ); } else directions = dirCont; unsigned int numDir = directions->size(); if (m_MaxNumDirections>0 && numDir>m_MaxNumDirections) numDir = m_MaxNumDirections; int count = 0; for (unsigned int i=0; i container = vtkSmartPointer::New(); itk::ContinuousIndex center; center[0] = index[0]; center[1] = index[1]; center[2] = index[2]; itk::Point worldCenter; m_MaskImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); DirectionType dir = directions->at(i); if (dir.magnitude()size()) { ItkDirectionImageType::Pointer directionImage = ItkDirectionImageType::New(); directionImage->SetSpacing( spacing ); directionImage->SetOrigin( origin ); directionImage->SetDirection( direction ); directionImage->SetRegions( imageRegion ); directionImage->Allocate(); Vector< float, 3 > nullVec; nullVec.Fill(0.0); directionImage->FillBuffer(nullVec); m_DirectionImageContainer->InsertElement(i, directionImage); } // set direction image pixel ItkDirectionImageType::Pointer directionImage = m_DirectionImageContainer->GetElement(i); Vector< float, 3 > pixel; pixel.SetElement(0, dir[0]); pixel.SetElement(1, dir[1]); pixel.SetElement(2, dir[2]); directionImage->SetPixel(index, pixel); } // add direction to vector field (with spacing compensation) itk::Point worldStart; worldStart[0] = worldCenter[0]-dir[0]/2*minSpacing; worldStart[1] = worldCenter[1]-dir[1]/2*minSpacing; worldStart[2] = worldCenter[2]-dir[2]/2*minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+dir[0]/2*minSpacing; worldEnd[1] = worldCenter[1]+dir[1]/2*minSpacing; worldEnd[2] = worldCenter[2]+dir[2]/2*minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); } dirIt.Set(count); ++dirIt; } vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); directionsPolyData->SetPoints(m_VtkPoints); directionsPolyData->SetLines(m_VtkCellArray); m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData); } template< class PixelType > TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths) { DirectionContainerType::Pointer outDirs = DirectionContainerType::New(); if (inDirs->size()<2) return inDirs; DirectionType oldMean, currentMean; std::vector< int > touched; // initialize touched.resize(inDirs->size(), 0); bool free = true; currentMean = inDirs->at(0); // initialize first seed currentMean.normalize(); double length = lengths.at(0); touched[0] = 1; std::vector< double > newLengths; bool meanChanged = false; double max = 0; while (free) { oldMean.fill(0.0); // start mean-shift clustering double angle = 0; while (fabs(dot_product(currentMean, oldMean))<0.99) { oldMean = currentMean; currentMean.fill(0.0); for (unsigned int i=0; isize(); i++) { angle = dot_product(oldMean, inDirs->at(i)); if (angle>=m_AngularThreshold) { currentMean += inDirs->at(i); if (meanChanged) length += lengths.at(i); touched[i] = 1; meanChanged = true; } else if (-angle>=m_AngularThreshold) { currentMean -= inDirs->at(i); if (meanChanged) length += lengths.at(i); touched[i] = 1; meanChanged = true; } } if(!meanChanged) currentMean = oldMean; else currentMean.normalize(); } // found stable mean outDirs->push_back(currentMean); newLengths.push_back(length); if (length>max) max = length; // find next unused seed free = false; for (unsigned int i=0; iat(i); free = true; meanChanged = false; length = lengths.at(i); touched[i] = 1; break; } } if (inDirs->size()==outDirs->size()) { if (max>0) for (unsigned int i=0; isize(); i++) outDirs->SetElement(i, outDirs->at(i)*newLengths.at(i)/max); return outDirs; } else return FastClustering(outDirs, newLengths); } //template< class PixelType > //std::vector< DirectionType > TractsToVectorImageFilter< PixelType >::Clustering(std::vector< DirectionType >& inDirs) //{ // std::vector< DirectionType > outDirs; // if (inDirs.empty()) // return outDirs; // DirectionType oldMean, currentMean, workingMean; // std::vector< DirectionType > normalizedDirs; // std::vector< int > touched; // for (std::size_t i=0; i0.0001) // { // counter = 0; // oldMean = currentMean; // workingMean = oldMean; // workingMean.normalize(); // currentMean.fill(0.0); // for (std::size_t i=0; i=m_AngularThreshold) // { // currentMean += inDirs[i]; // counter++; // } // else if (-angle>=m_AngularThreshold) // { // currentMean -= inDirs[i]; // counter++; // } // } // } // // found stable mean // if (counter>0) // { // bool add = true; // DirectionType normMean = currentMean; // normMean.normalize(); // for (std::size_t i=0; i0) // { // if (mag>max) // max = mag; // outDirs.push_back(currentMean); // } // } // } // } // if (m_NormalizeVectors) // for (std::size_t i=0; i0) // for (std::size_t i=0; i //TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::MeanShiftClustering(DirectionContainerType::Pointer dirCont) //{ // DirectionContainerType::Pointer container = DirectionContainerType::New(); // double max = 0; // for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it) // { // vnl_vector_fixed mean = ClusterStep(dirCont, it.Value()); // if (mean.is_zero()) // continue; // bool addMean = true; // for (DirectionContainerType::ConstIterator it2 = container->Begin(); it2!=container->End(); ++it2) // { // vnl_vector_fixed dir = it2.Value(); // double angle = fabs(dot_product(mean, dir)/(mean.magnitude()*dir.magnitude())); // if (angle>=m_Epsilon) // { // addMean = false; // break; // } // } // if (addMean) // { // if (m_NormalizeVectors) // mean.normalize(); // else if (mean.magnitude()>max) // max = mean.magnitude(); // container->InsertElement(container->Size(), mean); // } // } // // max normalize voxel directions // if (max>0 && !m_NormalizeVectors) // for (std::size_t i=0; iSize(); i++) // container->ElementAt(i) /= max; // if (container->Size()Size()) // return MeanShiftClustering(container); // else // return container; //} //template< class PixelType > //vnl_vector_fixed TractsToVectorImageFilter< PixelType >::ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed currentMean) //{ // vnl_vector_fixed newMean; newMean.fill(0); // for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it) // { // vnl_vector_fixed dir = it.Value(); // double angle = dot_product(currentMean, dir)/(currentMean.magnitude()*dir.magnitude()); // if (angle>=m_AngularThreshold) // newMean += dir; // else if (-angle>=m_AngularThreshold) // newMean -= dir; // } // if (fabs(dot_product(currentMean, newMean)/(currentMean.magnitude()*newMean.magnitude()))>=m_Epsilon || newMean.is_zero()) // return newMean; // else // return ClusterStep(dirCont, newMean); //} } diff --git a/Modules/DiffusionImaging/MiniApps/DFTracking.cpp b/Modules/DiffusionImaging/MiniApps/DFTracking.cpp index 93ab11c30f..093bed8d30 100755 --- a/Modules/DiffusionImaging/MiniApps/DFTracking.cpp +++ b/Modules/DiffusionImaging/MiniApps/DFTracking.cpp @@ -1,199 +1,194 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include //#include #include #include #include #include #include #include //#include #include #include #include #include #define _USE_MATH_DEFINES #include const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Machine Learning Based Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("image", "i", mitkCommandLineParser::String, "DWIs:", "input diffusion-weighted image", us::Any(), false); parser.addArgument("forest", "f", mitkCommandLineParser::String, "Forest:", "input forest", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle", us::Any(), false); parser.addArgument("stop", "st", mitkCommandLineParser::String, "Stop image:", "stop image", us::Any()); parser.addArgument("mask", "m", mitkCommandLineParser::String, "Mask image:", "mask image", us::Any()); parser.addArgument("seed", "s", mitkCommandLineParser::String, "Seed image:", "seed image", us::Any()); parser.addArgument("athres", "a", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold (in radians)", us::Any()); parser.addArgument("stepsize", "se", mitkCommandLineParser::Float, "Stepsize:", "stepsize", us::Any()); parser.addArgument("samples", "ns", mitkCommandLineParser::Int, "Samples:", "samples", us::Any()); parser.addArgument("samplingdist", "sd", mitkCommandLineParser::Float, "Sampling distance:", "sampling distance (in voxels)", us::Any()); parser.addArgument("seeds", "nse", mitkCommandLineParser::Int, "Seeds per voxel:", "seeds per voxel", us::Any()); - parser.addArgument("usedirection", "ud", mitkCommandLineParser::Bool, "Use previous direction:", "use previous direction as feature", us::Any()); parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Verbose:", "output additional images", us::Any()); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string imageFile = us::any_cast(parsedArgs["image"]); string forestFile = us::any_cast(parsedArgs["forest"]); string outFile = us::any_cast(parsedArgs["out"]); string maskFile = ""; if (parsedArgs.count("mask")) maskFile = us::any_cast(parsedArgs["mask"]); string seedFile = ""; if (parsedArgs.count("seed")) seedFile = us::any_cast(parsedArgs["seed"]); string stopFile = ""; if (parsedArgs.count("stop")) stopFile = us::any_cast(parsedArgs["stop"]); float stepsize = -1; if (parsedArgs.count("stepsize")) stepsize = us::any_cast(parsedArgs["stepsize"]); float athres = 0.7; if (parsedArgs.count("athres")) athres = us::any_cast(parsedArgs["athres"]); float samplingdist = 0.25; if (parsedArgs.count("samplingdist")) samplingdist = us::any_cast(parsedArgs["samplingdist"]); - bool useDirection = false; - if (parsedArgs.count("usedirection")) - useDirection = true; - bool verbose = false; if (parsedArgs.count("verbose")) verbose = true; int samples = 10; if (parsedArgs.count("samples")) samples = us::any_cast(parsedArgs["samples"]); int seeds = 1; if (parsedArgs.count("seeds")) seeds = us::any_cast(parsedArgs["seeds"]); typedef itk::Image ItkUcharImgType; MITK_INFO << "loading diffusion-weighted image"; mitk::Image::Pointer dwi = dynamic_cast(mitk::IOUtil::LoadImage(imageFile).GetPointer()); ItkUcharImgType::Pointer mask; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(maskFile).GetPointer()); mask = ItkUcharImgType::New(); mitk::CastToItkImage(img, mask); } ItkUcharImgType::Pointer seed; if (!seedFile.empty()) { MITK_INFO << "loading seed image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(seedFile).GetPointer()); seed = ItkUcharImgType::New(); mitk::CastToItkImage(img, seed); } ItkUcharImgType::Pointer stop; if (!stopFile.empty()) { MITK_INFO << "loading stop image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(stopFile).GetPointer()); stop = ItkUcharImgType::New(); mitk::CastToItkImage(img, stop); } MITK_INFO << "loading forest"; vigra::RandomForest rf; vigra::rf_import_HDF5(rf, forestFile); typedef itk::MLBSTrackingFilter<100> TrackerType; TrackerType::Pointer tracker = TrackerType::New(); tracker->SetInput(0, mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi)); tracker->SetGradientDirections( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi) ); tracker->SetB_Value( mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi) ); tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetSeedsPerVoxel(seeds); - tracker->SetUseDirection(useDirection); tracker->SetStepSize(stepsize); tracker->SetAngularThreshold(athres); tracker->SetDecisionForest(&rf); tracker->SetSamplingDistance(samplingdist); tracker->SetNumberOfSamples(samples); //tracker->SetAvoidStop(false); + tracker->SetVerbose(verbose); tracker->SetAposterioriCurvCheck(false); tracker->SetRemoveWmEndFibers(false); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); mitk::IOUtil::SaveBaseData(outFib, outFile); if (verbose) { MITK_INFO << "Writing images..."; string outName = itksys::SystemTools::GetFilenamePath(outFile)+"/"+itksys::SystemTools::GetFilenameWithoutLastExtension(outFile); itk::ImageFileWriter< TrackerType::ItkDoubleImgType >::Pointer writer = itk::ImageFileWriter< TrackerType::ItkDoubleImgType >::New(); writer->SetFileName(outName+"_WhiteMatter.nrrd"); writer->SetInput(tracker->GetWmImage()); writer->Update(); writer->SetFileName(outName+"_NotWhiteMatter.nrrd"); writer->SetInput(tracker->GetNotWmImage()); writer->Update(); writer->SetFileName(outName+"_AvoidStop.nrrd"); writer->SetInput(tracker->GetAvoidStopImage()); writer->Update(); } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/MiniApps/DFTraining.cpp b/Modules/DiffusionImaging/MiniApps/DFTraining.cpp index d4c2bc3d40..62b0b18593 100755 --- a/Modules/DiffusionImaging/MiniApps/DFTraining.cpp +++ b/Modules/DiffusionImaging/MiniApps/DFTraining.cpp @@ -1,529 +1,148 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include -#include -#include -#include -#include -#include - -#include -#include -#include +#include #define _USE_MATH_DEFINES #include -const int numOdfSamples = 200; // ODF is sampled in 200 directions but actuyll only 100 are used (symmetric) -typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; - -void TrainForest( vigra::RandomForest &rf, vigra::MultiArray<2, double> &labelData, vigra::MultiArray<2, double> &featureData, int numTrees, int max_tree_depth, double sample_fraction ) -{ - MITK_INFO << "Maximum tree depths: " << max_tree_depth; - MITK_INFO << "Sample fraction per tree: " << sample_fraction; - MITK_INFO << "Number of trees: " << numTrees; - vigra::rf::visitors::OOB_Error oob_v; - -// rf.set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal -// rf.set_options().sample_with_replacement(true); // if sampled with replacement or not -// rf.set_options().samples_per_tree(sample_fraction); // Fraction of samples that are used to train a tree -// rf.set_options().tree_count(1); // Number of trees that are calculated; -// rf.set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node -// rf.ext_param_.max_tree_depth = max_tree_depth; -// // rf.set_options().features_per_node(10); -// rf.learn(featureData, labelData, vigra::rf::visitors::create_visitor(oob_v)); - - std::vector< vigra::RandomForest > trees; - int count = 0; -#pragma omp parallel for - for (int i = 0; i < numTrees; ++i) - { - vigra::RandomForest lrf; - vigra::rf::visitors::OOB_Error loob_v; - - lrf.set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal - lrf.set_options().sample_with_replacement(true); // if sampled with replacement or not - lrf.set_options().samples_per_tree(sample_fraction); // Fraction of samples that are used to train a tree - lrf.set_options().tree_count(1); // Number of trees that are calculated; - lrf.set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node - lrf.ext_param_.max_tree_depth = max_tree_depth; - // lrf.set_options().features_per_node(10); - - lrf.learn(featureData, labelData);//, vigra::rf::visitors::create_visitor(loob_v)); -#pragma omp critical - { - count++; - MITK_INFO << "Tree " << count << " finished training."; - trees.push_back(lrf); - //rf.trees_.push_back(lrf.trees_[0]); - } - } - - for (int i = 1; i < numTrees; ++i) - trees.at(0).trees_.push_back(trees.at(i).trees_[0]); - - rf = trees.at(0); - rf.options_.tree_count_ = numTrees; - MITK_INFO << "Training finsihed"; - //MITK_INFO << "The out-of-bag error is: " << oob_v.oob_breiman << std::endl; -} - -SampledShImageType::PixelType GetImageValues(itk::Point itkP, SampledShImageType::Pointer image) -{ - itk::Index<3> idx; - itk::ContinuousIndex< double, 3> cIdx; - image->TransformPhysicalPointToIndex(itkP, idx); - image->TransformPhysicalPointToContinuousIndex(itkP, cIdx); - - SampledShImageType::PixelType pix; pix.Fill(0.0); - if ( image->GetLargestPossibleRegion().IsInside(idx) ) - pix = image->GetPixel(idx); - else - return pix; - - double frac_x = cIdx[0] - idx[0]; - double frac_y = cIdx[1] - idx[1]; - double frac_z = cIdx[2] - idx[2]; - if (frac_x<0) - { - idx[0] -= 1; - frac_x += 1; - } - if (frac_y<0) - { - idx[1] -= 1; - frac_y += 1; - } - if (frac_z<0) - { - idx[2] -= 1; - frac_z += 1; - } - frac_x = 1-frac_x; - frac_y = 1-frac_y; - frac_z = 1-frac_z; - - // int coordinates inside image? - if (idx[0] >= 0 && idx[0] < image->GetLargestPossibleRegion().GetSize(0)-1 && - idx[1] >= 0 && idx[1] < image->GetLargestPossibleRegion().GetSize(1)-1 && - idx[2] >= 0 && idx[2] < image->GetLargestPossibleRegion().GetSize(2)-1) - { - vnl_vector_fixed interpWeights; - interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); - interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); - interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); - interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); - interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); - interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); - interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); - interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); - - - pix = image->GetPixel(idx) * interpWeights[0]; - SampledShImageType::IndexType tmpIdx = idx; tmpIdx[0]++; - pix += image->GetPixel(tmpIdx) * interpWeights[1]; - tmpIdx = idx; tmpIdx[1]++; - pix += image->GetPixel(tmpIdx) * interpWeights[2]; - tmpIdx = idx; tmpIdx[2]++; - pix += image->GetPixel(tmpIdx) * interpWeights[3]; - tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; - pix += image->GetPixel(tmpIdx) * interpWeights[4]; - tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; - pix += image->GetPixel(tmpIdx) * interpWeights[5]; - tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; - pix += image->GetPixel(tmpIdx) * interpWeights[6]; - tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; - pix += image->GetPixel(tmpIdx) * interpWeights[7]; - } - - return pix; -} - int main(int argc, char* argv[]) { MITK_INFO << "DFTraining"; mitkCommandLineParser parser; parser.setTitle("Machine Learning Based Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("images", "i", mitkCommandLineParser::StringList, "DWIs:", "input diffusion-weighted images", us::Any(), false); - parser.addArgument("wmmasks", "w", mitkCommandLineParser::StringList, "WM-Masks:", "white matter mask images", us::Any(), false); + parser.addArgument("wmmasks", "w", mitkCommandLineParser::StringList, "WM-Masks:", "white matter mask images", us::Any()); parser.addArgument("tractograms", "t", mitkCommandLineParser::StringList, "Tractograms:", "input tractograms (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("masks", "m", mitkCommandLineParser::StringList, "Masks:", "mask images", us::Any()); parser.addArgument("forest", "f", mitkCommandLineParser::OutputFile, "Forest:", "output forest", us::Any(), false); parser.addArgument("stepsize", "s", mitkCommandLineParser::Float, "Stepsize:", "stepsize", us::Any()); parser.addArgument("gmsamples", "g", mitkCommandLineParser::Int, "Number of gray matter samples per voxel:", "Number of gray matter samples per voxel", us::Any()); parser.addArgument("numtrees", "n", mitkCommandLineParser::Int, "Number of trees:", "number of trees", us::Any()); parser.addArgument("max_tree_depth", "d", mitkCommandLineParser::Int, "Max. tree depth:", "maximum tree depth", us::Any()); parser.addArgument("sample_fraction", "sf", mitkCommandLineParser::Float, "Sample fraction:", "fraction of samples used per tree", us::Any()); - parser.addArgument("usedirection", "ud", mitkCommandLineParser::Bool, "bla:", "bla", us::Any()); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType imageFiles = us::any_cast(parsedArgs["images"]); - mitkCommandLineParser::StringContainerType wmMaskFiles = us::any_cast(parsedArgs["wmmasks"]); + mitkCommandLineParser::StringContainerType wmMaskFiles; + if (parsedArgs.count("wmmasks")) + wmMaskFiles = us::any_cast(parsedArgs["wmmasks"]); mitkCommandLineParser::StringContainerType maskFiles; if (parsedArgs.count("masks")) maskFiles = us::any_cast(parsedArgs["masks"]); string forestFile = us::any_cast(parsedArgs["forest"]); mitkCommandLineParser::StringContainerType tractogramFiles; if (parsedArgs.count("tractograms")) tractogramFiles = us::any_cast(parsedArgs["tractograms"]); int numTrees = 30; if (parsedArgs.count("numtrees")) numTrees = us::any_cast(parsedArgs["numtrees"]); int gmsamples = 50; if (parsedArgs.count("gmsamples")) gmsamples = us::any_cast(parsedArgs["gmsamples"]); float stepsize = -1; if (parsedArgs.count("stepsize")) stepsize = us::any_cast(parsedArgs["stepsize"]); int max_tree_depth = 50; if (parsedArgs.count("max_tree_depth")) max_tree_depth = us::any_cast(parsedArgs["max_tree_depth"]); double sample_fraction = 1.0; if (parsedArgs.count("sample_fraction")) sample_fraction = us::any_cast(parsedArgs["sample_fraction"]); - // load DWI images - if (imageFiles.size() QballFilterType; - MITK_INFO << "loading diffusion-weighted images and reconstructing feature images"; - std::vector< SampledShImageType::Pointer > sampledShImages; + MITK_INFO << "loading diffusion-weighted images"; + std::vector< mitk::Image::Pointer > rawData; for (unsigned int i=0; i(mitk::IOUtil::LoadImage(imageFiles.at(i)).GetPointer()); - - QballFilterType::Pointer qballfilter = QballFilterType::New(); - qballfilter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi) ); - qballfilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi)); - qballfilter->SetLambda(0.006); - qballfilter->SetNormalizationMethod(QballFilterType::QBAR_RAW_SIGNAL); - qballfilter->Update(); - // FeatureImageType::Pointer itkFeatureImage = qballfilter->GetCoefficientImage(); - // featureImageVector.push_back(itkFeatureImage); - sampledShImages.push_back(qballfilter->GetOutput()); + rawData.push_back(dwi); } - typedef itk::Image ItkUcharImgType; + MITK_INFO << "loading mask images"; std::vector< ItkUcharImgType::Pointer > maskImageVector; + for (unsigned int i=0; i(mitk::IOUtil::LoadImage(maskFiles.at(i)).GetPointer()); + ItkUcharImgType::Pointer mask = ItkUcharImgType::New(); + mitk::CastToItkImage(img, mask); + maskImageVector.push_back(mask); + } + + MITK_INFO << "loading white matter mask images"; std::vector< ItkUcharImgType::Pointer > wmMaskImageVector; + for (unsigned int i=0; i(mitk::IOUtil::LoadImage(wmMaskFiles.at(i)).GetPointer()); + ItkUcharImgType::Pointer wmmask = ItkUcharImgType::New(); + mitk::CastToItkImage(img, wmmask); + wmMaskImageVector.push_back(wmmask); + } MITK_INFO << "loading tractograms"; - int numSamples = 0; std::vector< mitk::FiberBundle::Pointer > tractograms; for (unsigned int t=0; t(mitk::IOUtil::LoadImage(maskFiles.at(t)).GetPointer()); - mask = ItkUcharImgType::New(); - mitk::CastToItkImage(img, mask); - maskImageVector.push_back(mask); - } - mitk::Image::Pointer img2 = dynamic_cast(mitk::IOUtil::LoadImage(wmMaskFiles.at(t)).GetPointer()); - ItkUcharImgType::Pointer wmmask = ItkUcharImgType::New(); - mitk::CastToItkImage(img2, wmmask); - wmMaskImageVector.push_back(wmmask); - - itk::ImageRegionConstIterator it(wmmask, wmmask->GetLargestPossibleRegion()); - int OUTOFWM = 0; // count voxels outside of the white matter mask - while(!it.IsAtEnd()) - { - if (it.Get()==0) - if (mask.IsNull() || (mask.IsNotNull() && mask->GetPixel(it.GetIndex())>0)) - OUTOFWM++; - ++it; - } - numSamples += gmsamples*OUTOFWM; // for each of the non-white matter voxels we add a certain number of sampling points. these sampling points are used to tell the classifier where to recognize non-WM tissue - - MITK_INFO << "Samples outside of WM: " << numSamples << " (" << gmsamples << " per non-WM voxel)"; - - // load and resample training tractograms mitk::FiberBundle::Pointer fib = dynamic_cast(mitk::IOUtil::Load(tractogramFiles.at(t)).at(0).GetPointer()); - if (stepsize<0) - { - SampledShImageType::Pointer image = sampledShImages.at(t); - float minSpacing = 1; - if(image->GetSpacing()[0]GetSpacing()[1] && image->GetSpacing()[0]GetSpacing()[2]) - minSpacing = image->GetSpacing()[0]; - else if (image->GetSpacing()[1] < image->GetSpacing()[2]) - minSpacing = image->GetSpacing()[1]; - else - minSpacing = image->GetSpacing()[2]; - stepsize = minSpacing*0.5; - } - fib->ResampleSpline(stepsize); tractograms.push_back(fib); - numSamples += fib->GetNumberOfPoints(); // each point of the fiber gives us a training direction - numSamples -= 2*fib->GetNumFibers(); // we don't use the first and last point because there we do not have a previous direction, which is needed as feature - } - MITK_INFO << "Number of samples: " << numSamples; - - // get ODF directions and number of features - vnl_vector_fixed ref; ref.fill(0); ref[0]=1; - itk::OrientationDistributionFunction< double, numOdfSamples > odf; - std::vector< int > directionIndices; - for (unsigned int f=0; f0) // we only use directions on one hemisphere (symmetric) - directionIndices.push_back(f); // remember indices that are on the desired hemisphere - } - const int numSignalFeatures = numOdfSamples/2; - int numDirectionFeatures = 0; - if (useDirection) - numDirectionFeatures = 3; - - vigra::MultiArray<2, double> featureData( vigra::Shape2(numSamples,numSignalFeatures+numDirectionFeatures) ); - MITK_INFO << "Number of features: " << featureData.shape(1); - vigra::MultiArray<2, double> labelData( vigra::Shape2(numSamples,1) ); - - itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); - m_RandGen->SetSeed(); - MITK_INFO << "Creating training data from tractograms and feature images"; - int sampleCounter = 0; - for (unsigned int t=0; t it(wmMask, wmMask->GetLargestPossibleRegion()); - while(!it.IsAtEnd()) - { - if (it.Get()==0 && (mask.IsNull() || (mask.IsNotNull() && mask->GetPixel(it.GetIndex())>0))) - { - SampledShImageType::PixelType pix = image->GetPixel(it.GetIndex()); - if (useDirection) - { - // null direction - for (unsigned int f=0; f probe; - probe[0] = m_RandGen->GetVariate()*2-1; - probe[1] = m_RandGen->GetVariate()*2-1; - probe[2] = m_RandGen->GetVariate()*2-1; - probe.normalize(); - if (dot_product(ref, probe)<0) - probe *= -1; - for (unsigned int f=numSignalFeatures; f idx; - idx[0] = it.GetIndex()[0]; - idx[1] = it.GetIndex()[1]; - idx[2] = it.GetIndex()[2]; - itk::Point itkP1; - image->TransformContinuousIndexToPhysicalPoint(idx, itkP1); - SampledShImageType::PixelType pix = GetImageValues(itkP1, image);; - for (unsigned int f=0; f idx; - idx[0] = it.GetIndex()[0] + m_RandGen->GetVariate()-0.5; - idx[1] = it.GetIndex()[1] + m_RandGen->GetVariate()-0.5; - idx[2] = it.GetIndex()[2] + m_RandGen->GetVariate()-0.5; - itk::Point itkP1; - image->TransformContinuousIndexToPhysicalPoint(idx, itkP1); - SampledShImageType::PixelType pix = GetImageValues(itkP1, image);; - for (unsigned int f=0; f polyData = fib->GetFiberPolyData(); - for (int i=0; iGetNumFibers(); i++) - { - vtkCell* cell = polyData->GetCell(i); - int numPoints = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - - vnl_vector_fixed dirOld; dirOld.fill(0.0); - - for (int j=0; jGetPoint(j); - itk::Point itkP1; - itkP1[0] = p1[0]; itkP1[1] = p1[1]; itkP1[2] = p1[2]; - - vnl_vector_fixed dir; dir.fill(0.0); - - itk::Point itkP2; - double* p2 = points->GetPoint(j+1); - itkP2[0] = p2[0]; itkP2[1] = p2[1]; itkP2[2] = p2[2]; - dir[0]=itkP2[0]-itkP1[0]; - dir[1]=itkP2[1]-itkP1[1]; - dir[2]=itkP2[2]-itkP1[2]; - - if (dir.magnitude()<0.0001) - { - MITK_INFO << "streamline error!"; - continue; - } - dir.normalize(); - if (dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2]) - { - MITK_INFO << "ERROR: NaN direction!"; - continue; - } - - if (j==0) - { - dirOld = dir; - continue; - } - - // get voxel values - SampledShImageType::PixelType pix = GetImageValues(itkP1, image); - for (unsigned int f=0; f0.0001) - { - int label = 0; - for (unsigned int f=0; fangle) - { - labelData(sampleCounter,0) = f; - angle = a; - label = f; - } - } - } - - dirOld = dir; - sampleCounter++; - } - } } - MITK_INFO << "Training forest"; - vigra::RandomForest rf; - TrainForest( rf, labelData, featureData, numTrees, max_tree_depth, sample_fraction ); - MITK_INFO << "Writing forest"; - vigra::rf_export_HDF5( rf, forestFile, "" ); - MITK_INFO << "Finished training"; + mitk::TrackingForestHandler<> forestHandler; + forestHandler.SetRawData(rawData); + forestHandler.SetTractograms(tractograms); + forestHandler.SetNumTrees(numTrees); + forestHandler.SetMaxTreeDepth(max_tree_depth); + forestHandler.SetGrayMatterSamplesPerVoxel(gmsamples); + forestHandler.SetSampleFraction(sample_fraction); + forestHandler.SetStepSize(stepsize); + forestHandler.StartTraining(); + forestHandler.SaveForest(forestFile); return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp b/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp index 326838cef7..64d7707814 100755 --- a/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp +++ b/Modules/DiffusionImaging/MiniApps/LocalDirectionalFiberPlausibility.cpp @@ -1,300 +1,319 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Local Directional Fiber Plausibility"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription(" "); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("reference", "r", mitkCommandLineParser::StringList, "Reference images:", "reference direction images", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::StringList, "Masks:", "mask images"); parser.addArgument("athresh", "a", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); + parser.addArgument("sthresh", "s", mitkCommandLineParser::Float, "Size threshold:", "Relative peak size threshold per voxel.", 0.0, true); + parser.addArgument("maxdirs", "md", mitkCommandLineParser::Int, "Max. Clusters:", "Maximum number of fiber clusters.", 0, true); parser.addArgument("verbose", "v", mitkCommandLineParser::Bool, "Verbose:", "output optional and intermediate calculation results"); parser.addArgument("ignore", "n", mitkCommandLineParser::Bool, "Ignore:", "don't increase error for missing or too many directions"); + parser.addArgument("empty", "e", mitkCommandLineParser::Bool, "Empty Voxels:", "don't increase error for empty voxels"); parser.addArgument("fileID", "id", mitkCommandLineParser::String, "ID:", "optional ID field"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType referenceImages = us::any_cast(parsedArgs["reference"]); mitkCommandLineParser::StringContainerType maskImages; if (parsedArgs.count("mask")) maskImages = us::any_cast(parsedArgs["mask"]); string fibFile = us::any_cast(parsedArgs["input"]); float angularThreshold = 25; if (parsedArgs.count("athresh")) angularThreshold = us::any_cast(parsedArgs["athresh"]); + float sizeThreshold = 0; + if (parsedArgs.count("sthresh")) + sizeThreshold = us::any_cast(parsedArgs["sthresh"]); + + int maxDirs = 0; + if (parsedArgs.count("maxdirs")) + maxDirs = us::any_cast(parsedArgs["maxdirs"]); + string outRoot = us::any_cast(parsedArgs["out"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); - bool ignore = false; + bool ignoreMissing = false; if (parsedArgs.count("ignore")) - ignore = us::any_cast(parsedArgs["ignore"]); + ignoreMissing = us::any_cast(parsedArgs["ignore"]); + + bool ignoreEmpty = false; + if (parsedArgs.count("empty")) + ignoreEmpty = us::any_cast(parsedArgs["empty"]); string fileID = ""; if (parsedArgs.count("fileID")) fileID = us::any_cast(parsedArgs["fileID"]); try { typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< unsigned int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::LoadDataNode(fibFile)->GetData()); // load reference directions ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData()); typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); } catch(...){ std::cout << "could not load: " << referenceImages.at(i); } } ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); itkMaskImage->SetSpacing( dirImg->GetSpacing() ); itkMaskImage->SetOrigin( dirImg->GetOrigin() ); itkMaskImage->SetDirection( dirImg->GetDirection() ); itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->Allocate(); itkMaskImage->FillBuffer(1); // extract directions from fiber bundle itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetMaskImage(itkMaskImage); fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); fOdfFilter->SetNormalizeVectors(true); fOdfFilter->SetUseWorkingCopy(false); + fOdfFilter->SetSizeThreshold(sizeThreshold); + fOdfFilter->SetMaxNumDirections(maxDirs); fOdfFilter->Update(); ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); if (verbose) { // write vector field mitk::FiberBundle::Pointer directions = fOdfFilter->GetOutputFiberBundle(); string outfilename = outRoot; outfilename.append("_VECTOR_FIELD.fib"); mitk::IOUtil::SaveBaseData(directions.GetPointer(), outfilename ); // write direction images for (unsigned int i=0; iSize(); i++) { itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_DIRECTION_"); outfilename.append(boost::lexical_cast(i)); outfilename.append(".nrrd"); writer->SetFileName(outfilename.c_str()); writer->SetInput(itkImg); writer->Update(); } // write num direction image { ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_NUM_DIRECTIONS.nrrd"); writer->SetFileName(outfilename.c_str()); writer->SetInput(numDirImage); writer->Update(); } } string logFile = outRoot; logFile.append("_ANGULAR_ERROR.csv"); ofstream file; file.open (logFile.c_str()); if (maskImages.size()>0) { for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(maskImages.at(i))->GetData()); mitk::CastToItkImage(mitkMaskImage, itkMaskImage); // evaluate directions EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); evaluationFilter->SetImageSet(directionImageContainer); evaluationFilter->SetReferenceImageSet(referenceImageContainer); evaluationFilter->SetMaskImage(itkMaskImage); - evaluationFilter->SetIgnoreMissingDirections(ignore); + evaluationFilter->SetIgnoreMissingDirections(ignoreMissing); + evaluationFilter->SetIgnoreEmptyVoxels(ignoreEmpty); evaluationFilter->Update(); if (verbose) { EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_ERROR_IMAGE.nrrd"); writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } string maskFileName = itksys::SystemTools::GetFilenameWithoutExtension(maskImages.at(i)); unsigned found = maskFileName.find_last_of("_"); string sens = itksys::SystemTools::GetFilenameWithoutLastExtension(fibFile); if (!fileID.empty()) sens = fileID; sens.append(","); sens.append(maskFileName.substr(found+1)); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append(";\n"); file << sens; } } else { // evaluate directions EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); evaluationFilter->SetImageSet(directionImageContainer); evaluationFilter->SetReferenceImageSet(referenceImageContainer); evaluationFilter->SetMaskImage(itkMaskImage); - evaluationFilter->SetIgnoreMissingDirections(ignore); + evaluationFilter->SetIgnoreMissingDirections(ignoreMissing); + evaluationFilter->SetIgnoreEmptyVoxels(ignoreEmpty); evaluationFilter->Update(); if (verbose) { EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_ERROR_IMAGE.nrrd"); writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } string sens = itksys::SystemTools::GetFilenameWithoutLastExtension(fibFile); if (!fileID.empty()) sens = fileID; sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append(";\n"); file << sens; } file.close(); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/LegacyGL/mitkPointSetGLMapper2D.cpp b/Modules/LegacyGL/mitkPointSetGLMapper2D.cpp index 7c71d1530e..a475a393ca 100644 --- a/Modules/LegacyGL/mitkPointSetGLMapper2D.cpp +++ b/Modules/LegacyGL/mitkPointSetGLMapper2D.cpp @@ -1,513 +1,513 @@ /*=================================================================== 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 "mitkPointSetGLMapper2D.h" #include "mitkPointSet.h" #include "mitkPlaneGeometry.h" #include "mitkColorProperty.h" #include "mitkProperties.h" #include "vtkLinearTransform.h" #include "mitkStringProperty.h" #include "mitkPointSet.h" #include "mitkVtkPropRenderer.h" #include "mitkGL.h" //const float selectedColor[]={1.0,0.0,0.6}; //for selected! mitk::PointSetGLMapper2D::PointSetGLMapper2D() : m_Polygon(false), m_ShowPoints(true), m_ShowDistances(false), m_DistancesDecimalDigits(1), m_ShowAngles(false), m_ShowDistantLines(true), m_LineWidth(1) { } mitk::PointSetGLMapper2D::~PointSetGLMapper2D() { } const mitk::PointSet *mitk::PointSetGLMapper2D::GetInput(void) { return static_cast ( GetDataNode()->GetData() ); } void mitk::PointSetGLMapper2D::ApplyAllProperties(mitk::BaseRenderer* renderer) { GLMapper::ApplyColorAndOpacityProperties( renderer ); const mitk::DataNode* node=GetDataNode(); if( node == NULL ) return; node->GetBoolProperty("show contour", m_Polygon); node->GetBoolProperty("close contour", m_PolygonClosed); node->GetBoolProperty("show points", m_ShowPoints); node->GetBoolProperty("show distances", m_ShowDistances); node->GetIntProperty("distance decimal digits", m_DistancesDecimalDigits); node->GetBoolProperty("show angles", m_ShowAngles); node->GetBoolProperty("show distant lines", m_ShowDistantLines); node->GetIntProperty("line width", m_LineWidth); node->GetIntProperty("point line width", m_PointLineWidth); - node->GetIntProperty("point 2D size", m_Point2DSize); + node->GetFloatProperty("point 2D size", m_Point2DSize); } static bool makePerpendicularVector2D(const mitk::Vector2D& in, mitk::Vector2D& out) { if((fabs(in[0])>0) && ( (fabs(in[0])>fabs(in[1])) || (in[1] == 0) ) ) { out[0]=-in[1]/in[0]; out[1]=1; out.Normalize(); return true; } else if(fabs(in[1])>0) { out[0]=1; out[1]=-in[0]/in[1]; out.Normalize(); return true; } else return false; } void mitk::PointSetGLMapper2D::Paint( mitk::BaseRenderer *renderer ) { const mitk::DataNode* node=GetDataNode(); if( node == NULL ) return; const int text2dDistance = 10; bool visible = true; GetDataNode()->GetVisibility(visible, renderer, "visible"); if ( !visible) return; // @FIXME: Logik fuer update bool updateNeccesary=true; if (updateNeccesary) { // ok, das ist aus GenerateData kopiert mitk::PointSet::Pointer input = const_cast(this->GetInput()); // Get the TimeGeometry of the input object const TimeGeometry* inputTimeGeometry = input->GetTimeGeometry(); if (( inputTimeGeometry == NULL ) || ( inputTimeGeometry->CountTimeSteps() == 0 ) ) { return; } // // get the world time // ScalarType time = renderer->GetTime(); // // convert the world time in time steps of the input object // int timeStep=0; if ( time > itk::NumericTraits::NonpositiveMin() ) timeStep = inputTimeGeometry->TimePointToTimeStep( time ); if ( inputTimeGeometry->IsValidTimeStep( timeStep ) == false ) { return; } mitk::PointSet::DataType::Pointer itkPointSet = input->GetPointSet( timeStep ); if ( itkPointSet.GetPointer() == NULL) { return; } mitk::DisplayGeometry::Pointer displayGeometry = renderer->GetDisplayGeometry(); assert(displayGeometry.IsNotNull()); //apply color and opacity read from the PropertyList this->ApplyAllProperties(renderer); vtkLinearTransform* transform = GetDataNode()->GetVtkTransform(); //List of the Points PointSet::DataType::PointsContainerConstIterator it, end; it = itkPointSet->GetPoints()->Begin(); end = itkPointSet->GetPoints()->End(); //iterator on the additional data of each point PointSet::DataType::PointDataContainerIterator selIt, selEnd; bool pointDataBroken = (itkPointSet->GetPointData()->Size() != itkPointSet->GetPoints()->Size()); selIt = itkPointSet->GetPointData()->Begin(); selEnd = itkPointSet->GetPointData()->End(); int counter = 0; //for writing text int j = 0; //for switching back to old color after using selected color float recallColor[4]; glGetFloatv(GL_CURRENT_COLOR,recallColor); //get the properties for coloring the points float unselectedColor[4] = {1.0, 1.0, 0.0, 1.0};//yellow //check if there is an unselected property if (dynamic_cast(node->GetPropertyList(renderer)->GetProperty("unselectedcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("unselectedcolor"))->GetValue(); unselectedColor[0] = tmpColor[0]; unselectedColor[1] = tmpColor[1]; unselectedColor[2] = tmpColor[2]; unselectedColor[3] = 1.0f; //!!define a new ColorProp to be able to pass alpha value } else if (dynamic_cast(node->GetPropertyList(NULL)->GetProperty("unselectedcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("unselectedcolor"))->GetValue(); unselectedColor[0] = tmpColor[0]; unselectedColor[1] = tmpColor[1]; unselectedColor[2] = tmpColor[2]; unselectedColor[3] = 1.0f; //!!define a new ColorProp to be able to pass alpha value } else { //get the color from the dataNode node->GetColor(unselectedColor, NULL); } //get selected property float selectedColor[4] = {1.0, 0.0, 0.6, 1.0}; if (dynamic_cast(node->GetPropertyList(renderer)->GetProperty("selectedcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("selectedcolor"))->GetValue(); selectedColor[0] = tmpColor[0]; selectedColor[1] = tmpColor[1]; selectedColor[2] = tmpColor[2]; selectedColor[3] = 1.0f; } else if (dynamic_cast(node->GetPropertyList(NULL)->GetProperty("selectedcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("selectedcolor"))->GetValue(); selectedColor[0] = tmpColor[0]; selectedColor[1] = tmpColor[1]; selectedColor[2] = tmpColor[2]; selectedColor[3] = 1.0f; } //check if there is an pointLineWidth property if (dynamic_cast(node->GetPropertyList(renderer)->GetProperty("point line width")) != NULL) { m_PointLineWidth = dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("point line width"))->GetValue(); } else if (dynamic_cast(node->GetPropertyList(NULL)->GetProperty("point line width")) != NULL) { m_PointLineWidth = dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("point line width"))->GetValue(); } //check if there is an point 2D size property - if (dynamic_cast(node->GetPropertyList(renderer)->GetProperty("point 2D size")) != NULL) + if (dynamic_cast(node->GetPropertyList(renderer)->GetProperty("point 2D size")) != NULL) { - m_Point2DSize = dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("point 2D size"))->GetValue(); + m_Point2DSize = dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("point 2D size"))->GetValue(); } - else if (dynamic_cast(node->GetPropertyList(NULL)->GetProperty("point 2D size")) != NULL) + else if (dynamic_cast(node->GetPropertyList(NULL)->GetProperty("point 2D size")) != NULL) { - m_Point2DSize = dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("point 2D size"))->GetValue(); + m_Point2DSize = dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("point 2D size"))->GetValue(); } Point3D p; // currently visited point Point3D lastP; // last visited point Vector3D vec; // p - lastP Vector3D lastVec; // lastP - point before lastP vec.Fill(0); mitk::Point3D projected_p; // p projected on viewplane Point2D pt2d; // projected_p in display coordinates Point2D lastPt2d; // last projected_p in display coordinates Point2D preLastPt2d;// projected_p in display coordinates before lastPt2d Point2D lastPt2DInPointSet; // The last point in the pointset in display coordinates mitk::PointSet::DataType::PointType plob; plob.Fill(0); itkPointSet->GetPoint( itkPointSet->GetNumberOfPoints()-1, &plob); //map lastPt2DInPointSet to display coordinates float vtkp[3]; itk2vtk(plob, vtkp); transform->TransformPoint(vtkp, vtkp); vtk2itk(vtkp,p); displayGeometry->Project(p, projected_p); displayGeometry->Map(projected_p, lastPt2DInPointSet); displayGeometry->WorldToDisplay(lastPt2DInPointSet, lastPt2DInPointSet); while(it!=end) // iterate over all points { lastP = p; // valid only for counter > 0 lastVec = vec; // valid only for counter > 1 preLastPt2d = lastPt2d; // valid only for counter > 1 lastPt2d = pt2d; // valid only for counter > 0 itk2vtk(it->Value(), vtkp); transform->TransformPoint(vtkp, vtkp); vtk2itk(vtkp,p); vec = p-lastP; // valid only for counter > 0 displayGeometry->Project(p, projected_p); Vector3D diff=p-projected_p; ScalarType scalardiff = diff.GetSquaredNorm(); //MouseOrientation bool isInputDevice=false; bool isRendererSlice = scalardiff < 0.00001; //cause roundoff error if(this->GetDataNode()->GetBoolProperty("inputdevice",isInputDevice) && isInputDevice && !isRendererSlice ) { displayGeometry->Map(projected_p, pt2d); displayGeometry->WorldToDisplay(pt2d, pt2d); //Point size depending of distance to slice /*float p_size = (1/scalardiff)*10*m_Point2DSize; if(p_size < m_Point2DSize * 0.6 ) p_size = m_Point2DSize * 0.6 ; else if ( p_size > m_Point2DSize ) p_size = m_Point2DSize;*/ float p_size = (1/scalardiff)*100.0; if(p_size < 6.0 ) p_size = 6.0 ; else if ( p_size > 10.0 ) p_size = 10.0; //draw Point float opacity = (p_size<8)?0.3:1.0;//don't get the opacity from the node? Feature not a bug! Otehrwise the 2D cross is hardly seen. glColor4f(unselectedColor[0],unselectedColor[1],unselectedColor[2],opacity); glPointSize(p_size); //glShadeModel(GL_FLAT); glBegin (GL_POINTS); glVertex2dv(&pt2d[0]); glEnd (); } //for point set if(!isInputDevice && ( (scalardiff<4.0) || (m_Polygon))) { Point2D tmp; displayGeometry->Map(projected_p, pt2d); displayGeometry->WorldToDisplay(pt2d, pt2d); Vector2D horz,vert; horz[0]=(float)m_Point2DSize-scalardiff*2; horz[1]=0; vert[0]=0; vert[1]=(float)m_Point2DSize-scalardiff*2; // now paint text if available if (dynamic_cast(this->GetDataNode() ->GetProperty("label")) != NULL) { const char * pointLabel = dynamic_cast( this->GetDataNode()->GetProperty("label"))->GetValue(); std::string l = pointLabel; if (input->GetSize()>1) { // char buffer[20]; // sprintf(buffer,"%d",it->Index()); std::stringstream ss; ss << it->Index(); l.append(ss.str()); } mitk::VtkPropRenderer* OpenGLrenderer = dynamic_cast( renderer ); float rgb[3];//yellow rgb[0] = unselectedColor[0]; rgb[1] = unselectedColor[1]; rgb[2] = unselectedColor[2]; OpenGLrenderer->WriteSimpleText(l, pt2d[0] + text2dDistance, pt2d[1] + text2dDistance,rgb[0], rgb[1],rgb[2]); } if((m_ShowPoints) && (scalardiff<4.0)) { //check if the point is to be marked as selected if(selIt != selEnd || pointDataBroken) { bool addAsSelected = false; if (pointDataBroken) addAsSelected = false; else if (selIt->Value().selected) addAsSelected = true; else addAsSelected = false; if (addAsSelected) { horz[0]=(float)m_Point2DSize; vert[1]=(float)m_Point2DSize; glColor3f(selectedColor[0],selectedColor[1],selectedColor[2]); glLineWidth(m_PointLineWidth); //a diamond around the point with the selected color glBegin (GL_LINE_LOOP); tmp=pt2d-horz; glVertex2dv(&tmp[0]); tmp=pt2d+vert; glVertex2dv(&tmp[0]); tmp=pt2d+horz; glVertex2dv(&tmp[0]); tmp=pt2d-vert; glVertex2dv(&tmp[0]); glEnd (); glLineWidth(1); //the actual point in the specified color to see the usual color of the point glColor3f(unselectedColor[0],unselectedColor[1],unselectedColor[2]); glPointSize(1); glBegin (GL_POINTS); tmp=pt2d; glVertex2dv(&tmp[0]); glEnd (); } else //if not selected { glColor3f(unselectedColor[0],unselectedColor[1],unselectedColor[2]); glLineWidth(m_PointLineWidth); //drawing crosses glBegin (GL_LINES); tmp=pt2d-horz; glVertex2dv(&tmp[0]); tmp=pt2d+horz; glVertex2dv(&tmp[0]); tmp=pt2d-vert; glVertex2dv(&tmp[0]); tmp=pt2d+vert; glVertex2dv(&tmp[0]); glEnd (); glLineWidth(1); } } } bool drawLinesEtc = true; if (!m_ShowDistantLines && counter > 0) // check, whether this line should be drawn { ScalarType currentDistance = displayGeometry->GetWorldGeometry()->SignedDistance(p); ScalarType lastDistance = displayGeometry->GetWorldGeometry()->SignedDistance(lastP); if ( currentDistance * lastDistance > 0.5 ) // points on same side of plane drawLinesEtc = false; } // draw a line if ((m_Polygon && counter>0 && drawLinesEtc) || (m_Polygon && m_PolygonClosed && drawLinesEtc)) { if ((counter == 0) && ( m_PolygonClosed)) { lastPt2d = lastPt2DInPointSet; } //get contour color property float contourColor[4] = {unselectedColor[0], unselectedColor[1], unselectedColor[2], unselectedColor[3]};//so if no property set, then use unselected color if (dynamic_cast(node->GetPropertyList(renderer)->GetProperty("contourcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(renderer)->GetProperty("contourcolor"))->GetValue(); contourColor[0] = tmpColor[0]; contourColor[1] = tmpColor[1]; contourColor[2] = tmpColor[2]; contourColor[3] = 1.0f; } else if (dynamic_cast(node->GetPropertyList(NULL)->GetProperty("contourcolor")) != NULL) { mitk::Color tmpColor = dynamic_cast(this->GetDataNode()->GetPropertyList(NULL)->GetProperty("contourcolor"))->GetValue(); contourColor[0] = tmpColor[0]; contourColor[1] = tmpColor[1]; contourColor[2] = tmpColor[2]; contourColor[3] = 1.0f; } //set this color glColor3f(contourColor[0],contourColor[1],contourColor[2]); glLineWidth( m_LineWidth ); glBegin (GL_LINES); glVertex2dv(&pt2d[0]); glVertex2dv(&lastPt2d[0]); glEnd (); glLineWidth(1.0); if(m_ShowDistances) // calculate and print a distance { std::stringstream buffer; float distance = vec.GetNorm(); buffer<( renderer ); OpenGLrenderer->WriteSimpleText(buffer.str(), pos2d[0], pos2d[1]); //this->WriteTextXY(pos2d[0], pos2d[1], buffer.str(),renderer); } if(m_ShowAngles && counter > 1 ) // calculate and print the angle btw. two lines { std::stringstream buffer; //buffer << angle(vec.Get_vnl_vector(), -lastVec.Get_vnl_vector())*180/vnl_math::pi << "�"; buffer << angle(vec.GetVnlVector(), -lastVec.GetVnlVector())*180/vnl_math::pi << (char)176; Vector2D vec2d = pt2d-lastPt2d; vec2d.Normalize(); Vector2D lastVec2d = lastPt2d-preLastPt2d; lastVec2d.Normalize(); vec2d=vec2d-lastVec2d; vec2d.Normalize(); Vector2D pos2d = lastPt2d.GetVectorFromOrigin()+vec2d*text2dDistance*text2dDistance; mitk::VtkPropRenderer* OpenGLrenderer = dynamic_cast( renderer ); OpenGLrenderer->WriteSimpleText(buffer.str(), pos2d[0], pos2d[1]); //this->WriteTextXY(pos2d[0], pos2d[1], buffer.str(),renderer); } } counter++; } ++it; if(selIt != selEnd && !pointDataBroken) ++selIt; j++; } //recall the color to the same color before this drawing glColor3f(recallColor[0],recallColor[1],recallColor[2]); } } void mitk::PointSetGLMapper2D::SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer, bool overwrite) { node->AddProperty( "line width", mitk::IntProperty::New(2), renderer, overwrite ); // width of the line from one point to another node->AddProperty( "point line width", mitk::IntProperty::New(1), renderer, overwrite ); //width of the cross marking a point - node->AddProperty( "point 2D size", mitk::IntProperty::New(8), renderer, overwrite ); // length of the cross marking a point // length of an edge of the box marking a point + node->AddProperty( "point 2D size", mitk::FloatProperty::New(8), renderer, overwrite ); // length of the cross marking a point // length of an edge of the box marking a point node->AddProperty( "show contour", mitk::BoolProperty::New(false), renderer, overwrite ); // contour of the line between points node->AddProperty( "close contour", mitk::BoolProperty::New(false), renderer, overwrite ); node->AddProperty( "show points", mitk::BoolProperty::New(true), renderer, overwrite ); //show or hide points node->AddProperty( "show distances", mitk::BoolProperty::New(false), renderer, overwrite ); //show or hide distance measure (not always available) node->AddProperty( "distance decimal digits", mitk::IntProperty::New(2), renderer, overwrite ); //set the number of decimal digits to be shown node->AddProperty( "show angles", mitk::BoolProperty::New(false), renderer, overwrite ); //show or hide angle measurement (not always available) node->AddProperty( "show distant lines", mitk::BoolProperty::New(false), renderer, overwrite ); //show the line between to points from a distant view (equals "always on top" option) node->AddProperty( "layer", mitk::IntProperty::New(1), renderer, overwrite ); // default to draw pointset above images (they have a default layer of 0) Superclass::SetDefaultProperties(node, renderer, overwrite); } diff --git a/Modules/LegacyGL/mitkPointSetGLMapper2D.h b/Modules/LegacyGL/mitkPointSetGLMapper2D.h index 9a2dc8f7a5..bf71f47102 100644 --- a/Modules/LegacyGL/mitkPointSetGLMapper2D.h +++ b/Modules/LegacyGL/mitkPointSetGLMapper2D.h @@ -1,99 +1,99 @@ /*=================================================================== 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 MITKPointSetMAPPER2D_H_HEADER_INCLUDED #define MITKPointSetMAPPER2D_H_HEADER_INCLUDED #include #include "mitkGLMapper.h" namespace mitk { class BaseRenderer; class PointSet; /** * @brief OpenGL-based mapper to display a mitk::PointSet in a 2D window. * * This mapper can actually more than just draw a number of points of a * mitk::PointSet. If you set the right properties of the mitk::DataNode, * which contains the point set, then this mapper will also draw lines * connecting the points, and calculate and display distances and angles * between adjacent points. Here is a complete list of boolean properties, * which might be of interest: * * - \b "show contour": Draw not only the points but also the connections between * them (default false) * - \b "line width": IntProperty which gives the width of the contour lines * - \b "show points": Wheter or not to draw the actual points (default true) * - \b "show distances": Wheter or not to calculate and print the distance * between adjacent points (default false) * - \b "show angles": Wheter or not to calculate and print the angle between * adjacent points (default false) * - \b "show distant lines": When true, the mapper will also draw contour * lines that are far away form the current slice (default true) * - \b "label": StringProperty with a label for this point set * * BUG 1321 - possible new features: * point-2d-size (length of lines in cross/diamond) * point-linewidth * * @ingroup Mapper */ /** \deprecatedSince{2013_06} This mapper is replaced by PointSetVtkMapper2D. The child classes of this class are deprecated. * To further ensure their functionality PointSetGLMapper2D cannot be removed and is set deprecated too. */ class MITKLEGACYGL_EXPORT PointSetGLMapper2D : public GLMapper { public: mitkClassMacro(PointSetGLMapper2D, GLMapper); itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** @brief Get the PointDataList to map */ virtual const mitk::PointSet * GetInput(void); virtual void Paint(mitk::BaseRenderer * renderer) override; virtual void ApplyAllProperties(mitk::BaseRenderer* renderer); static void SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer = NULL, bool overwrite = false); protected: PointSetGLMapper2D(); virtual ~PointSetGLMapper2D(); bool m_Polygon; bool m_PolygonClosed; bool m_ShowPoints; bool m_ShowDistances; int m_DistancesDecimalDigits; bool m_ShowAngles; bool m_ShowDistantLines; int m_LineWidth; int m_PointLineWidth; - int m_Point2DSize; + float m_Point2DSize; }; } // namespace mitk #endif /* MITKPointSetMapper2D_H_HEADER_INCLUDED */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.cpp index 992e0c577e..7aa099b1d7 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.cpp @@ -1,329 +1,364 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkMLBTView.h" #include "QmitkStdMultiWidget.h" // Qt #include #include #include #include "mitkNodePredicateDataType.h" #include #include #include #include #include +#include +#include +#include #define _USE_MATH_DEFINES #include const std::string QmitkMLBTView::VIEW_ID = "org.mitk.views.mlbtview"; using namespace berry; QmitkMLBTView::QmitkMLBTView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) { m_TrackingTimer = new QTimer(this); } // Destructor QmitkMLBTView::~QmitkMLBTView() { } void QmitkMLBTView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkMLBTViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_StartTrainingButton, SIGNAL ( clicked() ), this, SLOT( StartTrainingThread() ) ); connect( &m_TrainingWatcher, SIGNAL ( finished() ), this, SLOT( OnTrainingThreadStop() ) ); connect( m_Controls->m_StartTrackingButton, SIGNAL ( clicked() ), this, SLOT( StartTrackingThread() ) ); connect( &m_TrackingWatcher, SIGNAL ( finished() ), this, SLOT( OnTrackingThreadStop() ) ); connect( m_Controls->m_SaveForestButton, SIGNAL ( clicked() ), this, SLOT( SaveForest() ) ); connect( m_Controls->m_LoadForestButton, SIGNAL ( clicked() ), this, SLOT( LoadForest() ) ); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(BuildFibers()) ); connect( m_Controls->m_TimerIntervalBox, SIGNAL(valueChanged(int)), this, SLOT( ChangeTimerInterval(int) )); connect( m_Controls->m_DemoModeBox, SIGNAL(stateChanged(int)), this, SLOT( ToggleDemoMode(int) )); connect( m_Controls->m_PauseTrackingButton, SIGNAL ( clicked() ), this, SLOT( PauseTracking() ) ); connect( m_Controls->m_AbortTrackingButton, SIGNAL ( clicked() ), this, SLOT( AbortTracking() ) ); int numThread = itk::MultiThreader::GetGlobalDefaultNumberOfThreads(); m_Controls->m_NumberOfThreadsBox->setMaximum(numThread); m_Controls->m_NumberOfThreadsBox->setValue(numThread); m_Controls->m_TrackingMaskImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TrackingSeedImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TrackingStopImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TrackingRawImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isQbi = mitk::NodePredicateDataType::New("QBallImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isQbi); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); finalPredicate = mitk::NodePredicateAnd::New(finalPredicate, isBinaryPredicate); m_Controls->m_TrackingMaskImageBox->SetPredicate(finalPredicate); m_Controls->m_TrackingSeedImageBox->SetPredicate(finalPredicate); m_Controls->m_TrackingStopImageBox->SetPredicate(finalPredicate); - m_Controls->m_TrackingRawImageBox->SetPredicate(isDwi); + m_Controls->m_TrackingRawImageBox->SetPredicate(isMitkImage); } } void QmitkMLBTView::AbortTracking() { if (tracker.IsNotNull()) { tracker->m_AbortTracking = true; } } void QmitkMLBTView::PauseTracking() { if (tracker.IsNotNull()) { tracker->m_PauseTracking = !tracker->m_PauseTracking; } } void QmitkMLBTView::ChangeTimerInterval(int value) { m_TrackingTimer->setInterval(value); } void QmitkMLBTView::ToggleDemoMode(int state) { if (tracker.IsNotNull()) { tracker->SetDemoMode(m_Controls->m_DemoModeBox->isChecked()); tracker->m_Stop = false; } } void QmitkMLBTView::BuildFibers() { if (m_Controls->m_DemoModeBox->isChecked() && tracker.IsNotNull() && tracker->m_BuildFibersFinished) { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); - outFib->SetFiberColors(255,0,0); + outFib->SetFiberColors(255,255,255); m_TractogramNode->SetData(outFib); + + m_SamplingPointsNode->SetData(tracker->m_SamplingPointset); + m_AlternativePointsNode->SetData(tracker->m_AlternativePointset); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); tracker->m_BuildFibersFinished = false; - tracker->m_Stop = false; tracker->m_BuildFibersReady = 0; + tracker->m_Stop = false; } } void QmitkMLBTView::LoadForest() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Forest"), QDir::currentPath(), tr("HDF5 random forest file (*.rf)") ); if(filename.isEmpty() || filename.isNull()) return; m_ForestHandler.LoadForest( filename.toStdString() ); } void QmitkMLBTView::StartTrackingThread() { m_TractogramNode = mitk::DataNode::New(); m_TractogramNode->SetName("MLBT Result"); m_TractogramNode->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(20)); m_TractogramNode->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); m_TractogramNode->SetProperty("LineWidth", mitk::IntProperty::New(2)); m_TractogramNode->SetProperty("color",mitk::ColorProperty::New(0, 1, 1)); this->GetDataStorage()->Add(m_TractogramNode); + m_SamplingPointsNode = mitk::DataNode::New(); + m_SamplingPointsNode->SetName("SamplingPoints"); + m_SamplingPointsNode->SetProperty("pointsize", mitk::FloatProperty::New(0.2)); + m_SamplingPointsNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); + mitk::PointSetShapeProperty::Pointer bla = mitk::PointSetShapeProperty::New(); bla->SetValue(8); + m_SamplingPointsNode->SetProperty("Pointset.2D.shape", bla); + m_SamplingPointsNode->SetProperty("Pointset.2D.distance to plane", mitk::FloatProperty::New(1.5)); + m_SamplingPointsNode->SetProperty("point 2D size", mitk::FloatProperty::New(0.1)); + m_SamplingPointsNode->SetProperty("Pointset.2D.fill shape", mitk::BoolProperty::New(true)); + this->GetDataStorage()->Add(m_SamplingPointsNode); + + m_AlternativePointsNode = mitk::DataNode::New(); + m_AlternativePointsNode->SetName("AlternativePoints"); + m_AlternativePointsNode->SetProperty("pointsize", mitk::FloatProperty::New(0.2)); + m_AlternativePointsNode->SetProperty("color", mitk::ColorProperty::New(1,0,0)); + m_AlternativePointsNode->SetProperty("Pointset.2D.shape", bla); + m_AlternativePointsNode->SetProperty("Pointset.2D.distance to plane", mitk::FloatProperty::New(1.5)); + m_AlternativePointsNode->SetProperty("point 2D size", mitk::FloatProperty::New(0.1)); + m_AlternativePointsNode->SetProperty("Pointset.2D.fill shape", mitk::BoolProperty::New(true)); + this->GetDataStorage()->Add(m_AlternativePointsNode); + QFuture future = QtConcurrent::run( this, &QmitkMLBTView::StartTracking ); m_TrackingWatcher.setFuture(future); m_TrackingThreadIsRunning = true; m_Controls->m_StartTrackingButton->setEnabled(false); - m_TrackingTimer->start(10); + m_TrackingTimer->start(m_Controls->m_TimerIntervalBox->value()); } void QmitkMLBTView::OnTrackingThreadStop() { m_TrackingThreadIsRunning = false; m_Controls->m_StartTrackingButton->setEnabled(true); - vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); - outFib->SetFiberColors(255,0,0); + outFib->SetFiberColors(255,255,255); // mitk::DataNode::Pointer node = mitk::DataNode::New(); m_TractogramNode->SetData(outFib); + m_SamplingPointsNode->SetData(tracker->m_SamplingPointset); + m_AlternativePointsNode->SetData(tracker->m_AlternativePointset); tracker = NULL; m_TrackingTimer->stop(); + + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkMLBTView::StartTracking() { if ( m_Controls->m_TrackingRawImageBox->GetSelectedNode().IsNull() ) return; mitk::Image::Pointer dwi = dynamic_cast(m_Controls->m_TrackingRawImageBox->GetSelectedNode()->GetData()); tracker = TrackerType::New(); tracker->SetNumberOfThreads(m_Controls->m_NumberOfThreadsBox->value()); tracker->SetInput(0, mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi) ); tracker->SetGradientDirections( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi) ); tracker->SetB_Value( mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi) ); tracker->SetDemoMode(m_Controls->m_DemoModeBox->isChecked()); + if (m_Controls->m_DemoModeBox->isChecked()) + tracker->SetNumberOfThreads(1); if (m_Controls->m_TrackingUseMaskImageBox->isChecked() && m_Controls->m_TrackingMaskImageBox->GetSelectedNode().IsNotNull()) { mitk::Image::Pointer mask = dynamic_cast(m_Controls->m_TrackingMaskImageBox->GetSelectedNode()->GetData()); ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); mitk::CastToItkImage(mask, itkMask); tracker->SetMaskImage(itkMask); } if (m_Controls->m_TrackingUseSeedImageBox->isChecked() && m_Controls->m_TrackingSeedImageBox->GetSelectedNode().IsNotNull()) { mitk::Image::Pointer img = dynamic_cast(m_Controls->m_TrackingSeedImageBox->GetSelectedNode()->GetData()); ItkUcharImgType::Pointer itkImg = ItkUcharImgType::New(); mitk::CastToItkImage(img, itkImg); tracker->SetSeedImage(itkImg); } if (m_Controls->m_TrackingUseStopImageBox->isChecked() && m_Controls->m_TrackingStopImageBox->GetSelectedNode().IsNotNull()) { mitk::Image::Pointer img = dynamic_cast(m_Controls->m_TrackingStopImageBox->GetSelectedNode()->GetData()); ItkUcharImgType::Pointer itkImg = ItkUcharImgType::New(); mitk::CastToItkImage(img, itkImg); tracker->SetStoppingRegions(itkImg); } tracker->SetSeedsPerVoxel(m_Controls->m_NumberOfSeedsBox->value()); - tracker->SetUseDirection(true); tracker->SetStepSize(m_Controls->m_TrackingStepSizeBox->value()); tracker->SetAngularThreshold(cos((float)m_Controls->m_AngularThresholdBox->value()*M_PI/180)); tracker->SetMinTractLength(m_Controls->m_MinLengthBox->value()); tracker->SetMaxTractLength(m_Controls->m_MaxLengthBox->value()); + tracker->SetAposterioriCurvCheck(m_Controls->m_Curvcheck2->isChecked()); + tracker->SetRemoveWmEndFibers(false); + tracker->SetAvoidStop(m_Controls->m_AvoidStop->isChecked()); vigra::RandomForest forest = m_ForestHandler.GetForest(); tracker->SetDecisionForest(&forest); tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); + tracker->SetRandomSampling(m_Controls->m_RandomSampling->isChecked()); tracker->Update(); // vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); // mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); // outFib->SetColorCoding(mitk::FiberBundle::COLORCODING_CUSTOM); // mitk::DataNode::Pointer node = mitk::DataNode::New(); // m_TractogramNode->SetData(outFib); // node->SetData(outFib); // node->SetName("MLBT Result"); // this->GetDataStorage()->Add(node); // mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkMLBTView::SaveForest() { QString filename = QFileDialog::getSaveFileName(0, tr("Save Forest"), QDir::currentPath()+"/forest.rf", tr("HDF5 random forest file (*.rf)") ); if(filename.isEmpty() || filename.isNull()) return; if(!filename.endsWith(".rf")) filename += ".rf"; m_ForestHandler.SaveForest( filename.toStdString() ); } void QmitkMLBTView::StartTrainingThread() { QFuture future = QtConcurrent::run( this, &QmitkMLBTView::StartTraining ); m_TrainingWatcher.setFuture(future); m_Controls->m_StartTrainingButton->setEnabled(false); m_Controls->m_SaveForestButton->setEnabled(false); m_Controls->m_LoadForestButton->setEnabled(false); } void QmitkMLBTView::OnTrainingThreadStop() { m_Controls->m_StartTrainingButton->setEnabled(true); m_Controls->m_SaveForestButton->setEnabled(true); m_Controls->m_LoadForestButton->setEnabled(true); } void QmitkMLBTView::StartTraining() { m_ForestHandler.SetRawData(m_SelectedDiffImages); m_ForestHandler.SetTractograms(m_SelectedFB); m_ForestHandler.SetNumTrees(m_Controls->m_NumTreesBox->value()); m_ForestHandler.SetMaxTreeDepth(m_Controls->m_MaxDepthBox->value()); m_ForestHandler.SetGrayMatterSamplesPerVoxel(m_Controls->m_GmSamplingBox->value()); m_ForestHandler.SetSampleFraction(m_Controls->m_SampleFractionBox->value()); m_ForestHandler.SetStepSize(m_Controls->m_TrainingStepSizeBox->value()); - m_ForestHandler.SetUsePreviousDirection(m_Controls->m_TrackingUsePreviousDirectionBox->isChecked()); m_ForestHandler.StartTraining(); } void QmitkMLBTView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkMLBTView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkMLBTView::OnSelectionChanged( std::vector nodes ) { if ( !this->IsVisible() ) { // do nothing if nobody wants to see me :-( return; } m_SelectedFB.clear(); m_SelectedDiffImages.clear(); m_MaskImages.clear(); m_WhiteMatterImages.clear(); for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( dynamic_cast(node->GetData()) ) m_SelectedFB.push_back( dynamic_cast(node->GetData()) ); else if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(node)) m_SelectedDiffImages.push_back( dynamic_cast(node->GetData()) ); } } void QmitkMLBTView::Activated() { } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.h index fa566e1d40..f63b256870 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTView.h @@ -1,105 +1,107 @@ /*=================================================================== 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 "ui_QmitkMLBTViewControls.h" #ifndef Q_MOC_RUN #include "mitkDataStorage.h" #include "mitkDataStorageSelection.h" #include #include #endif #include #include #include /*! \brief */ // Forward Qt class declarations class QmitkMLBTView : public QmitkFunctionality { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; typedef itk::Image ItkUcharImgType; typedef itk::MLBSTrackingFilter<100> TrackerType; QmitkMLBTView(); virtual ~QmitkMLBTView(); virtual void CreateQtPartControl(QWidget *parent) override; virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) override; virtual void StdMultiWidgetNotAvailable() override; virtual void Activated() override; protected slots: void StartTrackingThread(); void OnTrackingThreadStop(); void StartTrainingThread(); void OnTrainingThreadStop(); void SaveForest(); void LoadForest(); void BuildFibers(); void ChangeTimerInterval(int value); void ToggleDemoMode(int state); void PauseTracking(); void AbortTracking(); protected: void StartTracking(); void StartTraining(); /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ) override; Ui::QmitkMLBTViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; mitk::TrackingForestHandler<> m_ForestHandler; std::vector< mitk::Image::Pointer > m_SelectedDiffImages; std::vector< mitk::FiberBundle::Pointer > m_SelectedFB; std::vector< ItkUcharImgType::Pointer > m_MaskImages; std::vector< ItkUcharImgType::Pointer > m_WhiteMatterImages; QFutureWatcher m_TrainingWatcher; QFutureWatcher m_TrackingWatcher; bool m_TrackingThreadIsRunning; TrackerType::Pointer tracker; QTimer* m_TrackingTimer; mitk::DataNode::Pointer m_TractogramNode; + mitk::DataNode::Pointer m_SamplingPointsNode; + mitk::DataNode::Pointer m_AlternativePointsNode; private: }; diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTViewControls.ui index cf19c52391..bf1799b4e8 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkMLBTViewControls.ui @@ -1,620 +1,623 @@ QmitkMLBTViewControls 0 0 359 1127 Form 0 0 0 0 1 0 0 359 1065 Training Load Forest Qt::Vertical 20 40 QFrame::NoFrame QFrame::Raised 0 0 0 0 - - - - Step size: - - - - - - - Num. Trees: - - - - - + + - Max. Depth: + GM Sampling Points: - - + + 3 - - -1.000000000000000 - - 999.000000000000000 + 1.000000000000000 0.100000000000000 - -1.000000000000000 + 1.000000000000000 1 999999999 - - + + - Sample Fraction: + Step size: 1 999999999 10 + + + + Num. Trees: + + + + + + + Max. Depth: + + + + + + + Sample Fraction: + + + 1 999999999 10 - - + + 3 + + -1.000000000000000 + - 1.000000000000000 + 999.000000000000000 0.100000000000000 - 1.000000000000000 - - - - - - - GM Sampling Points: - - - - - - - Use Previous Direction: - - - - - - - - - - true + -1.000000000000000 Start Training Save Forest 0 0 359 1065 Tractography - - + + - Use Previous Direction + Avoid premature termination true + + + + Qt::Vertical + + + + 20 + 40 + + + + - + + + Secondary curvature check + + + true + + + + + QFrame::NoFrame QFrame::Raised - + 0 0 0 0 - - - - ... - - - - :/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-pause.svg:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-pause.svg - - - - + - ... - - - - :/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-start.svg:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-start.svg + Demo Mode - - - - ... + + + + 1 - - - :/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-stop.svg:/org.mbi.gui.qt.diffusionimaginginternal/resources/Media-playback-stop.svg + + 1000 + + + 10 - - + + QFrame::NoFrame QFrame::Raised - + 0 0 0 0 - - + + - Demo Mode + Use Stop Image: - - - - 1 + + + + Use Seed Image: - - 1000 + + false - - 10 + + + + + + Use Mask Image: + + + false + + + + + + + + + QFrame::NoFrame QFrame::Raised 0 0 0 0 0.500000000000000 Max. Length 999999999 50 Num. Samples: Input DWI: Step Size: Sampling Distance: 1 999 Min. Length 0.100000000000000 0.500000000000000 999999999.000000000000000 1.000000000000000 20.000000000000000 Angular Threshold: Num. Seeds: 999999999.000000000000000 1.000000000000000 400.000000000000000 1 90.000000000000000 45.000000000000000 Num. Threads: 1 30 - - - - Qt::Vertical - - - - 20 - 40 - - - - - - + + QFrame::NoFrame QFrame::Raised - + 0 0 0 0 - - + + - Use Stop Image: + ... + + + + :/org_mitk_icons/icons/tango/scalable/actions/media-playback-pause.svg:/org_mitk_icons/icons/tango/scalable/actions/media-playback-pause.svg - - + + - Use Seed Image: + ... - - false + + + :/org_mitk_icons/icons/tango/scalable/actions/media-playback-start.svg:/org_mitk_icons/icons/tango/scalable/actions/media-playback-start.svg - - + + - Use Mask Image: + ... - - true + + + :/org_mitk_icons/icons/tango/scalable/actions/media-playback-stop.svg:/org_mitk_icons/icons/tango/scalable/actions/media-playback-stop.svg - - - - - - - - - + + + + Random sampling + + + false + + + QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 - +