diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp index ded085ffcc..79287c83ae 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,1643 +1,1722 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #define _USE_MATH_DEFINES #include "mitkFiberBundleX.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundleX::COLORCODING_ORIENTATION_BASED = "Color_Orient"; //const char* mitk::FiberBundleX::COLORCODING_FA_AS_OPACITY = "Color_Orient_FA_Opacity"; const char* mitk::FiberBundleX::COLORCODING_FA_BASED = "FA_Values"; const char* mitk::FiberBundleX::COLORCODING_CUSTOM = "custom"; const char* mitk::FiberBundleX::FIBER_ID_ARRAY = "Fiber_IDs"; using namespace std; mitk::FiberBundleX::FiberBundleX( vtkPolyData* fiberPolyData ) : m_CurrentColorCoding(NULL) , m_NumFibers(0) , m_FiberSampling(0) { m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != NULL) { m_FiberPolyData = fiberPolyData; //m_FiberPolyData->DeepCopy(fiberPolyData); this->DoColorCodingOrientationBased(); } m_NumFibers = m_FiberPolyData->GetNumberOfLines(); this->UpdateFiberGeometry(); this->SetColorCoding(COLORCODING_ORIENTATION_BASED); this->GenerateFiberIds(); } mitk::FiberBundleX::~FiberBundleX() { } mitk::FiberBundleX::Pointer mitk::FiberBundleX::GetDeepCopy() { mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(m_FiberPolyData); newFib->SetColorCoding(m_CurrentColorCoding); return newFib; } vtkSmartPointer mitk::FiberBundleX::GeneratePolyDataByIds(std::vector fiberIds) { MITK_DEBUG << "\n=====FINAL RESULT: fib_id ======\n"; MITK_DEBUG << "Number of new Fibers: " << fiberIds.size(); // iterate through the vectorcontainer hosting all desired fiber Ids vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); vtkSmartPointer newLineSet = vtkSmartPointer::New(); vtkSmartPointer newPointSet = vtkSmartPointer::New(); // if FA array available, initialize fa double array // if color orient array is available init color array vtkSmartPointer faValueArray; vtkSmartPointer colorsT; //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = sizeof(rgba); if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_FA_BASED)){ MITK_DEBUG << "FA VALUES AVAILABLE, init array for new fiberbundle"; faValueArray = vtkSmartPointer::New(); } if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)){ MITK_DEBUG << "colorValues available, init array for new fiberbundle"; colorsT = vtkUnsignedCharArray::New(); colorsT->SetNumberOfComponents(componentSize); colorsT->SetName(COLORCODING_ORIENTATION_BASED); } std::vector::iterator finIt = fiberIds.begin(); while ( finIt != fiberIds.end() ) { if (*finIt < 0 || *finIt>GetNumFibers()){ MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; break; } vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); vtkSmartPointer fibPoints = fiber->GetPoints(); vtkSmartPointer newFiber = vtkSmartPointer::New(); newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); for(int i=0; iGetNumberOfPoints(); i++) { // MITK_DEBUG << "id: " << fiber->GetPointId(i); // MITK_DEBUG << fibPoints->GetPoint(i)[0] << " | " << fibPoints->GetPoint(i)[1] << " | " << fibPoints->GetPoint(i)[2]; newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_FA_BASED)){ // MITK_DEBUG << m_FiberIdDataSet->GetPointData()->GetArray(FA_VALUE_ARRAY)->GetTuple(fiber->GetPointId(i)); } if (m_FiberIdDataSet->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED)){ // MITK_DEBUG << "ColorValue: " << m_FiberIdDataSet->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetTuple(fiber->GetPointId(i))[0]; } } newLineSet->InsertNextCell(newFiber); ++finIt; } newFiberPolyData->SetPoints(newPointSet); newFiberPolyData->SetLines(newLineSet); MITK_DEBUG << "new fiberbundle polydata points: " << newFiberPolyData->GetNumberOfPoints(); MITK_DEBUG << "new fiberbundle polydata lines: " << newFiberPolyData->GetNumberOfLines(); MITK_DEBUG << "=====================\n"; // mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(newFiberPolyData); return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundleX::Pointer mitk::FiberBundleX::AddBundle(mitk::FiberBundleX* fib) { if (fib==NULL) { MITK_WARN << "trying to call AddBundle with NULL argument"; return NULL; } vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); // add current fiber bundle int numFibers = GetNumFibers(); for( int i=0; iGetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } vLines = fib->m_FiberPolyData->GetLines(); vLines->InitTraversal(); // add new fiber bundle numFibers = fib->GetNumFibers(); for( int i=0; iGetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(fib->m_FiberPolyData->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); return newFib; } // subtract two fiber bundles mitk::FiberBundleX::Pointer mitk::FiberBundleX::SubtractBundle(mitk::FiberBundleX* fib) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); // iterate over current fibers int numFibers = GetNumFibers(); for( int i=0; iGetNextCell ( numPoints, points ); if (points==NULL) continue; vtkSmartPointer vLines2 = fib->m_FiberPolyData->GetLines(); vLines2->InitTraversal(); int numFibers2 = fib->GetNumFibers(); bool contained = false; for( int i2=0; i2GetNextCell ( numPoints2, points2 ); if (points2==NULL) continue; // check endpoints itk::Point point_start = GetItkPoint(m_FiberPolyData->GetPoint(points[0])); itk::Point point_end = GetItkPoint(m_FiberPolyData->GetPoint(points[numPoints-1])); itk::Point point2_start = GetItkPoint(fib->m_FiberPolyData->GetPoint(points2[0])); itk::Point point2_end = GetItkPoint(fib->m_FiberPolyData->GetPoint(points2[numPoints2-1])); if (point_start.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps || point_start.SquaredEuclideanDistanceTo(point2_end)<=mitk::eps && point_end.SquaredEuclideanDistanceTo(point2_start)<=mitk::eps) { // further checking ??? if (numPoints2==numPoints) contained = true; } } // add to result because fiber is not subtracted if (!contained) { vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } } if(vNewLines->GetNumberOfCells()==0) return NULL; // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundleX::Pointer newFib = mitk::FiberBundleX::New(vNewPolyData); return newFib; } itk::Point mitk::FiberBundleX::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } /* * set polydata (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundleX::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { if (fiberPD == NULL) this->m_FiberPolyData = vtkSmartPointer::New(); else { m_FiberPolyData->DeepCopy(fiberPD); DoColorCodingOrientationBased(); } m_NumFibers = m_FiberPolyData->GetNumberOfLines(); if (updateGeometry) UpdateFiberGeometry(); SetColorCoding(COLORCODING_ORIENTATION_BASED); GenerateFiberIds(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundleX::GetFiberPolyData() { return m_FiberPolyData; } void mitk::FiberBundleX::DoColorCodingOrientationBased() { //===== FOR WRITING A TEST ======================== // colorT size == tupelComponents * tupelElements // compare color results // to cover this code 100% also polydata needed, where colorarray already exists // + one fiber with exactly 1 point // + one fiber with 0 points //================================================= /* make sure that processing colorcoding is only called when necessary */ if ( m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) && m_FiberPolyData->GetNumberOfPoints() == m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetNumberOfTuples() ) { // fiberstructure is already colorcoded MITK_DEBUG << " NO NEED TO REGENERATE COLORCODING! " ; this->ResetFiberOpacity(); this->SetColorCoding(COLORCODING_ORIENTATION_BASED); return; } /* Finally, execute color calculation */ vtkPoints* extrPoints = NULL; extrPoints = m_FiberPolyData->GetPoints(); int numOfPoints = 0; if (extrPoints!=NULL) numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; // int componentSize = sizeof(rgba); int componentSize = 4; vtkSmartPointer colorsT = vtkSmartPointer::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 = vtkSmartPointer::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(ItkUcharImgType* mask, bool anyPoint) +{ + vtkSmartPointer polyData = this->GetFiberPolyData(); + if (anyPoint) + { + float minSpacing = 1; + if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) + minSpacing = mask->GetSpacing()[0]; + else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) + minSpacing = mask->GetSpacing()[1]; + else + minSpacing = mask->GetSpacing()[2]; + + mitk::FiberBundleX::Pointer fibCopy = this->GetDeepCopy(); + fibCopy->ResampleFibers(minSpacing/2); + polyData = fibCopy->GetFiberPolyData(); + } + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + vtkSmartPointer vLines = polyData->GetLines(); + vLines->InitTraversal(); + + MITK_INFO << "Extracting fibers"; + boost::progress_display disp(m_NumFibers); + for (int i=0; iGetNextCell ( numPoints, pointIds ); + + vtkSmartPointer container = vtkSmartPointer::New(); + + if (numPoints>1) + { + if (anyPoint) + { + for (int j=0; jGetPoint(pointIds[j]); + + itk::Point itkP; + itkP[0] = p[0]; itkP[1] = p[1]; p[2] = p[2]; + itk::Index<3> idx; + mask->TransformPhysicalPointToIndex(itkP, idx); + if ( mask->GetPixel(idx)>0 ) + { + for (int j=0; jGetPoint(pointIds[j]); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + break; + } + } + } + else + { + double* start = polyData->GetPoint(pointIds[0]); + itk::Point itkStart; + itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2]; + itk::Index<3> idxStart; + mask->TransformPhysicalPointToIndex(itkStart, idxStart); + + double* end = polyData->GetPoint(pointIds[numPoints-1]); + itk::Point itkEnd; + itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2]; + itk::Index<3> idxEnd; + mask->TransformPhysicalPointToIndex(itkEnd, idxEnd); + + if ( mask->GetPixel(idxStart)>0 && mask->GetPixel(idxEnd)>0 ) + { + for (int j=0; jGetPoint(pointIds[j]); + vtkIdType id = vtkNewPoints->InsertNextPoint(p); + container->GetPointIds()->InsertNextId(id); + } + } + } + } + + vtkNewCells->InsertNextCell(container); + } + + if (vtkNewCells->GetNumberOfCells()<=0) + return NULL; + + vtkSmartPointer newPolyData = vtkSmartPointer::New(); + newPolyData->SetPoints(vtkNewPoints); + newPolyData->SetLines(vtkNewCells); + return mitk::FiberBundleX::New(newPolyData); +} 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]); /* 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"; 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); } 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 = vtkSmartPointer::New(); //get the control points from pf and insert them to vtkPolygon unsigned int nrCtrlPnts = pf->GetNumberOfControlPoints(); for (int i=0; iGetPoints()->InsertNextPoint((double)pf->GetWorldControlPoint(i)[0], (double)pf->GetWorldControlPoint(i)[1], (double)pf->GetWorldControlPoint(i)[2] ); } //prepare everything for using pointInPolygon function double n[3]; polygonVtk->ComputeNormal(polygonVtk->GetPoints()->GetNumberOfPoints(), static_cast(polygonVtk->GetPoints()->GetData()->GetVoidPointer(0)), n); double bounds[6]; polygonVtk->GetPoints()->GetBounds(bounds); for (int i=0; iGetPoint(PointsOnPlane[i])[0], clipperout->GetPoint(PointsOnPlane[i])[1], clipperout->GetPoint(PointsOnPlane[i])[2]}; int isInPolygon = polygonVtk->PointInPolygon(checkIn, polygonVtk->GetPoints()->GetNumberOfPoints() , static_cast(polygonVtk->GetPoints()->GetData()->GetVoidPointer(0)), bounds, n); if( isInPolygon ) PointsInROI.push_back(PointsOnPlane[i]); } } MITK_DEBUG << "Step3: Identify fibers"; // we need to access the fiberId Array, so make sure that this array is available if (!clipperout->GetCellData()->HasArray(FIBER_ID_ARRAY)) { MITK_DEBUG << "ERROR: FiberID array does not exist, no correlation between points and fiberIds possible! Make sure calling GenerateFiberIds()"; return FibersInROI; // FibersInRoi is empty then } if (PointsInROI.size()<=0) return FibersInROI; // prepare a structure where each point id is represented as an indexId. // vector looks like: | pntId | fiberIdx | std::vector< long > pointindexFiberMap; // walk through the whole subline section and create an vector sorted by point index vtkCellArray *clipperlines = clipperout->GetLines(); clipperlines->InitTraversal(); long numOfLineCells = clipperlines->GetNumberOfCells(); long numofClippedPoints = clipperout->GetNumberOfPoints(); pointindexFiberMap.resize(numofClippedPoints); //prepare resulting vector FibersInROI.reserve(PointsInROI.size()); MITK_DEBUG << "\n===== Pointindex based structure initialized ======\n"; // go through resulting "sub"lines which are stored as cells, "i" corresponds to current line id. for (int i=0, ic=0 ; iGetCell(ic, npts, pts); // go through point ids in hosting subline, "j" corresponds to current pointindex in current line i. eg. idx[0]=45; idx[1]=46 for (long j=0; jGetCellData()->GetArray(FIBER_ID_ARRAY)->GetTuple(i)[0] << " to pointId: " << pts[j]; pointindexFiberMap[ pts[j] ] = clipperout->GetCellData()->GetArray(FIBER_ID_ARRAY)->GetTuple(i)[0]; // MITK_DEBUG << "in array: " << pointindexFiberMap[ pts[j] ]; } } MITK_DEBUG << "\n===== Pointindex based structure finalized ======\n"; // get all Points in ROI with according fiberID for (long k = 0; k < PointsInROI.size(); k++) { //MITK_DEBUG << "point " << PointsInROI[k] << " belongs to fiber " << pointindexFiberMap[ PointsInROI[k] ]; if (pointindexFiberMap[ PointsInROI[k] ]<=GetNumFibers() && pointindexFiberMap[ PointsInROI[k] ]>=0) FibersInROI.push_back(pointindexFiberMap[ PointsInROI[k] ]); else MITK_INFO << "ERROR in ExtractFiberIdSubset; impossible fiber id detected"; } m_PointsRoi = PointsInROI; } // detecting fiberId duplicates MITK_DEBUG << "check for duplicates"; sort(FibersInROI.begin(), FibersInROI.end()); bool hasDuplicats = false; for(long i=0; i::iterator it; it = unique (FibersInROI.begin(), FibersInROI.end()); FibersInROI.resize( it - FibersInROI.begin() ); } return FibersInROI; } void mitk::FiberBundleX::UpdateFiberGeometry() { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInput(m_FiberPolyData); cleaner->PointMergingOff(); cleaner->Update(); m_FiberPolyData = cleaner->GetOutput(); m_FiberLengths.clear(); m_MeanFiberLength = 0; m_MedianFiberLength = 0; m_LengthStDev = 0; m_NumFibers = m_FiberPolyData->GetNumberOfLines(); if (m_NumFibers<=0) // no fibers present; apply default geometry { m_MinFiberLength = 0; m_MaxFiberLength = 0; mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetImageGeometry(true); float b[] = {0, 1, 0, 1, 0, 1}; geometry->SetFloatBounds(b); SetGeometry(geometry); return; } float min = itk::NumericTraits::NonpositiveMin(); float max = itk::NumericTraits::max(); float b[] = {max, min, max, min, max, min}; vtkCellArray* cells = m_FiberPolyData->GetLines(); cells->InitTraversal(); for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float length = 0; for (int j=0; jGetPoint(j, p1); if (p1[0]b[1]) b[1]=p1[0]; if (p1[1]b[3]) b[3]=p1[1]; if (p1[2]b[5]) b[5]=p1[2]; // calculate statistics if (jGetPoint(j+1, p2); float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); length += dist; } } m_FiberLengths.push_back(length); m_MeanFiberLength += length; if (i==0) { m_MinFiberLength = length; m_MaxFiberLength = length; } else { if (lengthm_MaxFiberLength) m_MaxFiberLength = length; } } m_MeanFiberLength /= m_NumFibers; std::vector< float > sortedLengths = m_FiberLengths; std::sort(sortedLengths.begin(), sortedLengths.end()); for (int i=0; i1) m_LengthStDev /= (m_NumFibers-1); else m_LengthStDev = 0; m_LengthStDev = std::sqrt(m_LengthStDev); m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); // provide some border margin for(int i=0; i<=4; i+=2) b[i] -=10; for(int i=1; i<=5; i+=2) b[i] +=10; mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetFloatBounds(b); this->SetGeometry(geometry); } QStringList mitk::FiberBundleX::GetAvailableColorCodings() { QStringList availableColorCodings; int numColors = m_FiberPolyData->GetPointData()->GetNumberOfArrays(); for(int i=0; iGetPointData()->GetArrayName(i)); } //this controlstructure shall be implemented by the calling method if (availableColorCodings.isEmpty()) MITK_DEBUG << "no colorcodings available in fiberbundleX"; return availableColorCodings; } char* mitk::FiberBundleX::GetCurrentColorCoding() { return m_CurrentColorCoding; } void mitk::FiberBundleX::SetColorCoding(const char* requestedColorCoding) { if (requestedColorCoding==NULL) return; MITK_DEBUG << "SetColorCoding:" << requestedColorCoding; if( strcmp (COLORCODING_ORIENTATION_BASED,requestedColorCoding) == 0 ) { this->m_CurrentColorCoding = (char*) COLORCODING_ORIENTATION_BASED; } else if( strcmp (COLORCODING_FA_BASED,requestedColorCoding) == 0 ) { this->m_CurrentColorCoding = (char*) COLORCODING_FA_BASED; } else if( strcmp (COLORCODING_CUSTOM,requestedColorCoding) == 0 ) { this->m_CurrentColorCoding = (char*) COLORCODING_CUSTOM; } else { MITK_DEBUG << "FIBERBUNDLE X: UNKNOWN COLORCODING in FIBERBUNDLEX Datastructure"; this->m_CurrentColorCoding = (char*) COLORCODING_CUSTOM; //will cause blank colorcoding of fibers } } void mitk::FiberBundleX::RotateAroundAxis(double x, double y, double z) { MITK_INFO << "Rotating fibers"; x = x*M_PI/180; y = y*M_PI/180; z = z*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; mitk::Geometry3D::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::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]); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rotZ*rotY*rotX*dir; dir[0] += center[0]; dir[1] += center[1]; dir[2] += center[2]; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); } void mitk::FiberBundleX::ScaleFibers(double x, double y, double z) { MITK_INFO << "Scaling fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::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[0] *= x; p[1] *= y; p[2] *= z; 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(); } void mitk::FiberBundleX::TranslateFibers(double x, double y, double z) { MITK_INFO << "Translating fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::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[0] += x; p[1] += y; p[2] += z; 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(); } void mitk::FiberBundleX::MirrorFibers(unsigned int axis) { if (axis>2) return; MITK_INFO << "Mirroring fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); for (int i=0; iGetNextCell ( numPoints, pointIds ); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(pointIds[j]); p[axis] = -p[axis]; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); } bool mitk::FiberBundleX::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { if (minRadius<0) return true; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer vtkOldCells = m_FiberPolyData->GetLines(); vtkOldCells->InitTraversal(); MITK_INFO << "Applying curvature threshold"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkIdType numPoints(0); vtkIdType* points(NULL); vtkOldCells->GetNextCell ( numPoints, points ); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(points[j], p1); double p2[3]; m_FiberPolyData->GetPoint(points[j+1], p2); double p3[3]; m_FiberPolyData->GetPoint(points[j+2], p3); vnl_vector_fixed< float, 3 > v1, v2, v3; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; v2[0] = p3[0]-p2[0]; v2[1] = p3[1]-p2[1]; v2[2] = p3[2]-p2[2]; v3[0] = p1[0]-p3[0]; v3[1] = p1[1]-p3[1]; v3[2] = p1[2]-p3[2]; float a = v1.magnitude(); float b = v2.magnitude(); float c = v3.magnitude(); float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); if (deleteFibers && rInsertNextPoint(p2); // container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); container = vtkSmartPointer::New(); } else if (j==numPoints-3) { id = vtkNewPoints->InsertNextPoint(p2); container->GetPointIds()->InsertNextId(id); id = vtkNewPoints->InsertNextPoint(p3); container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); } } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } bool mitk::FiberBundleX::RemoveShortFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers return false; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); float min = m_MaxFiberLength; MITK_INFO << "Removing short fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetNextCell ( numPoints, pointIds ); if (m_FiberLengths.at(i)>=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(pointIds[j]); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); if (m_FiberLengths.at(i)GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } bool mitk::FiberBundleX::RemoveLongFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) return true; if (lengthInMM vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); MITK_INFO << "Removing long fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetNextCell ( numPoints, pointIds ); if (m_FiberLengths.at(i)<=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(pointIds[j]); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); UpdateColorCoding(); UpdateFiberGeometry(); return true; } void mitk::FiberBundleX::DoFiberSmoothing(int pointsPerCm, double tension, double continuity, double bias ) { vtkSmartPointer vtkSmoothPoints = vtkSmartPointer::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 = vtkSmartPointer::New(); //cellcontainer for smoothed lines vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); vtkIdType pointHelperCnt = 0; MITK_INFO << "Resampling fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetNextCell ( numPoints, pointIds ); vtkSmartPointer points = vtkSmartPointer::New(); for (int j=0; jInsertNextPoint(m_FiberPolyData->GetPoint(pointIds[j])); float length = m_FiberLengths.at(i); length /=10; int sampling = pointsPerCm*length; vtkSmartPointer xSpline = vtkSmartPointer::New(); vtkSmartPointer ySpline = vtkSmartPointer::New(); vtkSmartPointer zSpline = vtkSmartPointer::New(); xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); vtkSmartPointer spline = vtkSmartPointer::New(); spline->SetXSpline(xSpline); spline->SetYSpline(ySpline); spline->SetZSpline(zSpline); spline->SetPoints(points); vtkSmartPointer functionSource = vtkSmartPointer::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 = vtkSmartPointer::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(); m_FiberSampling = pointsPerCm; } void mitk::FiberBundleX::DoFiberSmoothing(int pointsPerCm) { DoFiberSmoothing(pointsPerCm, 0, 0, 0 ); } // Resample fiber to get equidistant points void mitk::FiberBundleX::ResampleFibers(float pointDistance) { if (pointDistance<=0.00001) return; vtkSmartPointer newPoly = vtkSmartPointer::New(); vtkSmartPointer newCellArray = vtkSmartPointer::New(); vtkSmartPointer newPoints = vtkSmartPointer::New(); vtkSmartPointer vLines = m_FiberPolyData->GetLines(); vLines->InitTraversal(); int numberOfLines = m_NumFibers; MITK_INFO << "Resampling fibers"; boost::progress_display disp(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(); m_FiberSampling = 10/pointDistance; } // 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/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.h index c369576fa4..3588003fdd 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleX.h @@ -1,155 +1,157 @@ /*=================================================================== 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 "FiberTrackingExports.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 FiberTracking_EXPORT FiberBundleX : public BaseData { public: + typedef itk::Image ItkUcharImgType; + // 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); void DoFiberSmoothing(int pointsPerCm, double tension, double continuity, double bias ); bool RemoveShortFibers(float lengthInMM); bool RemoveLongFibers(float lengthInMM); bool ApplyCurvatureThreshold(float minRadius, bool deleteFibers); void MirrorFibers(unsigned int axis); void RotateAroundAxis(double x, double y, double z); void TranslateFibers(double x, double y, double z); void ScaleFibers(double x, double y, double z); // 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); + FiberBundleX::Pointer ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint); + 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) itkGetMacro( FiberSampling, int) itkGetMacro( MinFiberLength, float ) itkGetMacro( MaxFiberLength, float ) itkGetMacro( MeanFiberLength, float ) itkGetMacro( MedianFiberLength, float ) itkGetMacro( LengthStDev, float ) std::vector GetPointsRoi() { return m_PointsRoi; } // copy fiber bundle mitk::FiberBundleX::Pointer GetDeepCopy(); // compare fiber bundles bool Equals(FiberBundleX* fib); 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; std::vector< float > m_FiberLengths; float m_MinFiberLength; float m_MaxFiberLength; float m_MeanFiberLength; float m_MedianFiberLength; float m_LengthStDev; int m_FiberSampling; std::vector m_PointsRoi; // this global variable needs to be refactored }; } // namespace mitk #endif /* _MITK_FiberBundleX_H */ 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 0869c112ea..dfa9bd2b38 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.cpp @@ -1,1751 +1,1703 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "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_CircleCounter(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(DoImageColorCoding())); connect( m_Controls->m_PruneFibersButton, SIGNAL(clicked()), this, SLOT(PruneBundle()) ); connect( m_Controls->m_CurvatureThresholdButton, SIGNAL(clicked()), this, SLOT(ApplyCurvatureThreshold()) ); connect( m_Controls->m_MirrorFibersButton, SIGNAL(clicked()), this, SLOT(MirrorFibers()) ); + connect( m_Controls->m_ExtractMask, SIGNAL(clicked()), this, SLOT(ExtractMask()) ); } } - -void QmitkFiberProcessingView::Extract3d() +void QmitkFiberProcessingView::ExtractMask() { - std::vector nodes = this->GetDataManagerSelection(); - if (nodes.empty()) + if (m_MaskImageNode.IsNull()) 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(m_MaskImageNode->GetData()); + 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(); + mitk::FiberBundleX::Pointer fib = dynamic_cast(m_SelectedFB.at(i)->GetData()); + QString name(m_SelectedFB.at(i)->GetName().c_str()); - vtkSmartPointer newPolyComplement = vtkSmartPointer::New(); - vtkSmartPointer newCellArrayComplement = vtkSmartPointer::New(); - vtkSmartPointer newPointsComplement = vtkSmartPointer::New(); + itkUCharImageType::Pointer mask = itkUCharImageType::New(); + mitk::CastToItkImage(mitkMask, mask); + mitk::FiberBundleX::Pointer newFib = fib->ExtractFiberSubset(mask, false); + DataNode::Pointer newNode = DataNode::New(); + newNode->SetData(newFib); + name += "_ending-in-mask"; + newNode->SetName(name.toStdString()); + GetDefaultDataStorage()->Add(newNode); + } +} - vtkSmartPointer vLines = polyFib->GetLines(); +void QmitkFiberProcessingView::Extract3d() +{ + if (m_MaskImageNode.IsNull()) + return; - vLines->InitTraversal(); - int numberOfLines = vLines->GetNumberOfCells(); - // each line - for (int j=0; j(m_MaskImageNode->GetData()); + for (int i=0; iGetNextCell ( numPoints, points ); - bool isPassing = false; + mitk::FiberBundleX::Pointer fib = dynamic_cast(m_SelectedFB.at(i)->GetData()); + QString name(m_SelectedFB.at(i)->GetName().c_str()); - // 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); - } + itkUCharImageType::Pointer mask = itkUCharImageType::New(); + mitk::CastToItkImage(mitkMask, mask); + mitk::FiberBundleX::Pointer newFib = fib->ExtractFiberSubset(mask, true); + DataNode::Pointer newNode = DataNode::New(); + newNode->SetData(newFib); + name += "_passing-mask"; + newNode->SetName(name.toStdString()); + GetDefaultDataStorage()->Add(newNode); } - - 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_DEBUG << "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_DEBUG << "Resampling requested image plane ... "; resampler->Update(); MITK_DEBUG << " ... 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_DEBUG << "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() { + m_Controls->m_Extract3dButton->setEnabled(false); + m_Controls->m_ExtractMask->setEnabled(false); + // are fiber bundles selected? if ( m_SelectedFB.empty() ) { m_Controls->m_InputData->setTitle("Please Select Input Data"); 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_CurvatureThresholdButton->setEnabled(false); if (m_SelectedSurfaces.size()>0) m_Controls->m_MirrorFibersButton->setEnabled(true); else m_Controls->m_MirrorFibersButton->setEnabled(false); } else { m_Controls->m_InputData->setTitle("Input Data"); 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_CurvatureThresholdButton->setEnabled(true); m_Controls->m_MirrorFibersButton->setEnabled(true); - if (m_SelectedSurfaces.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); + + if (m_MaskImageNode.IsNotNull()) + { + m_Controls->m_Extract3dButton->setEnabled(true); + m_Controls->m_ExtractMask->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 ) { //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection m_SelectedFB.clear(); m_SelectedPF.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = NULL; + m_MaskImageNode = NULL; m_Controls->m_FibLabel->setText("mandatory"); m_Controls->m_PfLabel->setText("needed for extraction"); for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( dynamic_cast(node->GetData()) ) { m_Controls->m_FibLabel->setText(node->GetName().c_str()); m_SelectedFB.push_back(node); } else if (dynamic_cast(node->GetData())) { m_Controls->m_PfLabel->setText(node->GetName().c_str()); m_SelectedPF.push_back(node); } else if (dynamic_cast(node->GetData())) + { m_SelectedImage = dynamic_cast(node->GetData()); + + bool isBinary = false; + node->GetPropertyValue("binary", isBinary); + if (isBinary) + m_MaskImageNode = node; + } else if (dynamic_cast(node->GetData())) { m_Controls->m_PfLabel->setText(node->GetName().c_str()); m_SelectedSurfaces.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_DEBUG << "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() { mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); this->GetDataStorage()->Modified(); 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() { } 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_DEBUG << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_DEBUG << "[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_DEBUG << 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_DEBUG << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; } else { MITK_DEBUG << "[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_DEBUG << savedPFchildNode->GetName() << " still exists"; } MITK_DEBUG << "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_DEBUG << "adding " << newPFCNode->GetName() << " to " << parentDataNode->GetName() ; GetDataStorage()->Add(newPFCNode, parentDataNode); } else { MITK_DEBUG << "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_DEBUG << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_DEBUG << "[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_DEBUG << 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_DEBUG << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_DEBUG << "[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_DEBUG << savedPFchildNode->GetName() << " still exists"; } MITK_DEBUG << "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_DEBUG << "adding " << newPFCNode->GetName() << " to " << parentDataNode->GetName() ; GetDataStorage()->Add(newPFCNode, parentDataNode); } else { MITK_DEBUG << "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_DEBUG << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_DEBUG << "[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_DEBUG << 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_DEBUG << savedPFchildNode->GetName() << " exists in DS...trying to remove it"; }else{ MITK_DEBUG << "[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_DEBUG << savedPFchildNode->GetName() << " still exists"; } MITK_DEBUG << "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_DEBUG << "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(); int maxLength = this->m_Controls->m_MaxPruneFibersSpinBox->value(); 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 if (!fib->RemoveLongFibers(maxLength)) QMessageBox::information(NULL, "No output generated:", "The resulting fiber bundle contains no fibers."); } GenerateStats(); RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyCurvatureThreshold() { int mm = this->m_Controls->m_MinCurvatureRadiusBox->value(); for (int i=0; i(m_SelectedFB.at(i)->GetData()); if (!fib->ApplyCurvatureThreshold(mm, this->m_Controls->m_RemoveFiberDueToCurvatureCheckbox->isChecked())) QMessageBox::information(NULL, "No output generated:", "The resulting fiber bundle contains no fibers."); } 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()); stats += "Number of fibers: "+ QString::number(fib->GetNumFibers()) + "\n"; stats += "Min. length: "+ QString::number(fib->GetMinFiberLength(),'f',1) + " mm\n"; stats += "Max. length: "+ QString::number(fib->GetMaxFiberLength(),'f',1) + " mm\n"; stats += "Mean length: "+ QString::number(fib->GetMeanFiberLength(),'f',1) + " mm\n"; stats += "Median length: "+ QString::number(fib->GetMedianFiberLength(),'f',1) + " mm\n"; stats += "Standard deviation: "+ QString::number(fib->GetLengthStDev(),'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, true); name += "_TDI"; break; case 1: newNode = GenerateTractDensityImage(fib, false, false); name += "_TDI"; break; case 2: newNode = GenerateTractDensityImage(fib, true, false); name += "_envelope"; break; case 3: newNode = GenerateColorHeatmap(fib); break; case 4: newNode = GenerateFiberEndingsImage(fib); name += "_fiber_endings"; break; case 5: 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->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, bool absolute) { const mitk::Geometry2D* bla = GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderer()->GetCurrentWorldGeometry2D(); typedef float OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(binary); generator->SetOutputAbsoluteValues(absolute); 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); } if (m_SelectedFB.size()>0) GenerateStats(); if (m_SelectedSurfaces.size()>0) { for (int i=0; i poly = surf->GetVtkPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); for (int i=0; iGetNumberOfPoints(); i++) { double* point = poly->GetPoint(i); point[axis] *= -1; vtkNewPoints->InsertNextPoint(point); } poly->SetPoints(vtkNewPoints); surf->CalculateBoundingBox(); } } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoImageColorCoding() { if (m_SelectedImage.IsNull()) 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 c3e6a43d98..806782e4d0 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingView.h @@ -1,176 +1,178 @@ /*=================================================================== 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 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 /*! \brief 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; 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(); ///< add circle interactors etc. void OnDrawPolygon(); ///< add circle interactors etc. void DoFiberExtraction(); ///< Extract fibers from selected bundle void GenerateAndComposite(); void GenerateOrComposite(); void GenerateNotComposite(); void PruneBundle(); ///< remove too short/too long fibers void MirrorFibers(); ///< mirror bundle on the specified plane void JoinBundles(); ///< merge selected fiber bundles void SubstractBundles(); ///< subtract bundle A from bundle B. Not commutative! Defined by order of selection. void GenerateRoiImage(); ///< generate binary image of selected planar figures. void ProcessSelectedBundles(); ///< start selected operation on fiber bundle (e.g. tract density image generation) void ResampleSelectedBundles(); ///< smooth fiber bundle using the specified number of sampling points per cm. void DoImageColorCoding(); ///< color fibers by selected scalar image void Extract3d(); ///< extract all fibers passing the selected surface mesh void ApplyCurvatureThreshold(); ///< remove/split fibers with a too high curvature threshold + void ExtractMask(); ///< extract all fibers passing the selected surface mesh 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(); ///< generate statistics of selected fiber bundles void UpdateGui(); ///< update button activity etc. dpending on current datamanager selection int m_CircleCounter; ///< used for data node naming int m_PolygonCounter; ///< used for data node naming std::vector m_SelectedFB; ///< selected fiber bundle nodes std::vector m_SelectedPF; ///< selected planar figure nodes std::vector m_SelectedSurfaces; mitk::Image::Pointer m_SelectedImage; mitk::Image::Pointer m_InternalImage; mitk::PlanarFigure::Pointer m_PlanarFigure; itkUCharImageType::Pointer m_InternalImageMask3D; itkUCharImageType::Pointer m_PlanarFigureImage; float m_UpsamplingFactor; ///< upsampling factor for all image generations + mitk::DataNode::Pointer m_MaskImageNode; 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, bool absolute); 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 179ceb0366..f6417c6711 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui @@ -1,954 +1,984 @@ QmitkFiberProcessingViewControls 0 0 492 - 809 + 866 Form Please Select Input Data <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> <html><head/><body><p><span style=" color:#969696;">needed for extraction</span></p></body></html> Input DTI Fiber Bundle: Binary seed ROI. If not specified, the whole image area is seeded. ROI: Fiber Extraction 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. + Returns all fibers contained in bundle X that are not contained in bundle Y (not commutative!). Select at least two fiber bundles to execute. - Extract + Substract - + 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. + Extract fibers starting/ending inside of the selected binary mask. - Substract + Ending in Mask - + + + + Qt::Horizontal + + + + 40 + 20 + + + + + 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 fibers passing through selected ROI or composite ROI. Select ROI and fiber bundle to execute. - Extract 3D + Extract - + 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 + + + + false + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + Extract fibers passing through the selected binary mask. + + + Passing Mask + + + 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 Processing QFormLayout::AllNonFixedFieldsGrow false 0 0 200 16777215 11 Perform selected operation on all selected fiber bundles. Generate Image QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Tract Density Image (TDI) Normalized TDI Binary Envelope Fiber Bundle Image Fiber Endings Image Fiber Endings Pointset Upsampling factor 1 0.100000000000000 10.000000000000000 0.100000000000000 1.000000000000000 false 0 0 200 16777215 11 Resample fibers using a Kochanek spline interpolation. Smooth Fibers Points per cm 1 50 5 false 0 0 200 16777215 11 Remove fibers shorter/longer than the specified length (in mm). Length Threshold QFrame::NoFrame QFrame::Raised 0 0 Minimum fiber length in mm 0 1000 20 Maximum fiber length in mm 0 10000 500 false 0 0 200 16777215 11 Remove fibers with a too high curvature Curvature Threshold QFrame::NoFrame QFrame::Raised 6 0 0 Minimum radius of circle created by three consecutive points of a fiber 100.000000000000000 0.100000000000000 2.000000000000000 Remove whole fiber if it is exceeding the curvature threshold, otherwise remove only high curvature part. Remove Fiber true false 0 0 200 16777215 11 Mirror fibers around specified axis. Mirror Fibers 0 3 3 Sagittal Coronal Axial false 0 0 200 16777215 11 Apply float image values (0-1) as color coding to the selected fiber bundle. Color By Scalar Map Fiber Statistics Courier 10 Pitch false true Qt::Vertical 20 40 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp index 41ce65dcc8..2b0ccd2ff1 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberfoxView.cpp @@ -1,1321 +1,1324 @@ /*=================================================================== 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. ===================================================================*/ //misc #define _USE_MATH_DEFINES #include // Blueberry #include #include // Qmitk #include "QmitkFiberfoxView.h" // MITK #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 #define _USE_MATH_DEFINES #include const std::string QmitkFiberfoxView::VIEW_ID = "org.mitk.views.fiberfoxview"; QmitkFiberfoxView::QmitkFiberfoxView() : QmitkAbstractView() , m_Controls( 0 ) , m_SelectedImage( NULL ) , m_SelectedBundle( NULL ) { } // Destructor QmitkFiberfoxView::~QmitkFiberfoxView() { } void QmitkFiberfoxView::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::QmitkFiberfoxViewControls; m_Controls->setupUi( parent ); m_Controls->m_VarianceBox->setVisible(false); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_T2bluringParamFrame->setVisible(false); m_Controls->m_KspaceParamFrame->setVisible(false); m_Controls->m_StickModelFrame->setVisible(false); m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); connect((QObject*) m_Controls->m_GenerateImageButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateImage())); connect((QObject*) m_Controls->m_GenerateFibersButton, SIGNAL(clicked()), (QObject*) this, SLOT(GenerateFibers())); connect((QObject*) m_Controls->m_CircleButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnDrawROI())); connect((QObject*) m_Controls->m_FlipButton, SIGNAL(clicked()), (QObject*) this, SLOT(OnFlipButton())); connect((QObject*) m_Controls->m_JoinBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(JoinBundles())); connect((QObject*) m_Controls->m_VarianceBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnVarianceChanged(double))); connect((QObject*) m_Controls->m_DistributionBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(OnDistributionChanged(int))); connect((QObject*) m_Controls->m_FiberDensityBox, SIGNAL(valueChanged(int)), (QObject*) this, SLOT(OnFiberDensityChanged(int))); connect((QObject*) m_Controls->m_FiberSamplingBox, SIGNAL(valueChanged(int)), (QObject*) this, SLOT(OnFiberSamplingChanged(int))); connect((QObject*) m_Controls->m_TensionBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnTensionChanged(double))); connect((QObject*) m_Controls->m_ContinuityBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnContinuityChanged(double))); connect((QObject*) m_Controls->m_BiasBox, SIGNAL(valueChanged(double)), (QObject*) this, SLOT(OnBiasChanged(double))); connect((QObject*) m_Controls->m_AddT2Smearing, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddT2Smearing(int))); connect((QObject*) m_Controls->m_AddGibbsRinging, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnAddGibbsRinging(int))); connect((QObject*) m_Controls->m_ConstantRadiusBox, SIGNAL(stateChanged(int)), (QObject*) this, SLOT(OnConstantRadius(int))); connect((QObject*) m_Controls->m_CopyBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(CopyBundles())); connect((QObject*) m_Controls->m_TransformBundlesButton, SIGNAL(clicked()), (QObject*) this, SLOT(ApplyTransform())); connect((QObject*) m_Controls->m_AlignOnGrid, SIGNAL(clicked()), (QObject*) this, SLOT(AlignOnGrid())); connect((QObject*) m_Controls->m_FiberCompartmentModelBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(FiberModelFrameVisibility(int))); connect((QObject*) m_Controls->m_NonFiberCompartmentModelBox, SIGNAL(currentIndexChanged(int)), (QObject*) this, SLOT(FiberModelFrameVisibility(int))); connect((QObject*) m_Controls->m_AdvancedFiberOptionsBox, SIGNAL( stateChanged(int)), (QObject*) this, SLOT(ShowAdvancedFiberOptions(int))); } } void QmitkFiberfoxView::ShowAdvancedFiberOptions(int state) { if (state) m_Controls->m_AdvancedFiberOptionsFrame->setVisible(true); else m_Controls->m_AdvancedFiberOptionsFrame->setVisible(false); } void QmitkFiberfoxView::FiberModelFrameVisibility(int index) { m_Controls->m_TensorModelFrame->setVisible(false); m_Controls->m_StickModelFrame->setVisible(false); switch (index) { case 0: m_Controls->m_TensorModelFrame->setVisible(true); break; case 1: m_Controls->m_StickModelFrame->setVisible(true); break; default: m_Controls->m_TensorModelFrame->setVisible(true); } } void QmitkFiberfoxView::NonFiberModelFrameVisibility(int index) { } void QmitkFiberfoxView::OnConstantRadius(int value) { if (value>0 && m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnAddT2Smearing(int value) { if (value>0) m_Controls->m_T2bluringParamFrame->setVisible(true); else m_Controls->m_T2bluringParamFrame->setVisible(false); } void QmitkFiberfoxView::OnAddGibbsRinging(int value) { if (value>0) m_Controls->m_KspaceParamFrame->setVisible(true); else m_Controls->m_KspaceParamFrame->setVisible(false); } void QmitkFiberfoxView::OnDistributionChanged(int value) { if (value==1) m_Controls->m_VarianceBox->setVisible(true); else m_Controls->m_VarianceBox->setVisible(false); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnVarianceChanged(double value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberDensityChanged(int value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFiberSamplingChanged(int value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnTensionChanged(double value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnContinuityChanged(double value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnBiasChanged(double value) { if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::AlignOnGrid() { for (int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::DataStorage::SetOfObjects::ConstPointer parentFibs = GetDataStorage()->GetSources(m_SelectedFiducials.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = parentFibs->begin(); it != parentFibs->end(); ++it ) { mitk::DataNode::Pointer pFibNode = *it; if ( pFibNode.IsNotNull() && dynamic_cast(pFibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(pFibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { mitk::Image::Pointer img = dynamic_cast(pImgNode->GetData()); mitk::Geometry3D::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); break; } } break; } } } for( int i=0; iGetSources(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it = sources->begin(); it != sources->end(); ++it ) { mitk::DataNode::Pointer imgNode = *it; if ( imgNode.IsNotNull() && dynamic_cast(imgNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::Image::Pointer img = dynamic_cast(imgNode->GetData()); mitk::Geometry3D::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } break; } } } for( int i=0; i(m_SelectedImages.at(i)->GetData()); mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) { mitk::DataStorage::SetOfObjects::ConstPointer derivations2 = GetDataStorage()->GetDerivations(fibNode); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations2->begin(); it2 != derivations2->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = dynamic_cast(fiducialNode->GetData()); mitk::Point3D wc0 = pe->GetWorldControlPoint(0); mitk::Geometry3D::Pointer geom = img->GetGeometry(); itk::Index<3> idx; geom->WorldToIndex(wc0, idx); mitk::Point3D cIdx; cIdx[0]=idx[0]; cIdx[1]=idx[1]; cIdx[2]=idx[2]; mitk::Point3D world; geom->IndexToWorld(cIdx,world); mitk::Vector3D trans = world - wc0; pe->GetGeometry()->Translate(trans); } } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::OnFlipButton() { if (m_SelectedFiducial.IsNull()) return; std::map::iterator it = m_DataNodeToPlanarFigureData.find(m_SelectedFiducial.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; data.m_Flipped += 1; data.m_Flipped %= 2; } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } QmitkFiberfoxView::GradientListType QmitkFiberfoxView::GenerateHalfShell(int NPoints) { NPoints *= 2; GradientListType pointshell; int numB0 = NPoints/10; if (numB0==0) numB0=1; GradientType g; g.Fill(0.0); for (int i=0; i theta; theta.set_size(NPoints); vnl_vector phi; phi.set_size(NPoints); double C = sqrt(4*M_PI); phi(0) = 0.0; phi(NPoints-1) = 0.0; for(int i=0; i0 && i std::vector > QmitkFiberfoxView::MakeGradientList() { std::vector > retval; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); // Add 0 vector for B0 int numB0 = ndirs/10; if (numB0==0) numB0=1; itk::Vector v; v.Fill(0.0); for (int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval.push_back(v); } return retval; } void QmitkFiberfoxView::OnAddBundle() { if (m_SelectedImage.IsNull()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedImage); mitk::FiberBundleX::Pointer bundle = mitk::FiberBundleX::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( bundle ); QString name = QString("Bundle_%1").arg(children->size()); node->SetName(name.toStdString()); m_SelectedBundle = node; m_SelectedBundles.push_back(node); UpdateGui(); GetDataStorage()->Add(node, m_SelectedImage); } void QmitkFiberfoxView::OnDrawROI() { if (m_SelectedBundle.IsNull()) OnAddBundle(); if (m_SelectedBundle.IsNull()) return; mitk::DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(m_SelectedBundle); mitk::PlanarEllipse::Pointer figure = mitk::PlanarEllipse::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( figure ); QList nodes = this->GetDataManagerSelection(); for( int i=0; iSetSelected(false); m_SelectedFiducial = node; QString name = QString("Fiducial_%1").arg(children->size()); node->SetName(name.toStdString()); node->SetSelected(true); GetDataStorage()->Add(node, m_SelectedBundle); this->DisableCrosshairNavigation(); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetInteractor()); if(figureInteractor.IsNull()) figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); UpdateGui(); } bool CompareLayer(mitk::DataNode::Pointer i,mitk::DataNode::Pointer j) { int li = -1; i->GetPropertyValue("layer", li); int lj = -1; j->GetPropertyValue("layer", lj); return liGetSources(m_SelectedFiducial); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) if(dynamic_cast((*it)->GetData())) m_SelectedBundles.push_back(*it); if (m_SelectedBundles.empty()) return; } vector< vector< mitk::PlanarEllipse::Pointer > > fiducials; vector< vector< unsigned int > > fliplist; for (int i=0; iGetDerivations(m_SelectedBundles.at(i)); std::vector< mitk::DataNode::Pointer > childVector; for( mitk::DataStorage::SetOfObjects::const_iterator it = children->begin(); it != children->end(); ++it ) childVector.push_back(*it); sort(childVector.begin(), childVector.end(), CompareLayer); vector< mitk::PlanarEllipse::Pointer > fib; vector< unsigned int > flip; float radius = 1; int count = 0; for( std::vector< mitk::DataNode::Pointer >::const_iterator it = childVector.begin(); it != childVector.end(); ++it ) { mitk::DataNode::Pointer node = *it; if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { mitk::PlanarEllipse* ellipse = dynamic_cast(node->GetData()); if (m_Controls->m_ConstantRadiusBox->isChecked()) { ellipse->SetTreatAsCircle(true); mitk::Point2D c = ellipse->GetControlPoint(0); mitk::Point2D p = ellipse->GetControlPoint(1); mitk::Vector2D v = p-c; if (count==0) { radius = v.GetVnlVector().magnitude(); ellipse->SetControlPoint(1, p); } else { v.Normalize(); v *= radius; ellipse->SetControlPoint(1, c+v); } } fib.push_back(ellipse); std::map::iterator it = m_DataNodeToPlanarFigureData.find(node.GetPointer()); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; flip.push_back(data.m_Flipped); } else flip.push_back(0); } count++; } if (fib.size()>1) { fiducials.push_back(fib); fliplist.push_back(flip); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); if (fib.size()<3) return; } itk::FibersFromPlanarFiguresFilter::Pointer filter = itk::FibersFromPlanarFiguresFilter::New(); filter->SetFiducials(fiducials); filter->SetFlipList(fliplist); switch(m_Controls->m_DistributionBox->currentIndex()){ case 0: filter->SetFiberDistribution(itk::FibersFromPlanarFiguresFilter::DISTRIBUTE_UNIFORM); break; case 1: filter->SetFiberDistribution(itk::FibersFromPlanarFiguresFilter::DISTRIBUTE_GAUSSIAN); filter->SetVariance(m_Controls->m_VarianceBox->value()); break; } filter->SetDensity(m_Controls->m_FiberDensityBox->value()); filter->SetTension(m_Controls->m_TensionBox->value()); filter->SetContinuity(m_Controls->m_ContinuityBox->value()); filter->SetBias(m_Controls->m_BiasBox->value()); filter->SetFiberSampling(m_Controls->m_FiberSamplingBox->value()); filter->Update(); vector< mitk::FiberBundleX::Pointer > fiberBundles = filter->GetFiberBundles(); for (int i=0; iSetData( fiberBundles.at(i) ); if (fiberBundles.at(i)->GetNumFibers()>50000) m_SelectedBundles.at(i)->SetVisibility(false); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::GenerateImage() { itk::ImageRegion<3> imageRegion; imageRegion.SetSize(0, m_Controls->m_SizeX->value()); imageRegion.SetSize(1, m_Controls->m_SizeY->value()); imageRegion.SetSize(2, m_Controls->m_SizeZ->value()); mitk::Vector3D spacing; spacing[0] = m_Controls->m_SpacingX->value(); spacing[1] = m_Controls->m_SpacingY->value(); spacing[2] = m_Controls->m_SpacingZ->value(); mitk::Point3D origin; origin[0] = spacing[0]/2; origin[1] = spacing[1]/2; origin[2] = spacing[2]/2; itk::Matrix directionMatrix; directionMatrix.SetIdentity(); if (m_SelectedBundle.IsNull()) { mitk::Image::Pointer image = mitk::ImageGenerator::GenerateGradientImage( m_Controls->m_SizeX->value(), m_Controls->m_SizeY->value(), m_Controls->m_SizeZ->value(), m_Controls->m_SpacingX->value(), m_Controls->m_SpacingY->value(), m_Controls->m_SpacingZ->value()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName("Dummy"); GetDataStorage()->Add(node); m_SelectedImage = node; mitk::BaseData::Pointer basedata = node->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeSlicedGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } UpdateGui(); return; } DiffusionSignalModel::GradientListType gradientList; double bVal = 1000; if (m_SelectedDWI.IsNull()) { gradientList = GenerateHalfShell(m_Controls->m_NumGradientsBox->value());; bVal = m_Controls->m_BvalueBox->value(); } else { mitk::DiffusionImage::Pointer dwi = dynamic_cast*>(m_SelectedDWI->GetData()); imageRegion = dwi->GetVectorImage()->GetLargestPossibleRegion(); spacing = dwi->GetVectorImage()->GetSpacing(); origin = dwi->GetVectorImage()->GetOrigin(); directionMatrix = dwi->GetVectorImage()->GetDirection(); bVal = dwi->GetB_Value(); mitk::DiffusionImage::GradientDirectionContainerType::Pointer dirs = dwi->GetDirections(); for (int i=0; iSize(); i++) { DiffusionSignalModel::GradientType g; g[0] = dirs->at(i)[0]; g[1] = dirs->at(i)[1]; g[2] = dirs->at(i)[2]; gradientList.push_back(g); } } // storage for generated phantom image mitk::DataNode::Pointer resultNode = mitk::DataNode::New(); // signal models QString signalModelString("Ball"); itk::TractsToDWIImageFilter::DiffusionModelList fiberModelList, nonFiberModelList; mitk::TensorModel tensorModel; mitk::StickModel stickModel; // free diffusion mitk::BallModel ballModel; ballModel.SetGradientList(gradientList); ballModel.SetBvalue(bVal); ballModel.SetDiffusivity(m_Controls->m_BallD->value()); ballModel.SetSignalScale(m_Controls->m_NonFiberS0Box->value()); ballModel.SetRelaxationT2(m_Controls->m_NonFiberRelaxationT2Box->value()); nonFiberModelList.push_back(&ballModel); resultNode->AddProperty("Fiberfox.Ball.Diffusivity", DoubleProperty::New(m_Controls->m_BallD->value())); resultNode->AddProperty("Fiberfox.Ball.Scaling", DoubleProperty::New(m_Controls->m_NonFiberS0Box->value())); if (m_Controls->m_AddT2Smearing->isChecked()) resultNode->AddProperty("Fiberfox.Ball.T2", DoubleProperty::New(m_Controls->m_NonFiberRelaxationT2Box->value())); // intra-axonal diffusion switch (m_Controls->m_FiberCompartmentModelBox->currentIndex()) { case 0: MITK_INFO << "Using zeppelin model"; tensorModel.SetGradientList(gradientList); tensorModel.SetBvalue(bVal); tensorModel.SetKernelFA(m_Controls->m_TensorFaBox->value()); tensorModel.SetSignalScale(m_Controls->m_FiberS0Box->value()); tensorModel.SetRelaxationT2(m_Controls->m_FiberRelaxationT2Box->value()); fiberModelList.push_back(&tensorModel); signalModelString += "-Zeppelin"; resultNode->AddProperty("Fiberfox.Zeppelin.FA", DoubleProperty::New(m_Controls->m_TensorFaBox->value())); resultNode->AddProperty("Fiberfox.Zeppelin.Scaling", DoubleProperty::New(m_Controls->m_FiberS0Box->value())); if (m_Controls->m_AddT2Smearing->isChecked()) resultNode->AddProperty("Fiberfox.Zeppelin.T2", DoubleProperty::New(m_Controls->m_FiberRelaxationT2Box->value())); break; case 1: MITK_INFO << "Using stick model"; stickModel.SetGradientList(gradientList); stickModel.SetDiffusivity(m_Controls->m_StickDiffusivityBox->value()); stickModel.SetSignalScale(m_Controls->m_FiberS0Box->value()); stickModel.SetRelaxationT2(m_Controls->m_FiberRelaxationT2Box->value()); fiberModelList.push_back(&stickModel); signalModelString += "-Stick"; resultNode->AddProperty("Fiberfox.Stick.Diffusivity", DoubleProperty::New(m_Controls->m_StickDiffusivityBox->value())); resultNode->AddProperty("Fiberfox.Stick.Scaling", DoubleProperty::New(m_Controls->m_FiberS0Box->value())); if (m_Controls->m_AddT2Smearing->isChecked()) resultNode->AddProperty("Fiberfox.Stick.T2", DoubleProperty::New(m_Controls->m_FiberRelaxationT2Box->value())); break; } itk::TractsToDWIImageFilter::KspaceArtifactList artifactList; // noise model double snr = m_Controls->m_NoiseLevel->value(); double noiseVariance = 0; if (snr <= 0) snr = 0.0001; if (snr<=99) { noiseVariance = (double)m_Controls->m_FiberS0Box->value()/snr; noiseVariance *= noiseVariance; } mitk::RicianNoiseModel noiseModel; noiseModel.SetNoiseVariance(noiseVariance); // artifact models + QString artifactModelString(""); mitk::GibbsRingingArtifact gibbsModel; if (m_Controls->m_AddGibbsRinging->isChecked()) { + artifactModelString += "_Gibbs-ringing"; resultNode->AddProperty("Fiberfox.k-Space-Undersampling", IntProperty::New(m_Controls->m_KspaceUndersamplingBox->currentText().toInt())); gibbsModel.SetKspaceCropping((double)m_Controls->m_KspaceUndersamplingBox->currentText().toInt()); artifactList.push_back(&gibbsModel); } mitk::T2SmearingArtifact t2Model; if (m_Controls->m_AddT2Smearing->isChecked()) { + artifactModelString += "_T2-blurring"; t2Model.SetReadoutPulseLength(1); artifactList.push_back(&t2Model); } for (int i=0; i(m_SelectedBundles.at(i)->GetData()); if (fiberBundle->GetNumFibers()<=0) continue; itk::TractsToDWIImageFilter::Pointer filter = itk::TractsToDWIImageFilter::New(); filter->SetImageRegion(imageRegion); filter->SetSpacing(spacing); filter->SetOrigin(origin); filter->SetDirectionMatrix(directionMatrix); filter->SetFiberBundle(fiberBundle); filter->SetFiberModels(fiberModelList); filter->SetNonFiberModels(nonFiberModelList); filter->SetNoiseModel(&noiseModel); filter->SetKspaceArtifacts(artifactList); filter->SetNumberOfRepetitions(m_Controls->m_RepetitionsBox->value()); filter->SetEnforcePureFiberVoxels(m_Controls->m_EnforcePureFiberVoxelsBox->isChecked()); if (m_TissueMask.IsNotNull()) { ItkUcharImgType::Pointer mask = ItkUcharImgType::New(); mitk::CastToItkImage(m_TissueMask, mask); filter->SetTissueMask(mask); } filter->Update(); mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(bVal); image->SetDirections(gradientList); image->InitializeFromVectorImage(); // image->GetGeometry()->SetImageGeometry(false); resultNode->SetData( image ); resultNode->SetName(m_SelectedBundle->GetName() +"_D"+QString::number(m_Controls->m_SizeX->value()).toStdString() +"-"+QString::number(m_Controls->m_SizeY->value()).toStdString() +"-"+QString::number(m_Controls->m_SizeZ->value()).toStdString() +"_S"+QString::number(spacing[0]).toStdString() +"-"+QString::number(spacing[1]).toStdString() +"-"+QString::number(spacing[2]).toStdString() +"_b"+QString::number(bVal).toStdString() +"_SNR"+QString::number(snr).toStdString() +"_"+signalModelString.toStdString()); GetDataStorage()->Add(resultNode, m_SelectedBundle); resultNode->AddProperty("Fiberfox.SNR", DoubleProperty::New(snr)); resultNode->AddProperty("Fiberfox.Repetitions", IntProperty::New(m_Controls->m_RepetitionsBox->value())); resultNode->AddProperty("Fiberfox.b-value", DoubleProperty::New(bVal)); resultNode->AddProperty("Fiberfox.Model", StringProperty::New(signalModelString.toStdString())); if (m_Controls->m_KspaceImageBox->isChecked()) { itk::Image::Pointer kspace = filter->GetKspaceImage(); mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(kspace.GetPointer()); image->SetVolume(kspace->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName(m_SelectedBundle->GetName()+"_k-space"); GetDataStorage()->Add(node, m_SelectedBundle); } mitk::BaseData::Pointer basedata = resultNode->GetData(); if (basedata.IsNotNull()) { mitk::RenderingManager::GetInstance()->InitializeViews( basedata->GetTimeSlicedGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } } void QmitkFiberfoxView::ApplyTransform() { vector< mitk::DataNode::Pointer > selectedBundles; for( int i=0; iGetDerivations(m_SelectedImages.at(i)); for( mitk::DataStorage::SetOfObjects::const_iterator it = derivations->begin(); it != derivations->end(); ++it ) { mitk::DataNode::Pointer fibNode = *it; if ( fibNode.IsNotNull() && dynamic_cast(fibNode->GetData()) ) selectedBundles.push_back(fibNode); } } if (selectedBundles.empty()) selectedBundles = m_SelectedBundles2; if (!selectedBundles.empty()) { std::vector::const_iterator it = selectedBundles.begin(); for (it; it!=selectedBundles.end(); ++it) { mitk::FiberBundleX::Pointer fib = dynamic_cast((*it)->GetData()); fib->RotateAroundAxis(m_Controls->m_XrotBox->value(), m_Controls->m_YrotBox->value(), m_Controls->m_ZrotBox->value()); fib->TranslateFibers(m_Controls->m_XtransBox->value(), m_Controls->m_YtransBox->value(), m_Controls->m_ZtransBox->value()); fib->ScaleFibers(m_Controls->m_XscaleBox->value(), m_Controls->m_YscaleBox->value(), m_Controls->m_ZscaleBox->value()); // handle child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse* pe = dynamic_cast(fiducialNode->GetData()); mitk::Geometry3D* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< float, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< float, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< float, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< float, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); // implicit translation mitk::Vector3D trans; trans[0] = geom->GetOrigin()[0]-fib->GetGeometry()->GetCenter()[0]; trans[1] = geom->GetOrigin()[1]-fib->GetGeometry()->GetCenter()[1]; trans[2] = geom->GetOrigin()[2]-fib->GetGeometry()->GetCenter()[2]; mitk::Vector3D newWc = rot*trans; newWc = newWc-trans; geom->Translate(newWc); } } } } } else { for (int i=0; i(m_SelectedFiducials.at(i)->GetData()); mitk::Geometry3D* geom = pe->GetGeometry(); // translate mitk::Vector3D world; world[0] = m_Controls->m_XtransBox->value(); world[1] = m_Controls->m_YtransBox->value(); world[2] = m_Controls->m_ZtransBox->value(); geom->Translate(world); // calculate rotation matrix double x = m_Controls->m_XrotBox->value()*M_PI/180; double y = m_Controls->m_YrotBox->value()*M_PI/180; double z = m_Controls->m_ZrotBox->value()*M_PI/180; itk::Matrix< float, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; itk::Matrix< float, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; itk::Matrix< float, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< float, 3, 3 > rot = rotZ*rotY*rotX; // transform control point coordinate into geometry translation geom->SetOrigin(pe->GetWorldControlPoint(0)); mitk::Point2D cp; cp.Fill(0.0); pe->SetControlPoint(0, cp); // rotate fiducial geom->GetIndexToWorldTransform()->SetMatrix(rot*geom->GetIndexToWorldTransform()->GetMatrix()); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::CopyBundles() { if ( m_SelectedBundles.size()<1 ){ QMessageBox::information( NULL, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } std::vector::const_iterator it = m_SelectedBundles.begin(); for (it; it!=m_SelectedBundles.end(); ++it) { // find parent image mitk::DataNode::Pointer parentNode; mitk::DataStorage::SetOfObjects::ConstPointer parentImgs = GetDataStorage()->GetSources(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = parentImgs->begin(); it2 != parentImgs->end(); ++it2 ) { mitk::DataNode::Pointer pImgNode = *it2; if ( pImgNode.IsNotNull() && dynamic_cast(pImgNode->GetData()) ) { parentNode = pImgNode; break; } } mitk::FiberBundleX::Pointer fib = dynamic_cast((*it)->GetData()); mitk::FiberBundleX::Pointer newBundle = fib->GetDeepCopy(); QString name((*it)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(fbNode, parentNode); else GetDataStorage()->Add(fbNode); // copy child fiducials if (m_Controls->m_IncludeFiducials->isChecked()) { mitk::DataStorage::SetOfObjects::ConstPointer derivations = GetDataStorage()->GetDerivations(*it); for( mitk::DataStorage::SetOfObjects::const_iterator it2 = derivations->begin(); it2 != derivations->end(); ++it2 ) { mitk::DataNode::Pointer fiducialNode = *it2; if ( fiducialNode.IsNotNull() && dynamic_cast(fiducialNode->GetData()) ) { mitk::PlanarEllipse::Pointer pe = mitk::PlanarEllipse::New(); pe->DeepCopy(dynamic_cast(fiducialNode->GetData())); mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetData(pe); newNode->SetName(fiducialNode->GetName()); GetDataStorage()->Add(newNode, fbNode); } } } } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::JoinBundles() { if ( m_SelectedBundles.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_SelectedBundles.begin(); mitk::FiberBundleX::Pointer newBundle = dynamic_cast((*it)->GetData()); QString name(""); name += QString((*it)->GetName().c_str()); ++it; for (it; it!=m_SelectedBundles.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); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberfoxView::UpdateGui() { m_Controls->m_FiberBundleLabel->setText("mandatory"); m_Controls->m_GeometryFrame->setEnabled(true); m_Controls->m_GeometryMessage->setVisible(false); m_Controls->m_DiffusionPropsMessage->setVisible(false); m_Controls->m_FiberGenMessage->setVisible(true); m_Controls->m_TransformBundlesButton->setEnabled(false); m_Controls->m_CopyBundlesButton->setEnabled(false); m_Controls->m_GenerateFibersButton->setEnabled(false); m_Controls->m_FlipButton->setEnabled(false); m_Controls->m_CircleButton->setEnabled(false); m_Controls->m_BvalueBox->setEnabled(true); m_Controls->m_NumGradientsBox->setEnabled(true); m_Controls->m_JoinBundlesButton->setEnabled(false); m_Controls->m_AlignOnGrid->setEnabled(false); if (m_SelectedFiducial.IsNotNull()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_FlipButton->setEnabled(true); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_SelectedImage.IsNotNull() || m_SelectedBundle.IsNotNull()) { m_Controls->m_TransformBundlesButton->setEnabled(true); m_Controls->m_CircleButton->setEnabled(true); m_Controls->m_FiberGenMessage->setVisible(false); m_Controls->m_AlignOnGrid->setEnabled(true); } if (m_TissueMask.IsNotNull()) { m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (m_SelectedDWI.IsNotNull()) { m_Controls->m_DiffusionPropsMessage->setVisible(true); m_Controls->m_BvalueBox->setEnabled(false); m_Controls->m_NumGradientsBox->setEnabled(false); m_Controls->m_GeometryMessage->setVisible(true); m_Controls->m_GeometryFrame->setEnabled(false); } if (m_SelectedBundle.IsNotNull()) { m_Controls->m_CopyBundlesButton->setEnabled(true); m_Controls->m_GenerateFibersButton->setEnabled(true); m_Controls->m_FiberBundleLabel->setText(m_SelectedBundle->GetName().c_str()); if (m_SelectedBundles.size()>1) m_Controls->m_JoinBundlesButton->setEnabled(true); } } void QmitkFiberfoxView::OnSelectionChanged( berry::IWorkbenchPart::Pointer, const QList& nodes ) { m_SelectedBundles2.clear(); m_SelectedImages.clear(); m_SelectedFiducials.clear(); m_SelectedFiducial = NULL; m_TissueMask = NULL; m_SelectedBundles.clear(); m_SelectedBundle = NULL; m_SelectedImage = NULL; m_SelectedDWI = NULL; m_Controls->m_TissueMaskLabel->setText("optional"); // iterate all selected objects, adjust warning visibility for( int i=0; i*>(node->GetData()) ) { m_SelectedDWI = node; m_SelectedImage = node; m_SelectedImages.push_back(node); } else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedImages.push_back(node); m_SelectedImage = node; bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) { m_TissueMask = dynamic_cast(node->GetData()); m_Controls->m_TissueMaskLabel->setText(node->GetName().c_str()); } } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedBundles2.push_back(node); if (m_Controls->m_RealTimeFibers->isChecked() && node!=m_SelectedBundle) { m_SelectedBundle = node; m_SelectedBundles.push_back(node); mitk::FiberBundleX::Pointer newFib = dynamic_cast(node->GetData()); if (newFib->GetNumFibers()!=m_Controls->m_FiberDensityBox->value()) GenerateFibers(); } else { m_SelectedBundle = node; m_SelectedBundles.push_back(node); } } else if ( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_SelectedFiducials.push_back(node); m_SelectedFiducial = node; m_SelectedBundles.clear(); mitk::DataStorage::SetOfObjects::ConstPointer parents = GetDataStorage()->GetSources(node); for( mitk::DataStorage::SetOfObjects::const_iterator it = parents->begin(); it != parents->end(); ++it ) { mitk::DataNode::Pointer pNode = *it; if ( pNode.IsNotNull() && dynamic_cast(pNode->GetData()) ) { m_SelectedBundle = pNode; m_SelectedBundles.push_back(pNode); } } } } UpdateGui(); } void QmitkFiberfoxView::EnableCrosshairNavigation() { MITK_DEBUG << "EnableCrosshairNavigation"; // enable the crosshair navigation if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "enabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(true); // linkedRenderWindow->EnableSlicingPlanes(true); } if (m_Controls->m_RealTimeFibers->isChecked()) GenerateFibers(); } void QmitkFiberfoxView::DisableCrosshairNavigation() { MITK_DEBUG << "DisableCrosshairNavigation"; // disable the crosshair navigation during the drawing if (mitk::ILinkedRenderWindowPart* linkedRenderWindow = dynamic_cast(this->GetRenderWindowPart())) { MITK_DEBUG << "disabling linked navigation"; linkedRenderWindow->EnableLinkedNavigation(false); // linkedRenderWindow->EnableSlicingPlanes(false); } } void QmitkFiberfoxView::NodeRemoved(const mitk::DataNode* node) { if (node == m_SelectedImage) m_SelectedImage = NULL; if (node == m_SelectedBundle) m_SelectedBundle = NULL; mitk::DataNode* nonConstNode = const_cast(node); std::map::iterator it = m_DataNodeToPlanarFigureData.find(nonConstNode); if( it != m_DataNodeToPlanarFigureData.end() ) { QmitkPlanarFigureData& data = it->second; // remove observers data.m_Figure->RemoveObserver( data.m_EndPlacementObserverTag ); data.m_Figure->RemoveObserver( data.m_SelectObserverTag ); data.m_Figure->RemoveObserver( data.m_StartInteractionObserverTag ); data.m_Figure->RemoveObserver( data.m_EndInteractionObserverTag ); m_DataNodeToPlanarFigureData.erase( it ); } } void QmitkFiberfoxView::NodeAdded( const mitk::DataNode* node ) { // add observer for selection in renderwindow mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); bool isPositionMarker (false); node->GetBoolProperty("isContourMarker", isPositionMarker); if( figure && !isPositionMarker ) { MITK_DEBUG << "figure added. will add interactor if needed."; mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetInteractor()); mitk::DataNode* nonConstNode = const_cast( node ); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", nonConstNode); } else { // just to be sure that the interactor is not added twice mitk::GlobalInteraction::GetInstance()->RemoveInteractor(figureInteractor); } MITK_DEBUG << "adding interactor to globalinteraction"; mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); MITK_DEBUG << "will now add observers for planarfigure"; QmitkPlanarFigureData data; data.m_Figure = figure; // // add observer for event when figure has been placed typedef itk::SimpleMemberCommand< QmitkFiberfoxView > SimpleCommandType; // SimpleCommandType::Pointer initializationCommand = SimpleCommandType::New(); // initializationCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureInitialized ); // data.m_EndPlacementObserverTag = figure->AddObserver( mitk::EndPlacementPlanarFigureEvent(), initializationCommand ); // add observer for event when figure is picked (selected) typedef itk::MemberCommand< QmitkFiberfoxView > MemberCommandType; MemberCommandType::Pointer selectCommand = MemberCommandType::New(); selectCommand->SetCallbackFunction( this, &QmitkFiberfoxView::PlanarFigureSelected ); data.m_SelectObserverTag = figure->AddObserver( mitk::SelectPlanarFigureEvent(), selectCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer startInteractionCommand = SimpleCommandType::New(); startInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::DisableCrosshairNavigation); data.m_StartInteractionObserverTag = figure->AddObserver( mitk::StartInteractionPlanarFigureEvent(), startInteractionCommand ); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberfoxView::EnableCrosshairNavigation); data.m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); m_DataNodeToPlanarFigureData[nonConstNode] = data; } } void QmitkFiberfoxView::PlanarFigureSelected( itk::Object* object, const itk::EventObject& ) { mitk::TNodePredicateDataType::Pointer isPf = mitk::TNodePredicateDataType::New(); mitk::DataStorage::SetOfObjects::ConstPointer allPfs = this->GetDataStorage()->GetSubset( isPf ); for ( mitk::DataStorage::SetOfObjects::const_iterator it = allPfs->begin(); it!=allPfs->end(); ++it) { mitk::DataNode* node = *it; if( node->GetData() == object ) { node->SetSelected(true); m_SelectedFiducial = node; } else node->SetSelected(false); } UpdateGui(); this->RequestRenderWindowUpdate(); } void QmitkFiberfoxView::SetFocus() { m_Controls->m_CircleButton->setFocus(); }