diff --git a/BlueBerry/Bundles/org.blueberry.ui.qt/files.cmake.user b/BlueBerry/Bundles/org.blueberry.ui.qt/files.cmake.user deleted file mode 100644 index 5280950531..0000000000 --- a/BlueBerry/Bundles/org.blueberry.ui.qt/files.cmake.user +++ /dev/null @@ -1,63 +0,0 @@ - - - - - - ProjectExplorer.Project.ActiveTarget - -1 - - - ProjectExplorer.Project.EditorSettings - - true - false - true - - Cpp - - CppGlobal - - - - QmlJS - - QmlJSGlobal - - - 2 - UTF-8 - false - 4 - false - true - 1 - true - 0 - true - 0 - 8 - true - 1 - true - true - true - true - - - - ProjectExplorer.Project.PluginSettings - - - - ProjectExplorer.Project.TargetCount - 0 - - - ProjectExplorer.Project.Updater.EnvironmentId - {3dea2ab9-fb51-4719-a902-a8416c2a1947} - - - ProjectExplorer.Project.Updater.FileVersion - 15 - - diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp index 6c437a400b..9f44d41d0a 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp @@ -1,556 +1,556 @@ #include "itkTractsToVectorImageFilter.h" // VTK #include #include #include // ITK #include #include // misc #define _USE_MATH_DEFINES #include #include namespace itk{ static bool CompareVectorLengths(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2) { return (v1.magnitude()>v2.magnitude()); } template< class PixelType > TractsToVectorImageFilter< PixelType >::TractsToVectorImageFilter(): m_AngularThreshold(0.7), m_Epsilon(0.999), m_MaskImage(NULL), m_NormalizeVectors(false), m_UseWorkingCopy(true), m_MaxNumDirections(3), m_SizeThreshold(0.2), m_NumDirectionsImage(NULL) { this->SetNumberOfRequiredOutputs(1); } template< class PixelType > TractsToVectorImageFilter< PixelType >::~TractsToVectorImageFilter() { } template< class PixelType > vnl_vector_fixed TractsToVectorImageFilter< PixelType >::GetVnlVector(double point[]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } template< class PixelType > itk::Point TractsToVectorImageFilter< PixelType >::GetItkPoint(double point[]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class PixelType > void TractsToVectorImageFilter< PixelType >::GenerateData() { mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry(); // calculate new image parameters itk::Vector spacing; itk::Point origin; itk::Matrix direction; ImageRegion<3> imageRegion; if (!m_MaskImage.IsNull()) { spacing = m_MaskImage->GetSpacing(); imageRegion = m_MaskImage->GetLargestPossibleRegion(); origin = m_MaskImage->GetOrigin(); direction = m_MaskImage->GetDirection(); } else { spacing = geometry->GetSpacing(); origin = geometry->GetOrigin(); mitk::BaseGeometry::BoundsArrayType bounds = geometry->GetBounds(); origin[0] += bounds.GetElement(0); origin[1] += bounds.GetElement(2); origin[2] += bounds.GetElement(4); for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]; imageRegion.SetSize(0, geometry->GetExtent(0)); imageRegion.SetSize(1, geometry->GetExtent(1)); imageRegion.SetSize(2, geometry->GetExtent(2)); m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( spacing ); m_MaskImage->SetOrigin( origin ); m_MaskImage->SetDirection( direction ); m_MaskImage->SetRegions( imageRegion ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } OutputImageType::RegionType::SizeType outImageSize = imageRegion.GetSize(); m_OutImageSpacing = m_MaskImage->GetSpacing(); m_ClusteredDirectionsContainer = ContainerType::New(); // initialize num directions image m_NumDirectionsImage = ItkUcharImgType::New(); m_NumDirectionsImage->SetSpacing( spacing ); m_NumDirectionsImage->SetOrigin( origin ); m_NumDirectionsImage->SetDirection( direction ); m_NumDirectionsImage->SetRegions( imageRegion ); m_NumDirectionsImage->Allocate(); m_NumDirectionsImage->FillBuffer(0); // initialize direction images m_DirectionImageContainer = DirectionImageContainerType::New(); // resample fiber bundle double minSpacing = 1; if(m_OutImageSpacing[0]GetDeepCopy(); // resample fiber bundle for sufficient voxel coverage m_FiberBundle->ResampleFibers(minSpacing/10); // iterate over all fibers vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); int numFibers = m_FiberBundle->GetNumFibers(); m_DirectionsContainer = ContainerType::New(); MITK_INFO << "Generating directions from tractogram"; boost::progress_display disp(numFibers); for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) continue; vnl_vector_fixed dir; itk::Point worldPos; vnl_vector v; for( int j=0; jGetPoint(j); worldPos = GetItkPoint(temp); itk::Index<3> index; m_MaskImage->TransformPhysicalPointToIndex(worldPos, index); if (!m_MaskImage->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0) continue; // get fiber tangent direction at this position v = GetVnlVector(temp); dir = GetVnlVector(points->GetPoint(j+1))-v; if (dir.is_zero()) continue; dir.normalize(); // add direction to container unsigned int idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]); DirectionContainerType::Pointer dirCont; if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->push_back(dir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->push_back(dir); } else { dirCont = DirectionContainerType::New(); dirCont->push_back(dir); m_DirectionsContainer->InsertElement(idx, dirCont); } } } vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); itk::ImageRegionIterator dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion()); MITK_INFO << "Clustering directions"; boost::progress_display disp2(outImageSize[0]*outImageSize[1]*outImageSize[2]); while(!dirIt.IsAtEnd()) { ++disp2; OutputImageType::IndexType index = dirIt.GetIndex(); int idx = index[0]+(index[1]+index[2]*outImageSize[1])*outImageSize[0]; if (!m_DirectionsContainer->IndexExists(idx)) { ++dirIt; continue; } DirectionContainerType::Pointer dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull() || dirCont->empty()) { ++dirIt; continue; } std::vector< double > lengths; lengths.resize(dirCont->size(), 1); // all peaks have size 1 DirectionContainerType::Pointer directions = FastClustering(dirCont, lengths); std::sort( directions->begin(), directions->end(), CompareVectorLengths ); unsigned int numDir = directions->size(); if (m_MaxNumDirections>0 && numDir>m_MaxNumDirections) numDir = m_MaxNumDirections; int count = 0; for (unsigned int i=0; i container = vtkSmartPointer::New(); itk::ContinuousIndex center; center[0] = index[0]; center[1] = index[1]; center[2] = index[2]; itk::Point worldCenter; m_MaskImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); DirectionType dir = directions->at(i); if (dir.magnitude()size()) { ItkDirectionImageType::Pointer directionImage = ItkDirectionImageType::New(); directionImage->SetSpacing( spacing ); directionImage->SetOrigin( origin ); directionImage->SetDirection( direction ); directionImage->SetRegions( imageRegion ); directionImage->Allocate(); Vector< float, 3 > nullVec; nullVec.Fill(0.0); directionImage->FillBuffer(nullVec); m_DirectionImageContainer->InsertElement(i, directionImage); } // set direction image pixel ItkDirectionImageType::Pointer directionImage = m_DirectionImageContainer->GetElement(i); Vector< float, 3 > pixel; pixel.SetElement(0, dir[0]); pixel.SetElement(1, dir[1]); pixel.SetElement(2, dir[2]); directionImage->SetPixel(index, pixel); // add direction to vector field (with spacing compensation) itk::Point worldStart; worldStart[0] = worldCenter[0]-dir[0]/2*minSpacing; worldStart[1] = worldCenter[1]-dir[1]/2*minSpacing; worldStart[2] = worldCenter[2]-dir[2]/2*minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+dir[0]/2*minSpacing; worldEnd[1] = worldCenter[1]+dir[1]/2*minSpacing; worldEnd[2] = worldCenter[2]+dir[2]/2*minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); } dirIt.Set(count); ++dirIt; } vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); directionsPolyData->SetPoints(m_VtkPoints); directionsPolyData->SetLines(m_VtkCellArray); m_OutputFiberBundle = mitk::FiberBundleX::New(directionsPolyData); } template< class PixelType > -typename TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths) +TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::FastClustering(DirectionContainerType::Pointer inDirs, std::vector< double > lengths) { DirectionContainerType::Pointer outDirs = DirectionContainerType::New(); if (inDirs->size()<2) return inDirs; DirectionType oldMean, currentMean; std::vector< int > touched; // initialize touched.resize(inDirs->size(), 0); bool free = true; currentMean = inDirs->at(0); // initialize first seed currentMean.normalize(); double length = lengths.at(0); touched[0] = 1; std::vector< double > newLengths; bool meanChanged = false; double max = 0; while (free) { oldMean.fill(0.0); // start mean-shift clustering double angle = 0; while (fabs(dot_product(currentMean, oldMean))<0.99) { oldMean = currentMean; currentMean.fill(0.0); for (unsigned int i=0; isize(); i++) { angle = dot_product(oldMean, inDirs->at(i)); if (angle>=m_AngularThreshold) { currentMean += inDirs->at(i); if (meanChanged) length += lengths.at(i); touched[i] = 1; meanChanged = true; } else if (-angle>=m_AngularThreshold) { currentMean -= inDirs->at(i); if (meanChanged) length += lengths.at(i); touched[i] = 1; meanChanged = true; } } if(!meanChanged) currentMean = oldMean; else currentMean.normalize(); } // found stable mean outDirs->push_back(currentMean); newLengths.push_back(length); if (length>max) max = length; // find next unused seed free = false; for (unsigned int i=0; iat(i); free = true; meanChanged = false; length = lengths.at(i); touched[i] = 1; break; } } if (inDirs->size()==outDirs->size()) { if (!m_NormalizeVectors && max>0) for (unsigned int i=0; isize(); i++) outDirs->SetElement(i, outDirs->at(i)*newLengths.at(i)/max); return outDirs; } else return FastClustering(outDirs, newLengths); } //template< class PixelType > //std::vector< DirectionType > TractsToVectorImageFilter< PixelType >::Clustering(std::vector< DirectionType >& inDirs) //{ // std::vector< DirectionType > outDirs; // if (inDirs.empty()) // return outDirs; // DirectionType oldMean, currentMean, workingMean; // std::vector< DirectionType > normalizedDirs; // std::vector< int > touched; // for (std::size_t i=0; i0.0001) // { // counter = 0; // oldMean = currentMean; // workingMean = oldMean; // workingMean.normalize(); // currentMean.fill(0.0); // for (std::size_t i=0; i=m_AngularThreshold) // { // currentMean += inDirs[i]; // counter++; // } // else if (-angle>=m_AngularThreshold) // { // currentMean -= inDirs[i]; // counter++; // } // } // } // // found stable mean // if (counter>0) // { // bool add = true; // DirectionType normMean = currentMean; // normMean.normalize(); // for (std::size_t i=0; i0) // { // if (mag>max) // max = mag; // outDirs.push_back(currentMean); // } // } // } // } // if (m_NormalizeVectors) // for (std::size_t i=0; i0) // for (std::size_t i=0; i //TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::MeanShiftClustering(DirectionContainerType::Pointer dirCont) //{ // DirectionContainerType::Pointer container = DirectionContainerType::New(); // double max = 0; // for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it) // { // vnl_vector_fixed mean = ClusterStep(dirCont, it.Value()); // if (mean.is_zero()) // continue; // bool addMean = true; // for (DirectionContainerType::ConstIterator it2 = container->Begin(); it2!=container->End(); ++it2) // { // vnl_vector_fixed dir = it2.Value(); // double angle = fabs(dot_product(mean, dir)/(mean.magnitude()*dir.magnitude())); // if (angle>=m_Epsilon) // { // addMean = false; // break; // } // } // if (addMean) // { // if (m_NormalizeVectors) // mean.normalize(); // else if (mean.magnitude()>max) // max = mean.magnitude(); // container->InsertElement(container->Size(), mean); // } // } // // max normalize voxel directions // if (max>0 && !m_NormalizeVectors) // for (std::size_t i=0; iSize(); i++) // container->ElementAt(i) /= max; // if (container->Size()Size()) // return MeanShiftClustering(container); // else // return container; //} //template< class PixelType > //vnl_vector_fixed TractsToVectorImageFilter< PixelType >::ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed currentMean) //{ // vnl_vector_fixed newMean; newMean.fill(0); // for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it) // { // vnl_vector_fixed dir = it.Value(); // double angle = dot_product(currentMean, dir)/(currentMean.magnitude()*dir.magnitude()); // if (angle>=m_AngularThreshold) // newMean += dir; // else if (-angle>=m_AngularThreshold) // newMean -= dir; // } // if (fabs(dot_product(currentMean, newMean)/(currentMean.magnitude()*newMean.magnitude()))>=m_Epsilon || newMean.is_zero()) // return newMean; // else // return ClusterStep(dirCont, newMean); //} }