diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp index 64ff3746b0..bc4961ed09 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,708 +1,704 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2010-03-31 16:40:27 +0200 (Mi, 31 Mrz 2010) $ Version: $Revision: 21975 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkFiberBundleX.h" #include //#include #include #include /* musthave */ //#include // without geometry, fibers are not rendered #include #include #include #include #include #include #include #include #include // baptize array names const char* mitk::FiberBundleX::COLORCODING_ORIENTATION_BASED = "Color_Orient"; const char* mitk::FiberBundleX::COLORCODING_FA_BASED = "Color_FA"; const char* mitk::FiberBundleX::FIBER_ID_ARRAY = "Fiber_IDs"; mitk::FiberBundleX::FiberBundleX(vtkSmartPointer fiberPolyData ) : m_currentColorCoding(NULL) , m_isModified(false) { //generate geometry of passed polydata if (fiberPolyData == NULL) this->m_FiberPolyData = vtkSmartPointer::New(); else this->m_FiberPolyData = fiberPolyData; this->UpdateFiberGeometry(); } mitk::FiberBundleX::~FiberBundleX() { } /* * set computed fibers from tractography algorithms */ void mitk::FiberBundleX::SetFiberPolyData(vtkSmartPointer fiberPD) { if (fiberPD == NULL) this->m_FiberPolyData = vtkSmartPointer::New(); else this->m_FiberPolyData = fiberPD; m_isModified = true; } /* * return fiberbundle as vtkPolyData * Depending on processing of input fibers, this method returns * the latest processed fibers. */ vtkSmartPointer mitk::FiberBundleX::GetFiberPolyData() { return m_FiberPolyData; } /*=================================== *++++ PROCESSING WITH FIBERS +++++++ ====================================*/ 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_INFO << " NO NEED TO REGENERATE COLORCODING! " ; return; } /* Finally, execute color calculation */ vtkPoints* extrPoints = m_FiberPolyData->GetPoints(); int 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); vtkUnsignedCharArray * 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_INFO << "\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_INFO << "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::abs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::abs(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::abs(diff1[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff1[1])); rgba[2] = (unsigned char) (255.0 * std::abs(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::abs(diff2[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff2[1])); rgba[2] = (unsigned char) (255.0 * std::abs(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_INFO << "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); m_isModified = true; // =========================== //mini test, shall be ported to MITK TESTINGS! if (colorsT->GetSize() != numOfPoints*componentSize) { MITK_INFO << "ALLOCATION ERROR IN INITIATING COLOR ARRAY"; } } void mitk::FiberBundleX::DoGenerateFiberIds() { if (m_FiberPolyData == NULL) return; // for (int i=0; i<10000000; ++i) // { // if(i%500 == 0) // MITK_INFO << i; // } // MITK_INFO << "Generating Fiber Ids"; 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_INFO << "Generating Fiber Ids...[done] | " << m_FiberIdDataSet->GetNumberOfCells(); } //temporarely include only #include //========================== std::vector mitk::FiberBundleX::DoExtractFiberIds(mitk::PlanarFigure::Pointer pf) { MITK_INFO << "Extracting fiber!"; /* Handle type of planarfigure */ // if incoming pf is a pfc mitk::PlanarFigureComposite::Pointer pfcomp= dynamic_cast(pf.GetPointer()); if (!pfcomp.IsNull()) { // process requested boolean operation of PFC - } else { + } 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_INFO << "planeOrigin: " << planeOrigin[0] << " | " << planeOrigin[1] << " | " << planeOrigin[2] << endl; MITK_INFO << "planeNormal: " << planeNormal[0] << " | " << planeNormal[1] << " | " << planeNormal[2] << endl; - - - - - /* init necessary vectors hosting pointIds and FiberIds */ // contains all pointIds which are crossing the cutting plane std::vector PointsOnPlane; // based on PointsOnPlane, all ROI relevant point IDs are stored here std::vector PointsInROI; // vector which is returned, contains all extracted FiberIds std::vector FibersInROI; /* 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(); planeR->SetOrigin(10.0,5.0,0.0); planeR->SetNormal(0.0,-1.0,0.0); /* get all points/fibers cutting the plane */ vtkSmartPointer clipper = vtkSmartPointer::New(); clipper->SetInput(m_FiberIdDataSet); clipper->SetClipFunction(plane); clipper->GenerateClipScalarsOn(); clipper->GenerateClippedOutputOn(); vtkSmartPointer clipperout1 = clipper->GetClippedOutput(); /* for some reason clipperoutput is not initialized for futher processing * so far only writing out clipped polydata provides requested */ vtkSmartPointer writerC = vtkSmartPointer::New(); writerC->SetInput(clipperout1); writerC->SetFileName("/vtkOutput/Cout1_FbId_clipLineId0+1+2-tests.vtk"); writerC->SetFileTypeToASCII(); writerC->Write(); vtkSmartPointer Rclipper = vtkSmartPointer::New(); Rclipper->SetInput(clipperout1); Rclipper->SetClipFunction(planeR); Rclipper->GenerateClipScalarsOn(); Rclipper->GenerateClippedOutputOn(); vtkSmartPointer clipperout = Rclipper->GetClippedOutput(); vtkSmartPointer writerC1 = vtkSmartPointer::New(); writerC1->SetInput(clipperout); writerC1->SetFileName("/vtkOutput/FbId_clipLineId0+1+2-tests.vtk"); writerC1->SetFileTypeToASCII(); writerC1->Write(); /*======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); std::cout << "distance of point " << i << " : " << distance[0] << std::endl; // 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) { std::cout << "adding " << i << endl; PointsOnPlane.push_back(i); //push back in combination with reserve is fastest way to fill vector with various values } } // 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++; } /*=======STEP 2===== * extract ROI relevant pointIds */ //ToDo mitk::PlanarCircle::Pointer circleName = mitk::PlanarCircle::New(); mitk::PlanarPolygon::Pointer polyName = mitk::PlanarPolygon::New(); if (pf->GetNameOfClass() == circleName->GetNameOfClass() ) { if( true /*point in ROI*/) { PointsInROI = PointsOnPlane; } } /*======STEP 3======= * identify fiberIds for points in ROI */ //prepare resulting vector FibersInROI.reserve(PointsInROI.size()); vtkCellArray *clipperlines = clipperout->GetLines(); clipperlines->InitTraversal(); long numOfLineCells = clipperlines->GetNumberOfCells(); // 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. for (long j=0; jGetCellData()->HasArray("FB_IDs")) { int originalFibId = clipperout->GetCellData()->GetArray("FB_IDs")->GetTuple(i)[0]; std::cout << "found pointid " << PointsInROI[k] << ": " << clipperout->GetPoint(PointsInROI[k])[0] << " | " << clipperout->GetPoint(PointsInROI[k])[1] << " | " << clipperout->GetPoint(PointsInROI[k])[2] << " in subline: " << i << " which belongs to fiber id: " << originalFibId << "\n" << endl; // do something to avoid duplicates int oldFibInRoiSize = FibersInROI.size(); if (oldFibInRoiSize != 0) { for (int f=0; f::iterator finIt = FibersInROI.begin(); - while ( finIt != FibersInROI.end() ) { + while ( finIt != FibersInROI.end() ) + { std::cout << *finIt << endl; ++finIt; } std::cout << "=====================\n"; } +} - void mitk::FiberBundleX::UpdateFiberGeometry() - { - +void mitk::FiberBundleX::UpdateFiberGeometry() +{ - float min = itk::NumericTraits::min(); - 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++) - { + float min = itk::NumericTraits::min(); + float max = itk::NumericTraits::max(); + float b[] = {max, min, max, min, max, min}; - vtkCell* cell = m_FiberPolyData->GetCell(i); - int p = cell->GetNumberOfPoints(); - vtkPoints* points = cell->GetPoints(); - for (int j=0; jGetPoint(j, p); + vtkCellArray* cells = m_FiberPolyData->GetLines(); + cells->InitTraversal(); + for (int i=0; iGetNumberOfCells(); i++) + { - if (p[0]b[1]) - b[1]=p[0]; + vtkCell* cell = m_FiberPolyData->GetCell(i); + int p = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + for (int j=0; jGetPoint(j, p); - if (p[1]b[3]) - b[3]=p[1]; + if (p[0]b[1]) + b[1]=p[0]; - if (p[2]b[5]) - b[5]=p[2]; + if (p[1]b[3]) + b[3]=p[1]; + if (p[2]b[5]) + b[5]=p[2]; - } } - // provide some buffer space at borders - - for(int i=0; i<=4; i+=2){ - b[i] -=10; - } + } - for(int i=1; i<=5; i+=2){ - b[i] +=10; - } + // provide some buffer space at borders - mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); - geometry->SetImageGeometry(true); - geometry->SetFloatBounds(b); - this->SetGeometry(geometry); + 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->SetImageGeometry(true); + geometry->SetFloatBounds(b); + this->SetGeometry(geometry); - } } /*============================== *++++ FIBER INFORMATION +++++++ ===============================*/ 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_INFO << "no colorcodings available in fiberbundleX"; // for(int i=0; im_currentColorCoding; } void mitk::FiberBundleX::SetColorCoding(const char* requestedColorCoding) { // MITK_INFO << "FbX try to set colorCoding: " << requestedColorCoding << " compare with: " << COLORCODING_ORIENTATION_BASED; if(strcmp (COLORCODING_ORIENTATION_BASED,requestedColorCoding) == 0 ) { this->m_currentColorCoding = (char*) COLORCODING_ORIENTATION_BASED; this->m_isModified = true; } else if(strcmp (COLORCODING_FA_BASED,requestedColorCoding) == 0 ) { this->m_currentColorCoding = (char*) COLORCODING_FA_BASED; this->m_isModified = true; } else { MITK_INFO << "FIBERBUNDLE X: UNKNOWN COLORCODING in FIBERBUNDLEX Datastructure"; this->m_currentColorCoding = "---"; //will cause blank colorcoding of fibers this->m_isModified = true; } } bool mitk::FiberBundleX::isFiberBundleXModified() { return m_isModified; } void mitk::FiberBundleX::setFBXModificationDone() { m_isModified = false; } // Resample fiber to get equidistant points void mitk::FiberBundleX::ResampleFibers(float len) { vtkSmartPointer newPoly = vtkSmartPointer::New(); vtkSmartPointer newCellArray = vtkSmartPointer::New(); vtkSmartPointer newPoints = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); int numberOfLines = vLines->GetNumberOfCells(); 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 <= len && 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 >= len) { 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-len)/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-len; } newCellArray->InsertNextCell(container); } newPoly->SetPoints(newPoints); newPoly->SetLines(newCellArray); m_FiberPolyData = newPoly; UpdateFiberGeometry(); } /* 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 ) { }