diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp index f254971e40..74ba6604c6 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,1288 +1,1374 @@ /*=================================================================== 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 "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"; mitk::FiberBundleX::FiberBundleX( vtkPolyData* fiberPolyData ) : m_CurrentColorCoding(NULL) , m_NumFibers(0) { m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != NULL) { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInput(fiberPolyData); cleaner->Update(); fiberPolyData = cleaner->GetOutput(); m_FiberPolyData->DeepCopy(fiberPolyData); this->DoColorCodingOrientationBased(); } if(m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)) MITK_DEBUG << "ok"; vtkUnsignedCharArray* tmpColors = (vtkUnsignedCharArray*) m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED); if (tmpColors!=NULL) { int tmpColorss = tmpColors->GetNumberOfTuples(); int tmpColorc = tmpColors->GetNumberOfComponents(); } 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 ); vtkSmartPointer vLines2 = fib->m_FiberPolyData->GetLines(); vLines2->InitTraversal(); int numFibers2 = fib->GetNumFibers(); bool contained = false; for( int i2=0; i2GetNextCell ( numPoints2, points2 ); // 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(); } 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"; } } // 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() { if (m_NumFibers<=0) // no fibers present; apply default geometry { 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(); for (int j=0; jGetPoint(j, p); if (p[0]b[1]) b[1]=p[0]; if (p[1]b[3]) b[3]=p[1]; if (p[2]b[5]) b[5]=p[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"; // for(int i=0; im_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::DoFiberSmoothing(int pointsPerCm) +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::RemoveShortFibers(float lengthInMM) +{ + if (lengthInMM<=0) + return false; + + vtkSmartPointer vtkNewPoints = vtkPoints::New(); + vtkSmartPointer vtkNewCells = vtkCellArray::New(); + + vtkSmartPointer vLines = m_FiberPolyData->GetLines(); + vLines->InitTraversal(); + for (int i=0; iGetNextCell ( numPoints, pointIds ); + // calculate fiber length + float length = 0; + itk::Point lastP; + for (int j=0; jGetPoint(pointIds[j]); + if (j>0) + length += sqrt(pow(p[0]-lastP[0], 2)+pow(p[1]-lastP[1], 2)+pow(p[2]-lastP[2], 2)); + lastP[0] = p[0]; + lastP[1] = p[1]; + lastP[2] = p[2]; + } + + if (length>=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(); float length = 0; itk::Point lastP; for (int j=0; jGetPoint(pointIds[j]); points->InsertNextPoint(p); if (j>0) length += sqrt(pow(p[0]-lastP[0], 2)+pow(p[1]-lastP[1], 2)+pow(p[2]-lastP[2], 2)); lastP[0] = p[0]; lastP[1] = p[1]; lastP[2] = p[2]; } length /=10; int sampling = pointsPerCm*length; /////PROCESS POLYLINE SMOOTHING///// 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 ) { } diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h index a6a6f5e76f..05ea707efb 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h @@ -1,122 +1,125 @@ /*=================================================================== 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 _MITK_FiberBundleX_H #define _MITK_FiberBundleX_H //includes for MITK datastructure #include #include "MitkDiffusionImagingExports.h" #include //includes storing fiberdata #include //may be replaced by class precompile argument #include // may be replaced by class #include // my be replaced by class #include #include #include namespace mitk { /** * \brief Base Class for Fiber Bundles; */ class MitkDiffusionImaging_EXPORT FiberBundleX : public BaseData { public: // fiber colorcodings static const char* COLORCODING_ORIENTATION_BASED; static const char* COLORCODING_FA_BASED; static const char* COLORCODING_CUSTOM; static const char* FIBER_ID_ARRAY; virtual void UpdateOutputInformation(); virtual void SetRequestedRegionToLargestPossibleRegion(); virtual bool RequestedRegionIsOutsideOfTheBufferedRegion(); virtual bool VerifyRequestedRegion(); virtual void SetRequestedRegion( itk::DataObject *data ); mitkClassMacro( FiberBundleX, BaseData ) itkNewMacro( Self ) mitkNewMacro1Param(Self, vtkSmartPointer) // custom constructor // colorcoding related methods void SetColorCoding(const char*); void SetFAMap(mitk::Image::Pointer); void DoColorCodingOrientationBased(); void DoColorCodingFaBased(); void DoUseFaFiberOpacity(); void ResetFiberOpacity(); // fiber smoothing/resampling void ResampleFibers(float pointDistance = 1); void DoFiberSmoothing(int pointsPerCm); + bool RemoveShortFibers(float lengthInMM); + void MirrorFibers(unsigned int axis); // add/subtract fibers FiberBundleX::Pointer AddBundle(FiberBundleX* fib); FiberBundleX::Pointer SubtractBundle(FiberBundleX* fib); // fiber subset extraction - FiberBundleX::Pointer ExtractFiberSubset(PlanarFigure *pf); - std::vector ExtractFiberIdSubset(PlanarFigure* pf); - vtkSmartPointer GeneratePolyDataByIds( std::vector ); + FiberBundleX::Pointer ExtractFiberSubset(PlanarFigure *pf); + std::vector ExtractFiberIdSubset(PlanarFigure* pf); + vtkSmartPointer GeneratePolyDataByIds( std::vector ); // TODO: make protected + void GenerateFiberIds(); // TODO: make protected // get/set data void SetFiberPolyData(vtkSmartPointer, bool updateGeometry = true); vtkSmartPointer GetFiberPolyData(); QStringList GetAvailableColorCodings(); char* GetCurrentColorCoding(); itkGetMacro(NumFibers, int) // copy fiber bundle mitk::FiberBundleX::Pointer GetDeepCopy(); + // compare fiber bundles bool Equals(FiberBundleX* fib); - void GenerateFiberIds(); protected: FiberBundleX( vtkPolyData* fiberPolyData = NULL ); virtual ~FiberBundleX(); itk::Point GetItkPoint(double point[3]); // calculate geometry from fiber extent void UpdateFiberGeometry(); // calculate colorcoding values according to m_CurrentColorCoding void UpdateColorCoding(); private: // actual fiber container vtkSmartPointer m_FiberPolyData; // contains fiber ids vtkSmartPointer m_FiberIdDataSet; char* m_CurrentColorCoding; int m_NumFibers; }; } // namespace mitk #endif /* _MITK_FiberBundleX_H */ diff --git a/Modules/DiffusionImaging/Testing/mitkFiberBundleXTest.cpp b/Modules/DiffusionImaging/Testing/mitkFiberBundleXTest.cpp index 52ca547666..8a3fdfda4b 100644 --- a/Modules/DiffusionImaging/Testing/mitkFiberBundleXTest.cpp +++ b/Modules/DiffusionImaging/Testing/mitkFiberBundleXTest.cpp @@ -1,97 +1,100 @@ /*=================================================================== 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 "mitkTestingMacros.h" #include #include #include #include #include #include #include #include /**Documentation * Test for fiber bundle reader and writer */ int mitkFiberBundleXTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkFiberBundleXTest"); MITK_TEST_CONDITION_REQUIRED(argc>1,"check for fielename") mitk::FiberBundleXReader::Pointer reader = mitk::FiberBundleXReader::New(); mitk::FiberBundleX::Pointer fib1, fib2; // first test: did this work? // using MITK_TEST_CONDITION_REQUIRED makes the test stop after failure, since // it makes no sense to continue without an object. MITK_TEST_CONDITION_REQUIRED(reader.IsNotNull(),"reader instantiation") try{ RegisterDiffusionImagingObjectFactory(); // test if fib1 can be read const std::string s1="", s2=""; std::vector fibInfile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[1], s1, s2, false ); mitk::BaseData::Pointer baseData = fibInfile.at(0); fib1 = dynamic_cast(baseData.GetPointer()); MITK_TEST_CONDITION_REQUIRED(fib1.IsNotNull(),"check if reader 1 returned null") fibInfile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[1], s1, s2, false ); baseData = fibInfile.at(0); fib2 = dynamic_cast(baseData.GetPointer()); MITK_TEST_CONDITION_REQUIRED(fib2.IsNotNull(),"check if reader 2 returned null") MITK_TEST_CONDITION_REQUIRED(fib1->Equals(fib2),"check if equals method is working"); int randNum = rand()%20; - MITK_INFO << "DoFiberSmoothing " << randNum; fib2->DoFiberSmoothing(randNum); + MITK_INFO << "DoFiberSmoothing(" << randNum << ")" << randNum; fib2->DoFiberSmoothing(randNum); MITK_TEST_CONDITION_REQUIRED(!fib1->Equals(fib2),"check if fiber resampling method does something"); mitk::FiberBundleX::Pointer fib3 = fib1->AddBundle(fib2); MITK_TEST_CONDITION_REQUIRED(!fib1->Equals(fib3),"check if A+B!=A"); fib3 = fib3->SubtractBundle(fib2); MITK_TEST_CONDITION_REQUIRED(fib1->Equals(fib3),"check if A+B-B==A"); fib1->AddBundle(NULL); MITK_INFO << "GenerateFiberIds"; fib1->GenerateFiberIds(); MITK_INFO << "GetFiberPolyData"; fib1->GetFiberPolyData(); MITK_INFO << "GetAvailableColorCodings"; fib1->GetAvailableColorCodings(); MITK_INFO << "GetCurrentColorCoding"; fib1->GetCurrentColorCoding(); MITK_INFO << "SetFiberPolyData"; fib1->SetFiberPolyData(NULL); MITK_INFO << "ExtractFiberSubset"; fib1->ExtractFiberSubset(NULL); MITK_INFO << "ExtractFiberIdSubset"; fib1->ExtractFiberIdSubset(NULL); std::vector< long > tmp; MITK_INFO << "GeneratePolyDataByIds"; fib1->GeneratePolyDataByIds(tmp); MITK_INFO << "SetColorCoding"; fib1->SetColorCoding(NULL); MITK_INFO << "SetFAMap"; fib1->SetFAMap(NULL); MITK_INFO << "DoColorCodingOrientationBased"; fib1->DoColorCodingOrientationBased(); MITK_INFO << "DoColorCodingFaBased"; fib1->DoColorCodingFaBased(); MITK_INFO << "DoUseFaFiberOpacity"; fib1->DoUseFaFiberOpacity(); MITK_INFO << "ResetFiberOpacity"; fib1->ResetFiberOpacity(); + + float randFloat = rand()%300; + MITK_INFO << "RemoveShortFibers(" << randFloat << ")"; fib1->RemoveShortFibers(randFloat); } catch(...) { //this means that a wrong exception (i.e. no itk:Exception) has been thrown std::cout << "Wrong exception (i.e. no itk:Exception) caught during write [FAILED]" << std::endl; return EXIT_FAILURE; } // always end with this! MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index b74a54a9bf..4ee60e05b9 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,653 +1,661 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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 "itkGibbsTrackingFilter.h" #include #include #include "itkPointShell.h" #include "GibbsTracking/BuildFibres.cpp" #pragma GCC visibility push(default) #include #pragma GCC visibility pop #include #include #include #include #include #include #include #include #include #include struct LessDereference { template bool operator()(const T * lhs, const T * rhs) const { return *lhs < *rhs; } }; namespace itk{ template< class TInputOdfImage, class TInputROIImage > GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GibbsTrackingFilter(): m_TempStart(0.1), m_TempEnd(0.001), m_NumIt(500000), m_ParticleWeight(0), m_ParticleWidth(0), m_ParticleLength(0), m_ChempotConnection(10), m_ChempotParticle(0), m_InexBalance(0), m_Chempot2(0.2), m_FiberLength(10), m_AbortTracking(false), m_NumConnections(0), m_NumParticles(0), m_NumAcceptedFibers(0), m_CurrentStep(0), m_SubtractMean(true), m_BuildFibers(false), m_Sampler(NULL), m_Steps(10), m_Memory(0), m_ProposalAcceptance(0), m_GfaImage(NULL), m_CurvatureHardThreshold(0.7), m_Meanval_sq(0.0) { //this->m_MeasurementFrame.set_identity(); this->SetNumberOfRequiredInputs(2); //Filter needs a DWI image + a Mask Image } template< class TInputOdfImage, class TInputROIImage > GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::~GibbsTrackingFilter(){ delete BESSEL_APPROXCOEFF; if (m_Sampler!=NULL) delete m_Sampler; } template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::ComputeFiberCorrelationOriginal(){ float bD = 15; vnl_matrix_fixed bDir = *itk::PointShell >::DistributePointShell(); const int N = QBALL_ODFSIZE; vnl_matrix_fixed C = bDir.transpose()*bDir; vnl_matrix_fixed Q = C; for(int i=0; i P = Q*Q; std::vector pointer; pointer.reserve(N*N); double * start = C.data_block(); double * end = start + N*N; for (double * iter = start; iter != end; ++iter) { pointer.push_back(iter); } std::sort(pointer.begin(), pointer.end(), LessDereference()); vnl_vector_fixed alpha; vnl_vector_fixed beta; for (int i=0; i alpha_0; vnl_vector_fixed alpha_2; vnl_vector_fixed alpha_4; vnl_vector_fixed alpha_6; for(int i=0; i T; T.set_column(0,alpha_0); T.set_column(1,alpha_2); T.set_column(2,alpha_4); T.set_column(3,alpha_6); vnl_vector_fixed coeff = vnl_matrix_inverse(T).pinverse()*beta; BESSEL_APPROXCOEFF = new float[4]; BESSEL_APPROXCOEFF[0] = coeff(0); BESSEL_APPROXCOEFF[1] = coeff(1); BESSEL_APPROXCOEFF[2] = coeff(2); BESSEL_APPROXCOEFF[3] = coeff(3); // // OLD // BESSEL_APPROXCOEFF[0] = 0,1982; // BESSEL_APPROXCOEFF[1] = 0.3415; // BESSEL_APPROXCOEFF[2] = -0.9515; // BESSEL_APPROXCOEFF[3] = 1.3423; } template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::ComputeFiberCorrelation(){ // float bD = 15; // vnl_matrix_fixed bDir = // *itk::PointShell >::DistributePointShell(); // const int N = QBALL_ODFSIZE; // vnl_matrix_fixed temp = bDir.transpose(); // vnl_matrix_fixed C = temp*bDir; // vnl_matrix_fixed Q = C; // vnl_vector_fixed mean; // for(int i=0; i repMean; // for (int i=0; i P = Q*Q; // std::vector pointer; // pointer.reserve(N*N); // double * start = C.data_block(); // double * end = start + N*N; // for (double * iter = start; iter != end; ++iter) // { // pointer.push_back(iter); // } // std::sort(pointer.begin(), pointer.end(), LessDereference()); // vnl_vector_fixed alpha; // vnl_vector_fixed beta; // for (int i=0; im_Meanval_sq = (sum*sum)/N; // vnl_vector_fixed alpha_0; // vnl_vector_fixed alpha_2; // vnl_vector_fixed alpha_4; // vnl_vector_fixed alpha_6; // for(int i=0; i T; // T.set_column(0,alpha_0); // T.set_column(1,alpha_2); // T.set_column(2,alpha_4); // T.set_column(3,alpha_6); // vnl_vector_fixed coeff = vnl_matrix_inverse(T).pinverse()*beta; // MITK_INFO << "itkGibbsTrackingFilter: Bessel oefficients: " << coeff; BESSEL_APPROXCOEFF = new float[4]; // BESSEL_APPROXCOEFF[0] = coeff(0); // BESSEL_APPROXCOEFF[1] = coeff(1); // BESSEL_APPROXCOEFF[2] = coeff(2); // BESSEL_APPROXCOEFF[3] = coeff(3); BESSEL_APPROXCOEFF[0] = -0.1714; BESSEL_APPROXCOEFF[1] = 0.5332; BESSEL_APPROXCOEFF[2] = -1.4889; BESSEL_APPROXCOEFF[3] = 2.0389; } // build fibers from tracking result template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::BuildFibers(float* points, int numPoints) { double spacing[3]; spacing[0] = m_ItkQBallImage->GetSpacing().GetElement(0); spacing[1] = m_ItkQBallImage->GetSpacing().GetElement(1); spacing[2] = m_ItkQBallImage->GetSpacing().GetElement(2); m_FiberPolyData = FiberPolyDataType::New(); // initialize array of particles FiberBuilder fiberBuilder(points, numPoints, spacing, m_ItkQBallImage); // label the particles according to fiber affiliation and return polydata m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); MITK_INFO << "itkGibbsTrackingFilter: " << m_NumAcceptedFibers << " accepted"; } // fill output fiber bundle datastructure template< class TInputOdfImage, class TInputROIImage > typename GibbsTrackingFilter< TInputOdfImage, TInputROIImage >::FiberPolyDataType GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GetFiberBundle() { if (!m_AbortTracking) { m_BuildFibers = true; while (m_BuildFibers){} } return m_FiberPolyData; } // get memory allocated for particle grid template< class TInputOdfImage, class TInputROIImage > float GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GetMemoryUsage() { if (m_Sampler!=NULL) return m_Sampler->m_ParticleGrid.GetMemoryUsage(); return 0; } template< class TInputOdfImage, class TInputROIImage > bool GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::EstimateParticleWeight() { MITK_INFO << "itkGibbsTrackingFilter: estimating particle weight"; typedef itk::DiffusionQballGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(m_ItkQBallImage); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); m_GfaImage = gfaFilter->GetOutput(); float samplingStart = 1.0; float samplingStop = 0.66; // copy GFA image (original should not be changed) typedef itk::ImageDuplicator< GfaImageType > DuplicateFilterType; DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); duplicator->SetInputImage( m_GfaImage ); duplicator->Update(); m_GfaImage = duplicator->GetOutput(); //// GFA iterator //// typedef ImageRegionIterator< GfaImageType > GfaIteratorType; GfaIteratorType gfaIt(m_GfaImage, m_GfaImage->GetLargestPossibleRegion() ); //// Mask iterator //// typedef ImageRegionConstIterator< MaskImageType > MaskIteratorType; MaskIteratorType maskIt(m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); // set unmasked region of gfa image to 0 gfaIt.GoToBegin(); maskIt.GoToBegin(); while( !gfaIt.IsAtEnd() ) { if(maskIt.Get()<=0) gfaIt.Set(0); ++gfaIt; ++maskIt; } // rescale gfa image to [0,1] typedef itk::RescaleIntensityImageFilter< GfaImageType, GfaImageType > RescaleFilterType; RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New(); rescaleFilter->SetInput( m_GfaImage ); rescaleFilter->SetOutputMaximum( samplingStart ); rescaleFilter->SetOutputMinimum( 0 ); rescaleFilter->Update(); m_GfaImage = rescaleFilter->GetOutput(); gfaIt = GfaIteratorType(m_GfaImage, m_GfaImage->GetLargestPossibleRegion() ); //// Input iterator //// typedef ImageRegionConstIterator< InputQBallImageType > InputIteratorType; InputIteratorType git(m_ItkQBallImage, m_ItkQBallImage->GetLargestPossibleRegion() ); float upper = 0; int count = 0; for(float thr=samplingStart; thr>samplingStop; thr-=0.01) { git.GoToBegin(); gfaIt.GoToBegin(); while( !gfaIt.IsAtEnd() ) { if(gfaIt.Get()>thr) { itk::OrientationDistributionFunction odf(git.Get().GetDataPointer()); upper += odf.GetMaxValue()-odf.GetMeanValue(); ++count; } ++gfaIt; ++git; } } if (count>0) upper /= count; else return false; m_ParticleWeight = upper/6; return true; } // perform global tracking template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GenerateData(){ // input qball image m_ItkQBallImage = dynamic_cast(this->GetInput(0)); m_NumAcceptedFibers = 0; // approximationscoeffizienten der // teilchenkorrelationen im orientierungsraum // 4er vektor //ComputeFiberCorrelationOriginal(); ComputeFiberCorrelation(); // image sizes and spacing int qBallImageSize[4] = {QBALL_ODFSIZE, m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(0), m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(1), m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(2)}; double qBallImageSpacing[3] = {m_ItkQBallImage->GetSpacing().GetElement(0),m_ItkQBallImage->GetSpacing().GetElement(1),m_ItkQBallImage->GetSpacing().GetElement(2)}; // make sure image has enough slices if (qBallImageSize[1]<3 || qBallImageSize[2]<3 || qBallImageSize[3]<3) { MITK_INFO << "itkGibbsTrackingFilter: image size < 3 not supported"; m_AbortTracking = true; } // calculate rotation matrix vnl_matrix_fixed directionMatrix = m_ItkQBallImage->GetDirection().GetVnlMatrix(); vnl_vector_fixed d0 = directionMatrix.get_column(0); d0.normalize(); vnl_vector_fixed d1 = directionMatrix.get_column(1); d1.normalize(); vnl_vector_fixed d2 = directionMatrix.get_column(2); d2.normalize(); directionMatrix.set_column(0, d0); directionMatrix.set_column(1, d1); directionMatrix.set_column(2, d2); vnl_matrix_fixed I = directionMatrix*directionMatrix.transpose(); if(!I.is_identity(mitk::eps)){ MITK_INFO << "itkGibbsTrackingFilter: image direction is not a rotation matrix. Tracking not possible!"; m_AbortTracking = true; } // generate local working copy of image buffer int bufferSize = qBallImageSize[0]*qBallImageSize[1]*qBallImageSize[2]*qBallImageSize[3]; float* qBallImageBuffer = (float*) m_ItkQBallImage->GetBufferPointer(); float* workingQballImage = new float[bufferSize]; for (int i=0; i0 && i%qBallImageSize[0] == 0 && i>0) { sum /= qBallImageSize[0]; for (int j=i-qBallImageSize[0]; jGetBufferPointer(); maskImageSize[0] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(0); maskImageSize[1] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(1); maskImageSize[2] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(2); } else { mask = 0; maskImageSize[0] = qBallImageSize[1]; maskImageSize[1] = qBallImageSize[2]; maskImageSize[2] = qBallImageSize[3]; } int mask_oversamp_mult = maskImageSize[0]/qBallImageSize[1]; // load lookuptable QString applicationDir = QCoreApplication::applicationDirPath(); applicationDir.append("/"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); applicationDir.append("../"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); applicationDir.append("../../"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); std::string lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTBaryCoords.bin"); ifstream BaryCoords; BaryCoords.open(lutPath.c_str(), ios::in | ios::binary); float* coords; if (BaryCoords.is_open()) { float tmp; coords = new float [1630818]; BaryCoords.seekg (0, ios::beg); for (int i=0; i<1630818; i++) { BaryCoords.read((char *)&tmp, sizeof(tmp)); coords[i] = tmp; } BaryCoords.close(); } else { MITK_INFO << "itkGibbsTrackingFilter: unable to open barycoords file"; m_AbortTracking = true; } ifstream Indices; lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTIndices.bin"); Indices.open(lutPath.c_str(), ios::in | ios::binary); int* ind; if (Indices.is_open()) { int tmp; ind = new int [1630818]; Indices.seekg (0, ios::beg); for (int i=0; i<1630818; i++) { Indices.read((char *)&tmp, 4); ind[i] = tmp; } Indices.close(); } else { MITK_INFO << "itkGibbsTrackingFilter: unable to open indices file"; m_AbortTracking = true; } // initialize sphere interpolator with lookuptables SphereInterpolator *sinterp = new SphereInterpolator(coords, ind, QBALL_ODFSIZE, 301, 0.5); // get paramters float minSpacing; if(qBallImageSpacing[0]m_NumIt) { MITK_INFO << "itkGibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } - MITK_INFO << "itkGibbsTrackingFilter: steps: " << m_Steps; if (m_CurvatureHardThreshold < mitk::eps) m_CurvatureHardThreshold = 0; - MITK_INFO << "itkGibbsTrackingFilter: curvature threshold: " << m_CurvatureHardThreshold; unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); // setup metropolis hastings sampler MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; if (m_Sampler!=NULL) delete m_Sampler; m_Sampler = new RJMCMC(NULL, 0, workingQballImage, qBallImageSize, qBallImageSpacing, cellsize); // setup energy computer MITK_INFO << "itkGibbsTrackingFilter: setting up Energy-computer"; EnergyComputer encomp(workingQballImage,qBallImageSize,qBallImageSpacing,sinterp,&(m_Sampler->m_ParticleGrid),mask,mask_oversamp_mult, directionMatrix); encomp.setParameters(m_ParticleWeight,m_ParticleWidth,m_ChempotConnection*m_ParticleLength*m_ParticleLength,m_ParticleLength,m_CurvatureHardThreshold,m_InexBalance,m_Chempot2, m_Meanval_sq); m_Sampler->SetEnergyComputer(&encomp); m_Sampler->SetParameters(m_TempStart,singleIts,m_ParticleLength,m_CurvatureHardThreshold,m_ChempotParticle); + + MITK_INFO << "itkGibbsTrackingFilter: Iterations: " << m_NumIt; + MITK_INFO << "itkGibbsTrackingFilter: steps: " << m_Steps; + MITK_INFO << "itkGibbsTrackingFilter: Particle weight: " << m_ParticleWeight; + MITK_INFO << "itkGibbsTrackingFilter: Particle width: " << m_ParticleWidth; + MITK_INFO << "itkGibbsTrackingFilter: Particle length: " << m_ParticleLength; + MITK_INFO << "itkGibbsTrackingFilter: Min. fiber length: " << m_ParticleLength; + MITK_INFO << "itkGibbsTrackingFilter: Start temperature: " << m_TempStart; + MITK_INFO << "itkGibbsTrackingFilter: End temperature: " << m_TempEnd; + MITK_INFO << "itkGibbsTrackingFilter: Energy balance: " << m_InexBalance; + MITK_INFO << "itkGibbsTrackingFilter: Curvature threshold: " << m_CurvatureHardThreshold; + // main loop for( int step = 0; step < m_Steps; step++ ) { if (m_AbortTracking) break; m_CurrentStep = step+1; float temperature = m_TempStart * exp(alpha*(((1.0)*step)/((1.0)*m_Steps))); m_Sampler->SetTemperature(temperature); m_Sampler->Iterate(&m_ProposalAcceptance, &m_NumConnections, &m_NumParticles, &m_AbortTracking); MITK_INFO << "itkGibbsTrackingFilter: proposal acceptance: " << 100*m_ProposalAcceptance << "%"; MITK_INFO << "itkGibbsTrackingFilter: particles: " << m_NumParticles; MITK_INFO << "itkGibbsTrackingFilter: connections: " << m_NumConnections; MITK_INFO << "itkGibbsTrackingFilter: progress: " << 100*(float)step/m_Steps << "%"; if (m_BuildFibers) { int numPoints = m_Sampler->m_ParticleGrid.pcnt; float* points = new float[numPoints*m_Sampler->m_NumAttributes]; m_Sampler->WriteOutParticles(points); BuildFibers(points, numPoints); delete points; m_BuildFibers = false; } } int numPoints = m_Sampler->m_ParticleGrid.pcnt; float* points = new float[numPoints*m_Sampler->m_NumAttributes]; m_Sampler->WriteOutParticles(points); BuildFibers(points, numPoints); delete points; delete sinterp; delete coords; delete ind; delete workingQballImage; m_AbortTracking = true; m_BuildFibers = false; MITK_INFO << "itkGibbsTrackingFilter: done generate data"; } } diff --git a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp new file mode 100644 index 0000000000..e39f3987dd --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp @@ -0,0 +1,348 @@ +/*=================================================================== + +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 __itkStreamlineTrackingFilter_txx +#define __itkStreamlineTrackingFilter_txx + +#include +#include +#include + +#include "itkStreamlineTrackingFilter.h" +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +namespace itk { + +//#define QBALL_RECON_PI M_PI + +template< class TTensorPixelType, class TPDPixelType> +StreamlineTrackingFilter< TTensorPixelType, +TPDPixelType> +::StreamlineTrackingFilter(): + m_FaThreshold(0.2), + m_StepSize(0.2), + m_MaxLength(10000), + m_SeedsPerVoxel(1) +{ + // At least 1 inputs is necessary for a vector image. + // For images added one at a time we need at least six + this->SetNumberOfRequiredInputs( 1 ); +} + +template< class TTensorPixelType, + class TPDPixelType> +void StreamlineTrackingFilter< TTensorPixelType, +TPDPixelType> +::BeforeThreadedGenerateData() +{ + m_FiberPolyData = FiberPolyDataType::New(); + m_Points = vtkPoints::New(); + m_Cells = vtkCellArray::New(); + + typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); + m_ImageSize[0] = inputImage->GetLargestPossibleRegion().GetSize().GetElement(0); + m_ImageSize[1] = inputImage->GetLargestPossibleRegion().GetSize().GetElement(1); + m_ImageSize[2] = inputImage->GetLargestPossibleRegion().GetSize().GetElement(2); + + m_PolyDataContainer = itk::VectorContainer< int, FiberPolyDataType >::New(); + for (int i=0; iGetNumberOfThreads(); i++) + { + FiberPolyDataType poly = FiberPolyDataType::New(); + m_PolyDataContainer->InsertElement(i, poly); + } + + if (m_MaskImage.IsNull()) + { + itk::Vector spacing = inputImage->GetSpacing(); + itk::Point origin = inputImage->GetOrigin(); + itk::Matrix direction = inputImage->GetDirection(); + ImageRegion<3> imageRegion = inputImage->GetLargestPossibleRegion(); + + // initialize crossings image + 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); + } + std::cout << "starting streamline tracking" << std::endl; +} + +template< class TTensorPixelType, + class TPDPixelType> +void StreamlineTrackingFilter< TTensorPixelType, +TPDPixelType> +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int threadId) +{ + FiberPolyDataType poly = m_PolyDataContainer->GetElement(threadId); + vtkSmartPointer Points = vtkPoints::New(); + vtkSmartPointer Cells = vtkCellArray::New(); + + typedef itk::DiffusionTensor3D TensorType; + typedef ImageRegionConstIterator< InputImageType > InputIteratorType; + typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; + typedef typename InputImageType::PixelType InputTensorType; + typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); + + InputIteratorType it(inputImage, outputRegionForThread ); + MaskIteratorType mit(m_MaskImage, outputRegionForThread ); + it.GoToBegin(); + mit.GoToBegin(); + while( !it.IsAtEnd() ) + { + if (mit.Value()==0) + { + ++mit; + ++it; + continue; + } + + typename TensorType::EigenValuesArrayType eigenvalues; + typename TensorType::EigenVectorsMatrixType eigenvectors; + + for (int s=0; s container = vtkSmartPointer::New(); + std::vector< vtkIdType > pointISs; + typename InputImageType::IndexType index = it.GetIndex(); + itk::ContinuousIndex pos; + itk::ContinuousIndex start; + + if (m_SeedsPerVoxel>1) + { + pos[0] = index[0]+(double)(rand()%99-49)/100; + pos[1] = index[1]+(double)(rand()%99-49)/100; + pos[2] = index[2]+(double)(rand()%99-49)/100; + } + else + { + pos[0] = index[0]; + pos[1] = index[1]; + pos[2] = index[2]; + } + start = pos; + + int step = 0; + vnl_vector_fixed dirOld; dirOld.fill(0.0); + // do forward tracking + while (step < m_MaxLength) + { + ++step; + + index[0] = round(pos[0]); + index[1] = round(pos[1]); + index[2] = round(pos[2]); + + if (index[0] < 0 || index[0]>=m_ImageSize[0]) + break; + if (index[1] < 0 || index[1]>=m_ImageSize[1]) + break; + if (index[2] < 0 || index[2]>=m_ImageSize[2]) + break; + + typename InputImageType::PixelType tensor = inputImage->GetPixel(index); + if(tensor.GetTrace()!=0 && tensor.GetFractionalAnisotropy()>m_FaThreshold) + { + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + + int eIndex = 2; + if( (eigenvalues[0] >= eigenvalues[1]) && (eigenvalues[0] >= eigenvalues[2]) ) + eIndex = 0; + else if(eigenvalues[1] >= eigenvalues[2]) + eIndex = 1; + + vnl_vector_fixed dir; + dir[0] = eigenvectors(eIndex, 0); + dir[1] = eigenvectors(eIndex, 1); + dir[2] = eigenvectors(eIndex, 2); + dir.normalize(); + + if (!dirOld.is_zero()) + { + float angle = dot_product(dirOld, dir); + if (angle<0) + dir *= -1; + angle = fabs(dot_product(dirOld, dir)); + if (angle<0.7) + break; + } + dirOld = dir; + + dir *= m_StepSize; + + itk::Point worldPos; + inputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); + + vtkIdType id = Points->InsertNextPoint(worldPos.GetDataPointer()); + pointISs.push_back(id); + counter++; + + pos[0] += dir[0]; + pos[1] += dir[1]; + pos[2] += dir[2]; + } + } + + // insert reverse IDs + while (!pointISs.empty()) + { + container->GetPointIds()->InsertNextId(pointISs.back()); + pointISs.pop_back(); + } + + // do backward tracking + index = it.GetIndex(); + pos = start; + dirOld.fill(0.0); + while (step < m_MaxLength) + { + ++step; + + index[0] = round(pos[0]); + index[1] = round(pos[1]); + index[2] = round(pos[2]); + + if (index[0] < 0 || index[0]>=m_ImageSize[0]) + break; + if (index[1] < 0 || index[1]>=m_ImageSize[1]) + break; + if (index[2] < 0 || index[2]>=m_ImageSize[2]) + break; + + typename InputImageType::PixelType tensor = inputImage->GetPixel(index); + if(tensor.GetTrace()!=0 && tensor.GetFractionalAnisotropy()>m_FaThreshold) + { + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + + int eIndex = 2; + if( (eigenvalues[0] >= eigenvalues[1]) && (eigenvalues[0] >= eigenvalues[2]) ) + eIndex = 0; + else if(eigenvalues[1] >= eigenvalues[2]) + eIndex = 1; + + vnl_vector_fixed dir; + dir[0] = eigenvectors(eIndex, 0); + dir[1] = eigenvectors(eIndex, 1); + dir[2] = eigenvectors(eIndex, 2); + dir.normalize(); + dir *= -1; // reverse direction + + if (!dirOld.is_zero()) + { + float angle = dot_product(dirOld, dir); + if (angle<0) + dir *= -1; + angle = fabs(dot_product(dirOld, dir)); + if (angle<0.7) + break; + } + dirOld = dir; + + dir *= m_StepSize; + + itk::Point worldPos; + inputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); + + vtkIdType id = Points->InsertNextPoint(worldPos.GetDataPointer()); + container->GetPointIds()->InsertNextId(id); + counter++; + + pos[0] += dir[0]; + pos[1] += dir[1]; + pos[2] += dir[2]; + } + } + + if (counter>0) + Cells->InsertNextCell(container); + } + ++mit; + ++it; + } + + poly->SetPoints(Points); + poly->SetLines(Cells); + + std::cout << "Thread " << threadId << " finished tracking" << std::endl; +} + +template< class TTensorPixelType, + class TPDPixelType> +vtkSmartPointer< vtkPolyData > StreamlineTrackingFilter< TTensorPixelType, +TPDPixelType> +::AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2) +{ + vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); + vtkSmartPointer vNewLines = poly1->GetLines(); + vtkSmartPointer vNewPoints = poly1->GetPoints(); + + vtkSmartPointer vLines = poly2->GetLines(); + vLines->InitTraversal(); + for( int i=0; iGetNumberOfCells(); i++ ) + { + vtkIdType numPoints(0); + vtkIdType* points(NULL); + vLines->GetNextCell ( numPoints, points ); + + vtkSmartPointer container = vtkSmartPointer::New(); + for( int j=0; jInsertNextPoint(poly2->GetPoint(points[j])); + container->GetPointIds()->InsertNextId(id); + } + vNewLines->InsertNextCell(container); + } + + // initialize polydata + vNewPolyData->SetPoints(vNewPoints); + vNewPolyData->SetLines(vNewLines); + + return vNewPolyData; +} +template< class TTensorPixelType, + class TPDPixelType> +void StreamlineTrackingFilter< TTensorPixelType, +TPDPixelType> +::AfterThreadedGenerateData() +{ + MITK_INFO << "Generating polydata "; + m_FiberPolyData = m_PolyDataContainer->GetElement(0); + for (int i=1; iGetNumberOfThreads(); i++) + { + m_FiberPolyData = AddPolyData(m_FiberPolyData, m_PolyDataContainer->GetElement(i)); + } + MITK_INFO << "done"; +} + +template< class TTensorPixelType, + class TPDPixelType> +void StreamlineTrackingFilter< TTensorPixelType, +TPDPixelType> +::PrintSelf(std::ostream& os, Indent indent) const +{ +} + +} +#endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx diff --git a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h new file mode 100644 index 0000000000..d6aaa4b759 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h @@ -0,0 +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. + +===================================================================*/ + +/*=================================================================== + +This file is based heavily on a corresponding ITK filter. + +===================================================================*/ +#ifndef __itkStreamlineTrackingFilter_h_ +#define __itkStreamlineTrackingFilter_h_ + +#include "MitkDiffusionImagingExports.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace itk{ + /** \class StreamlineTrackingFilter + */ + + template< class TTensorPixelType, class TPDPixelType=double> + class StreamlineTrackingFilter : + public ImageToImageFilter< Image< DiffusionTensor3D, 3 >, + Image< Vector< TPDPixelType, 3 >, 3 > > + { + + public: + + typedef StreamlineTrackingFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< Vector< TPDPixelType, 3 >, 3 > > Superclass; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(StreamlineTrackingFilter, ImageToImageFilter) + + typedef TTensorPixelType TensorComponentType; + typedef TPDPixelType DirectionPixelType; + typedef typename Superclass::InputImageType InputImageType; + typedef typename Superclass::OutputImageType OutputImageType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + typedef itk::Image ItkUcharImgType; + + typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; + + itkGetMacro( FiberPolyData, FiberPolyDataType ) + itkSetMacro( MaskImage, ItkUcharImgType::Pointer) + itkSetMacro( SeedsPerVoxel, int) + itkSetMacro( FaThreshold, float) + + protected: + StreamlineTrackingFilter(); + ~StreamlineTrackingFilter() {} + void PrintSelf(std::ostream& os, Indent indent) const; + + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int threadId); + void AfterThreadedGenerateData(); + + FiberPolyDataType AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2); + + FiberPolyDataType m_FiberPolyData; + vtkSmartPointer m_Points; + vtkSmartPointer m_Cells; + float m_FaThreshold; + float m_StepSize; + int m_MaxLength; + int m_SeedsPerVoxel; + int m_ImageSize[3]; + ItkUcharImgType::Pointer m_MaskImage; + + itk::VectorContainer< int, FiberPolyDataType >::Pointer m_PolyDataContainer; + + private: + + }; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkStreamlineTrackingFilter.cpp" +#endif + +#endif //__itkStreamlineTrackingFilter_h_ + diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index 73b96d2312..cfadef37e8 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,227 +1,228 @@ set(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionImagingObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures -> Tbss IODataStructures/TbssImages/mitkTbssImageSource.cpp IODataStructures/TbssImages/mitkTbssRoiImageSource.cpp IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssImageIOFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageIOFactory.cpp IODataStructures/TbssImages/mitkTbssImage.cpp IODataStructures/TbssImages/mitkTbssRoiImage.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriterFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriterFactory.cpp IODataStructures/TbssImages/mitkTbssImporter.cpp # DataStructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp Rendering/mitkTbssImageMapper.cpp Rendering/mitkPlanarCircleMapper3D.cpp Rendering/mitkPlanarPolygonMapper3D.cpp Rendering/mitkConnectomicsNetworkMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/mitkTractAnalyzer.cpp # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.cpp Algorithms/Connectomics/mitkConnectomicsHistogramBase.cpp Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.cpp Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.cpp Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.cpp Algorithms/Connectomics/mitkConnectomicsHistogramCache.cpp Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.cpp # Tractography Tractography/itkStochasticTractographyFilter.h # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection mitkDiffusionFunctionCollection.h - + # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkTbssImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h Rendering/mitkPlanarCircleMapper3D.h Rendering/mitkPlanarPolygonMapper3D.h Rendering/mitkConnectomicsNetworkMapper3D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h Reconstruction/itkRegularizedIVIMReconstructionFilter.h Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/TbssImages/mitkTbssImporter.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h # Datastructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.h # Tractography Tractography/itkGibbsTrackingFilter.h Tractography/itkStochasticTractographyFilter.h + Tractography/itkStreamlineTrackingFilter.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkB0ImageExtractionToSeparateImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkSkeletonizationFilter.h Algorithms/itkReduceDirectionGradientsFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.h Algorithms/Connectomics/mitkConnectomicsHistogramBase.h Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.h Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.h Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.h Algorithms/Connectomics/mitkConnectomicsHistogramCache.h Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.h Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h - + ) set( TOOL_FILES ) if(WIN32) endif(WIN32) #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/files.cmake b/Plugins/org.mitk.gui.qt.diffusionimaging/files.cmake index 68930fafcd..dcdcef089d 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/files.cmake +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/files.cmake @@ -1,121 +1,124 @@ set(SRC_CPP_FILES QmitkODFDetailsWidget.cpp QmitkODFRenderWidget.cpp QmitkPartialVolumeAnalysisWidget.cpp QmitkIVIMWidget.cpp QmitkTbssRoiAnalysisWidget.cpp QmitkResidualAnalysisWidget.cpp QmitkResidualViewWidget.cpp ) set(INTERNAL_CPP_FILES mitkPluginActivator.cpp QmitkQBallReconstructionView.cpp QmitkPreprocessingView.cpp QmitkDiffusionDicomImportView.cpp QmitkDiffusionQuantificationView.cpp QmitkTensorReconstructionView.cpp QmitkDiffusionImagingPublicPerspective.cpp QmitkControlVisualizationPropertiesView.cpp QmitkODFDetailsView.cpp QmitkGibbsTrackingView.cpp QmitkStochasticFiberTrackingView.cpp + QmitkStreamlineTrackingView.cpp QmitkFiberProcessingView.cpp QmitkFiberBundleDeveloperView.cpp QmitkPartialVolumeAnalysisView.cpp QmitkIVIMView.cpp QmitkTractbasedSpatialStatisticsView.cpp QmitkTbssTableModel.cpp QmitkTbssMetaTableModel.cpp QmitkTbssSkeletonizationView.cpp Connectomics/QmitkBrainNetworkAnalysisView.cpp Connectomics/QmitkNetworkHistogramCanvas.cpp ) set(UI_FILES src/internal/QmitkQBallReconstructionViewControls.ui src/internal/QmitkPreprocessingViewControls.ui src/internal/QmitkDiffusionDicomImportViewControls.ui src/internal/QmitkDiffusionQuantificationViewControls.ui src/internal/QmitkTensorReconstructionViewControls.ui src/internal/QmitkControlVisualizationPropertiesViewControls.ui src/internal/QmitkODFDetailsViewControls.ui src/internal/QmitkGibbsTrackingViewControls.ui src/internal/QmitkStochasticFiberTrackingViewControls.ui + src/internal/QmitkStreamlineTrackingViewControls.ui src/internal/QmitkFiberProcessingViewControls.ui src/internal/QmitkFiberBundleDeveloperViewControls.ui src/internal/QmitkPartialVolumeAnalysisViewControls.ui src/internal/QmitkIVIMViewControls.ui src/internal/QmitkTractbasedSpatialStatisticsViewControls.ui src/internal/QmitkTbssSkeletonizationViewControls.ui src/internal/Connectomics/QmitkBrainNetworkAnalysisViewControls.ui ) set(MOC_H_FILES src/internal/mitkPluginActivator.h src/internal/QmitkQBallReconstructionView.h src/internal/QmitkPreprocessingView.h src/internal/QmitkDiffusionDicomImportView.h src/internal/QmitkDiffusionImagingPublicPerspective.h src/internal/QmitkDiffusionQuantificationView.h src/internal/QmitkTensorReconstructionView.h src/internal/QmitkControlVisualizationPropertiesView.h src/internal/QmitkODFDetailsView.h src/QmitkODFRenderWidget.h src/QmitkODFDetailsWidget.h src/internal/QmitkGibbsTrackingView.h src/internal/QmitkStochasticFiberTrackingView.h + src/internal/QmitkStreamlineTrackingView.h src/internal/QmitkFiberProcessingView.h src/internal/QmitkFiberBundleDeveloperView.h src/internal/QmitkPartialVolumeAnalysisView.h src/QmitkPartialVolumeAnalysisWidget.h src/internal/QmitkIVIMView.h src/internal/QmitkTractbasedSpatialStatisticsView.h src/internal/QmitkTbssSkeletonizationView.h src/QmitkTbssRoiAnalysisWidget.h src/QmitkResidualAnalysisWidget.h src/QmitkResidualViewWidget.h src/internal/Connectomics/QmitkBrainNetworkAnalysisView.h src/internal/Connectomics/QmitkNetworkHistogramCanvas.h ) set(CACHED_RESOURCE_FILES # list of resource files which can be used by the plug-in # system without loading the plug-ins shared library, # for example the icon used in the menu and tabs for the # plug-in views in the workbench plugin.xml resources/preprocessing.png resources/dwiimport.png resources/quantification.png resources/reconodf.png resources/recontensor.png resources/vizControls.png resources/OdfDetails.png resources/GibbsTracking.png resources/FiberBundleOperations.png resources/PartialVolumeAnalysis_24.png resources/IVIM_48.png resources/stochFB.png resources/tbss.png resources/QmitkBrainNetworkAnalysisViewIcon_48.png resources/arrow.png ) set(QRC_FILES # uncomment the following line if you want to use Qt resources resources/QmitkDiffusionImaging.qrc #resources/QmitkTractbasedSpatialStatisticsView.qrc ) set(CPP_FILES ) foreach(file ${SRC_CPP_FILES}) set(CPP_FILES ${CPP_FILES} src/${file}) endforeach(file ${SRC_CPP_FILES}) foreach(file ${INTERNAL_CPP_FILES}) set(CPP_FILES ${CPP_FILES} src/internal/${file}) endforeach(file ${INTERNAL_CPP_FILES}) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/plugin.xml b/Plugins/org.mitk.gui.qt.diffusionimaging/plugin.xml index e8b9846c4b..76277778da 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/plugin.xml +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/plugin.xml @@ -1,125 +1,132 @@ + + + + diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp index 397b9cef81..2fb8a396fd 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp @@ -1,1766 +1,1797 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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 "QmitkFiberProcessingView.h" #include // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include +#include // ITK #include #include #include #include #include #include #include #include const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberProcessingView::QmitkFiberProcessingView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) , m_EllipseCounter(0) , m_PolygonCounter(0) , m_UpsamplingFactor(5) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { } void QmitkFiberProcessingView::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::QmitkFiberProcessingViewControls; m_Controls->setupUi( parent ); m_Controls->doExtractFibersButton->setDisabled(true); m_Controls->PFCompoANDButton->setDisabled(true); m_Controls->PFCompoORButton->setDisabled(true); m_Controls->PFCompoNOTButton->setDisabled(true); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->m_RectangleButton->setVisible(false); connect( m_Controls->doExtractFibersButton, SIGNAL(clicked()), this, SLOT(DoFiberExtraction()) ); connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) ); connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) ); connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) ); connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); - connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); - connect(m_Controls->m_Extract3dButton, SIGNAL(clicked()), this, SLOT(Extract3d())); connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) ); connect( m_Controls->m_ResampleFibersButton, SIGNAL(clicked()), this, SLOT(ResampleSelectedBundles()) ); connect(m_Controls->m_FaColorFibersButton, SIGNAL(clicked()), this, SLOT(DoFaColorCoding())); + connect( m_Controls->m_PruneFibersButton, SIGNAL(clicked()), this, SLOT(PruneBundle()) ); + connect( m_Controls->m_MirrorFibersButton, SIGNAL(clicked()), this, SLOT(MirrorFibers()) ); + } } void QmitkFiberProcessingView::Extract3d() { std::vector nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(); mitk::Surface::Pointer roi = mitk::Surface::New(); bool fibB = false; bool roiB = false; for (int i=0; i(nodes.at(i)->GetData())) { fib = dynamic_cast(nodes.at(i)->GetData()); fibB = true; } else if (dynamic_cast(nodes.at(i)->GetData())) { roi = dynamic_cast(nodes.at(i)->GetData()); roiB = true; } } if (!fibB) return; if (!roiB) return; vtkSmartPointer polyRoi = roi->GetVtkPolyData(); vtkSmartPointer polyFib = fib->GetFiberPolyData(); vtkSmartPointer selectEnclosedPoints = vtkSmartPointer::New(); selectEnclosedPoints->SetInput(polyFib); selectEnclosedPoints->SetSurface(polyRoi); selectEnclosedPoints->Update(); vtkSmartPointer newPoly = vtkSmartPointer::New(); vtkSmartPointer newCellArray = vtkSmartPointer::New(); vtkSmartPointer newPoints = vtkSmartPointer::New(); vtkSmartPointer newPolyComplement = vtkSmartPointer::New(); vtkSmartPointer newCellArrayComplement = vtkSmartPointer::New(); vtkSmartPointer newPointsComplement = vtkSmartPointer::New(); vtkSmartPointer vLines = polyFib->GetLines(); vLines->InitTraversal(); int numberOfLines = vLines->GetNumberOfCells(); // each line for (int j=0; jGetNextCell ( numPoints, points ); bool isPassing = false; // each point of this line for (int k=0; kIsInside(points[k])) { isPassing = true; // fill new polydata vtkSmartPointer container = vtkSmartPointer::New(); for (int k=0; kGetPoint(points[k]); vtkIdType pointId = newPoints->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); } newCellArray->InsertNextCell(container); break; } } if (!isPassing) { vtkSmartPointer container = vtkSmartPointer::New(); for (int k=0; kGetPoint(points[k]); vtkIdType pointId = newPointsComplement->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); } newCellArrayComplement->InsertNextCell(container); } } newPoly->SetPoints(newPoints); newPoly->SetLines(newCellArray); mitk::FiberBundleX::Pointer fb = mitk::FiberBundleX::New(newPoly); DataNode::Pointer newNode = DataNode::New(); newNode->SetData(fb); newNode->SetName("passing surface"); GetDefaultDataStorage()->Add(newNode); newPolyComplement->SetPoints(newPointsComplement); newPolyComplement->SetLines(newCellArrayComplement); mitk::FiberBundleX::Pointer fbComplement = mitk::FiberBundleX::New(newPolyComplement); DataNode::Pointer newNodeComplement = DataNode::New(); newNodeComplement->SetData(fbComplement); newNodeComplement->SetName("not passing surface"); GetDefaultDataStorage()->Add(newNodeComplement); } void QmitkFiberProcessingView::GenerateRoiImage(){ if (m_SelectedPF.empty()) return; mitk::Geometry3D::Pointer geometry; if (!m_SelectedFB.empty()) { mitk::FiberBundleX::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); geometry = fib->GetGeometry(); } else return; mitk::Vector3D spacing = geometry->GetSpacing(); spacing /= m_UpsamplingFactor; mitk::Point3D newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); itk::Matrix direction; itk::ImageRegion<3> imageRegion; for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j]; imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); m_PlanarFigureImage = itkUCharImageType::New(); m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin m_PlanarFigureImage->SetDirection( direction ); // Set the image direction m_PlanarFigureImage->SetRegions( imageRegion ); m_PlanarFigureImage->Allocate(); m_PlanarFigureImage->FillBuffer( 0 ); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); for (int i=0; iInitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); node->SetData(tmpImage); node->SetName("ROI Image"); this->GetDefaultDataStorage()->Add(node); } void QmitkFiberProcessingView::CompositeExtraction(mitk::DataNode::Pointer node, mitk::Image* image) { if (dynamic_cast(node.GetPointer()->GetData()) && !dynamic_cast(node.GetPointer()->GetData())) { m_PlanarFigure = dynamic_cast(node.GetPointer()->GetData()); AccessFixedDimensionByItk_2( image, InternalReorientImagePlane, 3, m_PlanarFigure->GetGeometry(), -1); // itk::Image< unsigned char, 3 >::Pointer outimage = itk::Image< unsigned char, 3 >::New(); // outimage->SetSpacing( m_PlanarFigure->GetGeometry()->GetSpacing()/m_UpsamplingFactor ); // Set the image spacing // mitk::Point3D origin = m_PlanarFigure->GetGeometry()->GetOrigin(); // mitk::Point3D indexOrigin; // m_PlanarFigure->GetGeometry()->WorldToIndex(origin, indexOrigin); // indexOrigin[0] = indexOrigin[0] - .5 * (1.0-1.0/m_UpsamplingFactor); // indexOrigin[1] = indexOrigin[1] - .5 * (1.0-1.0/m_UpsamplingFactor); // indexOrigin[2] = indexOrigin[2] - .5 * (1.0-1.0/m_UpsamplingFactor); // mitk::Point3D newOrigin; // m_PlanarFigure->GetGeometry()->IndexToWorld(indexOrigin, newOrigin); // outimage->SetOrigin( newOrigin ); // Set the image origin // itk::Matrix matrix; // for (int i=0; i<3; i++) // for (int j=0; j<3; j++) // matrix[j][i] = m_PlanarFigure->GetGeometry()->GetMatrixColumn(i)[j]/m_PlanarFigure->GetGeometry()->GetSpacing().GetElement(i); // outimage->SetDirection( matrix ); // Set the image direction // itk::ImageRegion<3> upsampledRegion; // upsampledRegion.SetSize(0, m_PlanarFigure->GetGeometry()->GetParametricExtentInMM(0)/m_PlanarFigure->GetGeometry()->GetSpacing()[0]); // upsampledRegion.SetSize(1, m_PlanarFigure->GetGeometry()->GetParametricExtentInMM(1)/m_PlanarFigure->GetGeometry()->GetSpacing()[1]); // upsampledRegion.SetSize(2, 1); // typename itk::Image< unsigned char, 3 >::RegionType::SizeType upsampledSize = upsampledRegion.GetSize(); // for (unsigned int n = 0; n < 2; n++) // { // upsampledSize[n] = upsampledSize[n] * m_UpsamplingFactor; // } // upsampledRegion.SetSize( upsampledSize ); // outimage->SetRegions( upsampledRegion ); // outimage->Allocate(); // this->m_InternalImage = mitk::Image::New(); // this->m_InternalImage->InitializeByItk( outimage.GetPointer() ); // this->m_InternalImage->SetVolume( outimage->GetBufferPointer() ); AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2, node->GetName() ); } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::Geometry3D* planegeo3D, int additionalIndex ) { MITK_INFO << "InternalReorientImagePlane() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > FloatImageType; typedef itk::ResampleImageFilter ResamplerType; typename ResamplerType::Pointer resampler = ResamplerType::New(); mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); float upsamp = m_UpsamplingFactor; float gausssigma = 0.5; // Spacing typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); spacing[0] = image->GetSpacing()[0] / upsamp; spacing[1] = image->GetSpacing()[1] / upsamp; spacing[2] = image->GetSpacing()[2]; resampler->SetOutputSpacing( spacing ); // Size typename ResamplerType::SizeType size; size[0] = planegeo->GetParametricExtentInMM(0) / spacing[0]; size[1] = planegeo->GetParametricExtentInMM(1) / spacing[1]; size[2] = 1; resampler->SetSize( size ); // Origin typename mitk::Point3D orig = planegeo->GetOrigin(); typename mitk::Point3D corrorig; planegeo3D->WorldToIndex(orig,corrorig); corrorig[0] += 0.5/upsamp; corrorig[1] += 0.5/upsamp; corrorig[2] += 0; planegeo3D->IndexToWorld(corrorig,corrorig); resampler->SetOutputOrigin(corrorig ); // Direction typename ResamplerType::DirectionType direction; typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); for(int c=0; cSetOutputDirection( direction ); // Gaussian interpolation if(gausssigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) { sigma[d] = gausssigma * image->GetSpacing()[d]; } double alpha = 2.0; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); interpolator->SetInputImage( image ); interpolator->SetParameters( sigma, alpha ); resampler->SetInterpolator( interpolator ); } else { // typedef typename itk::BSplineInterpolateImageFunction // InterpolatorType; typedef typename itk::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } // Other resampling options resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); MITK_INFO << "Resampling requested image plane ... "; resampler->Update(); MITK_INFO << " ... done"; if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string nodeName ) { MITK_INFO << "InternalCalculateMaskFromPlanarFigure() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, itkUCharImageType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". itkUCharImageType::Pointer newMaskImage = itkUCharImageType::New(); newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); newMaskImage->Allocate(); newMaskImage->FillBuffer( 1 ); // Generate VTK polygon from (closed) PlanarFigure polyline // (The polyline points are shifted by -0.5 in z-direction to make sure // that the extrusion filter, which afterwards elevates all points by +0.5 // in z-direction, creates a 3D object which is cut by the the plane z=0) const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const Geometry3D *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); vtkPolyData *polyline = vtkPolyData::New(); polyline->Allocate( 1, 1 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // Create VTK polydata object of polyline contour vtkPoints *points = vtkPoints::New(); PlanarFigure::PolyLineType::const_iterator it; std::vector indices; unsigned int numberOfPoints = 0; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image Point2D point2D = it->Point; planarFigureGeometry2D->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigureGeometry2D->IndexToWorld(point2D, point2D); planarFigureGeometry2D->Map( point2D, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { float bounds[2] = {0,0}; bounds[0] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0); bounds[1] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1); imageGeometry3D->WorldToIndex( point3D, point3D ); // if (point3D[i0]<0) // point3D[i0] = 0.5; // else if (point3D[i0]>bounds[0]) // point3D[i0] = bounds[0]-0.5; // if (point3D[i1]<0) // point3D[i1] = 0.5; // else if (point3D[i1]>bounds[1]) // point3D[i1] = bounds[1]-0.5; if (point3D[i0]<0) point3D[i0] = 0.0; else if (point3D[i0]>bounds[0]) point3D[i0] = bounds[0]-0.001; if (point3D[i1]<0) point3D[i1] = 0.0; else if (point3D[i1]>bounds[1]) point3D[i1] = bounds[1]-0.001; points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } else { imageGeometry3D->WorldToIndex( point3D, point3D ); // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } } polyline->SetPoints( points ); points->Delete(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( vtkIdType i = 0; i < numberOfPoints; ++i ) { ptIds[i] = i; } polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->SetInput( polyline ); extrudeFilter->SetScaleFactor( 1 ); extrudeFilter->SetExtrusionTypeToNormalExtrusion(); extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); // Make a stencil from the extruded polygon vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); polyDataToImageStencil->SetInput( extrudeFilter->GetOutput() ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< itkUCharImageType > ImageImportType; typedef itk::VTKImageExport< itkUCharImageType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( newMaskImage ); vtkImageImport *vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); vtkImporter->Update(); // Apply the generated image stencil to the input image vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencil( polyDataToImageStencil->GetOutput() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkImageExport *vtkExporter = vtkImageExport::New(); vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // calculate cropping bounding box m_InternalImageMask3D = itkImporter->GetOutput(); m_InternalImageMask3D->SetDirection(image->GetDirection()); itk::ImageRegionConstIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itk::ImageRegionIterator itimage(image, image->GetLargestPossibleRegion()); itmask = itmask.Begin(); itimage = itimage.Begin(); typename ImageType::SizeType lowersize = {{9999999999,9999999999,9999999999}}; typename ImageType::SizeType uppersize = {{0,0,0}}; while( !itmask.IsAtEnd() ) { if(itmask.Get() == 0) { itimage.Set(0); } else { typename ImageType::IndexType index = itimage.GetIndex(); typename ImageType::SizeType signedindex; signedindex[0] = index[0]; signedindex[1] = index[1]; signedindex[2] = index[2]; lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; } ++itmask; ++itimage; } typename ImageType::IndexType index; index[0] = lowersize[0]; index[1] = lowersize[1]; index[2] = lowersize[2]; typename ImageType::SizeType size; size[0] = uppersize[0] - lowersize[0] + 1; size[1] = uppersize[1] - lowersize[1] + 1; size[2] = uppersize[2] - lowersize[2] + 1; itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size); // crop internal mask typedef itk::RegionOfInterestImageFilter< itkUCharImageType, itkUCharImageType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(cropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); Image::Pointer tmpImage2 = Image::New(); tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); const Geometry3D *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); const Geometry3D *intImageGeometry3D = tmpImage->GetGeometry( 0 ); typedef itk::ImageRegionIteratorWithIndex IteratorType; IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion()); imageIterator.GoToBegin(); while ( !imageIterator.IsAtEnd() ) { unsigned char val = imageIterator.Value(); if (val>0) { itk::Index<3> index = imageIterator.GetIndex(); Point3D point; point[0] = index[0]; point[1] = index[1]; point[2] = index[2]; intImageGeometry3D->IndexToWorld(point, point); pfImageGeometry3D->WorldToIndex(point, point); point[i0] += 0.5; point[i1] += 0.5; index[0] = point[0]; index[1] = point[1]; index[2] = point[2]; if (pfImageGeometry3D->IsIndexInside(index)) m_PlanarFigureImage->SetPixel(index, 1); } ++imageIterator; } // Clean up VTK objects polyline->Delete(); extrudeFilter->Delete(); polyDataToImageStencil->Delete(); vtkImporter->Delete(); imageStencilFilter->Delete(); //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? delete[] ptIds; } void QmitkFiberProcessingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkFiberProcessingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } /* OnSelectionChanged is registered to SelectionService, therefore no need to implement SelectionService Listener explicitly */ void QmitkFiberProcessingView::UpdateGui() { // are fiber bundles selected? if ( m_SelectedFB.empty() ) { m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); m_Controls->m_ProcessFiberBundleButton->setEnabled(false); m_Controls->doExtractFibersButton->setEnabled(false); m_Controls->m_Extract3dButton->setEnabled(false); m_Controls->m_ResampleFibersButton->setEnabled(false); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->m_FaColorFibersButton->setEnabled(false); + m_Controls->m_PruneFibersButton->setEnabled(false); + m_Controls->m_MirrorFibersButton->setEnabled(false); } else { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_ProcessFiberBundleButton->setEnabled(true); m_Controls->m_ResampleFibersButton->setEnabled(true); + m_Controls->m_PruneFibersButton->setEnabled(true); + m_Controls->m_MirrorFibersButton->setEnabled(true); if (m_Surfaces.size()>0) m_Controls->m_Extract3dButton->setEnabled(true); // one bundle and one planar figure needed to extract fibers if (!m_SelectedPF.empty()) m_Controls->doExtractFibersButton->setEnabled(true); // more than two bundles needed to join/subtract if (m_SelectedFB.size() > 1) { m_Controls->m_JoinBundles->setEnabled(true); m_Controls->m_SubstractBundles->setEnabled(true); } else { m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); } if (m_SelectedImage.IsNotNull()) m_Controls->m_FaColorFibersButton->setEnabled(true); } // are planar figures selected? if ( m_SelectedPF.empty() ) { m_Controls->doExtractFibersButton->setEnabled(false); m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(false); m_Controls->m_GenerateRoiImage->setEnabled(false); } else { if ( !m_SelectedFB.empty() ) m_Controls->m_GenerateRoiImage->setEnabled(true); else m_Controls->m_GenerateRoiImage->setEnabled(false); if (m_SelectedPF.size() > 1) { m_Controls->PFCompoANDButton->setEnabled(true); m_Controls->PFCompoORButton->setEnabled(true); m_Controls->PFCompoNOTButton->setEnabled(false); } else { m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(true); } } } void QmitkFiberProcessingView::OnSelectionChanged( std::vector nodes ) { if ( !this->IsVisible() ) return; //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection m_SelectedFB.clear(); m_SelectedPF.clear(); m_Surfaces.clear(); m_SelectedImage = NULL; 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(node); else if (dynamic_cast(node->GetData())) m_SelectedPF.push_back(node); else if (dynamic_cast(node->GetData())) m_SelectedImage = dynamic_cast(node->GetData()); else if (dynamic_cast(node->GetData())) m_Surfaces.push_back(dynamic_cast(node->GetData())); } UpdateGui(); GenerateStats(); } void QmitkFiberProcessingView::OnDrawPolygon() { // bool checked = m_Controls->m_PolygonButton->isChecked(); // if(!this->AssertDrawingIsPossible(checked)) // return; mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->ClosedOn(); this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); MITK_INFO << "PlanarPolygon created ..."; mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDefaultDataStorage()->GetAll(); mitk::DataNode* node = 0; mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; mitk::PlanarFigure* figureP = 0; for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() ; it++) { node = const_cast(it->Value().GetPointer()); figureP = dynamic_cast(node->GetData()); if(figureP) { figureInteractor = dynamic_cast(node->GetInteractor()); if(figureInteractor.IsNull()) figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); } } } void QmitkFiberProcessingView::OnDrawCircle() { //bool checked = m_Controls->m_CircleButton->isChecked(); //if(!this->AssertDrawingIsPossible(checked)) // return; mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_EllipseCounter)); this->GetDataStorage()->Modified(); MITK_INFO << "PlanarCircle created ..."; //call mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDefaultDataStorage()->GetAll(); mitk::DataNode* node = 0; mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; mitk::PlanarFigure* figureP = 0; for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() ; it++) { node = const_cast(it->Value().GetPointer()); figureP = dynamic_cast(node->GetData()); if(figureP) { figureInteractor = dynamic_cast(node->GetInteractor()); if(figureInteractor.IsNull()) figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); } } } void QmitkFiberProcessingView::Activated() { MITK_INFO << "FB OPerations ACTIVATED()"; /* mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDefaultDataStorage()->GetAll(); mitk::DataNode* node = 0; mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; mitk::PlanarFigure* figure = 0; for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() ; it++) { node = const_cast(it->Value().GetPointer()); figure = dynamic_cast(node->GetData()); if(figure) { figureInteractor = dynamic_cast(node->GetInteractor()); if(figureInteractor.IsNull()) figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); } } */ } void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey, mitk::BaseProperty *property ) { // initialize figure's geometry with empty geometry mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); figure->SetGeometry2D( emptygeometry ); //set desired data to DataNode where Planarfigure is stored mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetName(name.toStdString()); newNode->SetData(figure); newNode->AddProperty( "planarfigure.default.line.color", mitk::ColorProperty::New(1.0,0.0,0.0)); newNode->AddProperty( "planarfigure.line.width", mitk::FloatProperty::New(2.0)); newNode->AddProperty( "planarfigure.drawshadow", mitk::BoolProperty::New(true)); newNode->AddProperty( "selected", mitk::BoolProperty::New(true) ); newNode->AddProperty( "planarfigure.ishovering", mitk::BoolProperty::New(true) ); newNode->AddProperty( "planarfigure.drawoutline", mitk::BoolProperty::New(true) ); newNode->AddProperty( "planarfigure.drawquantities", mitk::BoolProperty::New(false) ); newNode->AddProperty( "planarfigure.drawshadow", mitk::BoolProperty::New(true) ); newNode->AddProperty( "planarfigure.line.width", mitk::FloatProperty::New(3.0) ); newNode->AddProperty( "planarfigure.shadow.widthmodifier", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.outline.width", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.helperline.width", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.default.line.color", mitk::ColorProperty::New(1.0,1.0,1.0) ); newNode->AddProperty( "planarfigure.default.line.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.default.outline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.default.outline.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.default.helperline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.default.helperline.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.default.markerline.color", mitk::ColorProperty::New(0.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.default.markerline.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.default.marker.color", mitk::ColorProperty::New(1.0,1.0,1.0) ); newNode->AddProperty( "planarfigure.default.marker.opacity",mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.hover.line.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.hover.line.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.hover.outline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.hover.outline.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.hover.helperline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.hover.helperline.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.hover.markerline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.hover.markerline.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.hover.marker.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.hover.marker.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.selected.line.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.selected.line.opacity",mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.selected.outline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.selected.outline.opacity", mitk::FloatProperty::New(2.0)); newNode->AddProperty( "planarfigure.selected.helperline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.selected.helperline.opacity",mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.selected.markerline.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.selected.markerline.opacity", mitk::FloatProperty::New(2.0) ); newNode->AddProperty( "planarfigure.selected.marker.color", mitk::ColorProperty::New(1.0,0.0,0.0) ); newNode->AddProperty( "planarfigure.selected.marker.opacity",mitk::FloatProperty::New(2.0)); // figure drawn on the topmost layer / image newNode->SetColor(1.0,1.0,1.0); newNode->SetOpacity(0.8); GetDataStorage()->Add(newNode ); std::vector selectedNodes = GetDataManagerSelection(); for(unsigned int i = 0; i < selectedNodes.size(); i++) { selectedNodes[i]->SetSelected(false); } newNode->SetSelected(true); } void QmitkFiberProcessingView::DoFiberExtraction() { if ( m_SelectedFB.empty() ){ QMessageBox::information( NULL, "Warning", "No fibe bundle selected!"); MITK_WARN("QmitkFiberProcessingView") << "no fibe bundle selected"; return; } for (int i=0; i(m_SelectedFB.at(i)->GetData()); mitk::PlanarFigure::Pointer roi = dynamic_cast (m_SelectedPF.at(0)->GetData()); mitk::FiberBundleX::Pointer extFB = fib->ExtractFiberSubset(roi); if (extFB->GetNumFibers()<=0) continue; mitk::DataNode::Pointer node; node = mitk::DataNode::New(); node->SetData(extFB); QString name(m_SelectedFB.at(i)->GetName().c_str()); name += "_"; name += m_SelectedPF.at(0)->GetName().c_str(); node->SetName(name.toStdString()); GetDataStorage()->Add(node); m_SelectedFB.at(i)->SetVisibility(false); } } void QmitkFiberProcessingView::GenerateAndComposite() { mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); mitk::PlaneGeometry* currentGeometry2D = dynamic_cast( const_cast(GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderer()->GetCurrentWorldGeometry2D())); PFCAnd->SetGeometry2D(currentGeometry2D); PFCAnd->setOperationType(mitk::PFCOMPOSITION_AND_OPERATION); for( std::vector::iterator it = m_SelectedPF.begin(); it != m_SelectedPF.end(); ++it ) { mitk::DataNode::Pointer nodePF = *it; mitk::PlanarFigure::Pointer tmpPF = dynamic_cast( nodePF->GetData() ); PFCAnd->addPlanarFigure( tmpPF ); PFCAnd->addDataNode( nodePF ); PFCAnd->setDisplayName("AND_COMPO"); } AddCompositeToDatastorage(PFCAnd, NULL); } void QmitkFiberProcessingView::GenerateOrComposite() { mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); mitk::PlaneGeometry* currentGeometry2D = dynamic_cast( const_cast(GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderer()->GetCurrentWorldGeometry2D())); PFCOr->SetGeometry2D(currentGeometry2D); PFCOr->setOperationType(mitk::PFCOMPOSITION_OR_OPERATION); for( std::vector::iterator it = m_SelectedPF.begin(); it != m_SelectedPF.end(); ++it ) { mitk::DataNode::Pointer nodePF = *it; mitk::PlanarFigure::Pointer tmpPF = dynamic_cast( nodePF->GetData() ); PFCOr->addPlanarFigure( tmpPF ); PFCOr->addDataNode( nodePF ); PFCOr->setDisplayName("OR_COMPO"); } AddCompositeToDatastorage(PFCOr, NULL); } void QmitkFiberProcessingView::GenerateNotComposite() { mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); mitk::PlaneGeometry* currentGeometry2D = dynamic_cast( const_cast(GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderer()->GetCurrentWorldGeometry2D())); PFCNot->SetGeometry2D(currentGeometry2D); PFCNot->setOperationType(mitk::PFCOMPOSITION_NOT_OPERATION); for( std::vector::iterator it = m_SelectedPF.begin(); it != m_SelectedPF.end(); ++it ) { mitk::DataNode::Pointer nodePF = *it; mitk::PlanarFigure::Pointer tmpPF = dynamic_cast( nodePF->GetData() ); PFCNot->addPlanarFigure( tmpPF ); PFCNot->addDataNode( nodePF ); PFCNot->setDisplayName("NOT_COMPO"); } AddCompositeToDatastorage(PFCNot, NULL); } /* CLEANUP NEEDED */ void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::PlanarFigureComposite::Pointer pfcomp, mitk::DataNode::Pointer parentDataNode ) { mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName( pfcomp->getDisplayName() ); newPFCNode->SetData(pfcomp); newPFCNode->SetVisibility(true); switch (pfcomp->getOperationType()) { case 0: { if (!parentDataNode.IsNull()) { GetDataStorage()->Add(newPFCNode, parentDataNode); } else { GetDataStorage()->Add(newPFCNode); } //iterate through its childs for(int i=0; igetNumberOfChildren(); ++i) { mitk::PlanarFigure::Pointer tmpPFchild = pfcomp->getChildAt(i); mitk::DataNode::Pointer savedPFchildNode = pfcomp->getDataNodeAt(i); mitk::PlanarFigureComposite::Pointer pfcompcast= dynamic_cast(tmpPFchild.GetPointer()); if ( !pfcompcast.IsNull() ) { // child is of type planar Figure composite // make new node of the child, cuz later the child has to be removed of its old position in datamanager // feed new dataNode with information of the savedDataNode, which is gonna be removed soon mitk::DataNode::Pointer newChildPFCNode; newChildPFCNode = mitk::DataNode::New(); newChildPFCNode->SetData(tmpPFchild); newChildPFCNode->SetName( savedPFchildNode->GetName() ); pfcompcast->setDisplayName( savedPFchildNode->GetName() ); //name might be changed in DataManager by user //update inside vector the dataNodePointer pfcomp->replaceDataNodeAt(i, newChildPFCNode); AddCompositeToDatastorage(pfcompcast, newPFCNode); //the current PFCNode becomes the childs parent // remove savedNode here, cuz otherwise its children will change their position in the dataNodeManager // without having its parent anymore //GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); } // remove old child position in dataStorage GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " still exists"; } } else { // child is not of type PlanarFigureComposite, so its one of the planarFigures // create new dataNode containing the data of the old dataNode, but position in dataManager will be // modified cuz we re setting a (new) parent. mitk::DataNode::Pointer newPFchildNode = mitk::DataNode::New(); newPFchildNode->SetName(savedPFchildNode->GetName() ); newPFchildNode->SetData(tmpPFchild); newPFchildNode->SetVisibility(true); // replace the dataNode in PFComp DataNodeVector pfcomp->replaceDataNodeAt(i, newPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; } else { MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); } // remove old child position in dataStorage GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " still exists"; } MITK_INFO << "adding " << newPFchildNode->GetName() << " to " << newPFCNode->GetName(); //add new child to datamanager with its new position as child of newPFCNode parent GetDataStorage()->Add(newPFchildNode, newPFCNode); } } GetDataStorage()->Modified(); break; } case 1: { if (!parentDataNode.IsNull()) { MITK_INFO << "adding " << newPFCNode->GetName() << " to " << parentDataNode->GetName() ; GetDataStorage()->Add(newPFCNode, parentDataNode); } else { MITK_INFO << "adding " << newPFCNode->GetName(); GetDataStorage()->Add(newPFCNode); } for(int i=0; igetNumberOfChildren(); ++i) { mitk::PlanarFigure::Pointer tmpPFchild = pfcomp->getChildAt(i); mitk::DataNode::Pointer savedPFchildNode = pfcomp->getDataNodeAt(i); mitk::PlanarFigureComposite::Pointer pfcompcast= dynamic_cast(tmpPFchild.GetPointer()); if ( !pfcompcast.IsNull() ) { // child is of type planar Figure composite // make new node of the child, cuz later the child has to be removed of its old position in datamanager // feed new dataNode with information of the savedDataNode, which is gonna be removed soon mitk::DataNode::Pointer newChildPFCNode; newChildPFCNode = mitk::DataNode::New(); newChildPFCNode->SetData(tmpPFchild); newChildPFCNode->SetName( savedPFchildNode->GetName() ); pfcompcast->setDisplayName( savedPFchildNode->GetName() ); //name might be changed in DataManager by user //update inside vector the dataNodePointer pfcomp->replaceDataNodeAt(i, newChildPFCNode); AddCompositeToDatastorage(pfcompcast, newPFCNode); //the current PFCNode becomes the childs parent // remove savedNode here, cuz otherwise its children will change their position in the dataNodeManager // without having its parent anymore //GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); } // remove old child position in dataStorage GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " still exists"; } } else { // child is not of type PlanarFigureComposite, so its one of the planarFigures // create new dataNode containing the data of the old dataNode, but position in dataManager will be // modified cuz we re setting a (new) parent. mitk::DataNode::Pointer newPFchildNode = mitk::DataNode::New(); newPFchildNode->SetName(savedPFchildNode->GetName() ); newPFchildNode->SetData(tmpPFchild); newPFchildNode->SetVisibility(true); // replace the dataNode in PFComp DataNodeVector pfcomp->replaceDataNodeAt(i, newPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); } // remove old child position in dataStorage GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " still exists"; } MITK_INFO << "adding " << newPFchildNode->GetName() << " to " << newPFCNode->GetName(); //add new child to datamanager with its new position as child of newPFCNode parent GetDataStorage()->Add(newPFchildNode, newPFCNode); } } GetDataStorage()->Modified(); break; } case 2: { if (!parentDataNode.IsNull()) { MITK_INFO << "adding " << newPFCNode->GetName() << " to " << parentDataNode->GetName() ; GetDataStorage()->Add(newPFCNode, parentDataNode); } else { MITK_INFO << "adding " << newPFCNode->GetName(); GetDataStorage()->Add(newPFCNode); } //iterate through its childs for(int i=0; igetNumberOfChildren(); ++i) { mitk::PlanarFigure::Pointer tmpPFchild = pfcomp->getChildAt(i); mitk::DataNode::Pointer savedPFchildNode = pfcomp->getDataNodeAt(i); mitk::PlanarFigureComposite::Pointer pfcompcast= dynamic_cast(tmpPFchild.GetPointer()); if ( !pfcompcast.IsNull() ) { // child is of type planar Figure composite // makeRemoveBundle new node of the child, cuz later the child has to be removed of its old position in datamanager // feed new dataNode with information of the savedDataNode, which is gonna be removed soon mitk::DataNode::Pointer newChildPFCNode; newChildPFCNode = mitk::DataNode::New(); newChildPFCNode->SetData(tmpPFchild); newChildPFCNode->SetName( savedPFchildNode->GetName() ); pfcompcast->setDisplayName( savedPFchildNode->GetName() ); //name might be changed in DataManager by user //update inside vector the dataNodePointer pfcomp->replaceDataNodeAt(i, newChildPFCNode); AddCompositeToDatastorage(pfcompcast, newPFCNode); //the current PFCNode becomes the childs parent // remove savedNode here, cuz otherwise its children will change their position in the dataNodeManager // without having its parent anymore //GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); } // remove old child position in dataStorage GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " still exists"; } } else { // child is not of type PlanarFigureComposite, so its one of the planarFigures // create new dataNode containing the data of the old dataNode, but position in dataManager will be // modified cuz we re setting a (new) parent. mitk::DataNode::Pointer newPFchildNode = mitk::DataNode::New(); newPFchildNode->SetName(savedPFchildNode->GetName() ); newPFchildNode->SetData(tmpPFchild); newPFchildNode->SetVisibility(true); // replace the dataNode in PFComp DataNodeVector pfcomp->replaceDataNodeAt(i, newPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_INFO << "[ERROR] does NOT exist, but can I read its Name? " << savedPFchildNode->GetName(); } // remove old child position in dataStorage GetDataStorage()->Remove(savedPFchildNode); if ( GetDataStorage()->Exists(savedPFchildNode)) { MITK_INFO << savedPFchildNode->GetName() << " still exists"; } MITK_INFO << "adding " << newPFchildNode->GetName() << " to " << newPFCNode->GetName(); //add new child to datamanager with its new position as child of newPFCNode parent GetDataStorage()->Add(newPFchildNode, newPFCNode); } } GetDataStorage()->Modified(); break; } default: MITK_INFO << "we have an UNDEFINED composition... ERROR" ; break; } } void QmitkFiberProcessingView::JoinBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( NULL, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } std::vector::const_iterator it = m_SelectedFB.begin(); mitk::FiberBundleX::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (it; it!=m_SelectedFB.end(); ++it) { newBundle = newBundle->AddBundle(dynamic_cast((*it)->GetData())); name += "+"+QString((*it)->GetName().c_str()); } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } void QmitkFiberProcessingView::SubstractBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( NULL, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } std::vector::const_iterator it = m_SelectedFB.begin(); mitk::FiberBundleX::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (it; it!=m_SelectedFB.end(); ++it) { newBundle = newBundle->SubtractBundle(dynamic_cast((*it)->GetData())); if (newBundle.IsNull()) break; name += "-"+QString((*it)->GetName().c_str()); } if (newBundle.IsNull()) { QMessageBox::information(NULL, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!"); return; } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } +void QmitkFiberProcessingView::PruneBundle() +{ + int minLength = this->m_Controls->m_PruneFibersSpinBox->value(); + bool doneSomething = false; + for (int i=0; i(m_SelectedFB.at(i)->GetData()); + if (!fib->RemoveShortFibers(minLength)) + QMessageBox::information(NULL, "No output generated:", "The resulting fiber bundle contains no fibers."); + else + doneSomething = true; + } + + if (doneSomething) + { + GenerateStats(); + RenderingManager::GetInstance()->RequestUpdateAll(); + } +} + void QmitkFiberProcessingView::GenerateStats() { if ( m_SelectedFB.empty() ) return; QString stats(""); for( int i=0; i(node->GetData())) { if (i>0) stats += "\n-----------------------------\n"; stats += QString(node->GetName().c_str()) + "\n"; mitk::FiberBundleX::Pointer fib = dynamic_cast(node->GetData()); vtkSmartPointer fiberPolyData = fib->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int numberOfLines = vLines->GetNumberOfCells(); stats += "Number of fibers: "+ QString::number(numberOfLines) + "\n"; float length = 0; std::vector lengths; for (int i=0; iGetNextCell ( numPoints, points ); float l=0; for (unsigned int j=0; j p1; itk::Point p2; fiberPolyData->GetPoint(points[j], p1.GetDataPointer()); fiberPolyData->GetPoint(points[j+1], p2.GetDataPointer()); float dist = p1.EuclideanDistanceTo(p2); length += dist; l += dist; } itk::Point p2; fiberPolyData->GetPoint(points[numPoints-1], p2.GetDataPointer()); lengths.push_back(l); } std::sort(lengths.begin(), lengths.end()); if (numberOfLines>0) length /= numberOfLines; float dev=0; int count = 0; vLines->InitTraversal(); for (int i=0; iGetNextCell ( numPoints, points ); float l=0; for (unsigned int j=0; j p1; itk::Point p2; fiberPolyData->GetPoint(points[j], p1.GetDataPointer()); fiberPolyData->GetPoint(points[j+1], p2.GetDataPointer()); float dist = p1.EuclideanDistanceTo(p2); l += dist; } dev += (length-l)*(length-l); count++; } if (numberOfLines>1) dev /= (numberOfLines-1); else dev = 0; stats += "Min. length: "+ QString::number(lengths.front(),'f',1) + " mm\n"; stats += "Max. length: "+ QString::number(lengths.back(),'f',1) + " mm\n"; stats += "Mean length: "+ QString::number(length,'f',1) + " mm\n"; stats += "Median length: "+ QString::number(lengths.at(lengths.size()/2),'f',1) + " mm\n"; stats += "Standard deviation: "+ QString::number(sqrt(dev),'f',1) + " mm\n"; } } this->m_Controls->m_StatsTextEdit->setText(stats); } void QmitkFiberProcessingView::ProcessSelectedBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( NULL, "Warning", "No fibe bundle selected!"); MITK_WARN("QmitkFiberProcessingView") << "no fibe bundle selected"; return; } int generationMethod = m_Controls->m_GenerationBox->currentIndex(); for( int i=0; i(node->GetData())) { mitk::FiberBundleX::Pointer fib = dynamic_cast(node->GetData()); QString name(node->GetName().c_str()); DataNode::Pointer newNode = NULL; switch(generationMethod){ case 0: newNode = GenerateTractDensityImage(fib, false); name += "_TDI"; break; case 1: newNode = GenerateTractDensityImage(fib, true); name += "_envelope"; break; case 2: newNode = GenerateColorHeatmap(fib); break; case 3: newNode = GenerateFiberEndingsImage(fib); name += "_fiber_endings"; break; case 4: newNode = GenerateFiberEndingsPointSet(fib); name += "_fiber_endings"; break; } if (newNode.IsNotNull()) { newNode->SetName(name.toStdString()); GetDataStorage()->Add(newNode); } } } } // generate pointset displaying the fiber endings mitk::DataNode::Pointer QmitkFiberProcessingView::GenerateFiberEndingsPointSet(mitk::FiberBundleX::Pointer fib) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); vtkSmartPointer fiberPolyData = fib->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int count = 0; int numFibers = fib->GetNumFibers(); for( int i=0; iGetNextCell ( numPoints, points ); if (numPoints>0) { double* point = fiberPolyData->GetPoint(points[0]); itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; pointSet->InsertPoint(count, itkPoint); count++; } if (numPoints>2) { double* point = fiberPolyData->GetPoint(points[numPoints-1]); itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; pointSet->InsertPoint(count, itkPoint); count++; } } mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( pointSet ); return node; } // generate image displaying the fiber endings mitk::DataNode::Pointer QmitkFiberProcessingView::GenerateFiberEndingsImage(mitk::FiberBundleX::Pointer fib) { typedef unsigned char OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetInvertImage(m_Controls->m_InvertCheckbox->isChecked()); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image OutImageType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate rgba heatmap from fiber bundle mitk::DataNode::Pointer QmitkFiberProcessingView::GenerateColorHeatmap(mitk::FiberBundleX::Pointer fib) { typedef itk::RGBAPixel OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { itk::Image::Pointer itkImage = itk::Image::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate tract density image from fiber bundle mitk::DataNode::Pointer QmitkFiberProcessingView::GenerateTractDensityImage(mitk::FiberBundleX::Pointer fib, bool binary) { typedef float OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetInvertImage(m_Controls->m_InvertCheckbox->isChecked()); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } void QmitkFiberProcessingView::ResampleSelectedBundles() { int factor = this->m_Controls->m_ResampleFibersSpinBox->value(); for (int i=0; i(m_SelectedFB.at(i)->GetData()); fib->DoFiberSmoothing(factor); } + GenerateStats(); + RenderingManager::GetInstance()->RequestUpdateAll(); +} + +void QmitkFiberProcessingView::MirrorFibers() +{ + unsigned int axis = this->m_Controls->m_AxisSelectionBox->currentIndex(); + for (int i=0; i(m_SelectedFB.at(i)->GetData()); + fib->MirrorFibers(axis); + } + GenerateStats(); + RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoFaColorCoding() { if (m_SelectedImage.IsNull()) return; -// mitk::PixelType pType = mitk::MakeScalarPixelType(); -// if (m_SelectedImage->GetPixelType()!=pType) -// { -// //mitk::Image bla; bla.GetPixelType().GetNameOfClass() -// MITK_INFO << m_SelectedImage->GetPixelType().GetNameOfClass(); -// QMessageBox::warning(NULL, "Wrong Image Type", "FA/GFA image should be of type float"); -//// return; -// } for( int i=0; i(m_SelectedFB.at(i)->GetData()); fib->SetFAMap(m_SelectedImage); fib->SetColorCoding(mitk::FiberBundleX::COLORCODING_FA_BASED); fib->DoColorCodingFaBased(); } if(m_MultiWidget) m_MultiWidget->RequestUpdate(); } - diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.h index c1f606e896..db1daf65c2 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.h @@ -1,197 +1,199 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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 QmitkFiberProcessingView_h #define QmitkFiberProcessingView_h #include #include "ui_QmitkFiberProcessingViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include /*! \brief QmitkFiberProcessingView \warning View to process fiber bundles. Supplies methods to extract fibers from the bundle, join and subtract bundles, generate images from the selected bundle and much more. \sa QmitkFunctionality \ingroup Functionalities */ class QmitkFiberProcessingView : 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: typedef itk::Image< unsigned char, 3 > itkUCharImageType; typedef itk::Image< float, 3 > itkFloatImageType; static const std::string VIEW_ID; QmitkFiberProcessingView(); virtual ~QmitkFiberProcessingView(); virtual void CreateQtPartControl(QWidget *parent); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); virtual void Activated(); protected slots: void OnDrawCircle(); void OnDrawPolygon(); void DoFiberExtraction(); void GenerateAndComposite(); void GenerateOrComposite(); void GenerateNotComposite(); - + void PruneBundle(); + void MirrorFibers(); void JoinBundles(); void SubstractBundles(); void GenerateRoiImage(); void ProcessSelectedBundles(); void ResampleSelectedBundles(); void DoFaColorCoding(); void Extract3d(); virtual void AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey = NULL, mitk::BaseProperty *property = NULL ); protected: /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); Ui::QmitkFiberProcessingViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; /** Connection from VTK to ITK */ template void ConnectPipelines(VTK_Exporter* exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } template void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string nodeName ); template < typename TPixel, unsigned int VImageDimension > void InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::Geometry3D* planegeo3D, int additionalIndex ); void GenerateStats(); void UpdateGui(); berry::ISelectionListener::Pointer m_SelListener; berry::IStructuredSelection::ConstPointer m_CurrentSelection; private: int m_EllipseCounter; int m_PolygonCounter; //contains the selected FiberBundles std::vector m_SelectedFB; //contains the selected PlanarFigures std::vector m_SelectedPF; mitk::Image::Pointer m_SelectedImage; mitk::Image::Pointer m_InternalImage; mitk::PlanarFigure::Pointer m_PlanarFigure; float m_UpsamplingFactor; itkUCharImageType::Pointer m_InternalImageMask3D; itkUCharImageType::Pointer m_PlanarFigureImage; std::vector m_Surfaces; void AddCompositeToDatastorage(mitk::PlanarFigureComposite::Pointer, mitk::DataNode::Pointer); void debugPFComposition(mitk::PlanarFigureComposite::Pointer , int ); void CompositeExtraction(mitk::DataNode::Pointer node, mitk::Image* image); mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundleX::Pointer fib, bool binary); mitk::DataNode::Pointer GenerateColorHeatmap(mitk::FiberBundleX::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsImage(mitk::FiberBundleX::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsPointSet(mitk::FiberBundleX::Pointer fib); }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui index e907ff8063..34374ce3c9 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui @@ -1,683 +1,790 @@ QmitkFiberProcessingViewControls 0 0 665 - 587 + 671 Form 0 9 3 9 3 Fiber Bundle Modification 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 30 30 Draw circular ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/circle.png:/QmitkDiffusionImaging/circle.png 32 32 false true 30 30 Draw rectangular ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/rectangle.png:/QmitkDiffusionImaging/rectangle.png 32 32 true true 30 30 Draw polygonal ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/polygon.png:/QmitkDiffusionImaging/polygon.png 32 32 true true Qt::Horizontal 40 20 QFrame::NoFrame QFrame::Raised 0 false 0 0 200 16777215 11 Extract fibers passing through selected ROI or composite ROI. Select ROI and fiber bundle to execute. Extract false 0 0 200 16777215 11 Returns all fibers contained in bundle X that are not contained in bundle Y (not commutative!). Select at least two fiber bundles to execute. Substract false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Join Qt::Horizontal 40 20 false 0 0 200 16777215 11 Extract fibers passing through selected surface mesh. Select surface mesh and fiber bundle to execute. Extract 3D false 0 0 16777215 16777215 11 Generate a binary image containing all selected ROIs. Select at least one ROI (planar figure) and a reference fiber bundle or image. ROI Image 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 Qt::Horizontal 40 20 false 60 16777215 Create AND composition with selected ROIs. AND false 60 16777215 Create OR composition with selected ROIs. OR false 60 16777215 Create NOT composition from selected ROI. NOT Fiber Bundle Processing + + QFormLayout::AllNonFixedFieldsGrow + 0 0 Tract Density Image (TDI) Binary Envelope Fiber Bundle Image Fiber Endings Image Fiber Endings Pointset Upsampling factor 1 10 2 false 0 0 200 16777215 11 Perform selected operation on all selected fiber bundles. Generate If selected operation generates an image, the inverse image is returned. Invert false 0 0 200 16777215 11 Resample fibers using a Kochanek spline interpolation. Smooth Fibers Points per cm 1 50 10 + + + false + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + Remove fibers shorten than the specified length (in mm). + + + Prune Bundle + + + + + + + Minimum fiber length in mm + + + 0 + + + 1000 + + + 20 + + + + + + + false + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + Mirror fibers around specified axis. + + + Mirror Fibers + + + + false 0 0 200 16777215 11 Apply float image values (0-1) as color coding to the selected fiber bundle. Color By Scalar Map + + + + 0 + + + 3 + + + 3 + + + + x - Axis + + + + + y - Axis + + + + + z - Axis + + + + Fiber Bundle Statistics Courier 10 Pitch false true Qt::Vertical 20 40 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui index b99d3213ef..de4d981ce9 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui @@ -1,1045 +1,1045 @@ QmitkGibbsTrackingViewControls 0 0 463 1011 0 0 0 0 QmitkTemplate 0 9 3 9 3 Data Q-Ball Image: Mandatory input - Mask Image: Optional input to limit the algorithms search space. - Parameters 0 Iterations: 10^7 Specify number of iterations for the tracking algorithm. 10 6 Qt::Horizontal QSlider::TicksBelow true Activate continuous visualization of intermediate results. Visualize Tractography true Visualize intermediate result. :/QmitkDiffusionImaging/Refresh_48.png:/QmitkDiffusionImaging/Refresh_48.png true Advanced Settings Output File: QFrame::NoFrame QFrame::Plain 0 0 0 Select output file name and folder. ... N/A true true QFrame::StyledPanel QFrame::Raised 9 0 9 0 4 auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter auto = 0.5 * min. spacing; sigma 100 1 Qt::Horizontal QSlider::NoTicks auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter automatic estimation from gfa map and q-ball data. 0 1000 1 0 Qt::Horizontal true QSlider::NoTicks 0.1 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 1 100 1 10 Qt::Horizontal false false QSlider::NoTicks 0.001 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 1 99 1 10 Qt::Horizontal QSlider::NoTicks 0 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter IE Bias < 0 < EE Bias -50 50 1 Qt::Horizontal QSlider::NoTicks 40mm Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter Only fibers longer than specified are accepted. - 100 + 500 1 40 Qt::Horizontal QSlider::NoTicks Particle Length: Particle Width: Particle Weight: Start Temperature: End Temperature: Balance In/Ex Energy: Min. Fiber Length: true Use mean subtracted ODFs (recommended). Subtract ODF Mean true Qt::Horizontal QSizePolicy::Fixed 60 20 Curvature Threshold: 45° Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter Allow only fiber curvature values smaller than the selected threshold. 180 1 45 Qt::Horizontal QSlider::NoTicks auto = 1.5 * min. spacing; l 100 1 Qt::Horizontal QSlider::NoTicks QFrame::NoFrame QFrame::Plain 0 0 0 true Save current parameters as xml (.gtp) Qt::LeftToRight Save Parameters :/qmitk/btnMoveDown.png:/qmitk/btnMoveDown.png true Load parameters from xml file (.gtp) Qt::LeftToRight Load Parameters :/qmitk/btnMoveUp.png:/qmitk/btnMoveUp.png false No Q-Ball image selected. Qt::LeftToRight Start Tractography :/qmitk/play.xpm:/qmitk/play.xpm false Qt::LeftToRight Stop Tractography :/qmitk/stop.xpm:/qmitk/stop.xpm Monitor Progress: - Will only be updated if tracking is visualized Will only be updated if tracking is visualized Accepted Fibers: Connections: Particles: Proposal Acceptance Rate: Tracking Time: Will only be updated if tracking is visualized - - - - - Qt::Vertical QSizePolicy::Expanding 0 0 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.cpp new file mode 100644 index 0000000000..fe502112c6 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.cpp @@ -0,0 +1,197 @@ +/*=================================================================== + +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 +#include + +// Qmitk +#include "QmitkStreamlineTrackingView.h" +#include "QmitkStdMultiWidget.h" + +// Qt +#include + +// MITK +#include +#include + +// VTK +#include +#include +#include +#include +#include +#include + + +const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; +const std::string id_DataManager = "org.mitk.views.datamanager"; +using namespace berry; + +QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() + : QmitkFunctionality() + , m_Controls( 0 ) + , m_MultiWidget( NULL ) + , m_TensorImage( NULL ) + , m_SeedRoi( NULL ) +{ +} + +// Destructor +QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() +{ + +} + +void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) +{ + if ( !m_Controls ) + { + // create GUI widgets from the Qt Designer's .ui file + m_Controls = new Ui::QmitkStreamlineTrackingViewControls; + m_Controls->setupUi( parent ); + + connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); + connect( m_Controls->m_SeedsPerVoxelSlider, SIGNAL(valueChanged(int)), this, SLOT(OnSeedsPerVoxelChanged(int)) ); + connect( m_Controls->m_MinTractLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(OnMinTractLengthChanged(int)) ); + connect( m_Controls->m_FaThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(OnFaThresholdChanged(int)) ); + } +} + +void QmitkStreamlineTrackingView::OnSeedsPerVoxelChanged(int value) +{ + m_Controls->m_SeedsPerVoxelLabel->setText(QString("Seeds per Voxel: ")+QString::number(value)); +} + +void QmitkStreamlineTrackingView::OnMinTractLengthChanged(int value) +{ + m_Controls->m_MinTractLengthLabel->setText(QString("Min. Tract Length: ")+QString::number(value)+QString("mm")); +} + +void QmitkStreamlineTrackingView::OnFaThresholdChanged(int value) +{ + m_Controls->m_FaThresholdLabel->setText(QString("FA Threshold: ")+QString::number((float)value/100)); +} + +void QmitkStreamlineTrackingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) +{ + m_MultiWidget = &stdMultiWidget; +} + + +void QmitkStreamlineTrackingView::StdMultiWidgetNotAvailable() +{ + m_MultiWidget = NULL; +} + +void QmitkStreamlineTrackingView::OnSelectionChanged( std::vector nodes ) +{ + m_TensorImageNode = NULL; + m_TensorImage = NULL; + m_SeedRoi = NULL; + m_Controls->m_TensorImageLabel->setText("-"); + m_Controls->m_RoiImageLabel->setText("-"); + + if(nodes.empty()) + return; + + for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) + { + mitk::DataNode::Pointer node = *it; + + if( node.IsNotNull() && dynamic_cast(node->GetData()) ) + { + if( dynamic_cast(node->GetData()) ) + { + m_TensorImageNode = node; + m_TensorImage = dynamic_cast(node->GetData()); + m_Controls->m_TensorImageLabel->setText(node->GetName().c_str()); + } + else + { + bool isBinary = false; + node->GetPropertyValue("binary", isBinary); + if (isBinary) + { + m_SeedRoi = dynamic_cast(node->GetData()); + m_Controls->m_RoiImageLabel->setText(node->GetName().c_str()); + } + } + } + } + + if(m_TensorImage.IsNotNull()) + m_Controls->commandLinkButton->setEnabled(true); + else + m_Controls->commandLinkButton->setEnabled(false); +} + + + +void QmitkStreamlineTrackingView::DoFiberTracking() +{ + if (m_TensorImage.IsNull()) + return; + + typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; + typedef mitk::ImageToItk CastType; + typedef mitk::ImageToItk CastType2; + + CastType::Pointer caster = CastType::New(); + caster->SetInput(m_TensorImage); + caster->Update(); + TensorImageType::Pointer image = caster->GetOutput(); + + typedef itk::StreamlineTrackingFilter< float > FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(image); + filter->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelSlider->value()); + filter->SetFaThreshold((float)m_Controls->m_FaThresholdSlider->value()/100); + + if (m_SeedRoi.IsNotNull()) + { + CastType2::Pointer caster2 = CastType2::New(); + caster2->SetInput(m_SeedRoi); + caster2->Update(); + ItkUCharImageType::Pointer mask = caster2->GetOutput(); + filter->SetMaskImage(mask); + } + + filter->Update(); + + vtkSmartPointer fiberBundle = filter->GetFiberPolyData(); + if ( fiberBundle->GetNumberOfLines()==0 ) + return; + mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(fiberBundle); + + fib->RemoveShortFibers(m_Controls->m_MinTractLengthSlider->value()); + + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(fib); + + QString name(m_TensorImageNode->GetName().c_str()); + name += "_FiberBundle"; + node->SetName(name.toStdString()); + + node->SetVisibility(true); + + GetDataStorage()->Add(node); +} + + diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.h new file mode 100644 index 0000000000..3c436283e3 --- /dev/null +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.h @@ -0,0 +1,87 @@ +/*=================================================================== + +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 QmitkStreamlineTrackingView_h +#define QmitkStreamlineTrackingView_h + +#include + +#include "ui_QmitkStreamlineTrackingViewControls.h" + +#include +#include +#include +#include +#include + + +/*! +\brief QmitkStreamlineTrackingView + +\warning Implements standard streamline tracking as proposed by Mori et al. 1999 "Three-Dimensional Tracking of Axonal Projections in the Brain by Magnetic Resonance Imaging" + +\sa QmitkFunctionality +\ingroup Functionalities +*/ +class QmitkStreamlineTrackingView : 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< unsigned char, 3 > ItkUCharImageType; + + QmitkStreamlineTrackingView(); + virtual ~QmitkStreamlineTrackingView(); + + virtual void CreateQtPartControl(QWidget *parent); + + virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); + virtual void StdMultiWidgetNotAvailable(); + + protected slots: + + void DoFiberTracking(); + +protected: + + /// \brief called by QmitkFunctionality when DataManager's selection has changed + virtual void OnSelectionChanged( std::vector nodes ); + Ui::QmitkStreamlineTrackingViewControls* m_Controls; + QmitkStdMultiWidget* m_MultiWidget; + +protected slots: + + void OnSeedsPerVoxelChanged(int value); + void OnMinTractLengthChanged(int value); + void OnFaThresholdChanged(int value); + +private: + + mitk::Image::Pointer m_SeedRoi; + mitk::TensorImage::Pointer m_TensorImage; + mitk::DataNode::Pointer m_TensorImageNode; + +}; + + + +#endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED + diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui new file mode 100644 index 0000000000..5a90ff232e --- /dev/null +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui @@ -0,0 +1,212 @@ + + + QmitkStreamlineTrackingViewControls + + + + 0 + 0 + 480 + 553 + + + + + 0 + 0 + + + + QmitkTemplate + + + + 3 + + + 3 + + + 0 + + + + + Data + + + + + + Tensor Image: + + + + + + + - + + + + + + + Seed ROI Image: + + + + + + + - + + + + + + + + + + Parameters + + + + + + Qt::Horizontal + + + QSizePolicy::Fixed + + + + 200 + 0 + + + + + + + + Number of tracts started in each voxel of the seed ROI. + + + Seeds per Voxel: 1 + + + + + + + Minimum tract length in mm. + + + 0 + + + 500 + + + 40 + + + Qt::Horizontal + + + + + + + Number of tracts started in each voxel of the seed ROI. + + + 1 + + + 10 + + + Qt::Horizontal + + + + + + + Minimum tract length in mm. + + + Min. Tract Length: 40mm + + + + + + + Minimum tract length in mm. + + + FA Threshold: 0.2 + + + + + + + Fractional Anisotropy Threshold + + + 0 + + + 100 + + + 20 + + + Qt::Horizontal + + + + + + + + + + false + + + Start Tracking + + + + + + + Qt::Vertical + + + QSizePolicy::Expanding + + + + 20 + 220 + + + + + + + + + commandLinkButton + + + + diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp index 50f76f44a1..de283acb87 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp @@ -1,71 +1,73 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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 "mitkPluginActivator.h" #include #include "src/internal/QmitkDiffusionImagingPublicPerspective.h" #include "src/internal/QmitkQBallReconstructionView.h" #include "src/internal/QmitkPreprocessingView.h" #include "src/internal/QmitkDiffusionDicomImportView.h" #include "src/internal/QmitkDiffusionQuantificationView.h" #include "src/internal/QmitkTensorReconstructionView.h" #include "src/internal/QmitkControlVisualizationPropertiesView.h" #include "src/internal/QmitkODFDetailsView.h" #include "src/internal/QmitkGibbsTrackingView.h" #include "src/internal/QmitkStochasticFiberTrackingView.h" #include "src/internal/QmitkFiberProcessingView.h" #include "src/internal/QmitkFiberBundleDeveloperView.h" #include "src/internal/QmitkPartialVolumeAnalysisView.h" #include "src/internal/QmitkIVIMView.h" #include "src/internal/QmitkTractbasedSpatialStatisticsView.h" #include "src/internal/QmitkTbssSkeletonizationView.h" +#include "src/internal/QmitkStreamlineTrackingView.h" #include "src/internal/Connectomics/QmitkBrainNetworkAnalysisView.h" namespace mitk { void PluginActivator::start(ctkPluginContext* context) { BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionImagingPublicPerspective, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkQBallReconstructionView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkPreprocessingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionDicomImport, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionQuantificationView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkTensorReconstructionView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkControlVisualizationPropertiesView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkODFDetailsView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkGibbsTrackingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkStochasticFiberTrackingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkFiberProcessingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkFiberBundleDeveloperView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkPartialVolumeAnalysisView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkIVIMView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkTractbasedSpatialStatisticsView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkTbssSkeletonizationView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkBrainNetworkAnalysisView, context) + BERRY_REGISTER_EXTENSION_CLASS(QmitkStreamlineTrackingView, context) } void PluginActivator::stop(ctkPluginContext* context) { Q_UNUSED(context) } } Q_EXPORT_PLUGIN2(org_mitk_gui_qt_diffusionimaging, mitk::PluginActivator) diff --git a/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppFiberTractographyPerspective.cpp b/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppFiberTractographyPerspective.cpp index e03ad86c59..c476b26af0 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppFiberTractographyPerspective.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppFiberTractographyPerspective.cpp @@ -1,56 +1,60 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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 "QmitkDIAppFiberTractographyPerspective.h" #include "berryIViewLayout.h" void QmitkDIAppFiberTractographyPerspective::CreateInitialLayout(berry::IPageLayout::Pointer layout) { ///////////////////////////////////////////////////// // all di-app perspectives should have the following: ///////////////////////////////////////////////////// std::string editorArea = layout->GetEditorArea(); layout->AddStandaloneView("org.mitk.views.datamanager", false, berry::IPageLayout::LEFT, 0.3f, editorArea); layout->AddStandaloneView("org.mitk.views.controlvisualizationpropertiesview", false, berry::IPageLayout::BOTTOM, .2f, "org.mitk.views.datamanager"); berry::IFolderLayout::Pointer left = layout->CreateFolder("org.mbi.diffusionimaginginternal.leftcontrols", berry::IPageLayout::BOTTOM, 0.15f, "org.mitk.views.controlvisualizationpropertiesview"); layout->AddStandaloneView("org.mitk.views.imagenavigator", false, berry::IPageLayout::BOTTOM, .4f, "org.mbi.diffusionimaginginternal.leftcontrols"); ///////////////////////////////////////////// // here goes the perspective specific stuff ///////////////////////////////////////////// left->AddView("org.mitk.views.gibbstracking"); berry::IViewLayout::Pointer lo = layout->GetViewLayout("org.mitk.views.gibbstracking"); lo->SetCloseable(false); left->AddView("org.mitk.views.stochasticfibertracking"); lo = layout->GetViewLayout("org.mitk.views.stochasticfibertracking"); lo->SetCloseable(false); + left->AddView("org.mitk.views.streamlinetracking"); + lo = layout->GetViewLayout("org.mitk.views.streamlinetracking"); + lo->SetCloseable(false); + left->AddView("org.mitk.views.fiberprocessing"); lo = layout->GetViewLayout("org.mitk.views.fiberprocessing"); lo->SetCloseable(false); }