diff --git a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp index 713ecbf..04e9d90 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp +++ b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp @@ -1,2761 +1,2792 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBundle.h" #include #include #include #include "mitkImagePixelReadAccessor.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs"; mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData ) : m_NumFibers(0) { m_TrackVisHeader.hdr_size = 0; m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != nullptr) m_FiberPolyData = fiberPolyData; else { this->m_FiberPolyData->SetPoints(vtkSmartPointer::New()); this->m_FiberPolyData->SetLines(vtkSmartPointer::New()); } this->UpdateFiberGeometry(); this->GenerateFiberIds(); this->ColorFibersByOrientation(); } mitk::FiberBundle::~FiberBundle() { } mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy() { mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData); newFib->SetFiberColors(this->m_FiberColors); newFib->SetFiberWeights(this->m_FiberWeights); newFib->SetTrackVisHeader(this->GetTrackVisHeader()); return newFib; } vtkSmartPointer mitk::FiberBundle::GeneratePolyDataByIds(std::vector fiberIds, vtkSmartPointer weights) { vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); vtkSmartPointer newLineSet = vtkSmartPointer::New(); vtkSmartPointer newPointSet = vtkSmartPointer::New(); weights->SetNumberOfValues(fiberIds.size()); int counter = 0; auto finIt = fiberIds.begin(); while ( finIt != fiberIds.end() ) { if (*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++) { newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); } weights->InsertValue(counter, this->GetFiberWeight(*finIt)); newLineSet->InsertNextCell(newFiber); ++finIt; ++counter; } newFiberPolyData->SetPoints(newPointSet); newFiberPolyData->SetLines(newLineSet); return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundles(std::vector< mitk::FiberBundle::Pointer > fibs) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); auto num_weights = this->GetNumFibers(); for (auto fib : fibs) num_weights += fib->GetNumFibers(); weights->SetNumberOfValues(num_weights); unsigned int counter = 0; for (unsigned int i=0; iGetNumberOfCells(); ++i) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (unsigned int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } for (auto fib : fibs) { // add new fiber bundle for (unsigned int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (unsigned int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, fib->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib) { if (fib==nullptr) return this->GetDeepCopy(); MITK_INFO << "Adding fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers()); unsigned int counter = 0; for (unsigned int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (unsigned int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // add new fiber bundle for (unsigned int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (unsigned int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, fib->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // Only retain fibers with a weight larger than the specified threshold mitk::FiberBundle::Pointer mitk::FiberBundle::FilterByWeights(float weight_thr, bool invert) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector weights; for (unsigned int i=0; iGetNumFibers(); i++) { if ( (invert && this->GetFiberWeight(i)>weight_thr) || (!invert && this->GetFiberWeight(i)<=weight_thr)) continue; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); weights.push_back(this->GetFiberWeight(i)); } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); for (unsigned int i=0; iSetFiberWeight(i, weights.at(i)); newFib->SetTrackVisHeader(this->GetTrackVisHeader()); return newFib; } // Only retain a subsample of the fibers mitk::FiberBundle::Pointer mitk::FiberBundle::SubsampleFibers(float factor, bool random_seed) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); unsigned int new_num_fibs = static_cast(std::round(this->GetNumFibers()*factor)); MITK_INFO << "Subsampling fibers with factor " << factor << "(" << new_num_fibs << "/" << this->GetNumFibers() << ")"; // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(new_num_fibs); std::vector< unsigned int > ids; for (unsigned int i=0; iGetNumFibers(); i++) ids.push_back(i); if (random_seed) std::srand(static_cast(std::time(nullptr))); else std::srand(0); std::random_shuffle(ids.begin(), ids.end()); unsigned int counter = 0; for (unsigned int i=0; iGetCell(ids.at(i)); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(ids.at(i))); vNewLines->InsertNextCell(container); counter++; } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); newFib->SetTrackVisHeader(this->GetTrackVisHeader()); return newFib; } // subtract two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib) { if (fib==nullptr) return this->GetDeepCopy(); MITK_INFO << "Subtracting fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector< std::vector< itk::Point > > points1; for(unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start = mitk::imv::GetItkPoint(points->GetPoint(0)); itk::Point end = mitk::imv::GetItkPoint(points->GetPoint(numPoints-1)); points1.push_back( {start, end} ); } std::vector< std::vector< itk::Point > > points2; for(unsigned int i=0; iGetNumFibers(); i++ ) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start =mitk::imv::GetItkPoint(points->GetPoint(0)); itk::Point end =mitk::imv::GetItkPoint(points->GetPoint(numPoints-1)); points2.push_back( {start, end} ); } // int progress = 0; std::vector< int > ids; #pragma omp parallel for for (int i=0; i(points1.size()); i++) { bool match = false; for (unsigned int j=0; j(i)); auto v2 = points2.at(j); float dist=0; for (unsigned int c=0; cGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(points->GetPoint(j)); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if(vNewLines->GetNumberOfCells()==0) return mitk::FiberBundle::New(); // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle return mitk::FiberBundle::New(vNewPolyData); } /* * set PolyData (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { if (fiberPD == nullptr) this->m_FiberPolyData = vtkSmartPointer::New(); else m_FiberPolyData->DeepCopy(fiberPD); m_NumFibers = static_cast(m_FiberPolyData->GetNumberOfLines()); if (updateGeometry) UpdateFiberGeometry(); GenerateFiberIds(); ColorFibersByOrientation(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundle::GetFiberPolyData() const { return m_FiberPolyData; } -void mitk::FiberBundle::ColorFibersByLength(bool opacity, bool normalize) +void mitk::FiberBundle::ColorFibersByLength(bool opacity, bool normalize, bool weight_fibers) { if (m_MaxFiberLength<=0) return; auto numOfPoints = this->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(numOfPoints * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); auto numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) return; mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); unsigned int count = 0; for (unsigned int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); float l = m_FiberLengths.at(i)/m_MaxFiberLength; if (!normalize) { l = m_FiberLengths.at(i)/255.0f; if (l > 1.0f) l = 1.0; } for (int j=0; jGetColor(1.0 - static_cast(l), color); rgba[0] = static_cast(255.0 * color[0]); rgba[1] = static_cast(255.0 * color[1]); rgba[2] = static_cast(255.0 * color[2]); if (opacity) rgba[3] = static_cast(255.0f * l); else rgba[3] = static_cast(255.0); m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba); count++; } + + if (weight_fibers) + this->SetFiberWeight(i, m_FiberLengths.at(i)); } + m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByOrientation() { //===== 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 //================================================= vtkPoints* extrPoints = m_FiberPolyData->GetPoints(); vtkIdType numOfPoints = 0; if (extrPoints!=nullptr) numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(numOfPoints * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); auto numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) return; /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); /* 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] = static_cast(255.0 * std::fabs(diff[0])); rgba[1] = static_cast(255.0 * std::fabs(diff[1])); rgba[2] = static_cast(255.0 * std::fabs(diff[2])); rgba[3] = static_cast(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] = static_cast(255.0 * std::fabs(diff1[0])); rgba[1] = static_cast(255.0 * std::fabs(diff1[1])); rgba[2] = static_cast(255.0 * std::fabs(diff1[2])); rgba[3] = static_cast(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] = static_cast(255.0 * std::fabs(diff2[0])); rgba[1] = static_cast(255.0 * std::fabs(diff2[1])); rgba[2] = static_cast(255.0 * std::fabs(diff2[2])); rgba[3] = static_cast(255.0); } m_FiberColors->InsertTypedTuple(idList[i], rgba); } } else if (pointsPerFiber == 1) { /* a single point does not define a fiber (use vertex mechanisms instead */ continue; } else { MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } -void mitk::FiberBundle::ColorFibersByCurvature(bool, bool normalize) +void mitk::FiberBundle::ColorFibersByCurvature(bool, bool normalize, bool weight_fibers) { double window = 5; //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); std::vector< double > values; double min = 1; double max = 0; MITK_INFO << "Coloring fibers by curvature"; boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); + for (int i=0; iGetNumberOfCells(); i++) { ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); + double mean_curv = 0; // calculate curvatures for (int j=0; j > vectors; vnl_vector_fixed< double, 3 > meanV; meanV.fill(0.0); while(dist1) { double p1[3]; points->GetPoint(c-1, p1); double p2[3]; points->GetPoint(c, p2); vnl_vector_fixed< double, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); meanV += v; c--; } c = j; dist = 0; while(distGetPoint(c, p1); double p2[3]; points->GetPoint(c+1, p2); vnl_vector_fixed< double, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); meanV += v; c++; } meanV.normalize(); double dev = 0; for (unsigned int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/itk::Math::pi; } if (vectors.size()>0) dev /= vectors.size(); + if (weight_fibers) + mean_curv += dev; dev = 1.0-dev/180.0; values.push_back(dev); if (devmax) max = dev; } + + if (weight_fibers) + this->SetFiberWeight(i, mean_curv/numPoints); } + unsigned int count = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); for (int j=0; j1) dev = 1; lookupTable->GetColor(dev, color); rgba[0] = static_cast(255.0 * color[0]); rgba[1] = static_cast(255.0 * color[1]); rgba[2] = static_cast(255.0 * color[2]); rgba[3] = static_cast(255.0); m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba); count++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray) { for(long i=0; iGetNumberOfTuples(); i++) { double faValue = FAValArray->GetValue(i); faValue = faValue * 255.0; m_FiberColors->SetComponent(i,3, static_cast(faValue) ); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ResetFiberOpacity() { for(long i=0; iGetNumberOfTuples(); i++) m_FiberColors->SetComponent(i,3, 255.0 ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } -void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity, bool normalize) +void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity, bool normalize, bool weight_fibers) { - mitkPixelTypeMultiplex3( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity, normalize ); + mitkPixelTypeMultiplex4( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity, normalize, weight_fibers ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } template -void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize) +void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize, bool weight_fibers) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); unsigned char rgba[4] = {0,0,0,0}; vtkPoints* pointSet = m_FiberPolyData->GetPoints(); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); double min = 999999; double max = -999999; - for(long i=0; iGetNumberOfPoints(); ++i) + + for (unsigned int i=0; iGetPoint(i)[0]; - px[1] = pointSet->GetPoint(i)[1]; - px[2] = pointSet->GetPoint(i)[2]; - auto pixelValue = static_cast(readimage.GetPixelByWorldCoordinates(px)); - if (pixelValue>max) - max = pixelValue; - if (pixelValueGetCell(i); + auto numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + double mean_val = 0; + + for (int j=0; jGetPoint(j, p); + + Point3D px; + px[0] = p[0]; + px[1] = p[1]; + px[2] = p[2]; + auto pixelValue = static_cast(readimage.GetPixelByWorldCoordinates(px)); + + if (pixelValue>max) + max = pixelValue; + if (pixelValueSetFiberWeight(i, mean_val/numPoints); } 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]; auto pixelValue = static_cast(readimage.GetPixelByWorldCoordinates(px)); if (normalize) pixelValue = (pixelValue-min)/(max-min); else if (pixelValue>1) pixelValue = 1; double color[3]; lookupTable->GetColor(1-pixelValue, color); rgba[0] = static_cast(255.0 * color[0]); rgba[1] = static_cast(255.0 * color[1]); rgba[2] = static_cast(255.0 * color[2]); if (opacity) rgba[3] = static_cast(255.0 * pixelValue); else rgba[3] = static_cast(255.0); m_FiberColors->InsertTypedTuple(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByFiberWeights(bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); unsigned char rgba[4] = {0,0,0,0}; unsigned int counter = 0; float max = -999999; float min = 999999; for (unsigned int i=0; iGetFiberWeight(i); if (weight>max) max = weight; if (weightGetCell(i); auto numPoints = cell->GetNumberOfPoints(); auto weight = this->GetFiberWeight(i); for (int j=0; j1) v = 1; double color[3]; lookupTable->GetColor(static_cast(1-v), color); rgba[0] = static_cast(255.0 * color[0]); rgba[1] = static_cast(255.0 * color[1]); rgba[2] = static_cast(255.0 * color[2]); if (opacity) rgba[3] = static_cast(255.0f * v); else rgba[3] = static_cast(255.0); m_FiberColors->InsertTypedTuple(counter, rgba); counter++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); unsigned char rgba[4] = {0,0,0,0}; for(long i=0; iGetNumberOfPoints(); ++i) { rgba[0] = static_cast(r); rgba[1] = static_cast(g); rgba[2] = static_cast(b); rgba[3] = static_cast(alpha); m_FiberColors->InsertTypedTuple(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::GenerateFiberIds() { if (m_FiberPolyData == nullptr) return; vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); idFiberFilter->SetInputData(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(); } float mitk::FiberBundle::GetNumEpFractionInMask(ItkUcharImgType* mask, bool different_label) { vtkSmartPointer PolyData = m_FiberPolyData; MITK_INFO << "Calculating EP-Fraction"; boost::progress_display disp(m_NumFibers); unsigned int in_mask = 0; for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); itk::Point startVertex =mitk::imv::GetItkPoint(points->GetPoint(0)); itk::Index<3> startIndex; mask->TransformPhysicalPointToIndex(startVertex, startIndex); itk::Point endVertex =mitk::imv::GetItkPoint(points->GetPoint(numPoints-1)); itk::Index<3> endIndex; mask->TransformPhysicalPointToIndex(endVertex, endIndex); if (mask->GetLargestPossibleRegion().IsInside(startIndex) && mask->GetLargestPossibleRegion().IsInside(endIndex)) { float v1 = mask->GetPixel(startIndex); if (v1 < 0.5f) continue; float v2 = mask->GetPixel(startIndex); if (v2 < 0.5f) continue; if (!different_label) ++in_mask; else if (fabs(v1-v2)>0.00001f) ++in_mask; } } return float(in_mask)/m_NumFibers; } std::tuple mitk::FiberBundle::GetDirectionalOverlap(ItkUcharImgType* mask, mitk::PeakImage::ItkPeakImageType* peak_image) { vtkSmartPointer PolyData = m_FiberPolyData; MITK_INFO << "Calculating overlap"; auto spacing = mask->GetSpacing(); boost::progress_display disp(m_NumFibers); double length_sum = 0; double in_mask_length = 0; double aligned_length = 0; for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; j startVertex =mitk::imv::GetItkPoint(points->GetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; mask->TransformPhysicalPointToIndex(startVertex, startIndex); mask->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex =mitk::imv::GetItkPoint(points->GetPoint(j + 1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; mask->TransformPhysicalPointToIndex(endVertex, endIndex); mask->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); vnl_vector_fixed< float, 3 > fdir; fdir[0] = endVertex[0] - startVertex[0]; fdir[1] = endVertex[1] - startVertex[1]; fdir[2] = endVertex[2] - startVertex[2]; fdir.normalize(); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > segment : segments) { if ( mask->GetLargestPossibleRegion().IsInside(segment.first) && mask->GetPixel(segment.first) > 0 ) { in_mask_length += segment.second; mitk::PeakImage::ItkPeakImageType::IndexType idx4; idx4[0] = segment.first[0]; idx4[1] = segment.first[1]; idx4[2] = segment.first[2]; vnl_vector_fixed< float, 3 > peak; idx4[3] = 0; peak[0] = peak_image->GetPixel(idx4); idx4[3] = 1; peak[1] = peak_image->GetPixel(idx4); idx4[3] = 2; peak[2] = peak_image->GetPixel(idx4); if (std::isnan(peak[0]) || std::isnan(peak[1]) || std::isnan(peak[2]) || peak.magnitude()<0.0001f) continue; peak.normalize(); double f = 1.0 - std::acos(std::fabs(static_cast(dot_product(fdir, peak)))) * 2.0/itk::Math::pi; aligned_length += segment.second * f; } length_sum += segment.second; } } } if (length_sum<=0.0001) { MITK_INFO << "Fiber length sum is zero!"; return std::make_tuple(0,0); } return std::make_tuple(aligned_length/length_sum, in_mask_length/length_sum); } float mitk::FiberBundle::GetOverlap(ItkUcharImgType* mask) { vtkSmartPointer PolyData = m_FiberPolyData; MITK_INFO << "Calculating overlap"; auto spacing = mask->GetSpacing(); boost::progress_display disp(m_NumFibers); double length_sum = 0; double in_mask_length = 0; for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; j startVertex =mitk::imv::GetItkPoint(points->GetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; mask->TransformPhysicalPointToIndex(startVertex, startIndex); mask->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex =mitk::imv::GetItkPoint(points->GetPoint(j + 1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; mask->TransformPhysicalPointToIndex(endVertex, endIndex); mask->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > segment : segments) { if ( mask->GetLargestPossibleRegion().IsInside(segment.first) && mask->GetPixel(segment.first) > 0 ) in_mask_length += segment.second; length_sum += segment.second; } } } if (length_sum<=0.000001) { MITK_INFO << "Fiber length sum is zero!"; return 0; } return static_cast(in_mask_length/length_sum); } mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); std::vector< float > fib_weights; MITK_INFO << "Cutting fibers"; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); int newNumPoints = 0; if (numPoints>1) { for (int j=0; j itkP =mitk::imv::GetItkPoint(points->GetPoint(j)); itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); bool inside = false; if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx)!=0 ) inside = true; if (inside && !invert) { vtkIdType id = vtkNewPoints->InsertNextPoint(itkP.GetDataPointer()); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if ( !inside && invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(itkP.GetDataPointer()); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if (newNumPoints>1) { fib_weights.push_back(this->GetFiberWeight(i)); vtkNewCells->InsertNextCell(container); newNumPoints = 0; container = vtkSmartPointer::New(); } else { newNumPoints = 0; container = vtkSmartPointer::New(); } } if (newNumPoints>1) { fib_weights.push_back(this->GetFiberWeight(i)); vtkNewCells->InsertNextCell(container); } } } vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(static_cast(fib_weights.size())); if (vtkNewCells->GetNumberOfCells()<=0) return nullptr; for (unsigned int i=0; iGetNumberOfValues(); i++) newFiberWeights->SetValue(i, fib_weights.at(i)); // vtkSmartPointer newFiberColors = vtkSmartPointer::New(); // newFiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); // newFiberColors->SetNumberOfComponents(4); // newFiberColors->SetName("FIBER_COLORS"); // unsigned char rgba[4] = {0,0,0,0}; // for(long i=0; iGetNumberOfPoints(); ++i) // { // rgba[0] = (unsigned char) r; // rgba[1] = (unsigned char) g; // rgba[2] = (unsigned char) b; // rgba[3] = (unsigned char) alpha; // m_FiberColors->InsertTypedTuple(i, rgba); // } vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData); newFib->SetFiberWeights(newFiberWeights); // newFib->Compress(0.1); newFib->SetTrackVisHeader(this->GetTrackVisHeader()); return newFib; } mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage) { if (roi==nullptr || !(dynamic_cast(roi->GetData()) || dynamic_cast(roi->GetData())) ) return nullptr; std::vector tmp = ExtractFiberIdSubset(roi, storage); if (tmp.size()<=0) return mitk::FiberBundle::New(); vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); fib->SetFiberWeights(weights); fib->SetTrackVisHeader(this->GetTrackVisHeader()); return fib; } std::vector mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage) { std::vector result; if (roi==nullptr || roi->GetData()==nullptr) return result; mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast(roi->GetData()); if (!pfc.IsNull()) // handle composite { DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi); if (children->size()==0) return result; switch (pfc->getOperationType()) { case 0: // AND { MITK_INFO << "AND"; result = this->ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { std::vector inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(std::min(result.size(),inRoi.size())); it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( static_cast(it - rest.begin()) ); result = rest; } break; } case 1: // OR { MITK_INFO << "OR"; result = ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { it = result.end(); std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); result.insert(it, inRoi.begin(), inRoi.end()); } // remove duplicates sort(result.begin(), result.end()); it = unique(result.begin(), result.end()); result.resize( static_cast(it - result.begin()) ); break; } case 2: // NOT { MITK_INFO << "NOT"; for(unsigned int i=0; iGetNumFibers(); i++) result.push_back(i); std::vector::iterator it; for (unsigned int i=0; iSize(); ++i) { std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(result.size()-inRoi.size()); it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( static_cast(it - rest.begin()) ); result = rest; } break; } } } else if ( dynamic_cast(roi->GetData()) ) // actual extraction { if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarPoly = dynamic_cast(roi->GetData()); //create vtkPolygon using controlpoints from planarFigure polygon vtkSmartPointer polygonVtk = vtkSmartPointer::New(); for (unsigned int i=0; iGetNumberOfControlPoints(); ++i) { itk::Point p = planarPoly->GetWorldControlPoint(i); vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] ); polygonVtk->GetPointIds()->InsertNextId(id); } MITK_INFO << "Extracting with polygon"; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); double tolerance = 0.001; // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection double pcoords[3] = {0,0,0}; int subId = 0; int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId); if (iD!=0) { result.push_back(i); break; } } } } else if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(roi->GetData()); Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal(); planeNormal.Normalize(); //calculate circle radius mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint double radius = V1w.EuclideanDistanceTo(V2w); radius *= radius; MITK_INFO << "Extracting with circle"; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x); if (iD!=0) { double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]); if( dist <= radius) { result.push_back(i); break; } } } } } return result; } return result; } void mitk::FiberBundle::UpdateFiberGeometry() { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInputData(m_FiberPolyData); cleaner->PointMergingOff(); cleaner->Update(); m_FiberPolyData = cleaner->GetOutput(); m_FiberLengths.clear(); m_MeanFiberLength = 0; m_MedianFiberLength = 0; m_LengthStDev = 0; m_NumFibers = static_cast(m_FiberPolyData->GetNumberOfCells()); if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints()) this->ColorFibersByOrientation(); if (m_FiberWeights->GetNumberOfValues()!=m_NumFibers) { m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberWeights->SetNumberOfValues(m_NumFibers); this->SetFiberWeights(1); } 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(false); float b[] = {0, 1, 0, 1, 0, 1}; geometry->SetFloatBounds(b); SetGeometry(geometry); return; } double b[6]; m_FiberPolyData->GetBounds(b); // calculate statistics for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float length = 0; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); double 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 += static_cast(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 (unsigned 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); mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetFloatBounds(b); this->SetGeometry(geometry); GetTrackVisHeader(); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const { return m_FiberWeights->GetValue(fiber); } void mitk::FiberBundle::SetFiberWeights(float newWeight) { for (int i=0; iGetNumberOfValues(); i++) m_FiberWeights->SetValue(i, newWeight); } void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer weights) { if (m_NumFibers!=weights->GetNumberOfValues()) { MITK_INFO << "Weights array not equal to number of fibers! " << weights->GetNumberOfValues() << " vs " << m_NumFibers; return; } for (int i=0; iGetNumberOfValues(); i++) m_FiberWeights->SetValue(i, weights->GetValue(i)); m_FiberWeights->SetName("FIBER_WEIGHTS"); } void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight) { m_FiberWeights->SetValue(fiber, weight); } void mitk::FiberBundle::SetFiberColors(vtkSmartPointer fiberColors) { for(long i=0; iGetNumberOfPoints(); ++i) { unsigned char source[4] = {0,0,0,0}; fiberColors->GetTypedTuple(i, source); unsigned char target[4] = {0,0,0,0}; target[0] = source[0]; target[1] = source[1]; target[2] = source[2]; target[3] = source[3]; m_FiberColors->InsertTypedTuple(i, target); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz) { rx = rx*itk::Math::pi/180; ry = ry*itk::Math::pi/180; rz = rz*itk::Math::pi/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; m = rot*m; return m; } void mitk::FiberBundle::TransformFibers(itk::ScalableAffineTransform< mitk::ScalarType >::Pointer transform) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; j p =mitk::imv::GetItkPoint(points->GetPoint(j)); p = transform->TransformPoint(p); vtkIdType id = vtkNewPoints->InsertNextPoint(p.GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz) { vnl_matrix_fixed< double, 3, 3 > rot = mitk::imv::GetRotationMatrixVnl(rx, ry, rz); mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(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 = rot*dir; dir[0] += center[0]+tx; dir[1] += center[1]+ty; dir[2] += center[2]+tz; 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); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z) { x = x*itk::Math::pi/180; y = y*itk::Math::pi/180; z = z*itk::Math::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::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(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); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter) { MITK_INFO << "Scaling fibers"; boost::progress_display disp(m_NumFibers); mitk::BaseGeometry* geom = this->GetGeometry(); mitk::Point3D c = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); if (subtractCenter) { p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2]; } p[0] *= x; p[1] *= y; p[2] *= z; if (subtractCenter) { p[0] += c[0]; p[1] += c[1]; p[2] += c[2]; } vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TranslateFibers(double x, double y, double z) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(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); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::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(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(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); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RemoveDir(vnl_vector_fixed dir, double threshold) { dir.normalize(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); bool discard = false; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); vnl_vector_fixed< double, 3 > v1; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; if (v1.magnitude()>0.001) { v1.normalize(); if (fabs(dot_product(v1,dir))>threshold) { discard = true; break; } } } if (!discard) { for (int j=0; jGetPoint(j, p1); vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); // UpdateColorCoding(); // UpdateFiberGeometry(); } bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { if (minRadius<0) return true; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Applying curvature threshold"; boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); double p3[3]; points->GetPoint(j+2, p3); vnl_vector_fixed< float, 3 > v1, v2, v3; v1[0] = static_cast(p2[0]-p1[0]); v1[1] = static_cast(p2[1]-p1[1]); v1[2] = static_cast(p2[2]-p1[2]); v2[0] = static_cast(p3[0]-p2[0]); v2[1] = static_cast(p3[1]-p2[1]); v2[2] = static_cast(p3[2]-p2[2]); v3[0] = static_cast(p1[0]-p3[0]); v3[1] = static_cast(p1[1]-p3[1]); v3[2] = static_cast(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 && rInsertNextCell(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); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM) { MITK_INFO << "Removing short fibers"; if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers { MITK_WARN << "Process aborted. No fibers would be left!"; return false; } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); float min = m_MaxFiberLength; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)>=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(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); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) return true; if (lengthInMM vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Removing long fibers"; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)<=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(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); this->SetFiberPolyData(m_FiberPolyData, true); return true; } void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias ) { if (pointDistance<=0) return; 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 MITK_INFO << "Smoothing fibers"; vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); std::vector< vtkSmartPointer > resampled_streamlines; resampled_streamlines.resize(m_NumFibers); boost::progress_display disp(m_NumFibers); #pragma omp parallel for for (int i=0; i(m_NumFibers); i++) { vtkSmartPointer newPoints = vtkSmartPointer::New(); float length = 0; #pragma omp critical { length = m_FiberLengths.at(static_cast(i)); ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jInsertNextPoint(points->GetPoint(j)); } int sampling = static_cast(std::ceil(length/pointDistance)); 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(newPoints); 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(); #pragma omp critical { for (int j=0; jGetNumberOfPoints(); j++) { vtkIdType id = vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); smoothLine->GetPointIds()->InsertNextId(id); } resampled_streamlines[static_cast(i)] = smoothLine; } } for (auto container : resampled_streamlines) { vtkSmoothCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkSmoothPoints); m_FiberPolyData->SetLines(vtkSmoothCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ResampleSpline(float pointDistance) { ResampleSpline(pointDistance, 0, 0, 0 ); } unsigned int mitk::FiberBundle::GetNumberOfPoints() const { unsigned int points = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); points += cell->GetNumberOfPoints(); } return points; } void mitk::FiberBundle::Compress(float error) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Compressing fibers"; unsigned int numRemovedPoints = 0; boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for for (int i=0; i(m_FiberPolyData->GetNumberOfCells()); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; #pragma omp critical { ++disp; weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } // calculate curvatures auto numPoints = vertices.size(); std::vector< int > removedPoints; removedPoints.resize(numPoints, 0); removedPoints[0]=-1; removedPoints[numPoints-1]=-1; vtkSmartPointer container = vtkSmartPointer::New(); unsigned int remCounter = 0; bool pointFound = true; while (pointFound) { pointFound = false; double minError = static_cast(error); unsigned int removeIndex = 0; for (unsigned int j=0; j candV = vertices.at(j); int validP = -1; vnl_vector_fixed< double, 3 > pred; for (int k=static_cast(j)-1; k>=0; k--) if (removedPoints[static_cast(k)]<=0) { pred = vertices.at(static_cast(k)); validP = k; break; } int validS = -1; vnl_vector_fixed< double, 3 > succ; for (unsigned int k=j+1; k(k); break; } if (validP>=0 && validS>=0) { double a = (candV-pred).magnitude(); double b = (candV-succ).magnitude(); double c = (pred-succ).magnitude(); double s=0.5*(a+b+c); double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c))); if (hcInsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); numRemovedPoints += remCounter; vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()>0) { MITK_INFO << "Removed points: " << numRemovedPoints; SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } void mitk::FiberBundle::ResampleToNumPoints(unsigned int targetPoints) { if (targetPoints<2) mitkThrow() << "Minimum two points required for resampling!"; MITK_INFO << "Resampling fibers (number of points " << targetPoints << ")"; bool unequal_fibs = true; while (unequal_fibs) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); unequal_fibs = false; for (unsigned int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; double seg_len = 0; { weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); if (numPoints!=targetPoints) seg_len = static_cast(this->GetFiberLength(i)/(targetPoints-1)); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } vtkSmartPointer container = vtkSmartPointer::New(); vnl_vector_fixed< double, 3 > lastV = vertices.at(0); { vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); container->GetPointIds()->InsertNextId(id); } for (unsigned int j=1; j vec = vertices.at(j) - lastV; double new_dist = vec.magnitude(); if (new_dist >= seg_len && seg_len>0) { vnl_vector_fixed< double, 3 > newV = lastV; if ( new_dist-seg_len <= mitk::eps ) { vec.normalize(); newV += vec * seg_len; } else { // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') vnl_vector_fixed< double, 3 > p = vertices.at(j-1); vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - seg_len*seg_len; double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); if (v1>0) newV = p + d * v1; else if (v2>0) newV = p + d * v2; else MITK_INFO << "ERROR1 - linear resampling"; j--; } //#pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); container->GetPointIds()->InsertNextId(id); } lastV = newV; } else if ( (j==vertices.size()-1 && new_dist>0.0001) || seg_len<=0.0000001) { //#pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } //#pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); vtkNewCells->InsertNextCell(container); if (container->GetNumberOfPoints()!=targetPoints) unequal_fibs = true; } } if (vtkNewCells->GetNumberOfCells()>0) { SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } } void mitk::FiberBundle::ResampleLinear(double pointDistance) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Resampling fibers (linear)"; boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); std::vector< vtkSmartPointer > resampled_streamlines; resampled_streamlines.resize(static_cast(m_FiberPolyData->GetNumberOfCells())); #pragma omp parallel for for (int i=0; i(m_FiberPolyData->GetNumberOfCells()); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; #pragma omp critical { ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } vtkSmartPointer container = vtkSmartPointer::New(); vnl_vector_fixed< double, 3 > lastV = vertices.at(0); #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); container->GetPointIds()->InsertNextId(id); } for (unsigned int j=1; j vec = vertices.at(j) - lastV; double new_dist = vec.magnitude(); if (new_dist >= pointDistance) { vnl_vector_fixed< double, 3 > newV = lastV; if ( new_dist-pointDistance <= mitk::eps ) { vec.normalize(); newV += vec * pointDistance; } else { // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') vnl_vector_fixed< double, 3 > p = vertices.at(j-1); vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance; double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); if (v1>0) newV = p + d * v1; else if (v2>0) newV = p + d * v2; else MITK_INFO << "ERROR1 - linear resampling"; j--; } #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); container->GetPointIds()->InsertNextId(id); } lastV = newV; } else if (j==vertices.size()-1 && new_dist>0.0001) { #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { resampled_streamlines[static_cast(i)] = container; } } for (auto container : resampled_streamlines) { vtkNewCells->InsertNextCell(container); } if (vtkNewCells->GetNumberOfCells()>0) { m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } // reapply selected colorcoding in case PolyData structure has changed bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps) { if (fib==nullptr) { MITK_INFO << "Reference bundle is nullptr!"; return false; } if (m_NumFibers!=fib->GetNumFibers()) { MITK_INFO << "Unequal number of fibers!"; MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers(); return false; } for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i); auto numPoints2 = cell2->GetNumberOfPoints(); vtkPoints* points2 = cell2->GetPoints(); if (numPoints2!=numPoints) { MITK_INFO << "Unequal number of points in fiber " << i << "!"; MITK_INFO << numPoints2 << " vs. " << numPoints; return false; } for (int j=0; jGetPoint(j); double* p2 = points2->GetPoint(j); if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps) { MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!"; MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2]; MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2]; return false; } } } return true; } void mitk::FiberBundle::PrintSelf(std::ostream &os, itk::Indent indent) const { os << indent << "Number of fibers: " << this->GetNumFibers() << std::endl; os << indent << "Min. fiber length: " << this->GetMinFiberLength() << std::endl; os << indent << "Max. fiber length: " << this->GetMaxFiberLength() << std::endl; os << indent << "Mean fiber length: " << this->GetMeanFiberLength() << std::endl; os << indent << "Median fiber length: " << this->GetMedianFiberLength() << std::endl; os << indent << "STDEV fiber length: " << this->GetLengthStDev() << std::endl; os << indent << "Number of points: " << this->GetNumberOfPoints() << std::endl; os << indent << "Extent x: " << this->GetGeometry()->GetExtentInMM(0) << "mm" << std::endl; os << indent << "Extent y: " << this->GetGeometry()->GetExtentInMM(1) << "mm" << std::endl; os << indent << "Extent z: " << this->GetGeometry()->GetExtentInMM(2) << "mm" << std::endl; os << indent << "Diagonal: " << this->GetGeometry()->GetDiagonalLength() << "mm" << std::endl; os << "\nReference geometry:" << std::endl; os << indent << "Size: [" << std::defaultfloat << m_TrackVisHeader.dim[0] << " " << m_TrackVisHeader.dim[1] << " " << m_TrackVisHeader.dim[2] << "]" << std::endl; os << indent << "Voxel size: [" << m_TrackVisHeader.voxel_size[0] << " " << m_TrackVisHeader.voxel_size[1] << " " << m_TrackVisHeader.voxel_size[2] << "]" << std::endl; os << indent << "Origin: [" << m_TrackVisHeader.origin[0] << " " << m_TrackVisHeader.origin[1] << " " << m_TrackVisHeader.origin[2] << "]" << std::endl; os << indent << "Matrix: " << std::scientific << std::endl; os << indent << "[[" << m_TrackVisHeader.vox_to_ras[0][0] << ", " << m_TrackVisHeader.vox_to_ras[0][1] << ", " << m_TrackVisHeader.vox_to_ras[0][2] << ", " << m_TrackVisHeader.vox_to_ras[0][3] << "]" << std::endl; os << indent << " [" << m_TrackVisHeader.vox_to_ras[1][0] << ", " << m_TrackVisHeader.vox_to_ras[1][1] << ", " << m_TrackVisHeader.vox_to_ras[1][2] << ", " << m_TrackVisHeader.vox_to_ras[1][3] << "]" << std::endl; os << indent << " [" << m_TrackVisHeader.vox_to_ras[2][0] << ", " << m_TrackVisHeader.vox_to_ras[2][1] << ", " << m_TrackVisHeader.vox_to_ras[2][2] << ", " << m_TrackVisHeader.vox_to_ras[2][3] << "]" << std::endl; os << indent << " [" << m_TrackVisHeader.vox_to_ras[3][0] << ", " << m_TrackVisHeader.vox_to_ras[3][1] << ", " << m_TrackVisHeader.vox_to_ras[3][2] << ", " << m_TrackVisHeader.vox_to_ras[3][3] << "]]" << std::defaultfloat << std::endl; if (m_FiberWeights!=nullptr) { std::vector< float > weights; for (int i=0; iGetSize(); i++) weights.push_back(m_FiberWeights->GetValue(i)); std::sort(weights.begin(), weights.end()); os << "\nFiber weight statistics" << std::endl; os << indent << "Min: " << weights.front() << std::endl; os << indent << "1% quantile: " << weights.at(static_cast(weights.size()*0.01)) << std::endl; os << indent << "5% quantile: " << weights.at(static_cast(weights.size()*0.05)) << std::endl; os << indent << "25% quantile: " << weights.at(static_cast(weights.size()*0.25)) << std::endl; os << indent << "Median: " << weights.at(static_cast(weights.size()*0.5)) << std::endl; os << indent << "75% quantile: " << weights.at(static_cast(weights.size()*0.75)) << std::endl; os << indent << "95% quantile: " << weights.at(static_cast(weights.size()*0.95)) << std::endl; os << indent << "99% quantile: " << weights.at(static_cast(weights.size()*0.99)) << std::endl; os << indent << "Max: " << weights.back() << std::endl; } else os << indent << "\n\nNo fiber weight array found." << std::endl; Superclass::PrintSelf(os, 0); } mitk::FiberBundle::TrackVis_header mitk::FiberBundle::GetTrackVisHeader() { if (m_TrackVisHeader.hdr_size==0) { mitk::Geometry3D::Pointer geom = dynamic_cast(this->GetGeometry()); SetTrackVisHeader(geom); } return m_TrackVisHeader; } void mitk::FiberBundle::SetTrackVisHeader(const mitk::FiberBundle::TrackVis_header &TrackVisHeader) { m_TrackVisHeader = TrackVisHeader; } void mitk::FiberBundle::SetTrackVisHeader(mitk::BaseGeometry* geometry) { vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); matrix->Identity(); if (geometry==nullptr) return; for(int i=0; i<3 ;i++) { m_TrackVisHeader.dim[i] = geometry->GetExtent(i); m_TrackVisHeader.voxel_size[i] = geometry->GetSpacing()[i]; m_TrackVisHeader.origin[i] = geometry->GetOrigin()[i]; matrix = geometry->GetVtkMatrix(); } for (int i=0; i<4; ++i) for (int j=0; j<4; ++j) m_TrackVisHeader.vox_to_ras[i][j] = matrix->GetElement(i, j); m_TrackVisHeader.n_scalars = 0; m_TrackVisHeader.n_properties = 0; sprintf(m_TrackVisHeader.voxel_order,"LPS"); m_TrackVisHeader.image_orientation_patient[0] = 1.0; m_TrackVisHeader.image_orientation_patient[1] = 0.0; m_TrackVisHeader.image_orientation_patient[2] = 0.0; m_TrackVisHeader.image_orientation_patient[3] = 0.0; m_TrackVisHeader.image_orientation_patient[4] = 1.0; m_TrackVisHeader.image_orientation_patient[5] = 0.0; m_TrackVisHeader.pad1[0] = 0; m_TrackVisHeader.pad1[1] = 0; m_TrackVisHeader.pad2[0] = 0; m_TrackVisHeader.pad2[1] = 0; m_TrackVisHeader.invert_x = 0; m_TrackVisHeader.invert_y = 0; m_TrackVisHeader.invert_z = 0; m_TrackVisHeader.swap_xy = 0; m_TrackVisHeader.swap_yz = 0; m_TrackVisHeader.swap_zx = 0; m_TrackVisHeader.n_count = 0; m_TrackVisHeader.version = 2; m_TrackVisHeader.hdr_size = 1000; std::string id = "TRACK"; strcpy(m_TrackVisHeader.id_string, id.c_str()); } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundle::UpdateOutputInformation() { } void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundle::VerifyRequestedRegion() { return true; } void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* ) { } diff --git a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.h b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.h index 131b02e..70defee 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.h +++ b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.h @@ -1,244 +1,244 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. 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_FiberBundle_H #define _MITK_FiberBundle_H //includes for MITK datastructure #include #include #include #include #include #include #include #include //includes storing fiberdata #include #include #include #include #include #include #include #include namespace mitk { /** * \brief Base Class for Fiber Bundles; */ class MITKDIFFUSIONCORE_EXPORT FiberBundle : public BaseData { public: typedef itk::Image ItkUcharImgType; // fiber colorcodings static const char* FIBER_ID_ARRAY; void UpdateOutputInformation() override; void SetRequestedRegionToLargestPossibleRegion() override; bool RequestedRegionIsOutsideOfTheBufferedRegion() override; bool VerifyRequestedRegion() override; void SetRequestedRegion(const itk::DataObject*) override; mitkClassMacro( FiberBundle, BaseData ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(Self, vtkSmartPointer) // custom constructor // colorcoding related methods void ColorFibersByFiberWeights(bool opacity, bool normalize); - void ColorFibersByCurvature(bool opacity, bool normalize); - void ColorFibersByLength(bool opacity, bool normalize); - void ColorFibersByScalarMap(mitk::Image::Pointer, bool opacity, bool normalize); + void ColorFibersByCurvature(bool opacity, bool normalize, bool weight_fibers); + void ColorFibersByLength(bool opacity, bool normalize, bool weight_fibers); + void ColorFibersByScalarMap(mitk::Image::Pointer, bool opacity, bool normalize, bool weight_fibers); template - void ColorFibersByScalarMap(const mitk::PixelType pixelType, mitk::Image::Pointer, bool opacity, bool normalize); + void ColorFibersByScalarMap(const mitk::PixelType pixelType, mitk::Image::Pointer, bool opacity, bool normalize, bool weight_fibers); void ColorFibersByOrientation(); void SetFiberOpacity(vtkDoubleArray *FAValArray); void ResetFiberOpacity(); void SetFiberColors(vtkSmartPointer fiberColors); void SetFiberColors(float r, float g, float b, float alpha=255); vtkSmartPointer GetFiberColors() const { return m_FiberColors; } // fiber compression void Compress(float error = 0.0); // fiber resampling void ResampleSpline(float pointDistance=1); void ResampleSpline(float pointDistance, double tension, double continuity, double bias ); void ResampleLinear(double pointDistance=1); void ResampleToNumPoints(unsigned int targetPoints); mitk::FiberBundle::Pointer FilterByWeights(float weight_thr, bool invert=false); 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, bool subtractCenter=true); void TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz); void TransformFibers(itk::ScalableAffineTransform< mitk::ScalarType >::Pointer transform); void RemoveDir(vnl_vector_fixed dir, double threshold); template< class TType=float > void TransformPoint(itk::Point& point, itk::Matrix< TType, 3, 3>& rot, TType& tx, TType& ty, TType& tz) { mitk::Point3D center = this->GetGeometry()->GetCenter(); point[0] -= center[0]; point[1] -= center[1]; point[2] -= center[2]; point = rot*point; point[0] += center[0]+tx; point[1] += center[1]+ty; point[2] += center[2]+tz; } template< class TType=float > void TransformPoint(itk::Point& point, TType rx, TType ry, TType rz, TType tx, TType ty, TType tz) { auto rot = mitk::imv::GetRotationMatrixItk(rx, ry, rz); mitk::Point3D center = this->GetGeometry()->GetCenter(); point[0] -= center[0]; point[1] -= center[1]; point[2] -= center[2]; point = rot*point; point[0] += center[0]+tx; point[1] += center[1]+ty; point[2] += center[2]+tz; } itk::Matrix< double, 3, 3 > TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz); // add/subtract fibers FiberBundle::Pointer AddBundle(FiberBundle* fib); mitk::FiberBundle::Pointer AddBundles(std::vector< mitk::FiberBundle::Pointer > fibs); FiberBundle::Pointer SubtractBundle(FiberBundle* fib); // fiber subset extraction FiberBundle::Pointer ExtractFiberSubset(DataNode *roi, DataStorage* storage); std::vector ExtractFiberIdSubset(DataNode* roi, DataStorage* storage); FiberBundle::Pointer RemoveFibersOutside(ItkUcharImgType* mask, bool invert=false); float GetOverlap(ItkUcharImgType* mask); std::tuple GetDirectionalOverlap(ItkUcharImgType* mask, mitk::PeakImage::ItkPeakImageType* peak_image); float GetNumEpFractionInMask(ItkUcharImgType* mask, bool different_label); mitk::FiberBundle::Pointer SubsampleFibers(float factor, bool random_seed); // get/set data float GetFiberLength(unsigned int index) const { return m_FiberLengths.at(index); } vtkSmartPointer GetFiberWeights() const { return m_FiberWeights; } float GetFiberWeight(unsigned int fiber) const; void SetFiberWeights(float newWeight); void SetFiberWeight(unsigned int fiber, float weight); void SetFiberWeights(vtkSmartPointer weights); void SetFiberPolyData(vtkSmartPointer, bool updateGeometry = true); vtkSmartPointer GetFiberPolyData() const; itkGetConstMacro( NumFibers, unsigned int) //itkGetMacro( FiberSampling, int) itkGetConstMacro( MinFiberLength, float ) itkGetConstMacro( MaxFiberLength, float ) itkGetConstMacro( MeanFiberLength, float ) itkGetConstMacro( MedianFiberLength, float ) itkGetConstMacro( LengthStDev, float ) itkGetConstMacro( UpdateTime2D, itk::TimeStamp ) itkGetConstMacro( UpdateTime3D, itk::TimeStamp ) void RequestUpdate2D(){ m_UpdateTime2D.Modified(); } void RequestUpdate3D(){ m_UpdateTime3D.Modified(); } void RequestUpdate(){ m_UpdateTime2D.Modified(); m_UpdateTime3D.Modified(); } unsigned int GetNumberOfPoints() const; // copy fiber bundle mitk::FiberBundle::Pointer GetDeepCopy(); // compare fiber bundles bool Equals(FiberBundle* fib, double eps=0.01); vtkSmartPointer GeneratePolyDataByIds(std::vector fiberIds, vtkSmartPointer weights); // Structure to hold metadata of a TrackVis file struct TrackVis_header { char id_string[6]; short int dim[3]; float voxel_size[3]; float origin[3]; short int n_scalars; char scalar_name[10][20]; short int n_properties; char property_name[10][20]; float vox_to_ras[4][4]; char reserved[444]; char voxel_order[4]; char pad2[4]; float image_orientation_patient[6]; char pad1[2]; unsigned char invert_x; unsigned char invert_y; unsigned char invert_z; unsigned char swap_xy; unsigned char swap_yz; unsigned char swap_zx; int n_count; int version; int hdr_size; }; TrackVis_header GetTrackVisHeader(); void SetTrackVisHeader(const TrackVis_header &TrackVisHeader); void SetTrackVisHeader(BaseGeometry *geometry); protected: FiberBundle( vtkPolyData* fiberPolyData = nullptr ); ~FiberBundle() override; void GenerateFiberIds(); void UpdateFiberGeometry(); void PrintSelf(std::ostream &os, itk::Indent indent) const override; private: // actual fiber container vtkSmartPointer m_FiberPolyData; // contains fiber ids vtkSmartPointer m_FiberIdDataSet; unsigned int m_NumFibers; vtkSmartPointer m_FiberColors; vtkSmartPointer m_FiberWeights; std::vector< float > m_FiberLengths; float m_MinFiberLength; float m_MaxFiberLength; float m_MeanFiberLength; float m_MedianFiberLength; float m_LengthStDev; itk::TimeStamp m_UpdateTime2D; itk::TimeStamp m_UpdateTime3D; TrackVis_header m_TrackVisHeader; }; } // namespace mitk #endif /* _MITK_FiberBundle_H */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp index e2be804..371d68f 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp @@ -1,1728 +1,1732 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkFiberProcessingView.h" #include #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include "mitkNodePredicateDataType.h" #include #include #include #include #include #include #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() : QmitkAbstractView() , m_Controls( 0 ) , m_CircleCounter(0) , m_PolygonCounter(0) , m_UpsamplingFactor(1) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { RemoveObservers(); } 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 ); 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_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); - connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); + connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubtractBundles()) ); connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) ); connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract())); connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove())); connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify())); connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged())); m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isOdf = mitk::NodePredicateDataType::New("OdfImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isOdf); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_ColorMapBox->SetPredicate(finalPredicate); } UpdateGui(); OnMaskExtractionChanged(); } void QmitkFiberProcessingView::OnMaskExtractionChanged() { m_Controls->m_FiberExtractionFractionLabel->setVisible(false); m_Controls->m_FiberExtractionFractionBox->setVisible(false); m_Controls->m_FiberExtractionThresholdLabel->setVisible(false); m_Controls->m_FiberExtractionThresholdBox->setVisible(false); m_Controls->m_InterpolateRoiBox->setVisible(false); m_Controls->m_BothEnds->setVisible(false); m_Controls->m_LabelsBox->setVisible(false); m_Controls->m_LabelsLabel->setVisible(false); if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2 || m_Controls->m_ExtractionBoxMask->currentIndex() == 3) { m_Controls->m_FiberExtractionFractionLabel->setVisible(true); m_Controls->m_FiberExtractionFractionBox->setVisible(true); m_Controls->m_FiberExtractionThresholdLabel->setVisible(true); m_Controls->m_FiberExtractionThresholdBox->setVisible(true); m_Controls->m_InterpolateRoiBox->setVisible(true); } else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 0 || m_Controls->m_ExtractionBoxMask->currentIndex() == 1) { if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3) m_Controls->m_BothEnds->setVisible(true); m_Controls->m_InterpolateRoiBox->setVisible(true); m_Controls->m_FiberExtractionThresholdLabel->setVisible(true); m_Controls->m_FiberExtractionThresholdBox->setVisible(true); } else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 4) { m_Controls->m_BothEnds->setVisible(true); m_Controls->m_LabelsBox->setVisible(true); m_Controls->m_LabelsLabel->setVisible(true); } } void QmitkFiberProcessingView::SetFocus() { m_Controls->toolBoxx->setFocus(); } void QmitkFiberProcessingView::Modify() { switch (m_Controls->m_ModificationMethodBox->currentIndex()) { case 0: { ResampleSelectedBundlesSpline(); break; } case 1: { ResampleSelectedBundlesLinear(); break; } case 2: { CompressSelectedBundles(); break; } case 3: { DoImageColorCoding(); break; } case 4: { MirrorFibers(); break; } case 5: { WeightFibers(); break; } case 6: { DoCurvatureColorCoding(); break; } case 7: { DoWeightColorCoding(); break; } case 8: { DoLengthColorCoding(); break; } } } void QmitkFiberProcessingView::WeightFibers() { float weight = this->m_Controls->m_BundleWeightBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->SetFiberWeights(weight); } } void QmitkFiberProcessingView::Remove() { switch (m_Controls->m_RemovalMethodBox->currentIndex()) { case 0: { RemoveDir(); break; } case 1: { PruneBundle(); break; } case 2: { ApplyCurvatureThreshold(); break; } case 3: { RemoveWithMask(false); break; } case 4: { RemoveWithMask(true); break; } case 5: { ApplyWeightThreshold(); break; } case 6: { ApplyDensityThreshold(); break; } } } void QmitkFiberProcessingView::Extract() { switch (m_Controls->m_ExtractionMethodBox->currentIndex()) { case 0: { ExtractWithPlanarFigure(); break; } case 1: { switch (m_Controls->m_ExtractionBoxMask->currentIndex()) { { case 0: ExtractWithMask(true, false, false); break; } { case 1: ExtractWithMask(true, true, false); break; } { case 2: ExtractWithMask(false, false, false); break; } { case 3: ExtractWithMask(false, true, false); break; } { case 4: ExtractWithMask(true, false, true); break; } } break; } } } void QmitkFiberProcessingView::PruneBundle() { int minLength = this->m_Controls->m_PruneFibersMinBox->value(); int maxLength = this->m_Controls->m_PruneFibersMaxBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (!fib->RemoveShortFibers(minLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); else if (!fib->RemoveLongFibers(maxLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyWeightThreshold() { float thr = this->m_Controls->m_WeightThresholdBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->FilterByWeights(thr); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByFiberWeights(false, true); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyDensityThreshold() { float thr = this->m_Controls->m_DensityThresholdBox->value(); float ol = this->m_Controls->m_DensityOverlapBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::TractDensityImageFilter< ItkFloatImageType >::Pointer generator = itk::TractDensityImageFilter< ItkFloatImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(false); generator->SetOutputAbsoluteValues(false); generator->Update(); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetRoiImages({generator->GetOutput()}); extractor->SetInputFiberBundle(fib); extractor->SetOverlapFraction(ol); extractor->SetInterpolate(true); extractor->SetThreshold(thr); extractor->SetNoNegatives(true); extractor->Update(); if (extractor->GetPositives().empty()) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::FiberBundle::Pointer newFib = extractor->GetPositives().at(0); if (newFib->GetNumFibers()>0) node->SetData(newFib); else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyCurvatureThreshold() { int angle = this->m_Controls->m_CurvSpinBox->value(); int dist = this->m_Controls->m_CurvDistanceSpinBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(angle); filter->SetDistance(dist); filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked()); filter->Update(); mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle(); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByOrientation(); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveDir() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); vnl_vector_fixed dir; dir[0] = m_Controls->m_ExtractDirX->value(); dir[1] = m_Controls->m_ExtractDirY->value(); dir[2] = m_Controls->m_ExtractDirZ->value(); fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*itk::Math::pi/180)); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveWithMask(bool removeInside) { if (m_RoiImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkMask, mask); mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside); if (newFib->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } node->SetData(newFib); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert, bool labelmap) { if (m_RoiImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { std::string roi_name = m_RoiImageNode->GetName(); mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(mitkMask, mask); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(fib); extractor->SetRoiImages({mask}); extractor->SetRoiImageNames({roi_name}); extractor->SetThreshold(m_Controls->m_FiberExtractionThresholdBox->value()); extractor->SetOverlapFraction(m_Controls->m_FiberExtractionFractionBox->value()); extractor->SetBothEnds(m_Controls->m_BothEnds->isChecked()); extractor->SetInterpolate(m_Controls->m_InterpolateRoiBox->isChecked()); extractor->SetMinFibersPerTract(m_Controls->m_MinExtractedFibersBox->value()); extractor->SetSplitByRoi(true); if (invert) extractor->SetNoPositives(true); else extractor->SetNoNegatives(true); if (labelmap) { std::string labels_string = m_Controls->m_LabelsBox->text().toStdString(); if (labels_string!="ALL") { std::vector strs; boost::split(strs,labels_string,boost::is_any_of(" ,;\t")); std::vector< unsigned short > labels_vector; for (auto v : strs) { try{ unsigned short l = boost::lexical_cast(v); labels_vector.push_back(l); } catch(...) { } } extractor->SetLabels(labels_vector); } extractor->SetInterpolate(false); extractor->SetInputType(itk::FiberExtractionFilter::INPUT::LABEL_MAP); extractor->SetSplitLabels(true); onlyEnds = true; } if (onlyEnds) extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); extractor->Update(); std::vector< mitk::FiberBundle::Pointer > newFibs; if (invert) newFibs = extractor->GetNegatives(); else newFibs = extractor->GetPositives(); if (newFibs.empty()) { QMessageBox::information(nullptr, "No output generated:", "No fibers could be extracted."); continue; } auto labels = extractor->GetPositiveLabels(); for (unsigned int i=0; iSetData(fib); std::string name = roi_name; if (iSetName(name); GetDataStorage()->Add(newNode, node); } node->SetVisibility(false); } } void QmitkFiberProcessingView::GenerateRoiImage() { if (m_SelectedPF.empty()) return; mitk::BaseGeometry::Pointer geometry; if (!m_SelectedFB.empty()) { mitk::FiberBundle::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); geometry = fib->GetGeometry(); } else if (m_SelectedImage) geometry = m_SelectedImage->GetGeometry(); else return; itk::Vector 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()); std::string name = m_SelectedPF.at(0)->GetName(); WritePfToImage(m_SelectedPF.at(0), tmpImage); for (unsigned int i=1; iGetName(); WritePfToImage(m_SelectedPF.at(i), tmpImage); } DataNode::Pointer node = DataNode::New(); tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); node->SetData(tmpImage); node->SetName(name); this->GetDataStorage()->Add(node); } void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image) { if (dynamic_cast(node->GetData())) { m_PlanarFigure = dynamic_cast(node->GetData()); AccessFixedDimensionByItk_2( image, InternalReorientImagePlane, 3, m_PlanarFigure->GetGeometry(), -1); AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2, node->GetName() ); } else if (dynamic_cast(node->GetData())) { DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node); for (unsigned int i=0; iSize(); i++) { WritePfToImage(children->at(i), image); } } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ) { 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->GetExtentInMM(0) / spacing[0]; size[1] = planegeo->GetExtentInMM(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(unsigned 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::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); resampler->Update(); 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 ) { typedef itk::Image< TPixel, VImageDimension > ImageType; // 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 PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const BaseGeometry *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; 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; planarFigurePlaneGeometry->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigurePlaneGeometry->IndexToWorld(point2D, point2D); planarFigurePlaneGeometry->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.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->SetInputData( 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->SetInputConnection( extrudeFilter->GetOutputPort() ); // 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->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); 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.GoToBegin(); itimage.GoToBegin(); typename ImageType::SizeType lowersize = {{itk::NumericTraits::max(),itk::NumericTraits::max(),itk::NumericTraits::max()}}; 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 BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); const BaseGeometry *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::UpdateGui() { m_Controls->m_FibLabel->setText("mandatory"); m_Controls->m_PfLabel->setText("needed for extraction"); m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_RemoveButton->setEnabled(false); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(false); m_Controls->m_GenerateRoiImage->setEnabled(false); m_Controls->m_ExtractFibersButton->setEnabled(false); m_Controls->m_ModifyButton->setEnabled(false); m_Controls->m_CopyBundle->setEnabled(false); m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); // disable alle frames m_Controls->m_BundleWeightFrame->setVisible(false); m_Controls->m_ExtactionFramePF->setVisible(false); m_Controls->m_RemoveDirectionFrame->setVisible(false); m_Controls->m_RemoveLengthFrame->setVisible(false); m_Controls->m_RemoveCurvatureFrame->setVisible(false); m_Controls->m_RemoveByWeightFrame->setVisible(false); m_Controls->m_RemoveByDensityFrame->setVisible(false); m_Controls->m_SmoothFibersFrame->setVisible(false); m_Controls->m_CompressFibersFrame->setVisible(false); m_Controls->m_ColorFibersFrame->setVisible(false); m_Controls->m_MirrorFibersFrame->setVisible(false); m_Controls->m_MaskExtractionFrame->setVisible(false); m_Controls->m_ColorMapBox->setVisible(false); + m_Controls->m_ValueAsWeightBox->setVisible(false); bool pfSelected = !m_SelectedPF.empty(); bool fibSelected = !m_SelectedFB.empty(); bool multipleFibsSelected = (m_SelectedFB.size()>1); bool maskSelected = m_RoiImageNode.IsNotNull(); bool imageSelected = m_SelectedImage.IsNotNull(); // toggle visibility of elements according to selected method switch ( m_Controls->m_ExtractionMethodBox->currentIndex() ) { case 0: m_Controls->m_ExtactionFramePF->setVisible(true); break; case 1: m_Controls->m_MaskExtractionFrame->setVisible(true); break; } switch ( m_Controls->m_RemovalMethodBox->currentIndex() ) { case 0: m_Controls->m_RemoveDirectionFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 1: m_Controls->m_RemoveLengthFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 2: m_Controls->m_RemoveCurvatureFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 3: break; case 4: break; case 5: m_Controls->m_RemoveByWeightFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 6: m_Controls->m_RemoveByDensityFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; } switch ( m_Controls->m_ModificationMethodBox->currentIndex() ) { case 0: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 1: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 2: m_Controls->m_CompressFibersFrame->setVisible(true); break; case 3: + m_Controls->m_ValueAsWeightBox->setVisible(true); m_Controls->m_ColorFibersFrame->setVisible(true); m_Controls->m_ColorMapBox->setVisible(true); break; case 4: m_Controls->m_MirrorFibersFrame->setVisible(true); if (m_SelectedSurfaces.size()>0) m_Controls->m_ModifyButton->setEnabled(true); break; case 5: m_Controls->m_BundleWeightFrame->setVisible(true); break; case 6: + m_Controls->m_ValueAsWeightBox->setVisible(true); m_Controls->m_ColorFibersFrame->setVisible(true); break; case 7: m_Controls->m_ColorFibersFrame->setVisible(true); break; case 8: + m_Controls->m_ValueAsWeightBox->setVisible(true); m_Controls->m_ColorFibersFrame->setVisible(true); break; } // are fiber bundles selected? if ( fibSelected ) { m_Controls->m_CopyBundle->setEnabled(true); m_Controls->m_ModifyButton->setEnabled(true); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str())); // one bundle and one planar figure needed to extract fibers if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } // more than two bundles needed to join/subtract if (multipleFibsSelected) { m_Controls->m_FibLabel->setText("multiple bundles selected"); m_Controls->m_JoinBundles->setEnabled(true); m_Controls->m_SubstractBundles->setEnabled(true); } if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) ) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str())); m_Controls->m_RemoveButton->setEnabled(true); } } // are planar figures selected? if (pfSelected) { if ( fibSelected || m_SelectedImage.IsNotNull()) m_Controls->m_GenerateRoiImage->setEnabled(true); if (m_SelectedPF.size() > 1) { m_Controls->PFCompoANDButton->setEnabled(true); m_Controls->PFCompoORButton->setEnabled(true); } else m_Controls->PFCompoNOTButton->setEnabled(true); } // is image selected if (imageSelected || maskSelected) { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); } } void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* node ) { for (auto fnode: m_SelectedFB) if (node == fnode) { m_SelectedFB.clear(); break; } berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* ) { if (!m_Controls->m_InteractiveBox->isChecked()) { berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } } void QmitkFiberProcessingView::OnEndInteraction() { if (m_Controls->m_InteractiveBox->isChecked()) ExtractWithPlanarFigure(true); } void QmitkFiberProcessingView::AddObservers() { typedef itk::SimpleMemberCommand< QmitkFiberProcessingView > SimpleCommandType; for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) { figure->RemoveAllObservers(); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnEndInteraction); m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); } } } void QmitkFiberProcessingView::RemoveObservers() { for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) figure->RemoveAllObservers(); } } void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { RemoveObservers(); //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection std::vector lastSelectedFB = m_SelectedFB; m_SelectedFB.clear(); m_SelectedPF.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = nullptr; m_RoiImageNode = nullptr; for (auto node: nodes) { if ( dynamic_cast(node->GetData()) ) m_SelectedFB.push_back(node); else if (dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData())) m_SelectedPF.push_back(node); else if (dynamic_cast(node->GetData())) { m_SelectedImage = dynamic_cast(node->GetData()); if (m_SelectedImage->GetDimension()==3) m_RoiImageNode = node; } else if (dynamic_cast(node->GetData())) m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); } // if we perform interactive fiber extraction, we want to avoid auto-selection of the extracted bundle if (m_SelectedFB.empty() && m_Controls->m_InteractiveBox->isChecked()) m_SelectedFB = lastSelectedFB; // if no fibers or surfaces are selected, select topmost if (m_SelectedFB.empty() && m_SelectedSurfaces.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedFB.clear(); m_SelectedFB.push_back(nodes->at(i)); } } } // if no plar figure is selected, select topmost if (m_SelectedPF.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedPF.clear(); m_SelectedPF.push_back(nodes->at(i)); } } } AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::OnDrawPolygon() { mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->ClosedOn(); this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); } void QmitkFiberProcessingView::OnDrawCircle() { mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); } void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* ) { // initialize figure's geometry with empty geometry mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); figure->SetPlaneGeometry( 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->SetBoolProperty("planarfigure.3drendering", true); newNode->SetBoolProperty("planarfigure.3drendering.fill", true); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(newNode->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode(newNode); } // figure drawn on the topmost layer / image GetDataStorage()->Add(newNode ); RemoveObservers(); for(unsigned int i = 0; i < m_SelectedPF.size(); i++) m_SelectedPF[i]->SetSelected(false); newNode->SetSelected(true); m_SelectedPF.clear(); m_SelectedPF.push_back(newNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive) { if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); return; } try { std::vector fiberBundles = m_SelectedFB; mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0); for (unsigned int i=0; i(fiberBundles.at(i)->GetData()); mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage()); if (interactive && m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } float op = 5.0/sqrt(fib->GetNumFibers()); float currentOp = 0; fiberBundles.at(i)->GetFloatProperty("opacity", currentOp); if (currentOp!=op) { fib->SetFiberColors(255, 255, 255); fiberBundles.at(i)->SetFloatProperty("opacity", op); fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false); } m_InteractiveNode->SetData(extFB); } else { if (extFB->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::DataNode::Pointer node; node = mitk::DataNode::New(); node->SetData(extFB); QString name(fiberBundles.at(i)->GetName().c_str()); name += "*"; node->SetName(name.toStdString()); fiberBundles.at(i)->SetVisibility(false); GetDataStorage()->Add(node); } } } catch(const std::out_of_range& ) { QMessageBox::warning( nullptr, "Fiber extraction failed", "Did you only create the planar figure, using the circle or polygon button, but forgot to actually place it in the image afterwards? \nAfter creating a planar figure, simply left-click at the desired position in the image or on the tractogram to place it."); } } void QmitkFiberProcessingView::GenerateAndComposite() { mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("AND"); newPFCNode->SetData(PFCAnd); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); RemoveObservers(); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::GenerateOrComposite() { mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); PFCOr->setOperationType(mitk::PlanarFigureComposite::OR); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("OR"); newPFCNode->SetData(PFCOr); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); UpdateGui(); } void QmitkFiberProcessingView::GenerateNotComposite() { mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("NOT"); newPFCNode->SetData(PFCNot); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode ) { pfc->SetSelected(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(pfc, parentNode); else GetDataStorage()->Add(pfc); for (auto child : children) { if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); newChild->SetBoolProperty("planarfigure.3drendering", true); newChild->SetBoolProperty("planarfigure.3drendering.fill", true); GetDataStorage()->Add(newChild, pfc); GetDataStorage()->Remove(child); } else if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); std::vector< mitk::DataNode::Pointer > grandChildVector; mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child); for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it ) grandChildVector.push_back(*it); AddCompositeToDatastorage(newChild, grandChildVector, pfc); GetDataStorage()->Remove(child); } } UpdateGui(); } void QmitkFiberProcessingView::CopyBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy(); node->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newFib); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } UpdateGui(); } void QmitkFiberProcessingView::JoinBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } m_SelectedFB.at(0)->SetVisibility(false); mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); std::vector< mitk::FiberBundle::Pointer > tractograms; for (unsigned int i=1; iSetVisibility(false); tractograms.push_back(dynamic_cast(m_SelectedFB.at(i)->GetData())); } newBundle = newBundle->AddBundles(tractograms); mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName("Joined_Tractograms"); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } -void QmitkFiberProcessingView::SubstractBundles() +void QmitkFiberProcessingView::SubtractBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); m_SelectedFB.at(0)->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); for (unsigned int i=1; iSubtractBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); if (newBundle.IsNull()) break; name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str()); m_SelectedFB.at(i)->SetVisibility(false); } if (newBundle.IsNull()) { QMessageBox::information(nullptr, "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); UpdateGui(); } void QmitkFiberProcessingView::ResampleSelectedBundlesSpline() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleSpline(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ResampleSelectedBundlesLinear() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleLinear(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::CompressSelectedBundles() { double factor = this->m_Controls->m_ErrorThresholdBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->Compress(factor); fib->ColorFibersByOrientation(); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoImageColorCoding() { if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available."); return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); + fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked(), m_Controls->m_ValueAsWeightBox->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoCurvatureColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); + fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked(), m_Controls->m_ValueAsWeightBox->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoLengthColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - fib->ColorFibersByLength(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); + fib->ColorFibersByLength(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked(), m_Controls->m_ValueAsWeightBox->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoWeightColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::MirrorFibers() { unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->MirrorFibers(axis); } for (auto surf : m_SelectedSurfaces) { vtkSmartPointer 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(); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h index eba1582..a49c5d3 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h @@ -1,203 +1,203 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. 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, fiber resampling, mirroring, join and subtract bundles and much more. */ class QmitkFiberProcessingView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: typedef itk::Image< unsigned char, 3 > ItkUCharImageType; typedef itk::Image< float, 3 > ItkFloatImageType; static const std::string VIEW_ID; QmitkFiberProcessingView(); virtual ~QmitkFiberProcessingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; protected slots: void OnDrawCircle(); ///< add circle interactors etc. void OnDrawPolygon(); ///< add circle interactors etc. void GenerateAndComposite(); void GenerateOrComposite(); void GenerateNotComposite(); void CopyBundles(); ///< add copies of selected bundles to data storage void JoinBundles(); ///< merge selected fiber bundles - void SubstractBundles(); ///< subtract bundle A from bundle B. Not commutative! Defined by order of selection. + void SubtractBundles(); ///< subtract bundle A from bundle B. Not commutative! Defined by order of selection. void GenerateRoiImage(); ///< generate binary image of selected planar figures. void Remove(); void Extract(); void Modify(); void UpdateGui(); ///< update button activity etc. dpending on current datamanager selection void OnMaskExtractionChanged(); virtual void AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey = nullptr, mitk::BaseProperty *property = nullptr ); protected: void MirrorFibers(); ///< mirror bundle on the specified plane void ResampleSelectedBundlesSpline(); ///< void ResampleSelectedBundlesLinear(); ///< void DoImageColorCoding(); ///< color fibers by selected scalar image void DoWeightColorCoding(); ///< color fibers by their respective weights void DoLengthColorCoding(); ///< color fibers by length void DoCurvatureColorCoding(); ///< color fibers by curvature void CompressSelectedBundles(); ///< remove points below certain error threshold void WeightFibers(); void ApplyWeightThreshold(); void ApplyDensityThreshold(); void RemoveWithMask(bool removeInside); void RemoveDir(); void ApplyCurvatureThreshold(); ///< remove/split fibers with a too high curvature threshold void PruneBundle(); ///< remove too short/too long fibers void ExtractWithMask(bool onlyEnds, bool invert, bool labelmap); void ExtractWithPlanarFigure(bool interactive=false); void OnEndInteraction(); /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkFiberProcessingViewControls* m_Controls; /** 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::BaseGeometry* planegeo3D, int additionalIndex ); 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_RoiImageNode; unsigned int m_StartInteractionObserverTag; unsigned int m_EndInteractionObserverTag; mitk::DataNode::Pointer m_InteractiveNode; void AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode=nullptr); void debugPFComposition(mitk::PlanarFigureComposite::Pointer , int ); void WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image); mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute); mitk::DataNode::Pointer GenerateColorHeatmap(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib); void NodeAdded( const mitk::DataNode* node ) override; void NodeRemoved(const mitk::DataNode* node) override; void RemoveObservers(); void AddObservers(); }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui index 8708175..59fd8c3 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui @@ -1,1731 +1,1744 @@ QmitkFiberProcessingViewControls 0 0 385 711 Form QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 9 9 9 9 75 true 0 5 0 0 353 503 Fiber Extraction Extract a fiber subset from the selected fiber bundle using manually placed planar figures as waypoints or binary regions of interest. Qt::Vertical 20 40 QFrame::NoFrame QFrame::Raised 0 0 0 0 Min. overlap: Threshold: Min. num. fibers: Threshold on ROI image for positions to be considered as positive. 3 9999.000000000000000 0.100000000000000 0.500000000000000 0 0 Ending in ROI Not ending in ROI Passing ROI Not passing ROI Labelmap Both ends true Interpolate ROI true 1 999999999 Extract fibers: Minimum overlap of streamlines and ROI in terms of streamline length. Zero means that one streamline point inside the ROI is enough to be considered as "overlapping". 3 1.000000000000000 0.100000000000000 Labels: Label values seperated by whitespace. ALL false 0 0 200 16777215 11 Extract fibers passing through selected ROI or composite ROI. Select ROI and fiber bundle to execute. Extract 0 0 Extract using planar figures Extract using ROI image QFrame::NoFrame QFrame::Raised 0 0 0 0 6 Interactive Extraction 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 0 0 0 30 30 Draw circular ROI. Select reference fiber bundle to execute. - + :/org.mitk.gui.qt.diffusionimaging.fiberprocessing/resources/circle.png:/org.mitk.gui.qt.diffusionimaging.fiberprocessing/resources/circle.png 32 32 false true Qt::Horizontal 40 20 30 30 Draw polygonal ROI. Select reference fiber bundle to execute. - + :/org.mitk.gui.qt.diffusionimaging.fiberprocessing/resources/polygon.png:/org.mitk.gui.qt.diffusionimaging.fiberprocessing/resources/polygon.png 32 32 true true 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 0 0 0 false 60 16777215 Create NOT composition from selected ROI. NOT false 60 16777215 Create OR composition with selected ROIs. OR Qt::Horizontal 40 20 false 60 16777215 Create AND composition with selected ROIs. AND 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. Generate ROI Image 0 0 353 456 Fiber Removal Remove fibers that satisfy certain criteria from the selected bundle. QFrame::NoFrame QFrame::Raised 0 0 0 0 0 X: Y: Z: Angle: Angular deviation threshold in degree 1 90.000000000000000 1.000000000000000 25.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Weight threshold: Only fibers with weight larger than this threshold are kept. 5 99999.000000000000000 0.100000000000000 false 0 0 200 16777215 11 Remove 0 0 Remove fibers in direction Remove fibers by length Remove fibers by curvature Remove fiber parts outside mask Remove fiber parts inside mask Remove fibers by weight Remove fibers by tract density QFrame::NoFrame QFrame::Raised 0 0 0 0 If unchecked, the fiber exceeding the threshold will be split in two instead of removed. Remove Fiber false QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Max. Angular Deviation: Qt::Horizontal 40 20 Maximum angular deviation in degree 180.000000000000000 0.100000000000000 30.000000000000000 Distance: Distance in mm 1 999.000000000000000 1.000000000000000 10.000000000000000 Qt::Vertical 20 40 QFrame::NoFrame QFrame::Raised 0 0 0 0 Qt::Horizontal 40 20 Minimum fiber length in mm 0 999999999 20 Max. Length: Min. Length: Maximum fiber length in mm 0 999999999 300 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 5 0.000000000000000 10.000000000000000 0.100000000000000 0.010000000000000 Density threshold: Min. overlap: Tracts have to spend at least this fraction of their length inside the specified density region. 3 1.000000000000000 0.100000000000000 0.500000000000000 0 0 - 367 - 408 + 353 + 426 Bundle Modification Modify the selected bundle with operations such as fiber resampling, FA coloring, etc. QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 Error threshold in mm: 999999999.000000000000000 0.100000000000000 0.100000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 Sagittal Coronal Axial Select direction: QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 + + + + Scalar map: + + + If checked, the image values are not only used to color the fibers but are also used as opaxity values. Values as opacity false - - - - + + + The values used to color the fibers are min-max normalized. If not checked, the values should be between 0 and 1. + - Scalar map: + Normalize color values + + + true + + + - + - The values used to color the fibers are min-max normalized. If not checked, the values should be between 0 and 1. + If checked, the (mean) unnormalized values per fiber are also used as fiber weights, e.g. the fiber length in mm. - Normalize values + Values as weight - true + false 0 0 Resample fibers (spline) Resample fibers (linear) Compress fibers Color fibers by scalar map (e.g. FA) Mirror fibers Weight bundle Color fibers by curvature Color fibers by fiber weights Color fibers by length QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 0.010000000000000 999999999.000000000000000 0.100000000000000 1.000000000000000 Point distance in mm: Qt::Vertical 20 40 false 0 0 200 16777215 11 Execute QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Weight: 7 999999999.000000000000000 0.100000000000000 1.000000000000000 0 0 367 182 Bundle Operations Join, subtract or copy bundles. false 0 0 200 16777215 11 Returns all fibers contained in bundle X that are not contained in bundle Y (not commutative!). Select at least two fiber bundles to execute. Substract Qt::Vertical 20 40 false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Join false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Copy Please Select Input Data 6 6 6 6 <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true <html><head/><body><p><span style=" color:#969696;">needed for extraction</span></p></body></html> true Input DTI Fiber Bundle: Binary seed ROI. If not specified, the whole image area is seeded. ROI: Qt::Vertical 20 40 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
- +
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMView.cpp index aa5bb54..1b13de9 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMView.cpp @@ -1,1095 +1,1095 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. 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 "QmitkIVIMView.h" // qt #include "qmessagebox.h" #include "qclipboard.h" // mitk #include "mitkImage.h" #include "mitkImageCast.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include #include // itk #include "itkScalarImageToHistogramGenerator.h" #include "itkRegionOfInterestImageFilter.h" #include "itkImageRegionConstIteratorWithIndex.h" // itk/mitk #include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" #include "itkRegularizedIVIMReconstructionFilter.h" #include "mitkImageCast.h" #include #include #include #include #include #include #include const std::string QmitkIVIMView::VIEW_ID = "org.mitk.views.ivim"; QmitkIVIMView::QmitkIVIMView() : QmitkAbstractView() , m_Controls( 0 ) , m_Active(false) , m_Visible(false) , m_HoldUpdate(false) { } QmitkIVIMView::~QmitkIVIMView() { } void QmitkIVIMView::CreateQtPartControl( QWidget *parent ) { // hold update untill all elements are set this->m_HoldUpdate = true; // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkIVIMViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_ButtonStart, SIGNAL(clicked()), this, SLOT(FittIVIMStart()) ); connect( m_Controls->m_ButtonAutoThres, SIGNAL(clicked()), this, SLOT(AutoThreshold()) ); connect( m_Controls->m_MethodCombo, SIGNAL(currentIndexChanged(int)), this, SLOT(OnIvimFitChanged(int)) ); connect( m_Controls->m_DStarSlider, SIGNAL(valueChanged(int)), this, SLOT(DStarSlider(int)) ); connect( m_Controls->m_BThreshSlider, SIGNAL(valueChanged(int)), this, SLOT(BThreshSlider(int)) ); connect( m_Controls->m_S0ThreshSlider, SIGNAL(valueChanged(int)), this, SLOT(S0ThreshSlider(int)) ); connect( m_Controls->m_NumItSlider, SIGNAL(valueChanged(int)), this, SLOT(NumItsSlider(int)) ); connect( m_Controls->m_LambdaSlider, SIGNAL(valueChanged(int)), this, SLOT(LambdaSlider(int)) ); connect( m_Controls->m_CheckDStar, SIGNAL(clicked()), this, SLOT(Checkbox()) ); connect( m_Controls->m_CheckD, SIGNAL(clicked()), this, SLOT(Checkbox()) ); connect( m_Controls->m_Checkf, SIGNAL(clicked()), this, SLOT(Checkbox()) ); connect( m_Controls->m_CurveClipboard, SIGNAL(clicked()), this, SLOT(ClipboardCurveButtonClicked()) ); connect( m_Controls->m_ValuesClipboard, SIGNAL(clicked()), this, SLOT(ClipboardStatisticsButtonClicked()) ); connect( m_Controls->m_SavePlot, SIGNAL(clicked()), this, SLOT(SavePlotButtonClicked()) ); // connect all kurtosis actions to a recompute connect( m_Controls->m_KurtosisRangeWidget, SIGNAL( rangeChanged(double, double)), this, SLOT(OnKurtosisParamsChanged() ) ); connect( m_Controls->m_OmitBZeroCB, SIGNAL( stateChanged(int) ), this, SLOT( OnKurtosisParamsChanged() ) ); connect( m_Controls->m_KurtosisFitScale, SIGNAL( currentIndexChanged(int)), this, SLOT( OnKurtosisParamsChanged() ) ); connect( m_Controls->m_UseKurtosisBoundsCB, SIGNAL(clicked() ), this, SLOT( OnKurtosisParamsChanged() ) ); m_Controls->m_DwiBox->SetDataStorage(this->GetDataStorage()); mitk::NodePredicateIsDWI::Pointer isDwi = mitk::NodePredicateIsDWI::New(); m_Controls->m_DwiBox->SetPredicate( isDwi ); connect( (QObject*)(m_Controls->m_DwiBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); m_Controls->m_MaskBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskBox->SetZeroEntryText("--"); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateDimension::Pointer is3D = mitk::NodePredicateDimension::New(3); m_Controls->m_MaskBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, mitk::NodePredicateAnd::New(isImagePredicate, is3D)) ); connect( (QObject*)(m_Controls->m_MaskBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect( (QObject*)(m_Controls->m_ModelTabSelectionWidget), SIGNAL(currentChanged(int)), this, SLOT(OnModelTabChanged(int))); } QString dstar = QString::number(m_Controls->m_DStarSlider->value()/1000.0); m_Controls->m_DStarLabel->setText(dstar); QString bthresh = QString::number(m_Controls->m_BThreshSlider->value()*5.0); m_Controls->m_BThreshLabel->setText(bthresh); QString s0thresh = QString::number(m_Controls->m_S0ThreshSlider->value()*0.5); m_Controls->m_S0ThreshLabel->setText(s0thresh); QString numits = QString::number(m_Controls->m_NumItSlider->value()); m_Controls->m_NumItsLabel->setText(numits); QString lambda = QString::number(m_Controls->m_LambdaSlider->value()*.00001); m_Controls->m_LambdaLabel->setText(lambda); m_Controls->m_Warning->setVisible(false); OnIvimFitChanged(m_Controls->m_MethodCombo->currentIndex()); m_Controls->m_KurtosisRangeWidget->setSingleStep(0.1); m_Controls->m_KurtosisRangeWidget->setRange( 0.0, 10.0 ); m_Controls->m_KurtosisRangeWidget->setMaximumValue( 5.0 ); // LogScale not working yet, have to fix that first // m_Controls->m_KurtosisFitScale->setEnabled(false); //m_Controls->m_MaximalBValueWidget->setVisible( false ); // release update block after the UI-elements were all set this->m_HoldUpdate = false; QmitkIVIMView::InitChartIvim(); m_ListenerActive = false; if (this->GetRenderWindowPart()) { m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); m_ListenerActive = true; } } void QmitkIVIMView::OnModelTabChanged(int tab) { if (tab==0) InitChartIvim(); else if (tab==1) InitChartKurtosis(); UpdateGui(); } void QmitkIVIMView::AddSecondFitPlot() { if (m_Controls->m_ChartWidget->GetDataElementByLabel("signal values (used for second fit)") == nullptr) { std::vector< std::pair > init_data; m_Controls->m_ChartWidget->AddData2D(init_data, "signal values (used for second fit)", QmitkChartWidget::ChartType::scatter); m_Controls->m_ChartWidget->SetColor("signal values (used for second fit)", "white"); m_Controls->m_ChartWidget->SetMarkerSymbol("signal values (used for second fit)", QmitkChartWidget::MarkerSymbol::x_thin); m_Controls->m_ChartWidget->Show(true); m_Controls->m_ChartWidget->SetShowDataPoints(false); m_Controls->m_ChartWidget->SetShowSubchart(false); } } void QmitkIVIMView::RemoveSecondFitPlot() { if (m_Controls->m_ChartWidget->GetDataElementByLabel("signal values (used for second fit)") != nullptr) { m_Controls->m_ChartWidget->RemoveData("signal values (used for second fit)"); m_Controls->m_ChartWidget->Show(true); m_Controls->m_ChartWidget->SetShowDataPoints(false); m_Controls->m_ChartWidget->SetShowSubchart(false); } } void QmitkIVIMView::InitChartIvim() { m_Controls->m_ChartWidget->Clear(); std::vector< std::pair > init_data; m_Controls->m_ChartWidget->AddData2D(init_data, "D-part of fitted model", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData2D(init_data, "fitted model", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData2D(init_data, "signal values", QmitkChartWidget::ChartType::scatter); m_Controls->m_ChartWidget->SetColor("fitted model", "red"); m_Controls->m_ChartWidget->SetColor("signal values", "white"); m_Controls->m_ChartWidget->SetColor("D-part of fitted model", "cyan"); m_Controls->m_ChartWidget->SetYAxisScale(QmitkChartWidget::AxisScale::log); m_Controls->m_ChartWidget->SetYAxisLabel("S/S0"); m_Controls->m_ChartWidget->SetXAxisLabel("b-value"); m_Controls->m_ChartWidget->SetLineStyle("fitted model", QmitkChartWidget::LineStyle::solid); m_Controls->m_ChartWidget->SetLineStyle("D-part of fitted model", QmitkChartWidget::LineStyle::dashed); m_Controls->m_ChartWidget->SetMarkerSymbol("signal values", QmitkChartWidget::MarkerSymbol::diamond); m_Controls->m_ChartWidget->Show(true); m_Controls->m_ChartWidget->SetShowDataPoints(false); m_Controls->m_ChartWidget->SetShowSubchart(false); } void QmitkIVIMView::InitChartKurtosis() { m_Controls->m_ChartWidget->Clear(); std::vector< std::pair > init_data; m_Controls->m_ChartWidget->AddData2D(init_data, "D-part of fitted model", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData2D(init_data, "fitted model", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData2D(init_data, "signal values", QmitkChartWidget::ChartType::scatter); m_Controls->m_ChartWidget->SetColor("fitted model", "red"); m_Controls->m_ChartWidget->SetColor("signal values", "white"); m_Controls->m_ChartWidget->SetColor("D-part of fitted model", "cyan"); m_Controls->m_ChartWidget->SetYAxisScale(QmitkChartWidget::AxisScale::log); m_Controls->m_ChartWidget->SetYAxisLabel("S"); m_Controls->m_ChartWidget->SetXAxisLabel("b-value"); m_Controls->m_ChartWidget->SetLineStyle("fitted model", QmitkChartWidget::LineStyle::solid); m_Controls->m_ChartWidget->SetLineStyle("D-part of fitted model", QmitkChartWidget::LineStyle::dashed); m_Controls->m_ChartWidget->SetMarkerSymbol("signal values", QmitkChartWidget::MarkerSymbol::diamond); m_Controls->m_ChartWidget->Show(true); m_Controls->m_ChartWidget->SetShowDataPoints(false); m_Controls->m_ChartWidget->SetShowSubchart(false); } void QmitkIVIMView::SetFocus() { m_Controls->m_ButtonAutoThres->setFocus(); } void QmitkIVIMView::Checkbox() { OnSliceChanged(); } void QmitkIVIMView::OnIvimFitChanged(int val) { switch(val) { case 0: m_Controls->m_DstarFrame->setVisible(false); m_Controls->m_NeglSiFrame->setVisible(true); m_Controls->m_NeglBframe->setVisible(false); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; case 1: m_Controls->m_DstarFrame->setVisible(true); m_Controls->m_NeglSiFrame->setVisible(true); m_Controls->m_NeglBframe->setVisible(false); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; case 2: m_Controls->m_DstarFrame->setVisible(false); m_Controls->m_NeglSiFrame->setVisible(true); m_Controls->m_NeglBframe->setVisible(true); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; case 3: m_Controls->m_DstarFrame->setVisible(false); m_Controls->m_NeglSiFrame->setVisible(true); m_Controls->m_NeglBframe->setVisible(true); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; case 4: m_Controls->m_DstarFrame->setVisible(false); m_Controls->m_NeglSiFrame->setVisible(false); m_Controls->m_NeglBframe->setVisible(false); m_Controls->m_IterationsFrame->setVisible(false); m_Controls->m_LambdaFrame->setVisible(false); break; } OnSliceChanged(); } void QmitkIVIMView::DStarSlider (int val) { QString sval = QString::number(val/1000.0); m_Controls->m_DStarLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::BThreshSlider (int val) { QString sval = QString::number(val*5.0); m_Controls->m_BThreshLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::S0ThreshSlider (int val) { QString sval = QString::number(val*0.5); m_Controls->m_S0ThreshLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::NumItsSlider (int val) { QString sval = QString::number(val); m_Controls->m_NumItsLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::LambdaSlider (int val) { QString sval = QString::number(val*.00001); m_Controls->m_LambdaLabel->setText(sval); OnSliceChanged(); } void QmitkIVIMView::UpdateGui() { m_Controls->m_FittedParamsLabel->setText(""); if (m_Controls->m_DwiBox->GetSelectedNode().IsNotNull()) { m_Controls->m_ChartWidget->setVisible(true); m_HoldUpdate = false; } else { m_Controls->m_ChartWidget->setVisible(false); } m_Controls->m_ButtonStart->setEnabled( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ); m_Controls->m_ButtonAutoThres->setEnabled( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ); m_Controls->m_ControlsFrame->setEnabled( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ); m_Controls->m_BottomControlsFrame->setEnabled( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ); OnSliceChanged(); } void QmitkIVIMView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { // UpdateGui(); } void QmitkIVIMView::AutoThreshold() { if (m_Controls->m_DwiBox->GetSelectedNode().IsNull()) { // Nothing selected. Inform the user and return QMessageBox::information( nullptr, "Template", "Please load and select a diffusion image before starting image processing."); return; } mitk::Image* dimg = dynamic_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); if (!dimg) { // Nothing selected. Inform the user and return QMessageBox::information( nullptr, "Template", "No valid diffusion image was found."); return; } // find bzero index int index = -1; auto directions = mitk::DiffusionPropertyHelper::GetGradientContainer(dimg); for(DirContainerType::ConstIterator it = directions->Begin(); it != directions->End(); ++it) { index++; GradientDirectionType g = it.Value(); if(g[0] == 0 && g[1] == 0 && g[2] == 0 ) break; } VecImgType::Pointer vecimg = VecImgType::New(); mitk::CastToItkImage(dimg, vecimg); int vecLength = vecimg->GetVectorLength(); index = index > vecLength-1 ? vecLength-1 : index; MITK_INFO << "Performing Histogram Analysis on Channel" << index; typedef itk::Image ImgType; ImgType::Pointer img = ImgType::New(); mitk::CastToItkImage(dimg, img); itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); itw.GoToBegin(); itk::ImageRegionConstIterator itr (vecimg, vecimg->GetLargestPossibleRegion() ); itr.GoToBegin(); while(!itr.IsAtEnd()) { itw.Set(itr.Get().GetElement(index)); ++itr; ++itw; } typedef itk::Statistics::ScalarImageToHistogramGenerator< ImgType > HistogramGeneratorType; typedef HistogramGeneratorType::HistogramType HistogramType; HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( img ); histogramGenerator->SetMarginalScale( 10 ); // Defines y-margin width of histogram histogramGenerator->SetNumberOfBins( 100 ); // CT range [-1024, +2048] --> bin size 4 values histogramGenerator->SetHistogramMin( dimg->GetStatistics()->GetScalarValueMin() ); histogramGenerator->SetHistogramMax( dimg->GetStatistics()->GetScalarValueMax() * .5 ); histogramGenerator->Compute(); HistogramType::ConstIterator iter = histogramGenerator->GetOutput()->Begin(); float maxFreq = 0; float maxValue = 0; while ( iter != histogramGenerator->GetOutput()->End() ) { if(iter.GetFrequency() > maxFreq) { maxFreq = iter.GetFrequency(); maxValue = iter.GetMeasurementVector()[0]; } ++iter; } maxValue *= 2; int sliderPos = maxValue * 2; m_Controls->m_S0ThreshSlider->setValue(sliderPos); S0ThreshSlider(sliderPos); } void QmitkIVIMView::FittIVIMStart() { if (m_Controls->m_DwiBox->GetSelectedNode().IsNull()) { QMessageBox::information( nullptr, "Template", "No valid diffusion-weighted image selected."); return; } mitk::Image* img = dynamic_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); VecImgType::Pointer vecimg = VecImgType::New(); mitk::CastToItkImage(img, vecimg); OutImgType::IndexType dummy; if( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { // KURTOSIS KurtosisFilterType::Pointer filter = KurtosisFilterType::New(); filter->SetInput(vecimg); filter->SetReferenceBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(img)); filter->SetGradientDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(img)); filter->SetSmoothingSigma( m_Controls->m_SigmaSpinBox->value() ); if( m_Controls->m_UseKurtosisBoundsCB->isChecked() ) filter->SetBoundariesForKurtosis( m_Controls->m_KurtosisRangeWidget->minimumValue(), m_Controls->m_KurtosisRangeWidget->maximumValue() ); filter->SetFittingScale( static_cast(m_Controls->m_KurtosisFitScale->currentIndex() ) ); if( m_Controls->m_MaskBox->GetSelectedNode().IsNotNull() ) { mitk::Image::Pointer maskImg = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); typedef itk::Image MaskImgType; MaskImgType::Pointer maskItk; CastToItkImage( maskImg, maskItk ); filter->SetImageMask( maskItk ); } filter->Update(); mitk::LookupTable::Pointer kurt_map_lut = mitk::LookupTable::New(); kurt_map_lut->SetType( mitk::LookupTable::JET ); mitk::LookupTableProperty::Pointer kurt_lut_prop = mitk::LookupTableProperty::New(); kurt_lut_prop->SetLookupTable( kurt_map_lut ); mitk::Image::Pointer dimage = mitk::Image::New(); dimage->InitializeByItk( filter->GetOutput(0) ); dimage->SetVolume( filter->GetOutput(0)->GetBufferPointer()); mitk::Image::Pointer kimage = mitk::Image::New(); kimage->InitializeByItk( filter->GetOutput(1) ); kimage->SetVolume( filter->GetOutput(1)->GetBufferPointer()); QString new_dname = "Kurtosis_DMap"; new_dname.append("_Method-"+m_Controls->m_KurtosisFitScale->currentText()); QString new_kname = "Kurtosis_KMap"; new_kname.append("_Method-"+m_Controls->m_KurtosisFitScale->currentText()); if( m_Controls->m_CheckKurtD->isChecked() ) { mitk::DataNode::Pointer dnode = mitk::DataNode::New(); dnode->SetData( dimage ); dnode->SetName(new_dname.toLatin1()); dnode->SetProperty("LookupTable", kurt_lut_prop ); GetDataStorage()->Add(dnode, m_Controls->m_DwiBox->GetSelectedNode()); } if( m_Controls->m_CheckKurtK->isChecked() ) { mitk::DataNode::Pointer knode = mitk::DataNode::New(); knode->SetData( kimage ); knode->SetName(new_kname.toLatin1()); knode->SetProperty("LookupTable", kurt_lut_prop ); GetDataStorage()->Add(knode, m_Controls->m_DwiBox->GetSelectedNode()); } } else { FittIVIM(vecimg, mitk::DiffusionPropertyHelper::GetGradientContainer(img), mitk::DiffusionPropertyHelper::GetReferenceBValue(img), true, dummy); OutputToDatastorage(m_Controls->m_DwiBox->GetSelectedNode()); } } void QmitkIVIMView::OnKurtosisParamsChanged() { OnSliceChanged(); } void QmitkIVIMView::OnSliceChanged() { if(m_HoldUpdate || !m_Visible) return; m_Controls->m_Warning->setVisible(false); if(m_Controls->m_DwiBox->GetSelectedNode().IsNull()) return; mitk::Image::Pointer diffusionImg = dynamic_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); mitk::Image::Pointer maskImg = nullptr; if (m_Controls->m_MaskBox->GetSelectedNode().IsNotNull()) maskImg = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); if (!this->GetRenderWindowPart()) return; if (!m_ListenerActive) { m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); m_ListenerActive = true; } VecImgType::Pointer vecimg = VecImgType::New(); mitk::CastToItkImage(diffusionImg, vecimg); VecImgType::Pointer roiImage = VecImgType::New(); if(maskImg.IsNull()) { int roisize = 0; if(m_Controls->m_MethodCombo->currentIndex() == 4) roisize = 3; mitk::Point3D pos = this->GetRenderWindowPart()->GetSelectedPosition(); VecImgType::IndexType crosspos; diffusionImg->GetGeometry()->WorldToIndex(pos, crosspos); if (!vecimg->GetLargestPossibleRegion().IsInside(crosspos)) { m_Controls->m_Warning->setText(QString("Crosshair position not inside of selected diffusion weighted image. Reinit needed!")); m_Controls->m_Warning->setVisible(true); return; } else m_Controls->m_Warning->setVisible(false); VecImgType::IndexType index; index[0] = crosspos[0] - roisize; index[0] = index[0] < 0 ? 0 : index[0]; index[1] = crosspos[1] - roisize; index[1] = index[1] < 0 ? 0 : index[1]; index[2] = crosspos[2] - roisize; index[2] = index[2] < 0 ? 0 : index[2]; VecImgType::SizeType size; size[0] = roisize*2+1; size[1] = roisize*2+1; size[2] = roisize*2+1; VecImgType::SizeType maxSize = vecimg->GetLargestPossibleRegion().GetSize(); size[0] = index[0]+size[0] > maxSize[0] ? maxSize[0]-index[0] : size[0]; size[1] = index[1]+size[1] > maxSize[1] ? maxSize[1]-index[1] : size[1]; size[2] = index[2]+size[2] > maxSize[2] ? maxSize[2]-index[2] : size[2]; VecImgType::RegionType region; region.SetSize( size ); region.SetIndex( index ); vecimg->SetRequestedRegion( region ); VecImgType::IndexType newstart; newstart.Fill(0); VecImgType::RegionType newregion; newregion.SetSize( size ); newregion.SetIndex( newstart ); roiImage->CopyInformation( vecimg ); roiImage->SetRegions( newregion ); roiImage->SetOrigin( pos ); roiImage->Allocate(); // roiImage->SetPixel(newstart, vecimg->GetPixel(index)); typedef itk::ImageRegionIterator VectorIteratorType; VectorIteratorType vecit(vecimg, vecimg->GetRequestedRegion() ); vecit.GoToBegin(); typedef itk::ImageRegionIterator VectorIteratorType; VectorIteratorType vecit2(roiImage, roiImage->GetLargestPossibleRegion() ); vecit2.GoToBegin(); while( !vecit.IsAtEnd() ) { vecit2.Set(vecit.Get()); ++vecit; ++vecit2; } if( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { FitKurtosis(roiImage, mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImg), mitk::DiffusionPropertyHelper::GetReferenceBValue(diffusionImg), newstart); } else { FittIVIM(roiImage, mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImg), mitk::DiffusionPropertyHelper::GetReferenceBValue(diffusionImg), false, crosspos); } } else { typedef itk::Image MaskImgType; MaskImgType::Pointer maskItk; CastToItkImage( maskImg, maskItk ); mitk::Point3D pos; pos[0] = 0; pos[1] = 0; pos[2] = 0; VecImgType::IndexType index; index[0] = 0; index[1] = 0; index[2] = 0; VecImgType::SizeType size; size[0] = 1; size[1] = 1; size[2] = 1; VecImgType::RegionType region; region.SetSize( size ); region.SetIndex( index ); vecimg->SetRequestedRegion( region ); // iterators over output and input itk::ImageRegionConstIteratorWithIndex vecit(vecimg, vecimg->GetLargestPossibleRegion()); itk::VariableLengthVector avg(vecimg->GetVectorLength()); avg.Fill(0); float numPixels = 0; while ( ! vecit.IsAtEnd() ) { VecImgType::PointType point; vecimg->TransformIndexToPhysicalPoint(vecit.GetIndex(), point); MaskImgType::IndexType index; maskItk->TransformPhysicalPointToIndex(point, index); if(maskItk->GetPixel(index) != 0) { avg += vecit.Get(); numPixels += 1.0; } // update iterators ++vecit; } avg /= numPixels; m_Controls->m_Warning->setText(QString("Averaging ")+QString::number((int)numPixels)+QString(" voxels!")); m_Controls->m_Warning->setVisible(true); roiImage->CopyInformation( vecimg ); roiImage->SetRegions( region ); roiImage->SetOrigin( pos ); roiImage->Allocate(); roiImage->SetPixel(index, avg); if( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { FitKurtosis(roiImage, mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImg), mitk::DiffusionPropertyHelper::GetReferenceBValue(diffusionImg), index); } else { FittIVIM(roiImage, mitk::DiffusionPropertyHelper::GetGradientContainer(diffusionImg), mitk::DiffusionPropertyHelper::GetReferenceBValue(diffusionImg), false, index); } // do not update until selection changed, the values will remain the same as long as the mask is selected! m_HoldUpdate = true; } vecimg->SetRegions( vecimg->GetLargestPossibleRegion() ); } bool QmitkIVIMView::FitKurtosis( itk::VectorImage *vecimg, DirContainerType *dirs, float bval, OutImgType::IndexType &crosspos ) { KurtosisFilterType::Pointer filter = KurtosisFilterType::New(); itk::KurtosisFitConfiguration fit_config; fit_config.omit_bzero = m_Controls->m_OmitBZeroCB->isChecked(); if( m_Controls->m_UseKurtosisBoundsCB->isChecked() ) { fit_config.use_K_limits = true; vnl_vector_fixed k_limits; k_limits[0] = m_Controls->m_KurtosisRangeWidget->minimumValue(); k_limits[1] = m_Controls->m_KurtosisRangeWidget->maximumValue(); fit_config.K_limits = k_limits; } fit_config.fit_scale = static_cast(m_Controls->m_KurtosisFitScale->currentIndex() ); m_KurtosisSnap = filter->GetSnapshot( vecimg->GetPixel( crosspos ), dirs, bval, fit_config ); QString param_label_text("K=%2, D=%1"); param_label_text = param_label_text.arg( m_KurtosisSnap.m_K, 4); param_label_text = param_label_text.arg( m_KurtosisSnap.m_D, 4); m_Controls->m_FittedParamsLabel->setText(param_label_text); const double maxb = m_KurtosisSnap.bvalues.max_value(); double S0 = m_KurtosisSnap.m_Bzero; if (m_KurtosisSnap.m_fittedBZero) S0 = m_KurtosisSnap.m_BzeroFit; std::vector< std::pair > d_line; d_line.emplace_back(0, S0); d_line.emplace_back(maxb, S0*exp(-maxb * m_KurtosisSnap.m_D)); m_Controls->m_ChartWidget->UpdateData2D(d_line, "D-part of fitted model"); const unsigned int num_samples = 50; std::vector< std::pair > y; for( unsigned int i=0; i<=num_samples; ++i) { double b = (((1.0)*i)/(1.0*num_samples))*maxb; y.emplace_back(b, S0 * exp( -b * m_KurtosisSnap.m_D + b*b * m_KurtosisSnap.m_D * m_KurtosisSnap.m_D * m_KurtosisSnap.m_K / 6.0 )); } m_Controls->m_ChartWidget->UpdateData2D(y, "fitted model"); std::vector< std::pair > y_meas; for (unsigned int i=0; im_OmitBZeroCB->isChecked() || m_KurtosisSnap.bvalues[i] > 0.01) y_meas.emplace_back(m_KurtosisSnap.bvalues[i], m_KurtosisSnap.measurements[i]); } m_Controls->m_ChartWidget->UpdateData2D(y_meas, "signal values"); return true; } bool QmitkIVIMView::FittIVIM(itk::VectorImage* vecimg, DirContainerType* dirs, float bval, bool multivoxel, OutImgType::IndexType &crosspos) { IVIMFilterType::Pointer filter = IVIMFilterType::New(); filter->SetInput(vecimg); filter->SetGradientDirections(dirs); filter->SetBValue(bval); switch(m_Controls->m_MethodCombo->currentIndex()) { case 0: filter->SetMethod(IVIMFilterType::IVIM_FIT_ALL); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); break; case 1: filter->SetMethod(IVIMFilterType::IVIM_DSTAR_FIX); filter->SetDStar(m_Controls->m_DStarLabel->text().toDouble()); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); break; case 2: filter->SetMethod(IVIMFilterType::IVIM_D_THEN_DSTAR); filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); break; case 3: filter->SetMethod(IVIMFilterType::IVIM_LINEAR_D_THEN_F); filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); break; case 4: filter->SetMethod(IVIMFilterType::IVIM_REGULARIZED); filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); filter->SetNumberIterations(m_Controls->m_NumItsLabel->text().toInt()); filter->SetLambda(m_Controls->m_LambdaLabel->text().toDouble()); filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); break; } if(!multivoxel) { filter->SetFitDStar(true); } filter->SetNumberOfThreads(1); filter->SetVerbose(false); filter->SetCrossPosition(crosspos); try{ filter->Update(); m_IvimSnap = filter->GetSnapshot(); m_DStarMap = filter->GetOutput(2); m_DMap = filter->GetOutput(1); m_fMap = filter->GetOutput(); QString param_label_text("f=%1, D=%2, D*=%3"); param_label_text = param_label_text.arg(m_IvimSnap.currentF,4); param_label_text = param_label_text.arg(m_IvimSnap.currentD,4); param_label_text = param_label_text.arg(m_IvimSnap.currentDStar,4); m_Controls->m_FittedParamsLabel->setText(param_label_text); double maxb = m_IvimSnap.bvalues.max_value(); std::vector< std::pair > d_line; d_line.emplace_back(0, 1-m_IvimSnap.currentFunceiled); d_line.emplace_back(maxb, d_line[0].second*exp(-maxb * m_IvimSnap.currentD)); m_Controls->m_ChartWidget->UpdateData2D(d_line, "D-part of fitted model"); std::vector< std::pair > y; int nsampling = 50; double f = 1-m_IvimSnap.currentFunceiled; for(int i=0; i<=nsampling; i++) { double x = (((1.0)*i)/(1.0*nsampling))*maxb; y.emplace_back(x, f*exp(- x * m_IvimSnap.currentD) + (1-f)*exp(- x * (m_IvimSnap.currentD+m_IvimSnap.currentDStar))); } m_Controls->m_ChartWidget->UpdateData2D(y, "fitted model"); std::vector< std::pair > y_meas; for (unsigned int i=0; im_ChartWidget->UpdateData2D(y_meas, "signal values"); - if(m_IvimSnap.num_high) // num_high is 0 if there is only one fit + if(m_IvimSnap.bvals2.size() > 0) { AddSecondFitPlot(); std::vector< std::pair > additonal_meas; - for (int i=0; i(m_IvimSnap.bvals2.size()); ++i) additonal_meas.emplace_back(m_IvimSnap.bvals2[i], m_IvimSnap.meas2[i]); m_Controls->m_ChartWidget->UpdateData2D(additonal_meas, "signal values (used for second fit)"); } else { RemoveSecondFitPlot(); } } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; m_Controls->m_Warning->setText(QString("IVIM fit not possible: ")+ex.GetDescription()); m_Controls->m_Warning->setVisible(true); return false; } return true; } void QmitkIVIMView::OutputToDatastorage(mitk::DataNode::Pointer node) { mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType( mitk::LookupTable::JET ); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable( lut ); if(m_Controls->m_CheckDStar->isChecked()) { mitk::Image::Pointer dstarimage = mitk::Image::New(); dstarimage->InitializeByItk(m_DStarMap.GetPointer()); dstarimage->SetVolume(m_DStarMap->GetBufferPointer()); QString newname2 = ""; newname2 = newname2.append("D* Fit-%1").arg(m_Controls->m_MethodCombo->currentIndex()+1); mitk::DataNode::Pointer node2=mitk::DataNode::New(); node2->SetData( dstarimage ); node2->SetName(newname2.toLatin1()); node2->SetProperty("LookupTable", lut_prop ); mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); mitk::LevelWindow levelwindow; levelwindow.SetRangeMinMax(0, 0.2); levWinProp->SetLevelWindow(levelwindow); node2->SetProperty("levelwindow", levWinProp); GetDataStorage()->Add(node2, node); } if(m_Controls->m_CheckD->isChecked()) { mitk::Image::Pointer dimage = mitk::Image::New(); dimage->InitializeByItk(m_DMap.GetPointer()); dimage->SetVolume(m_DMap->GetBufferPointer()); QString newname1 = ""; newname1 = newname1.append("D Fit-%1").arg(m_Controls->m_MethodCombo->currentIndex()+1); mitk::DataNode::Pointer node1=mitk::DataNode::New(); node1->SetData( dimage ); node1->SetName(newname1.toLatin1()); node1->SetProperty("LookupTable", lut_prop ); mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); mitk::LevelWindow levelwindow; levelwindow.SetRangeMinMax(0, 0.003); levWinProp->SetLevelWindow(levelwindow); node1->SetProperty("levelwindow", levWinProp); GetDataStorage()->Add(node1, node); } if(m_Controls->m_Checkf->isChecked()) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(m_fMap.GetPointer()); image->SetVolume(m_fMap->GetBufferPointer()); QString newname0 = ""; newname0 = newname0.append("f Fit-%1").arg(m_Controls->m_MethodCombo->currentIndex()+1); mitk::DataNode::Pointer node3=mitk::DataNode::New(); node3->SetData( image ); node3->SetName(newname0.toLatin1()); node3->SetProperty("LookupTable", lut_prop ); mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); mitk::LevelWindow levelwindow; levelwindow.SetRangeMinMax(0, 1); levWinProp->SetLevelWindow(levelwindow); node3->SetProperty("levelwindow", levWinProp); GetDataStorage()->Add(node3, node); } this->GetRenderWindowPart()->RequestUpdate(); } void QmitkIVIMView::ClipboardCurveButtonClicked() { // Kurtosis if ( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { std::stringstream ss; QString clipboard("Measurement Points\n"); ss << m_KurtosisSnap.bvalues << "\n" << m_KurtosisSnap.measurements << "\n\n"; ss << "Fitted Values ( D K [b_0] ) \n" << m_KurtosisSnap.m_D << " " << m_KurtosisSnap.m_K; if( m_KurtosisSnap.m_fittedBZero ) ss << " " << m_KurtosisSnap.m_BzeroFit; ss << "\n\n"; clipboard.append( QString( ss.str().c_str() )); ss.str( std::string() ); ss.clear(); QApplication::clipboard()->setText(clipboard, QClipboard::Clipboard ); } else { std::stringstream ss; QString clipboard("Normalized Measurement Points\n"); ss << m_IvimSnap.bvalues << "\n" << m_IvimSnap.allmeas << "\n\n"; ss << "Fitted Values ( f D D* ) \n" << m_IvimSnap.currentF << " " << m_IvimSnap.currentD << " " << m_IvimSnap.currentDStar; ss << "\n\n"; clipboard.append( QString( ss.str().c_str() )); ss.str( std::string() ); ss.clear(); QApplication::clipboard()->setText(clipboard, QClipboard::Clipboard ); } } void QmitkIVIMView::SavePlotButtonClicked() { m_Controls->m_ChartWidget->SavePlotAsImage(); } void QmitkIVIMView::ClipboardStatisticsButtonClicked() { // Kurtosis if ( m_Controls->m_ModelTabSelectionWidget->currentIndex() ) { QString clipboard( "D \t K \n" ); clipboard = clipboard.append( "%L1 \t %L2" ) .arg( m_KurtosisSnap.m_D, 0, 'f', 10 ) .arg( m_KurtosisSnap.m_K, 0, 'f', 10 ) ; QApplication::clipboard()->setText(clipboard, QClipboard::Clipboard ); } else { QString clipboard( "f \t D \t D* \n" ); clipboard = clipboard.append( "%L1 \t %L2 \t %L3" ) .arg( m_IvimSnap.currentF, 0, 'f', 10 ) .arg( m_IvimSnap.currentD, 0, 'f', 10 ) .arg( m_IvimSnap.currentDStar, 0, 'f', 10 ) ; QApplication::clipboard()->setText(clipboard, QClipboard::Clipboard ); } } void QmitkIVIMView::Activated() { m_Active = true; } void QmitkIVIMView::Deactivated() { m_Active = false; } void QmitkIVIMView::Visible() { m_Visible = true; if (this->GetRenderWindowPart() && !m_ListenerActive) { m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); m_ListenerActive = true; } } void QmitkIVIMView::Hidden() { m_Visible = false; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMViewControls.ui index 175f5c7..cf4a641 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.ivim/src/internal/QmitkIVIMViewControls.ui @@ -1,1213 +1,1213 @@ QmitkIVIMViewControls 0 0 424 689 0 0 QmitkTemplate QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 9 9 9 9 9 QFrame::NoFrame QFrame::Raised 0 0 0 0 warning display Qt::RichText true 0 0 16 16 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 16 16 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 16 16 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 0 0 16777215 16777215 Generate Parameter Maps QFrame::StyledPanel QFrame::Raised 0 0 0 0 0 0 0 IVIM Parameters 9 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 51 16777215 200 Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter 80 16777215 D* 100 60 Qt::Horizontal QFrame::NoFrame QFrame::Plain 0 0 0 0 0 130 0 130 16777215 Signal threshold: 100 0 Qt::Horizontal 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 0 30 16777215 TextLabel Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter 15 16777215 Calculate threshold from histogram * QFrame::NoFrame QFrame::Raised 0 0 0 0 0 130 0 130 16777215 Ignore b< (first fit) 250 34 Qt::Horizontal 51 16777215 46.5 Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter QFrame::NoFrame QFrame::Raised 0 0 0 0 0 80 16777215 #iterations: 100 10 Qt::Horizontal 30 16777215 TextLabel Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter QFrame::NoFrame QFrame::Raised 0 0 0 0 0 80 16777215 lambda: 1000 10 Qt::Horizontal 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 0 30 16777215 TextLabel Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter 15 16777215 Calculate threshold from histogram * QFrame::NoFrame QFrame::Raised 0 0 0 0 Fit method: - 0 + 2 1. Jointly fit D, f and D* 2. Fit D & f with fixed D* value 3. Fit D & f (high b), then fit D* 4. Linearly fit D & f (high b), then fit D* QFrame::NoFrame QFrame::Raised 0 0 0 0 80 0 Output Maps f true D false D* false Kurtosis QFrame::StyledPanel QFrame::Raised 2 2 2 2 2 Smoothing sigma Select Fit Type Omit b=0 Measurements false 80 0 Output Maps Force the fitting of K to remain within the given boundaries Boundaries for K Select if the data is fitted directly (straight) or the logarithmic equation is used Straight Fit Logarithmic Fit 2 QLayout::SetMaximumSize D false K true Signa for gaussian smoothing applied prior to map computation 0.000000000000000 5.000000000000000 0.100000000000000 On Input Data 6 6 6 6 Optional ROI image ROI: DWI to analyze Raw DWI: Qt::Vertical 20 40 Copy the signal values in the current voxel to the clipboard. Save plot as image QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Copy the signal values in the current voxel to the clipboard. Signal values to clipboard Copy the parameters of the curve fit in the current voxel to the clipboard. Fit-parameters to clipboard Intra Voxel Incoherent Motion Estimation 6 9 6 6 Qt::AlignCenter QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
ctkRangeWidget QWidget
ctkRangeWidget.h
1
QmitkChartWidget QWidget
QmitkChartWidget.h
1