diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp index 23ce7e33ba..c9e01b2b45 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,1672 +1,1738 @@ /*=================================================================== 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. ===================================================================*/ #define _USE_MATH_DEFINES #include "mitkFiberBundleX.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundleX::COLORCODING_ORIENTATION_BASED = "Color_Orient"; //const char* mitk::FiberBundleX::COLORCODING_FA_AS_OPACITY = "Color_Orient_FA_Opacity"; const char* mitk::FiberBundleX::COLORCODING_FA_BASED = "FA_Values"; const char* mitk::FiberBundleX::COLORCODING_CUSTOM = "custom"; const char* mitk::FiberBundleX::FIBER_ID_ARRAY = "Fiber_IDs"; using namespace std; mitk::FiberBundleX::FiberBundleX( vtkPolyData* fiberPolyData ) : m_CurrentColorCoding(NULL) , m_NumFibers(0) { m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != NULL) { m_FiberPolyData = fiberPolyData; //m_FiberPolyData->DeepCopy(fiberPolyData); this->DoColorCodingOrientationBased(); } m_NumFibers = m_FiberPolyData->GetNumberOfLines(); this->UpdateFiberGeometry(); this->SetColorCoding(COLORCODING_ORIENTATION_BASED); this->GenerateFiberIds(); } mitk::FiberBundleX::~FiberBundleX() { } mitk::FiberBundleX::Pointer mitk::FiberBundleX::GetDeepCopy() { mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(m_FiberPolyData); if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)) MITK_DEBUG << "ok"; vtkUnsignedCharArray* tmpColors = (vtkUnsignedCharArray*) m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED); int tmpColorss = tmpColors->GetNumberOfTuples(); int tmpColorc = tmpColors->GetNumberOfComponents(); newFib->SetColorCoding(m_CurrentColorCoding); return newFib; } vtkSmartPointer mitk::FiberBundleX::GeneratePolyDataByIds(std::vector fiberIds) { MITK_DEBUG << "\n=====FINAL RESULT: fib_id ======\n"; MITK_DEBUG << "Number of new Fibers: " << fiberIds.size(); // iterate through the vectorcontainer hosting all desired fiber Ids vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); vtkSmartPointer newLineSet = vtkSmartPointer::New(); vtkSmartPointer newPointSet = vtkSmartPointer::New(); // if FA array available, initialize fa double array // if color orient array is available init color array vtkSmartPointer faValueArray; vtkSmartPointer colorsT; //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = sizeof(rgba); if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_FA_BASED)){ MITK_DEBUG << "FA VALUES AVAILABLE, init array for new fiberbundle"; faValueArray = vtkSmartPointer::New(); } if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)){ MITK_DEBUG << "colorValues available, init array for new fiberbundle"; colorsT = vtkUnsignedCharArray::New(); colorsT->SetNumberOfComponents(componentSize); colorsT->SetName(COLORCODING_ORIENTATION_BASED); } std::vector::iterator finIt = fiberIds.begin(); while ( finIt != fiberIds.end() ) { if (*finIt < 0 || *finIt>GetNumFibers()){ MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; break; } vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); vtkSmartPointer fibPoints = fiber->GetPoints(); vtkSmartPointer newFiber = vtkSmartPointer::New(); newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); for(int i=0; iGetNumberOfPoints(); i++) { // MITK_DEBUG << "id: " << fiber->GetPointId(i); // MITK_DEBUG << fibPoints->GetPoint(i)[0] << " | " << fibPoints->GetPoint(i)[1] << " | " << fibPoints->GetPoint(i)[2]; newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_FA_BASED)){ // MITK_DEBUG << m_FiberIdDataSet->GetPointData()->GetArray(FA_VALUE_ARRAY)->GetTuple(fiber->GetPointId(i)); } if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)){ // MITK_DEBUG << "ColorValue: " << m_FiberIdDataSet->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetTuple(fiber->GetPointId(i))[0]; } } newLineSet->InsertNextCell(newFiber); ++finIt; } newFiberPolyData->SetPoints(newPointSet); newFiberPolyData->SetLines(newLineSet); MITK_DEBUG << "new fiberbundle polydata points: " << newFiberPolyData->GetNumberOfPoints(); MITK_DEBUG << "new fiberbundle polydata lines: " << newFiberPolyData->GetNumberOfLines(); MITK_DEBUG << "=====================\n"; // mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(newFiberPolyData); return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundleX::Pointer mitk::FiberBundleX::AddBundle(mitk::FiberBundleX* fib) { if (fib==NULL) { MITK_WARN << "trying to call AddBundle with NULL argument"; return NULL; } vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); // add current fiber bundle int numFibers = GetNumFibers(); for( int i=0; iGetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } vLines = fib->m_FiberPolyData->GetLines(); vLines->InitTraversal(); // add new fiber bundle numFibers = fib->GetNumFibers(); for( int i=0; iGetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(fib->m_FiberPolyData->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); return newFib; } // subtract two fiber bundles mitk::FiberBundleX::Pointer mitk::FiberBundleX::SubtractBundle(mitk::FiberBundleX* fib) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); // iterate over current fibers int numFibers = GetNumFibers(); for( int i=0; iGetNextCell ( numPoints, points ); if (points==NULL) continue; vtkSmartPointer vLines2 = fib->m_FiberPolyData->GetLines(); vLines2->InitTraversal(); int numFibers2 = fib->GetNumFibers(); bool contained = false; for( int i2=0; i2GetNextCell ( numPoints2, points2 ); if (points2==NULL) continue; // check endpoints itk::Point point_start = GetItkPoint(m_FiberPolyData->GetPoint(points[0])); itk::Point point_end = GetItkPoint(m_FiberPolyData->GetPoint(points[numPoints-1])); itk::Point point2_start = GetItkPoint(fib->m_FiberPolyData->GetPoint(points2[0])); itk::Point point2_end = GetItkPoint(fib->m_FiberPolyData->GetPoint(points2[numPoints2-1])); if (point_start.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps || point_start.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps) { // further checking ??? if (numPoints2==numPoints) contained = true; } } // add to result because fiber is not subtracted if (!contained) { vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } } if(vNewLines->GetNumberOfCells()==0) return NULL; // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); return newFib; } itk::Point mitk::FiberBundleX::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } /* * set polydata (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundleX::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { if (fiberPD == NULL) this->m_FiberPolyData = vtkSmartPointer::New(); else { m_FiberPolyData->DeepCopy(fiberPD); DoColorCodingOrientationBased(); } m_NumFibers = m_FiberPolyData->GetNumberOfLines(); if (updateGeometry) UpdateFiberGeometry(); SetColorCoding(COLORCODING_ORIENTATION_BASED); GenerateFiberIds(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundleX::GetFiberPolyData() { return m_FiberPolyData; } void mitk::FiberBundleX::DoColorCodingOrientationBased() { //===== FOR WRITING A TEST ======================== // colorT size == tupelComponents * tupelElements // compare color results // to cover this code 100% also polydata needed, where colorarray already exists // + one fiber with exactly 1 point // + one fiber with 0 points //================================================= /* make sure that processing colorcoding is only called when necessary */ if ( m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) && m_FiberPolyData->GetNumberOfPoints() == m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetNumberOfTuples() ) { // fiberstructure is already colorcoded MITK_DEBUG << " NO NEED TO REGENERATE COLORCODING! " ; this->ResetFiberOpacity(); this->SetColorCoding(COLORCODING_ORIENTATION_BASED); return; } /* Finally, execute color calculation */ vtkPoints* extrPoints = NULL; extrPoints = m_FiberPolyData->GetPoints(); int numOfPoints = 0; if (extrPoints!=NULL) numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; // int componentSize = sizeof(rgba); int componentSize = 4; vtkSmartPointer colorsT = vtkUnsignedCharArray::New(); colorsT->Allocate(numOfPoints * componentSize); colorsT->SetNumberOfComponents(componentSize); colorsT->SetName(COLORCODING_ORIENTATION_BASED); /* checkpoint: does polydata contain any fibers */ int numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) { MITK_DEBUG << "\n ========= Number of Fibers is 0 and below ========= \n"; return; } /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); // MITK_DEBUG << "Fib#: " << fi << " of " << numOfFibers << " pnts in fiber: " << pointsPerFiber ; /* single fiber checkpoints: is number of points valid */ if (pointsPerFiber > 1) { /* operate on points of single fiber */ for (int i=0; i 0) { /* The color value of the current point is influenced by the previous point and next point. */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; vnl_vector_fixed< double, 3 > diff; diff = (diff1 - diff2) / 2.0; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if (i==0) { /* First point has no previous point, therefore only diff1 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; diff1.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); rgba[3] = (unsigned char) (255.0); } else if (i==pointsPerFiber-1) { /* Last point has no next point, therefore only diff2 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; diff2.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); rgba[3] = (unsigned char) (255.0); } colorsT->InsertTupleValue(idList[i], rgba); } //end for loop } else if (pointsPerFiber == 1) { /* a single point does not define a fiber (use vertex mechanisms instead */ continue; // colorsT->InsertTupleValue(0, rgba); } else { MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } }//end for loop m_FiberPolyData->GetPointData()->AddArray(colorsT); /*========================= - this is more relevant for renderer than for fiberbundleX datastructure - think about sourcing this to a explicit method which coordinates colorcoding */ this->SetColorCoding(COLORCODING_ORIENTATION_BASED); // =========================== //mini test, shall be ported to MITK TESTINGS! if (colorsT->GetSize() != numOfPoints*componentSize) MITK_DEBUG << "ALLOCATION ERROR IN INITIATING COLOR ARRAY"; } void mitk::FiberBundleX::DoColorCodingFaBased() { if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED) != 1 ) return; this->SetColorCoding(COLORCODING_FA_BASED); MITK_DEBUG << "FBX: done CC FA based"; this->GenerateFiberIds(); } void mitk::FiberBundleX::DoUseFaFiberOpacity() { if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED) != 1 ) return; if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) != 1 ) return; vtkDoubleArray* FAValArray = (vtkDoubleArray*) m_FiberPolyData->GetPointData()->GetArray(COLORCODING_FA_BASED); vtkUnsignedCharArray* ColorArray = dynamic_cast (m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)); for(long i=0; iGetNumberOfTuples(); i++) { double faValue = FAValArray->GetValue(i); faValue = faValue * 255.0; ColorArray->SetComponent(i,3, (unsigned char) faValue ); } this->SetColorCoding(COLORCODING_ORIENTATION_BASED); MITK_DEBUG << "FBX: done CC OPACITY"; this->GenerateFiberIds(); } void mitk::FiberBundleX::ResetFiberOpacity() { vtkUnsignedCharArray* ColorArray = dynamic_cast (m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)); if (ColorArray==NULL) return; for(long i=0; iGetNumberOfTuples(); i++) ColorArray->SetComponent(i,3, 255.0 ); } void mitk::FiberBundleX::SetFAMap(mitk::Image::Pointer FAimage) { MITK_DEBUG << "SetFAMap"; vtkSmartPointer faValues = vtkDoubleArray::New(); faValues->SetName(COLORCODING_FA_BASED); faValues->Allocate(m_FiberPolyData->GetNumberOfPoints()); // MITK_DEBUG << faValues->GetNumberOfTuples(); // MITK_DEBUG << faValues->GetSize(); faValues->SetNumberOfValues(m_FiberPolyData->GetNumberOfPoints()); // MITK_DEBUG << faValues->GetNumberOfTuples(); // MITK_DEBUG << faValues->GetSize(); vtkPoints* pointSet = m_FiberPolyData->GetPoints(); for(long i=0; iGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; double faPixelValue = 1-FAimage->GetPixelValueByWorldCoordinate(px); // faValues->InsertNextTuple1(faPixelValue); faValues->InsertValue(i, faPixelValue); // MITK_DEBUG << faPixelValue; // MITK_DEBUG << faValues->GetValue(i); } m_FiberPolyData->GetPointData()->AddArray(faValues); this->GenerateFiberIds(); if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_FA_BASED)) MITK_DEBUG << "FA VALUE ARRAY SET"; // vtkDoubleArray* valueArray = (vtkDoubleArray*) m_FiberPolyData->GetPointData()->GetArray(FA_VALUE_ARRAY); // for(long i=0; iGetNumberOfPoints(); i++) // { // MITK_DEBUG << "value at pos "<< i << ": " << valueArray->GetValue(i); // } } void mitk::FiberBundleX::GenerateFiberIds() { if (m_FiberPolyData == NULL) return; vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); idFiberFilter->SetInput(m_FiberPolyData); idFiberFilter->CellIdsOn(); // idFiberFilter->PointIdsOn(); // point id's are not needed idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); idFiberFilter->FieldDataOn(); idFiberFilter->Update(); m_FiberIdDataSet = idFiberFilter->GetOutput(); MITK_DEBUG << "Generating Fiber Ids...[done] | " << m_FiberIdDataSet->GetNumberOfCells(); } std::vector mitk::FiberBundleX::CutFiberBundle(mitk::PlanarFigure *pf) { std::vector bundles; // Get all fibers that do not pass the roi. these can be used directly. mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); mitk::PlaneGeometry* currentGeometry2D = dynamic_cast( const_cast(pf->GetGeometry2D()) ); PFCNot->SetGeometry2D(currentGeometry2D); PFCNot->setOperationType(mitk::PFCOMPOSITION_NOT_OPERATION); PFCNot->addPlanarFigure( pf ); mitk::FiberBundleX::Pointer notInRoi = this->ExtractFiberSubset(PFCNot); // Now find the fibers that pass the roi and need some more attention and get a list of all points that lie on the ROI mitk::FiberBundleX::Pointer inRoi = this->ExtractFiberSubset(pf); int num = inRoi->GetNumFibers(); vtkSmartPointer fiberPolyData = inRoi->GetFiberPolyData(); vtkCellArray* lines = fiberPolyData->GetLines(); lines->InitTraversal(); // New fiber polyData for clipped and residual bundles vtkSmartPointer clippedPolyData = vtkSmartPointer::New(); vtkSmartPointer residualPolyData = vtkSmartPointer::New(); // true when we are currently moving on points that should be thrown away - bool currentlyOnResidual = false; + bool currentlyOnPositive = false; + + + // plane equation: Ax+By+Cz+D=0 + Vector3D normal = currentGeometry2D->GetNormal(); + double A = normal[0]; + double B = normal[1]; + double C = normal[2]; + double D = - (A + (B + C) ); + + + + + + // Create a cell array to store the lines in and add the lines to it + vtkSmartPointer cells = + vtkSmartPointer::New(); + + + vtkSmartPointer points = vtkPoints::New(); + int lastId = 0; // iterate through all the lines for( int fiberID( 0 ); fiberID < num; fiberID++ ) { vtkIdType numPointsInCell(0); vtkIdType* pointsInCell(NULL); lines->GetNextCell ( numPointsInCell, pointsInCell ); + int pointsSinceLastCrossing = 0; + + // iterate trough all the points on the line and check on which side of the plane defined by the planar figure they are for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) { // push back point double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); + Point3D point; + point[0] = p[0]; + point[1] = p[1]; + point[2] = p[2]; + + if( (A*p[0]+B*p[1]+C*p[2]+D > 0) && currentlyOnPositive || + (A*p[0]+B*p[1]+C*p[2]+D) <= 0 && !currentlyOnPositive ) + { + // we remain on the same side of the plane + points->InsertNextPoint(point); + pointsSinceLastCrossing++; + } + else + { + // we change to the other side of the plane + currentlyOnPositive = !currentlyOnPositive; + + + // add the line we created thus far + vtkSmartPointer polyLine = vtkPolyLine::New(); + polyLine->GetPointIds()->SetNumberOfIds( pointsSinceLastCrossing ); + while(lastId < points->GetNumberOfPoints()) + { + polyLine->GetPointIds()->SetId(i,i); + lastId+++; + } + + cells->InsertNextCell(polyLine); + + pointsSinceLastCrossing=0; + + } + } - //tracts.push_back(singleTract); + // Done with all points, so add remaining line + + vtkSmartPointer polyLine = vtkPolyLine::New(); + polyLine->GetPointIds()->SetNumberOfIds( pointsSinceLastCrossing ); + while(lastId < points->GetNumberOfPoints()) + { + polyLine->GetPointIds()->SetId(i,i); + lastId+++; + } + + } /* if (pf==NULL) return bundles; // Obtain the geometry info of the PlanarFigure roi mitk::Geometry2D::ConstPointer pfgeometry = pf->GetGeometry2D(); const mitk::PlaneGeometry* planeGeometry = dynamic_cast (pfgeometry.GetPointer()); Vector3D planeNormal = planeGeometry->GetNormal(); planeNormal.Normalize(); Point3D planeOrigin = planeGeometry->GetOrigin(); //Define cutting plane by ROI (PlanarFigure) vtkSmartPointer plane = vtkSmartPointer::New(); plane->SetOrigin(planeOrigin[0],planeOrigin[1],planeOrigin[2]); plane->SetNormal(planeNormal[0],planeNormal[1],planeNormal[2]); // Use the plane that was defined by the planar figure and cut the fiber bundle in two parts: // One part on the min-side of the plane and one on the plus-side of the plane MITK_DEBUG << "start clipping"; vtkSmartPointer clipper = vtkSmartPointer::New(); clipper->SetInput(m_FiberIdDataSet); clipper->SetClipFunction(plane); clipper->GenerateClipScalarsOn(); clipper->GenerateClippedOutputOn(); vtkSmartPointer minusSide = clipper->GetClippedOutput(); vtkSmartPointer plusSide = clipper->GetOutput(); mitk::FiberBundleX::Pointer minusBundle = mitk::FiberBundleX::New(minusSide); mitk::FiberBundleX::Pointer plusBundle = mitk::FiberBundleX::New(plusSide); // Now from the bundles on the min-side of the plane mitk::PlaneGeometry* currentGeometry2D = dynamic_cast( const_cast(pf->GetGeometry2D()) ); PFCNot->SetGeometry2D(currentGeometry2D); PFCNot->setOperationType(mitk::PFCOMPOSITION_NOT_OPERATION); PFCNot->addPlanarFigure( pf ); mitk::FiberBundleX::Pointer notRoi = minusBundle->ExtractFiberSubset(PFCNot); mitk::FiberBundleX::Pointer plusNotRoi = plusBundle->ExtractFiberSubset(PFCNot); mitk::FiberBundleX::Pointer finalResult = notRoi->AddBundle(plusBundle); bundles.push_back(plusBundle); bundles.push_back(minusBundle); bundles.push_back(notRoi); bundles.push_back(plusNotRoi); bundles.push_back(finalResult); */ return bundles; } mitk::FiberBundleX::Pointer mitk::FiberBundleX::ExtractFiberSubset(mitk::PlanarFigure* pf) { if (pf==NULL) return NULL; std::vector tmp = ExtractFiberIdSubset(pf); if (tmp.size()<=0) return mitk::FiberBundleX::New(); vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp); return mitk::FiberBundleX::New(pTmp); } std::vector mitk::FiberBundleX::ExtractFiberIdSubset(mitk::PlanarFigure* pf) { MITK_DEBUG << "Extracting fibers!"; // vector which is returned, contains all extracted FiberIds std::vector FibersInROI; if (pf==NULL) return FibersInROI; /* Handle type of planarfigure */ // if incoming pf is a pfc mitk::PlanarFigureComposite::Pointer pfcomp= dynamic_cast(pf); if (!pfcomp.IsNull()) { // process requested boolean operation of PFC switch (pfcomp->getOperationType()) { case 0: { MITK_DEBUG << "AND PROCESSING"; //AND //temporarly store results of the child in this vector, we need that to accumulate the std::vector childResults = this->ExtractFiberIdSubset(pfcomp->getChildAt(0)); MITK_DEBUG << "first roi got fibers in ROI: " << childResults.size(); MITK_DEBUG << "sorting..."; std::sort(childResults.begin(), childResults.end()); MITK_DEBUG << "sorting done"; std::vector AND_Assamblage(childResults.size()); //std::vector AND_Assamblage; fill(AND_Assamblage.begin(), AND_Assamblage.end(), -1); //AND_Assamblage.reserve(childResults.size()); //max size AND can reach anyway std::vector::iterator it; for (int i=1; igetNumberOfChildren(); ++i) { std::vector tmpChild = this->ExtractFiberIdSubset(pfcomp->getChildAt(i)); MITK_DEBUG << "ROI " << i << " has fibers in ROI: " << tmpChild.size(); sort(tmpChild.begin(), tmpChild.end()); it = std::set_intersection(childResults.begin(), childResults.end(), tmpChild.begin(), tmpChild.end(), AND_Assamblage.begin() ); } MITK_DEBUG << "resize Vector"; long i=0; while (i < AND_Assamblage.size() && AND_Assamblage[i] != -1){ //-1 represents a placeholder in the array ++i; } AND_Assamblage.resize(i); MITK_DEBUG << "returning AND vector, size: " << AND_Assamblage.size(); return AND_Assamblage; // break; } case 1: { //OR std::vector OR_Assamblage = this->ExtractFiberIdSubset(pfcomp->getChildAt(0)); std::vector::iterator it; MITK_DEBUG << OR_Assamblage.size(); for (int i=1; igetNumberOfChildren(); ++i) { it = OR_Assamblage.end(); std::vector tmpChild = this->ExtractFiberIdSubset(pfcomp->getChildAt(i)); OR_Assamblage.insert(it, tmpChild.begin(), tmpChild.end()); MITK_DEBUG << "ROI " << i << " has fibers in ROI: " << tmpChild.size() << " OR Assamblage: " << OR_Assamblage.size(); } sort(OR_Assamblage.begin(), OR_Assamblage.end()); it = unique(OR_Assamblage.begin(), OR_Assamblage.end()); OR_Assamblage.resize( it - OR_Assamblage.begin() ); MITK_DEBUG << "returning OR vector, size: " << OR_Assamblage.size(); return OR_Assamblage; } case 2: { //NOT //get IDs of all fibers std::vector childResults; childResults.reserve(this->GetNumFibers()); vtkSmartPointer idSet = m_FiberIdDataSet->GetCellData()->GetArray(FIBER_ID_ARRAY); MITK_DEBUG << "m_NumOfFib: " << this->GetNumFibers() << " cellIdNum: " << idSet->GetNumberOfTuples(); for(long i=0; iGetNumFibers(); i++) { MITK_DEBUG << "i: " << i << " idset: " << idSet->GetTuple(i)[0]; childResults.push_back(idSet->GetTuple(i)[0]); } std::sort(childResults.begin(), childResults.end()); std::vector NOT_Assamblage(childResults.size()); //fill it with -1, otherwise 0 will be stored and 0 can also be an ID of fiber! fill(NOT_Assamblage.begin(), NOT_Assamblage.end(), -1); std::vector::iterator it; for (long i=0; igetNumberOfChildren(); ++i) { std::vector tmpChild = ExtractFiberIdSubset(pfcomp->getChildAt(i)); sort(tmpChild.begin(), tmpChild.end()); it = std::set_difference(childResults.begin(), childResults.end(), tmpChild.begin(), tmpChild.end(), NOT_Assamblage.begin() ); } MITK_DEBUG << "resize Vector"; long i=0; while (NOT_Assamblage[i] != -1){ //-1 represents a placeholder in the array ++i; } NOT_Assamblage.resize(i); return NOT_Assamblage; } default: MITK_DEBUG << "we have an UNDEFINED composition... ERROR" ; break; } } else { mitk::Geometry2D::ConstPointer pfgeometry = pf->GetGeometry2D(); const mitk::PlaneGeometry* planeGeometry = dynamic_cast (pfgeometry.GetPointer()); Vector3D planeNormal = planeGeometry->GetNormal(); planeNormal.Normalize(); Point3D planeOrigin = planeGeometry->GetOrigin(); MITK_DEBUG << "planeOrigin: " << planeOrigin[0] << " | " << planeOrigin[1] << " | " << planeOrigin[2] << endl; MITK_DEBUG << "planeNormal: " << planeNormal[0] << " | " << planeNormal[1] << " | " << planeNormal[2] << endl; std::vector PointsOnPlane; // contains all pointIds which are crossing the cutting plane std::vector PointsInROI; // based on PointsOnPlane, all ROI relevant point IDs are stored here /* Define cutting plane by ROI (PlanarFigure) */ vtkSmartPointer plane = vtkSmartPointer::New(); plane->SetOrigin(planeOrigin[0],planeOrigin[1],planeOrigin[2]); plane->SetNormal(planeNormal[0],planeNormal[1],planeNormal[2]); //same plane but opposite normal direction. so point cloud will be reduced -> better performance // vtkSmartPointer planeR = vtkSmartPointer::New(); //define new origin along the normal but close to the original one // OriginNew = OriginOld + 1*Normal // Vector3D extendedNormal; // int multiplyFactor = 1; // extendedNormal[0] = planeNormal[0] * multiplyFactor; // extendedNormal[1] = planeNormal[1] * multiplyFactor; // extendedNormal[2] = planeNormal[2] * multiplyFactor; // Point3D RplaneOrigin = planeOrigin - extendedNormal; // planeR->SetOrigin(RplaneOrigin[0],RplaneOrigin[1],RplaneOrigin[2]); // planeR->SetNormal(-planeNormal[0],-planeNormal[1],-planeNormal[2]); // MITK_DEBUG << "RPlaneOrigin: " << RplaneOrigin[0] << " | " << RplaneOrigin[1] // << " | " << RplaneOrigin[2]; /* get all points/fibers cutting the plane */ MITK_DEBUG << "start clipping"; vtkSmartPointer clipper = vtkSmartPointer::New(); clipper->SetInput(m_FiberIdDataSet); clipper->SetClipFunction(plane); clipper->GenerateClipScalarsOn(); clipper->GenerateClippedOutputOn(); vtkSmartPointer clipperout = clipper->GetClippedOutput(); MITK_DEBUG << "end clipping"; /* for some reason clipperoutput is not initialized for futher processing * so far only writing out clipped polydata provides requested */ // MITK_DEBUG << "writing clipper output"; // vtkSmartPointer writerC = vtkSmartPointer::New(); // writerC->SetInput(clipperout1); // writerC->SetFileName("/vtkOutput/Clipping.vtk"); // writerC->SetFileTypeToASCII(); // writerC->Write(); // MITK_DEBUG << "writing done"; MITK_DEBUG << "init and update clipperoutput"; clipperout->GetPointData()->Initialize(); clipperout->Update(); MITK_DEBUG << "init and update clipperoutput completed"; // MITK_DEBUG << "start clippingRecursive"; // vtkSmartPointer Rclipper = vtkSmartPointer::New(); // Rclipper->SetInput(clipperout1); // Rclipper->SetClipFunction(planeR); // Rclipper->GenerateClipScalarsOn(); // Rclipper->GenerateClippedOutputOn(); // vtkSmartPointer clipperout = Rclipper->GetClippedOutput(); // MITK_DEBUG << "end clipping recursive"; // MITK_DEBUG << "writing clipper output 2"; // vtkSmartPointer writerC1 = vtkSmartPointer::New(); // writerC1->SetInput(clipperout); // writerC1->SetFileName("/vtkOutput/RClipping.vtk"); // writerC1->SetFileTypeToASCII(); // writerC1->Write(); // MITK_DEBUG << "init and update clipperoutput"; // clipperout->GetPointData()->Initialize(); // clipperout->Update(); // MITK_DEBUG << "init and update clipperoutput completed"; MITK_DEBUG << "STEP 1: find all points which have distance 0 to the given plane"; /*======STEP 1====== * extract all points, which are crossing the plane */ // Scalar values describe the distance between each remaining point to the given plane. Values sorted by point index vtkSmartPointer distanceList = clipperout->GetPointData()->GetScalars(); vtkIdType sizeOfList = distanceList->GetNumberOfTuples(); PointsOnPlane.reserve(sizeOfList); /* use reserve for high-performant push_back, no hidden copy procedures are processed then! * size of list can be optimized by reducing allocation, but be aware of iterator and vector size*/ for (int i=0; iGetTuple(i); // check if point is on plane. // 0.01 due to some approximation errors when calculating distance if (distance[0] >= -0.01 && distance[0] <= 0.01) PointsOnPlane.push_back(i); } // DEBUG print out all interesting points, stop where array starts with value -1. after -1 no more interesting idx are set! // std::vector::iterator rit = PointsOnPlane.begin(); // while (rit != PointsOnPlane.end() ) { // std::cout << "interesting point: " << *rit << " coord: " << clipperout->GetPoint(*rit)[0] << " | " << clipperout->GetPoint(*rit)[1] << " | " << clipperout->GetPoint(*rit)[2] << endl; // rit++; // } MITK_DEBUG << "Num Of points on plane: " << PointsOnPlane.size(); MITK_DEBUG << "Step 2: extract Interesting points with respect to given extraction planarFigure"; PointsInROI.reserve(PointsOnPlane.size()); /*=======STEP 2===== * extract ROI relevant pointIds */ mitk::PlanarCircle::Pointer circleName = mitk::PlanarCircle::New(); mitk::PlanarPolygon::Pointer polyName = mitk::PlanarPolygon::New(); if ( pf->GetNameOfClass() == circleName->GetNameOfClass() ) { //calculate circle radius mitk::Point3D V1w = pf->GetWorldControlPoint(0); //centerPoint mitk::Point3D V2w = pf->GetWorldControlPoint(1); //radiusPoint double distPF = V1w.EuclideanDistanceTo(V2w); for (int i=0; iGetPoint(PointsOnPlane[i])[0] - V1w[0]) * (clipperout->GetPoint(PointsOnPlane[i])[0] - V1w[0]) + (clipperout->GetPoint(PointsOnPlane[i])[1] - V1w[1]) * (clipperout->GetPoint(PointsOnPlane[i])[1] - V1w[1]) + (clipperout->GetPoint(PointsOnPlane[i])[2] - V1w[2]) * (clipperout->GetPoint(PointsOnPlane[i])[2] - V1w[2])) ; if( XdistPnt <= distPF) PointsInROI.push_back(PointsOnPlane[i]); } } else if ( pf->GetNameOfClass() == polyName->GetNameOfClass() ) { //create vtkPolygon using controlpoints from planarFigure polygon vtkSmartPointer polygonVtk = vtkPolygon::New(); //get the control points from pf and insert them to vtkPolygon unsigned int nrCtrlPnts = pf->GetNumberOfControlPoints(); for (int i=0; iGetPoints()->InsertNextPoint((double)pf->GetWorldControlPoint(i)[0], (double)pf->GetWorldControlPoint(i)[1], (double)pf->GetWorldControlPoint(i)[2] ); } //prepare everything for using pointInPolygon function double n[3]; polygonVtk->ComputeNormal(polygonVtk->GetPoints()->GetNumberOfPoints(), static_cast(polygonVtk->GetPoints()->GetData()->GetVoidPointer(0)), n); double bounds[6]; polygonVtk->GetPoints()->GetBounds(bounds); for (int i=0; iGetPoint(PointsOnPlane[i])[0], clipperout->GetPoint(PointsOnPlane[i])[1], clipperout->GetPoint(PointsOnPlane[i])[2]}; int isInPolygon = polygonVtk->PointInPolygon(checkIn, polygonVtk->GetPoints()->GetNumberOfPoints() , static_cast(polygonVtk->GetPoints()->GetData()->GetVoidPointer(0)), bounds, n); if( isInPolygon ) PointsInROI.push_back(PointsOnPlane[i]); } } MITK_DEBUG << "Step3: Identify fibers"; // we need to access the fiberId Array, so make sure that this array is available if (!clipperout->GetCellData()->HasArray(FIBER_ID_ARRAY)) { MITK_DEBUG << "ERROR: FiberID array does not exist, no correlation between points and fiberIds possible! Make sure calling GenerateFiberIds()"; return FibersInROI; // FibersInRoi is empty then } if (PointsInROI.size()<=0) return FibersInROI; // prepare a structure where each point id is represented as an indexId. // vector looks like: | pntId | fiberIdx | std::vector< long > pointindexFiberMap; // walk through the whole subline section and create an vector sorted by point index vtkCellArray *clipperlines = clipperout->GetLines(); clipperlines->InitTraversal(); long numOfLineCells = clipperlines->GetNumberOfCells(); long numofClippedPoints = clipperout->GetNumberOfPoints(); pointindexFiberMap.resize(numofClippedPoints); //prepare resulting vector FibersInROI.reserve(PointsInROI.size()); MITK_DEBUG << "\n===== Pointindex based structure initialized ======\n"; // go through resulting "sub"lines which are stored as cells, "i" corresponds to current line id. for (int i=0, ic=0 ; iGetCell(ic, npts, pts); // go through point ids in hosting subline, "j" corresponds to current pointindex in current line i. eg. idx[0]=45; idx[1]=46 for (long j=0; jGetCellData()->GetArray(FIBER_ID_ARRAY)->GetTuple(i)[0] << " to pointId: " << pts[j]; pointindexFiberMap[ pts[j] ] = clipperout->GetCellData()->GetArray(FIBER_ID_ARRAY)->GetTuple(i)[0]; // MITK_DEBUG << "in array: " << pointindexFiberMap[ pts[j] ]; } } MITK_DEBUG << "\n===== Pointindex based structure finalized ======\n"; // get all Points in ROI with according fiberID for (long k = 0; k < PointsInROI.size(); k++) { //MITK_DEBUG << "point " << PointsInROI[k] << " belongs to fiber " << pointindexFiberMap[ PointsInROI[k] ]; if (pointindexFiberMap[ PointsInROI[k] ]<=GetNumFibers() && pointindexFiberMap[ PointsInROI[k] ]>=0) FibersInROI.push_back(pointindexFiberMap[ PointsInROI[k] ]); else MITK_INFO << "ERROR in ExtractFiberIdSubset; impossible fiber id detected"; } m_PointsRoi = FibersInROI; } // detecting fiberId duplicates MITK_DEBUG << "check for duplicates"; sort(FibersInROI.begin(), FibersInROI.end()); bool hasDuplicats = false; for(long i=0; i::iterator it; it = unique (FibersInROI.begin(), FibersInROI.end()); FibersInROI.resize( it - FibersInROI.begin() ); } return FibersInROI; } void mitk::FiberBundleX::UpdateFiberGeometry() { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInput(m_FiberPolyData); cleaner->PointMergingOff(); cleaner->Update(); m_FiberPolyData = cleaner->GetOutput(); m_FiberLengths.clear(); m_MeanFiberLength = 0; m_MedianFiberLength = 0; m_LengthStDev = 0; m_NumFibers = m_FiberPolyData->GetNumberOfLines(); if (m_NumFibers<=0) // no fibers present; apply default geometry { m_MinFiberLength = 0; m_MaxFiberLength = 0; mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetImageGeometry(true); float b[] = {0, 1, 0, 1, 0, 1}; geometry->SetFloatBounds(b); SetGeometry(geometry); return; } float min = itk::NumericTraits::NonpositiveMin(); float max = itk::NumericTraits::max(); float b[] = {max, min, max, min, max, min}; vtkCellArray* cells = m_FiberPolyData->GetLines(); cells->InitTraversal(); for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float length = 0; for (int j=0; jGetPoint(j, p1); if (p1[0]b[1]) b[1]=p1[0]; if (p1[1]b[3]) b[3]=p1[1]; if (p1[2]b[5]) b[5]=p1[2]; // calculate statistics if (jGetPoint(j+1, p2); float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); length += dist; } } m_FiberLengths.push_back(length); m_MeanFiberLength += length; if (i==0) { m_MinFiberLength = length; m_MaxFiberLength = length; } else { if (lengthm_MaxFiberLength) m_MaxFiberLength = length; } } m_MeanFiberLength /= m_NumFibers; std::vector< float > sortedLengths = m_FiberLengths; std::sort(sortedLengths.begin(), sortedLengths.end()); for (int i=0; i1) m_LengthStDev /= (m_NumFibers-1); else m_LengthStDev = 0; m_LengthStDev = std::sqrt(m_LengthStDev); m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); // provide some border margin for(int i=0; i<=4; i+=2) b[i] -=10; for(int i=1; i<=5; i+=2) b[i] +=10; mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetFloatBounds(b); this->SetGeometry(geometry); } QStringList mitk::FiberBundleX::GetAvailableColorCodings() { QStringList availableColorCodings; int numColors = m_FiberPolyData->GetPointData()->GetNumberOfArrays(); for(int i=0; iGetPointData()->GetArrayName(i)); } //this controlstructure shall be implemented by the calling method if (availableColorCodings.isEmpty()) MITK_DEBUG << "no colorcodings available in fiberbundleX"; return availableColorCodings; } char* mitk::FiberBundleX::GetCurrentColorCoding() { return m_CurrentColorCoding; } void mitk::FiberBundleX::SetColorCoding(const char* requestedColorCoding) { if (requestedColorCoding==NULL) return; MITK_DEBUG << "SetColorCoding:" << requestedColorCoding; if( strcmp (COLORCODING_ORIENTATION_BASED,requestedColorCoding) == 0 ) { this->m_CurrentColorCoding = (char*) COLORCODING_ORIENTATION_BASED; } else if( strcmp (COLORCODING_FA_BASED,requestedColorCoding) == 0 ) { this->m_CurrentColorCoding = (char*) COLORCODING_FA_BASED; } else if( strcmp (COLORCODING_CUSTOM,requestedColorCoding) == 0 ) { this->m_CurrentColorCoding = (char*) COLORCODING_CUSTOM; } else { MITK_DEBUG << "FIBERBUNDLE X: UNKNOWN COLORCODING in FIBERBUNDLEX Datastructure"; this->m_CurrentColorCoding = (char*) COLORCODING_CUSTOM; //will cause blank colorcoding of fibers } } void mitk::FiberBundleX::MirrorFibers(unsigned int axis) { if (axis>2) return; vtkSmartPointer vtkNewPoints = vtkPoints::New(); vtkSmartPointer vtkNewCells = vtkCellArray::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); for (int i=0; iGetNextCell ( numPoints, pointIds ); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(pointIds[j]); p[axis] = -p[axis]; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); } bool mitk::FiberBundleX::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { if (minRadius<0) return true; vtkSmartPointer vtkNewPoints = vtkPoints::New(); vtkSmartPointer vtkNewCells = vtkCellArray::New(); vtkSmartPointer vtkOldCells = m_FiberPolyData->GetLines(); vtkOldCells->InitTraversal(); for (int i=0; iGetNumberOfCells(); i++) { MITK_INFO << "Processing fiber " << i << "/" << m_NumFibers; vtkIdType numPoints(0); vtkIdType* points(NULL); vtkOldCells->GetNextCell ( numPoints, points ); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(points[j], p1); double p2[3]; m_FiberPolyData->GetPoint(points[j+1], p2); double p3[3]; m_FiberPolyData->GetPoint(points[j+2], p3); vnl_vector_fixed< float, 3 > v1, v2, v3; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; v2[0] = p3[0]-p2[0]; v2[1] = p3[1]-p2[1]; v2[2] = p3[2]-p2[2]; v3[0] = p1[0]-p3[0]; v3[1] = p1[1]-p3[1]; v3[2] = p1[2]-p3[2]; float a = v1.magnitude(); float b = v2.magnitude(); float c = v3.magnitude(); float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); if (deleteFibers && rInsertNextPoint(p2); // container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); container = vtkSmartPointer::New(); } else if (j==numPoints-3) { id = vtkNewPoints->InsertNextPoint(p2); container->GetPointIds()->InsertNextId(id); id = vtkNewPoints->InsertNextPoint(p3); container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); } } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } bool mitk::FiberBundleX::RemoveShortFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers return false; vtkSmartPointer vtkNewPoints = vtkPoints::New(); vtkSmartPointer vtkNewCells = vtkCellArray::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); float min = m_MaxFiberLength; for (int i=0; iGetNextCell ( numPoints, pointIds ); if (m_FiberLengths.at(i)>=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(pointIds[j]); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); if (m_FiberLengths.at(i)GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } bool mitk::FiberBundleX::RemoveLongFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) return true; if (lengthInMM vtkNewPoints = vtkPoints::New(); vtkSmartPointer vtkNewCells = vtkCellArray::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); for (int i=0; iGetNextCell ( numPoints, pointIds ); if (m_FiberLengths.at(i)<=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(pointIds[j]); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } void mitk::FiberBundleX::DoFiberSmoothing(int pointsPerCm) { vtkSmartPointer vtkSmoothPoints = vtkPoints::New(); //in smoothpoints the interpolated points representing a fiber are stored. //in vtkcells all polylines are stored, actually all id's of them are stored vtkSmartPointer vtkSmoothCells = vtkCellArray::New(); //cellcontainer for smoothed lines vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); vtkIdType pointHelperCnt = 0; for (int i=0; iGetNextCell ( numPoints, pointIds ); vtkSmartPointer points = vtkSmartPointer::New(); for (int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(pointIds[j])); float length = m_FiberLengths.at(i); length /=10; int sampling = pointsPerCm*length; vtkSmartPointer xSpline = vtkKochanekSpline::New(); vtkSmartPointer ySpline = vtkKochanekSpline::New(); vtkSmartPointer zSpline = vtkKochanekSpline::New(); vtkSmartPointer spline = vtkParametricSpline::New(); spline->SetXSpline(xSpline); spline->SetYSpline(ySpline); spline->SetZSpline(zSpline); spline->SetPoints(points); vtkSmartPointer functionSource = vtkParametricFunctionSource::New(); functionSource->SetParametricFunction(spline); functionSource->SetUResolution(sampling); functionSource->SetVResolution(sampling); functionSource->SetWResolution(sampling); functionSource->Update(); vtkPolyData* outputFunction = functionSource->GetOutput(); vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber vtkSmartPointer smoothLine = vtkPolyLine::New(); smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints()); for (int j=0; jGetNumberOfPoints(); j++) { smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt); vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); } vtkSmoothCells->InsertNextCell(smoothLine); pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints(); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkSmoothPoints); m_FiberPolyData->SetLines(vtkSmoothCells); UpdateColorCoding(); UpdateFiberGeometry(); } // Resample fiber to get equidistant points void mitk::FiberBundleX::ResampleFibers(float pointDistance) { vtkSmartPointer newPoly = vtkSmartPointer::New(); vtkSmartPointer newCellArray = vtkSmartPointer::New(); vtkSmartPointer newPoints = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); int numberOfLines = m_NumFibers; for (int i=0; iGetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); double* point = m_FiberPolyData->GetPoint(points[0]); vtkIdType pointId = newPoints->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); float dtau = 0; int cur_p = 1; itk::Vector dR; float normdR = 0; for (;;) { while (dtau <= pointDistance && cur_p < numPoints) { itk::Vector v1; point = m_FiberPolyData->GetPoint(points[cur_p-1]); v1[0] = point[0]; v1[1] = point[1]; v1[2] = point[2]; itk::Vector v2; point = m_FiberPolyData->GetPoint(points[cur_p]); v2[0] = point[0]; v2[1] = point[1]; v2[2] = point[2]; dR = v2 - v1; normdR = std::sqrt(dR.GetSquaredNorm()); dtau += normdR; cur_p++; } if (dtau >= pointDistance) { itk::Vector v1; point = m_FiberPolyData->GetPoint(points[cur_p-1]); v1[0] = point[0]; v1[1] = point[1]; v1[2] = point[2]; itk::Vector v2 = v1 - dR*( (dtau-pointDistance)/normdR ); pointId = newPoints->InsertNextPoint(v2.GetDataPointer()); container->GetPointIds()->InsertNextId(pointId); } else { point = m_FiberPolyData->GetPoint(points[numPoints-1]); pointId = newPoints->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); break; } dtau = dtau-pointDistance; } newCellArray->InsertNextCell(container); } newPoly->SetPoints(newPoints); newPoly->SetLines(newCellArray); m_FiberPolyData = newPoly; UpdateFiberGeometry(); UpdateColorCoding(); } // reapply selected colorcoding in case polydata structure has changed void mitk::FiberBundleX::UpdateColorCoding() { char* cc = GetCurrentColorCoding(); if( strcmp (COLORCODING_ORIENTATION_BASED,cc) == 0 ) DoColorCodingOrientationBased(); else if( strcmp (COLORCODING_FA_BASED,cc) == 0 ) DoColorCodingFaBased(); } // reapply selected colorcoding in case polydata structure has changed bool mitk::FiberBundleX::Equals(mitk::FiberBundleX* fib) { if (fib==NULL) return false; mitk::FiberBundleX::Pointer tempFib = this->SubtractBundle(fib); mitk::FiberBundleX::Pointer tempFib2 = fib->SubtractBundle(this); if (tempFib.IsNull() && tempFib2.IsNull()) return true; return false; } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundleX::UpdateOutputInformation() { } void mitk::FiberBundleX::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundleX::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundleX::VerifyRequestedRegion() { return true; } void mitk::FiberBundleX::SetRequestedRegion( itk::DataObject *data ) { }