diff --git a/Modules/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx b/Modules/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx index 885a635..9e39a9d 100644 --- a/Modules/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx +++ b/Modules/DiffusionCore/Algorithms/Reconstruction/itkOrientationDistributionFunction.txx @@ -1,1366 +1,1366 @@ /*=================================================================== 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 _itkOrientationDistributionFunction_txx #define _itkOrientationDistributionFunction_txx #include #include #include #include #include #include #include "itkPointShell.h" #ifdef _MSC_VER #if _MSC_VER <= 1700 #define fmin(a,b) ((a<=b)?(a):(b)) #define fmax(a,b) ((a>=b)?(a):(b)) #define isnan(c) (c!=c) #endif #endif #include #include #include #include #include #include namespace itk { template vtkPolyData* itk::OrientationDistributionFunction::m_BaseMesh = nullptr; template double itk::OrientationDistributionFunction::m_MaxChordLength = -1.0; template vnl_matrix_fixed* itk::OrientationDistributionFunction::m_Directions = itk::PointShell >::DistributePointShell(); template std::vector< std::vector* >* itk::OrientationDistributionFunction::m_NeighborIdxs = nullptr; template std::vector< std::vector* >* itk::OrientationDistributionFunction::m_AngularRangeIdxs = nullptr; template std::vector* itk::OrientationDistributionFunction::m_HalfSphereIdxs = nullptr; template itk::SimpleFastMutexLock itk::OrientationDistributionFunction::m_MutexBaseMesh; template itk::SimpleFastMutexLock itk::OrientationDistributionFunction::m_MutexHalfSphereIdxs; template itk::SimpleFastMutexLock itk::OrientationDistributionFunction::m_MutexNeighbors; template itk::SimpleFastMutexLock itk::OrientationDistributionFunction::m_MutexAngularRange; /** * Assignment Operator */ template OrientationDistributionFunction& OrientationDistributionFunction ::operator= (const Self& r) { BaseArray::operator=(r); return *this; } /** * Assignment Operator from a scalar constant */ template OrientationDistributionFunction& OrientationDistributionFunction ::operator= (const ComponentType & r) { BaseArray::operator=(&r); return *this; } /** * Assigment from a plain array */ template OrientationDistributionFunction& OrientationDistributionFunction ::operator= (const ComponentArrayType r ) { BaseArray::operator=(r); return *this; } /** * Returns a temporary copy of a vector */ template OrientationDistributionFunction OrientationDistributionFunction ::operator+(const Self & r) const { Self result; for( unsigned int i=0; i OrientationDistributionFunction OrientationDistributionFunction ::operator-(const Self & r) const { Self result; for( unsigned int i=0; i const OrientationDistributionFunction & OrientationDistributionFunction ::operator+=(const Self & r) { for( unsigned int i=0; i const OrientationDistributionFunction & OrientationDistributionFunction ::operator-=(const Self & r) { for( unsigned int i=0; i const OrientationDistributionFunction & OrientationDistributionFunction ::operator*=(const RealValueType & r) { for( unsigned int i=0; i const OrientationDistributionFunction & OrientationDistributionFunction ::operator/=(const RealValueType & r) { for( unsigned int i=0; i OrientationDistributionFunction OrientationDistributionFunction ::operator*(const RealValueType & r) const { Self result; for( unsigned int i=0; i OrientationDistributionFunction OrientationDistributionFunction ::operator/(const RealValueType & r) const { Self result; for( unsigned int i=0; i const typename OrientationDistributionFunction::ValueType & OrientationDistributionFunction ::operator()(unsigned int row, unsigned int col) const { unsigned int k; if( row < col ) { k = row * InternalDimension + col - row * ( row + 1 ) / 2; } else { k = col * InternalDimension + row - col * ( col + 1 ) / 2; } if( k >= InternalDimension ) { k = 0; } return (*this)[k]; } /** * Matrix notation access to elements */ template typename OrientationDistributionFunction::ValueType & OrientationDistributionFunction ::operator()(unsigned int row, unsigned int col) { unsigned int k; if( row < col ) { k = row * InternalDimension + col - row * ( row + 1 ) / 2; } else { k = col * InternalDimension + row - col * ( col + 1 ) / 2; } if( k >= InternalDimension ) { k = 0; } return (*this)[k]; } /** * Set the Tensor to an Identity. * Set ones in the diagonal and zeroes every where else. */ template void OrientationDistributionFunction ::SetIsotropic() { this->Fill(NumericTraits< T >::One / NOdfDirections); } /** * InitFromTensor() */ template void OrientationDistributionFunction ::InitFromTensor(itk::DiffusionTensor3D tensor) { m_EigenAnalysisCalculated = false; for(unsigned int i=0; i void OrientationDistributionFunction:: InitFromEllipsoid( itk::DiffusionTensor3D tensor ) { m_EigenAnalysisCalculated = false; FixedArray nulltensor; nulltensor.Fill(0.0); if( tensor == nulltensor ) { for ( unsigned int it=0; it < NOdfDirections; ++it ){ (*this)[it] = (T)0; } MITK_DEBUG << "OrientationDistributionFunction<" << typeid(T).name() << ", " << NOdfDirections << ">::InitFromEllipsoid(" << typeid(tensor).name() << ") encountered a nulltensor as dti input point and ignorend it."; return; } tensor.ComputeEigenAnalysis( m_EigenValues, this->m_EigenVectors ); // gives normalized eigenvectors as lines i.e. rows. m_EigenAnalysisCalculated = true; double a = m_EigenValues[0]; // those eigenvalues are the 3 |axes of the ellipsoid|, double b = m_EigenValues[1]; // ComputeEigenAnalysis gives eigenValues in ascending < order per default, double c = m_EigenValues[2]; // therefor the third eigenVector is the main direction of diffusion. if( a <= 0.0 || b <= 0.0 || c <= 0.0 ) { for ( unsigned int it=0; it < NOdfDirections; ++it ){ (*this)[it] = (T)0; } MITK_DEBUG << "OrientationDistributionFunction<" << typeid(T).name() << ", " << NOdfDirections << ">::InitFromEllipsoid(" << typeid(tensor).name() << ") encountered an eigenvalue <= 0 and ignored this input point."; return; } // check magnitude and scale towards 1 to minimize numerical condition kappa: #ifdef _MSC_VER #if _MSC_VER <= 1700 int exponent_a = floor(std::log(a)/std::log(2)); int exponent_b = floor(std::log(b)/std::log(2)); int exponent_c = floor(std::log(c)/std::log(2)); #else int exponent_a = std::ilogb(a); int exponent_b = std::ilogb(b); int exponent_c = std::ilogb(c); #endif #else int exponent_a = std::ilogb(a); int exponent_b = std::ilogb(b); int exponent_c = std::ilogb(c); #endif T min_exponent= fmin(exponent_a, fmin(exponent_b, exponent_c) ); T max_exponent= fmax(exponent_c, fmax(exponent_b, exponent_a) ); int scale_exponent = floor(0.5 * (min_exponent + max_exponent)); double scaling = pow(2, scale_exponent); a= a/scaling; b= b/scaling; c= c/scaling; vnl_matrix_fixed eigenBase; // for change of base system. for (int row = 0 ; row < 3; ++row) // Transposing because ComputeEigenAnalysis(,) gave us _row_ vectors. { for (int col = 0; col < 3; ++col) { eigenBase(row, col) = this->m_EigenVectors(col, row); } } eigenBase= vnl_inverse(eigenBase); // Assuming canonical orthonormal system x=[1;0;0];y=[0;1;0];z=[0;0;1] for original DT. eigenBase.assert_finite(); #ifndef NDEBUG double kappa=1.0; { // calculate numerical condition kappa= ||f(x)-f(x~)|| approximately: double gxaa = pow( a, -2.0); double gybb = pow( b, -2.0); double gzcc = pow( c, -2.0); kappa = sqrt( pow( a, 2.0)+ pow( b, 2.0)+ pow( c, 2.0) + 1 ) / (gxaa + gybb + gzcc) * sqrt( pow( a, -6.0)+ pow( b, -6.0)+ pow( c, -6.0) + pow( a, -4.0)+ pow( b, -4.0)+ pow( c, -4.0) ); MITK_DEBUG <<"kappa= "<< kappa << ", eigenvalues= [" << a <<", "<< b <<", "<< c <<"], eigenbase= [" < g( (*m_Directions)(0,i), (*m_Directions)(1,i), (*m_Directions)(2,i) ); g = eigenBase*g; // passive change of base system of g. g = g.normalize(); // unit vectors necessary. (*this)[i] = scaling / sqrt( (g[0]/a)*(g[0]/a) + (g[1]/b)*(g[1]/b) + (g[2]/c)*(g[2]/c) ); #ifndef NDEBUG { // boundary check for numerical stability, assuming sigma=6, ||f(x~)-f~(x~)|| <= eps*kappa*sigma. T min_ev= fmin(a, fmin(b, c) ); T max_ev= fmax(c, fmax(b, a) ); double eps= std::numeric_limits::epsilon(); assert( scaling*min_ev <= ((*this)[i] + eps*kappa*6.0) ); // we should be between smallest and assert( (*this)[i] <= (scaling*max_ev + eps*kappa*6.0) ); // biggest eigenvalue. } #endif if ( (*this)[i] < T(0) || (*this)[i] > T(1) || std::isnan((*this)[i]) ) // P∈[0;1] sanity check. { // C: NaN != NaN, C++11: isnan((*this)[i]). MITK_DEBUG << "OrientationDistributionFunction<" << typeid(T).name() << ", " << NOdfDirections << ">::InitFromEllipsoid(" << typeid(tensor).name() << ") encountered a probability value out of range [0;1] and set it to zero: (*this)[" << i <<"]= " << (*this)[i]; (*this)[i] = T(0); } } } /** * L2-Normalization */ template void OrientationDistributionFunction ::L2Normalize() { T sum = 0; for( unsigned int i=0; i void OrientationDistributionFunction ::Normalize() { T sum = 0; for( unsigned int i=0; i0) { for( unsigned int i=0; i OrientationDistributionFunction OrientationDistributionFunction ::MinMaxNormalize() const { T max = NumericTraits::NonpositiveMin(); T min = NumericTraits::max(); for( unsigned int i=0; i max ? (*this)[i] : max; min = (*this)[i] < min ? (*this)[i] : min; } Self retval; for( unsigned int i=0; i OrientationDistributionFunction OrientationDistributionFunction ::MaxNormalize() const { T max = NumericTraits::NonpositiveMin(); for( unsigned int i=0; i max ? (*this)[i] : max; } Self retval; for( unsigned int i=0; i T OrientationDistributionFunction ::GetMaxValue() const { T max = NumericTraits::NonpositiveMin(); for( unsigned int i=0; i= max ) { max = (*this)[i]; } } return max; } template T OrientationDistributionFunction ::GetMinValue() const { T min = NumericTraits::max(); for( unsigned int i=0; i= min ) { min = (*this)[i]; } } return min; } template T OrientationDistributionFunction ::GetMeanValue() const { T sum = 0; for( unsigned int i=0; i double OrientationDistributionFunction ::GetMaxChordLength() { if(m_MaxChordLength<0.0) { ComputeBaseMesh(); double max_dist = -1; vtkPoints* points = m_BaseMesh->GetPoints(); for(int i=0; iGetPoint(i,p); std::vector neighbors = GetNeighbors(i); for(std::size_t j=0; jGetPoint(neighbors[j],n); double d = sqrt( (p[0]-n[0])*(p[0]-n[0]) + (p[1]-n[1])*(p[1]-n[1]) + (p[2]-n[2])*(p[2]-n[2])); max_dist = d>max_dist ? d : max_dist; } } m_MaxChordLength = max_dist; } return m_MaxChordLength; } template void OrientationDistributionFunction ::ComputeBaseMesh() { m_MutexBaseMesh.Lock(); if(m_BaseMesh == nullptr) { vtkPoints* points = vtkPoints::New(); for(unsigned int j=0; jInsertNextPoint(az,elev,r); } vtkPolyData* polydata = vtkPolyData::New(); polydata->SetPoints( points ); vtkDelaunay2D *delaunay = vtkDelaunay2D::New(); delaunay->SetInputData( polydata ); delaunay->Update(); vtkCellArray* vtkpolys = delaunay->GetOutput()->GetPolys(); vtkCellArray* vtknewpolys = vtkCellArray::New(); - vtkIdType npts; vtkIdType *pts; + vtkIdType npts; vtkIdType const *pts; while(vtkpolys->GetNextCell(npts,pts)) { bool insert = true; for(int i=0; iGetPoint(pts[i]); double az = tmpPoint[0]; double elev = tmpPoint[1]; if((std::abs(az)>itk::Math::pi-0.5) || (std::abs(elev)>itk::Math::pi/2-0.5)) insert = false; } if(insert) vtknewpolys->InsertNextCell(npts, pts); } vtkPoints* points2 = vtkPoints::New(); for(unsigned int j=0; jInsertNextPoint(az,elev,r); } vtkPolyData* polydata2 = vtkPolyData::New(); polydata2->SetPoints( points2 ); vtkDelaunay2D *delaunay2 = vtkDelaunay2D::New(); delaunay2->SetInputData( polydata2 ); delaunay2->Update(); vtkpolys = delaunay2->GetOutput()->GetPolys(); while(vtkpolys->GetNextCell(npts,pts)) { bool insert = true; for(int i=0; iGetPoint(pts[i]); double az = tmpPoint[0]; double elev = tmpPoint[1]; if((std::abs(az)>itk::Math::pi-0.5) || (std::abs(elev)>itk::Math::pi/2-0.5)) insert = false; } if(insert) vtknewpolys->InsertNextCell(npts, pts); } polydata->SetPolys(vtknewpolys); for (unsigned int p = 0; p < NOdfDirections; p++) { points->SetPoint(p,m_Directions->get_column(p).data_block()); } polydata->SetPoints( points ); m_BaseMesh = polydata; } m_MutexBaseMesh.Unlock(); } /** * Extract the index of the principal diffusion direction */ template int OrientationDistributionFunction::GetPrincipalDiffusionDirectionIndex() const { T max = NumericTraits::NonpositiveMin(); int maxidx = -1; for( unsigned int i=0; i= max ) { max = (*this)[i]; maxidx = i; } } return maxidx; } /** * Extract the principal diffusion direction */ template vnl_vector_fixed OrientationDistributionFunction::GetPrincipalDiffusionDirection() const { if (m_EigenAnalysisCalculated) { vnl_vector_fixed vec; vec[0] = this->m_EigenVectors(2,0); vec[1] = this->m_EigenVectors(2,1); vec[2] = this->m_EigenVectors(2,2); vec.normalize(); return vec; } else { int idx = GetPrincipalDiffusionDirectionIndex(); if (idx>0 && idx<(int)NOdfDirections) return OrientationDistributionFunction::GetDirection(idx); vnl_vector_fixed vec; vec.fill(0); return vec; } } template std::vector OrientationDistributionFunction ::GetNeighbors(int idx) { ComputeBaseMesh(); m_MutexNeighbors.Lock(); if(m_NeighborIdxs == nullptr) { m_NeighborIdxs = new std::vector< std::vector* >(); vtkCellArray* polys = m_BaseMesh->GetPolys(); for(unsigned int i=0; i(); polys->InitTraversal(); - vtkIdType npts; vtkIdType *pts; + vtkIdType npts; vtkIdType const *pts; while(polys->GetNextCell(npts,pts)) { if( pts[0] == i ) { idxs->push_back(pts[1]); idxs->push_back(pts[2]); } else if( pts[1] == i ) { idxs->push_back(pts[0]); idxs->push_back(pts[2]); } else if( pts[2] == i ) { idxs->push_back(pts[0]); idxs->push_back(pts[1]); } } std::sort(idxs->begin(), idxs->end()); std::vector< int >::iterator endLocation; endLocation = std::unique( idxs->begin(), idxs->end() ); idxs->erase(endLocation, idxs->end()); m_NeighborIdxs->push_back(idxs); } } m_MutexNeighbors.Unlock(); return *m_NeighborIdxs->at(idx); } /** * Extract the n-th diffusion direction */ template int OrientationDistributionFunction ::GetNthDiffusionDirection(int n, vnl_vector_fixed rndVec) const { if( n == 0 ) return GetPrincipalDiffusionDirectionIndex(); m_MutexHalfSphereIdxs.Lock(); if( !m_HalfSphereIdxs ) { m_HalfSphereIdxs = new std::vector(); for( unsigned int i=0; iget_column(i),rndVec) > 0.0) { m_HalfSphereIdxs->push_back(i); } } } m_MutexHalfSphereIdxs.Unlock(); // collect indices of directions // that are local maxima std::vector localMaxima; std::vector::iterator it; for( it=m_HalfSphereIdxs->begin(); it!=m_HalfSphereIdxs->end(); it++) { std::vector nbs = GetNeighbors(*it); std::vector::iterator it2; bool max = true; for(it2 = nbs.begin(); it2 != nbs.end(); it2++) { if((*this)[*it2] > (*this)[*it]) { max = false; break; } } if(max) localMaxima.push_back(*it); } // delete n highest local maxima from list // and return remaining highest int maxidx = -1; for( int i=0; i<=n; i++ ) { maxidx = -1; T max = NumericTraits::NonpositiveMin(); for(it = localMaxima.begin(); it != localMaxima.end(); it++) { if((*this)[*it]>max) { max = (*this)[*it]; maxidx = *it; } } it = find(localMaxima.begin(), localMaxima.end(), maxidx); if(it!=localMaxima.end()) localMaxima.erase(it); } return maxidx; } template < typename TComponent, unsigned int NOdfDirections > vnl_vector_fixed itk::OrientationDistributionFunction ::GetDirection( int i ) { return m_Directions->get_column(i); } /** * Interpolate a position between sampled directions */ template T OrientationDistributionFunction ::GetInterpolatedComponent(vnl_vector_fixed dir, InterpolationMethods method) const { ComputeBaseMesh(); double retval = -1.0; switch(method) { case ODF_NEAREST_NEIGHBOR_INTERP: { vtkPoints* points = m_BaseMesh->GetPoints(); double current_min = NumericTraits::max(); int current_min_idx = -1; for(int i=0; i P(points->GetPoint(i)); double dist = (dir-P).two_norm(); current_min_idx = distGetNthComponent(current_min_idx); break; } case ODF_TRILINEAR_BARYCENTRIC_INTERP: { double maxChordLength = GetMaxChordLength(); vtkCellArray* polys = m_BaseMesh->GetPolys(); vtkPoints* points = m_BaseMesh->GetPoints(); - vtkIdType npts; vtkIdType *pts; + vtkIdType npts; vtkIdType const *pts; double current_min = NumericTraits::max(); polys->InitTraversal(); while(polys->GetNextCell(npts,pts)) { vnl_vector_fixed A(points->GetPoint(pts[0])); vnl_vector_fixed B(points->GetPoint(pts[1])); vnl_vector_fixed C(points->GetPoint(pts[2])); vnl_vector_fixed d1; d1.put(0,(dir-A).two_norm()); d1.put(1,(dir-B).two_norm()); d1.put(2,(dir-C).two_norm()); double maxval = d1.max_value(); if(maxval>maxChordLength) { continue; } // Compute vectors vnl_vector_fixed v0 = C - A; vnl_vector_fixed v1 = B - A; // Project direction to plane ABC vnl_vector_fixed v6 = dir; vnl_vector_fixed cross = vnl_cross_3d(v0, v1); cross = cross.normalize(); vtkPlane::ProjectPoint(v6.data_block(),A.data_block(),cross.data_block(),v6.data_block()); v6 = v6-A; // Calculate barycentric coords vnl_matrix_fixed mat; mat.set_column(0, v0); mat.set_column(1, v1); vnl_matrix_inverse inv(mat); vnl_matrix_fixed inver = inv.pinverse(); vnl_vector uv = inv.pinverse()*v6; // Check if point is in triangle double eps = 0.01; if( (uv(0) >= 0-eps) && (uv(1) >= 0-eps) && (uv(0) + uv(1) <= 1+eps) ) { // check if minimum angle is the max so far if(d1.two_norm() < current_min) { current_min = d1.two_norm(); vnl_vector barycentricCoords(3); barycentricCoords[2] = uv[0]<0 ? 0 : (uv[0]>1?1:uv[0]); barycentricCoords[1] = uv[1]<0 ? 0 : (uv[1]>1?1:uv[1]); barycentricCoords[0] = 1-(barycentricCoords[1]+barycentricCoords[2]); retval = barycentricCoords[0]*this->GetNthComponent(pts[0]) + barycentricCoords[1]*this->GetNthComponent(pts[1]) + barycentricCoords[2]*this->GetNthComponent(pts[2]); } } } break; } case ODF_SPHERICAL_GAUSSIAN_BASIS_FUNCTIONS: { double maxChordLength = GetMaxChordLength(); double sigma = asin(maxChordLength/2); // this is the contribution of each kernel to each sampling point on the // equator vnl_vector contrib; contrib.set_size(NOdfDirections); vtkPoints* points = m_BaseMesh->GetPoints(); double sum = 0; for(int i=0; i P(points->GetPoint(i)); double stv = dir[0]*P[0] + dir[1]*P[1] + dir[2]*P[2]; stv = (stv<-1.0) ? -1.0 : ( (stv>1.0) ? 1.0 : stv); double x = acos(stv); contrib[i] = (1.0/(sigma*sqrt(2.0*itk::Math::pi))) *exp((-x*x)/(2*sigma*sigma)); sum += contrib[i]; } retval = 0; for(int i=0; iGetNthComponent(i); } break; } } if(retval==-1) { std::cout << "Interpolation failed" << std::endl; return 0; } return retval; } /** * Calculate Generalized Fractional Anisotropy */ template T OrientationDistributionFunction ::GetGeneralizedFractionalAnisotropy() const { double mean = 0; double std = 0; double rms = 0; for( unsigned int i=0; i T itk::OrientationDistributionFunction ::GetGeneralizedGFA( int k, int p ) const { double mean = 0; double std = 0; double rms = 0; double max = NumericTraits::NonpositiveMin(); for( unsigned int i=0; i max ? val : max; } max = pow(max,(double)p); mean /= N; for( unsigned int i=0; i0) { rms += pow(val,(double)(p*k)); } } std /= N - 1; std = sqrt(std); if(k>0) { rms /= N; rms = pow(rms,(double)(1.0/k)); } else if(k<0) // lim k->inf gives us the maximum { rms = max; } else // k==0 undefined, we define zeros root from 1 as 1 { rms = 1; } if(rms == 0) { return 0; } else { return (T)(std/rms); } } /** * Calculate Nematic Order Parameter */ template < typename T, unsigned int N > T itk::OrientationDistributionFunction ::GetNematicOrderParameter() const { // not yet implemented return 0; } /** * Calculate StdDev by MaxValue */ template < typename T, unsigned int N > T itk::OrientationDistributionFunction ::GetStdDevByMaxValue() const { double mean = 0; double std = 0; T max = NumericTraits::NonpositiveMin(); for( unsigned int i=0; i max ? (*this)[i] : max; } mean /= InternalDimension; for( unsigned int i=0; i T itk::OrientationDistributionFunction ::GetPrincipleCurvature(double alphaMinDegree, double alphaMaxDegree, int invert) const { // following loop only performed once // (computing indices of each angular range) m_MutexAngularRange.Lock(); if(m_AngularRangeIdxs == nullptr) { m_AngularRangeIdxs = new std::vector< std::vector* >(); for(unsigned int i=0; i pDir = GetDirection(i); auto idxs = new std::vector(); for(unsigned int j=0; j cDir = GetDirection(j); double angle = ( 180 / itk::Math::pi ) * acos( dot_product(pDir, cDir) ); if( (angle < alphaMaxDegree) && (angle > alphaMinDegree) ) { idxs->push_back(j); } } m_AngularRangeIdxs->push_back(idxs); } } m_MutexAngularRange.Unlock(); // find the maximum (or minimum) direction (remember index and value) T mode; int pIdx = -1; if(invert == 0) { pIdx = GetPrincipalDiffusionDirectionIndex(); mode = (*this)[pIdx]; } else { mode = NumericTraits::max(); for( unsigned int i=0; i nbs = GetNeighbors(pIdx); //////std::vector modeAndNeighborVals; //////modeAndNeighborVals.push_back(mode); //////int numNeighbors = nbs.size(); //////for(int i=0; i odfValuesInAngularRange; int numInRange = m_AngularRangeIdxs->at(pIdx)->size(); for(int i=0; iat(pIdx))[i] ]); } // sort them by value std::sort( odfValuesInAngularRange.begin(), odfValuesInAngularRange.end() ); // median of angular range T median = odfValuesInAngularRange[floor(quantile*(double)numInRange+0.5)]; // compute and return final value if(mode > median) { return mode/median - 1.0; } else { return median/mode - 1.0; } } /** * Calculate Normalized Entropy */ template < typename T, unsigned int N > T itk::OrientationDistributionFunction ::GetNormalizedEntropy() const { double mean = 0; for( unsigned int i=0; i OrientationDistributionFunction OrientationDistributionFunction ::PreMultiply( const MatrixType & m ) const { Self result; typedef typename NumericTraits::AccumulateType AccumulateType; for(unsigned int r=0; r::ZeroValue(); for(unsigned int t=0; t( sum ); } } return result; } /** * Post-multiply the Tensor by a Matrix */ template OrientationDistributionFunction OrientationDistributionFunction ::PostMultiply( const MatrixType & m ) const { Self result; typedef typename NumericTraits::AccumulateType AccumulateType; for(unsigned int r=0; r::ZeroValue(); for(unsigned int t=0; t( sum ); } } return result; } /** * Print content to an ostream */ template std::ostream & operator<<(std::ostream& os,const OrientationDistributionFunction & c ) { for(unsigned int i=0; i::PrintType>(c[i]) << " "; } return os; } /** * Read content from an istream */ template std::istream & operator>>(std::istream& is, OrientationDistributionFunction & dt ) { for(unsigned int i=0; i < dt.GetNumberOfComponents(); i++) { is >> dt[i]; } return is; } } // end namespace itk #endif diff --git a/Modules/DiffusionCore/IODataStructures/Properties/mitkBValueMapPropertySerializer.cpp b/Modules/DiffusionCore/IODataStructures/Properties/mitkBValueMapPropertySerializer.cpp index c5e1654..1c12e40 100644 --- a/Modules/DiffusionCore/IODataStructures/Properties/mitkBValueMapPropertySerializer.cpp +++ b/Modules/DiffusionCore/IODataStructures/Properties/mitkBValueMapPropertySerializer.cpp @@ -1,131 +1,130 @@ /*=================================================================== 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 mitkBValueMapPropertySerializer_h_included #define mitkBValueMapPropertySerializer_h_included #include "mitkBasePropertySerializer.h" - #include "mitkBValueMapProperty.h" - #include +#include namespace mitk { class MITKDIFFUSIONCORE_EXPORT BValueMapPropertySerializer : public BasePropertySerializer { protected: void split(const std::string &s, char delim, std::vector &elems) { std::stringstream ss(s); std::string item; while (std::getline(ss, item, delim)) { elems.push_back(std::atoi(item.c_str())); } } std::vector split(const std::string &s, char delim) { std::vector elems; split(s, delim, elems); return elems; } public: mitkClassMacro( BValueMapPropertySerializer, BasePropertySerializer ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) - TiXmlElement* Serialize() override + tinyxml2::XMLElement* Serialize(tinyxml2::XMLDocument& doc) override { if (const BValueMapProperty* prop = dynamic_cast(m_Property.GetPointer())) { BValueMapProperty::BValueMap map = prop->GetBValueMap(); if(map.empty()) return nullptr; BValueMapProperty::BValueMap::const_iterator it = map.begin(); BValueMapProperty::BValueMap::const_iterator end = map.end(); - auto element = new TiXmlElement("bvaluemap"); + auto element = doc.NewElement("bvaluemap"); while (it != end) { - auto child = new TiXmlElement("entry"); + auto child = doc.NewElement("entry"); { std::stringstream ss; ss << it->first; - child->SetAttribute("key", ss.str()); + child->SetAttribute("key", ss.str().c_str()); } { std::stringstream ss; for(unsigned int i = 0 ; i < it->second.size(); i++) { ss << it->second[i] << ","; } - child->SetAttribute("value", ss.str()); + child->SetAttribute("value", ss.str().c_str()); } - element->InsertEndChild(*child); + element->InsertEndChild(child); ++it; } return element; } else return nullptr; } - BaseProperty::Pointer Deserialize(TiXmlElement* element) override + BaseProperty::Pointer Deserialize(const tinyxml2::XMLElement* element) override { if (!element) return nullptr; BValueMapProperty::BValueMap map; - TiXmlElement* entry = element->FirstChildElement( "entry" )->ToElement(); + const tinyxml2::XMLElement* entry = element->FirstChildElement( "entry" )->ToElement(); while(entry != nullptr){ std::string key, value; - entry->QueryStringAttribute("key",&key); - entry->QueryStringAttribute("value",&value); + key = std::string(entry->Attribute("key")); + value = std::string(entry->Attribute("value")); std::vector indices = split(value.c_str(), ','); map[std::atoi(key.c_str())] = indices; entry = entry->NextSiblingElement( "entry" ); } return BValueMapProperty::New(map).GetPointer(); } protected: BValueMapPropertySerializer(){} ~BValueMapPropertySerializer() override {} }; } // namespace // important to put this into the GLOBAL namespace (because it starts with 'namespace mitk') MITK_REGISTER_SERIALIZER(BValueMapPropertySerializer) #endif diff --git a/Modules/DiffusionCore/IODataStructures/Properties/mitkGradientDirectionsPropertySerializer.cpp b/Modules/DiffusionCore/IODataStructures/Properties/mitkGradientDirectionsPropertySerializer.cpp index 3c3b047..f3e0890 100644 --- a/Modules/DiffusionCore/IODataStructures/Properties/mitkGradientDirectionsPropertySerializer.cpp +++ b/Modules/DiffusionCore/IODataStructures/Properties/mitkGradientDirectionsPropertySerializer.cpp @@ -1,106 +1,105 @@ /*=================================================================== 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 mitkGradientDirectionPropertySerializer_h_included #define mitkGradientDirectionPropertySerializer_h_included #include "mitkBasePropertySerializer.h" - #include "mitkGradientDirectionsProperty.h" - #include "MitkDiffusionCoreExports.h" +#include namespace mitk { class MITKDIFFUSIONCORE_EXPORT GradientDirectionsPropertySerializer : public BasePropertySerializer { public: mitkClassMacro( GradientDirectionsPropertySerializer, BasePropertySerializer ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) - TiXmlElement* Serialize() override + tinyxml2::XMLElement* Serialize(tinyxml2::XMLDocument& doc) override { if (const GradientDirectionsProperty* prop = dynamic_cast(m_Property.GetPointer())) { typedef mitk::GradientDirectionsProperty::GradientDirectionsContainerType GradientDirectionsContainerType; GradientDirectionsContainerType::ConstPointer gdc = prop->GetGradientDirectionsContainer().GetPointer(); if(gdc.IsNull() || gdc->Size() == 0) return nullptr; GradientDirectionsContainerType::ConstIterator it = gdc->Begin(); GradientDirectionsContainerType::ConstIterator end = gdc->End(); - auto element = new TiXmlElement("gradientdirections"); + auto element = doc.NewElement("gradientdirections"); while (it != end) { - auto child = new TiXmlElement("entry"); + auto child = doc.NewElement("entry"); std::stringstream ss; ss << it.Value(); - child->SetAttribute("value", ss.str()); - element->InsertEndChild(*child); + child->SetAttribute("value", ss.str().c_str()); + element->InsertEndChild(child); ++it; } return element; } else return nullptr; } - BaseProperty::Pointer Deserialize(TiXmlElement* element) override + BaseProperty::Pointer Deserialize(const tinyxml2::XMLElement* element) override { if (!element) return nullptr; mitk::GradientDirectionsProperty::GradientDirectionsContainerType::Pointer gdc; gdc = mitk::GradientDirectionsProperty::GradientDirectionsContainerType::New(); - TiXmlElement* entry = element->FirstChildElement( "entry" )->ToElement(); + const tinyxml2::XMLElement* entry = element->FirstChildElement( "entry" )->ToElement(); while(entry != nullptr){ std::stringstream ss; std::string value; - entry->QueryStringAttribute("value",&value); + value = std::string(entry->Attribute("value")); ss << value; vnl_vector_fixed vector; vector.read_ascii(ss); gdc->push_back(vector); entry = entry->NextSiblingElement( "entry" ); } return GradientDirectionsProperty::New(gdc).GetPointer(); } protected: GradientDirectionsPropertySerializer() {} ~GradientDirectionsPropertySerializer() override {} }; } // namespace // important to put this into the GLOBAL namespace (because it starts with 'namespace mitk') MITK_REGISTER_SERIALIZER(GradientDirectionsPropertySerializer) #endif diff --git a/Modules/DiffusionCore/IODataStructures/Properties/mitkMeasurementFramePropertySerializer.cpp b/Modules/DiffusionCore/IODataStructures/Properties/mitkMeasurementFramePropertySerializer.cpp index 3251898..3638e25 100644 --- a/Modules/DiffusionCore/IODataStructures/Properties/mitkMeasurementFramePropertySerializer.cpp +++ b/Modules/DiffusionCore/IODataStructures/Properties/mitkMeasurementFramePropertySerializer.cpp @@ -1,89 +1,88 @@ /*=================================================================== 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 mitkMeasurementFramePropertySerializer_h_included #define mitkMeasurementFramePropertySerializer_h_included #include "mitkBasePropertySerializer.h" - #include "mitkMeasurementFrameProperty.h" - #include +#include namespace mitk { class MITKDIFFUSIONCORE_EXPORT MeasurementFramePropertySerializer : public BasePropertySerializer { public: mitkClassMacro( MeasurementFramePropertySerializer, BasePropertySerializer ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) - TiXmlElement* Serialize() override + tinyxml2::XMLElement* Serialize(tinyxml2::XMLDocument& doc) override { if (const MeasurementFrameProperty* prop = dynamic_cast(m_Property.GetPointer())) { typedef mitk::MeasurementFrameProperty::MeasurementFrameType MeasurementFrameType; const MeasurementFrameType & mft = prop->GetMeasurementFrame(); if(mft.is_zero()) return nullptr; - auto element = new TiXmlElement("measurementframe"); + auto element = doc.NewElement("measurementframe"); - auto child = new TiXmlElement("entry"); + auto child = doc.NewElement("entry"); std::stringstream ss; ss << mft; - child->SetAttribute("value", ss.str()); - element->InsertEndChild(*child); + child->SetAttribute("value", ss.str().c_str()); + element->InsertEndChild(child); return element; } else return nullptr; } - BaseProperty::Pointer Deserialize(TiXmlElement* element) override + BaseProperty::Pointer Deserialize(const tinyxml2::XMLElement* element) override { if (!element) return nullptr; - TiXmlElement* entry = element->FirstChildElement( "entry" )->ToElement(); + const tinyxml2::XMLElement* entry = element->FirstChildElement( "entry" )->ToElement(); std::stringstream ss; std::string value; - entry->QueryStringAttribute("value",&value); + value = std::string(entry->Attribute("value")); ss << value; MeasurementFrameProperty::MeasurementFrameType matrix; matrix.read_ascii(ss); return MeasurementFrameProperty::New(matrix).GetPointer(); } protected: MeasurementFramePropertySerializer() {} ~MeasurementFramePropertySerializer() override {} }; } // namespace // important to put this into the GLOBAL namespace (because it starts with 'namespace mitk') MITK_REGISTER_SERIALIZER(MeasurementFramePropertySerializer) #endif diff --git a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp index 74e41d0..befe565 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp +++ b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp @@ -1,2840 +1,2840 @@ /*=================================================================== 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, 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::ColorSinglePoint(int f_idx, int p_idx, double rgb[3]) { // 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(num_points, idList); // fiberList->GetCell(f_idx, num_points, idList); vtkCell* cell = m_FiberPolyData->GetCell(f_idx); // /* single fiber checkpoints: is number of points valid */ // if (p_idx < num_points) // { rgba[0] = static_cast(255.0 * rgb[0]); rgba[1] = static_cast(255.0 * rgb[1]); rgba[2] = static_cast(255.0 * rgb[2]); rgba[3] = 255; m_FiberColors->InsertTypedTuple(cell->GetPointId(p_idx), rgba); // } // } 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, 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, bool weight_fibers) { 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, 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 (unsigned int i=0; iGetCell(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->SetCellIdsArrayName(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 with max. error " << error << "mm"; 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/mitkPeakImage.cpp b/Modules/DiffusionCore/IODataStructures/mitkPeakImage.cpp index a4d0bdd..6da36f4 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkPeakImage.cpp +++ b/Modules/DiffusionCore/IODataStructures/mitkPeakImage.cpp @@ -1,221 +1,222 @@ /*=================================================================== 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 "mitkPeakImage.h" #include "mitkImageDataItem.h" #include "mitkImageCast.h" #include #include #include #include #include +#include mitk::PeakImage::PeakImage() : Image() { m_PolyData = nullptr; } mitk::PeakImage::~PeakImage() { } void mitk::PeakImage::ConstructPolydata() { MITK_INFO << "PeakImage constructing polydata"; if (this->GetDimensions()[3]%3!=0) mitkThrow() << "Fourth dimension needs to be a multiple of 3"; typedef mitk::ImageToItk< ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(this); caster->Update(); ItkPeakImageType::ConstPointer itkImg = caster->GetOutput(); int sz_x = static_cast(this->GetDimensions()[0]); int sz_y = static_cast(this->GetDimensions()[1]); int sz_z = static_cast(this->GetDimensions()[2]); unsigned int num_dirs = this->GetDimensions()[3]/3; double minSpacing = 1; ItkPeakImageType::SpacingType spacing = itkImg->GetSpacing(); if(spacing[0] m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); #ifdef WIN32 #pragma omp parallel for #else #pragma omp parallel for collapse(3) #endif for (int x=0; x idx4; idx4[0] = x; idx4[1] = y; idx4[2] = z; for (unsigned int i=0; i center; center[0] = idx4[0]; center[1] = idx4[1]; center[2] = idx4[2]; center[3] = 0; itk::Point worldCenter; itkImg->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); idx4[3] = 0; vnl_vector_fixed< float, 3 > dir; idx4[3] = i*3; dir[0] = itkImg->GetPixel(idx4); idx4[3] += 1; dir[1] = itkImg->GetPixel(idx4); idx4[3] += 1; dir[2] = itkImg->GetPixel(idx4); if (dir.magnitude()<0.001f) continue; #pragma omp critical { vtkSmartPointer container = vtkSmartPointer::New(); // add direction to vector field (with spacing compensation) itk::Point worldStart; worldStart[0] = worldCenter[0]-static_cast(dir[0])/2*minSpacing; worldStart[1] = worldCenter[1]-static_cast(dir[1])/2*minSpacing; worldStart[2] = worldCenter[2]-static_cast(dir[2])/2*minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+static_cast(dir[0])/2*minSpacing; worldEnd[1] = worldCenter[1]+static_cast(dir[1])/2*minSpacing; worldEnd[2] = worldCenter[2]+static_cast(dir[2])/2*minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); } } } m_PolyData = vtkSmartPointer::New(); m_PolyData->SetPoints(m_VtkPoints); m_PolyData->SetLines(m_VtkCellArray); MITK_INFO << "PeakImage finished constructing polydata"; ColorByOrientation(); this->Modified(); } void mitk::PeakImage::SetCustomColor(float r, float g, float b) { vtkPoints* extrPoints = nullptr; extrPoints = m_PolyData->GetPoints(); int numOfPoints = 0; if (extrPoints!=nullptr) numOfPoints = static_cast(extrPoints->GetNumberOfPoints()); int componentSize = 4; vtkSmartPointer colors = vtkSmartPointer::New(); colors->Allocate(numOfPoints * componentSize); colors->SetNumberOfComponents(componentSize); colors->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] = 255; colors->InsertTypedTuple(i, rgba); } m_PolyData->GetPointData()->AddArray(colors); this->Modified(); } void mitk::PeakImage::ColorByOrientation() { //===== 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 = nullptr; extrPoints = m_PolyData->GetPoints(); int numOfPoints = 0; if (extrPoints!=nullptr) numOfPoints = static_cast(extrPoints->GetNumberOfPoints()); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = 4; vtkSmartPointer colors = vtkSmartPointer::New(); colors->Allocate(numOfPoints * componentSize); colors->SetNumberOfComponents(componentSize); colors->SetName("FIBER_COLORS"); int num_peaks = static_cast(m_PolyData->GetNumberOfLines()); if (num_peaks < 1) return; vtkCellArray* fiberList = m_PolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(num_points, idList); if (num_points == 2) { vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[0])[0], extrPoints->GetPoint(idList[0])[1],extrPoints->GetPoint(idList[0])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[1])[0], extrPoints->GetPoint(idList[1])[1], extrPoints->GetPoint(idList[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] = 255; colors->InsertTypedTuple(idList[0], rgba); colors->InsertTypedTuple(idList[1], rgba); } else { MITK_INFO << "Peak with point number!=2 detected!" ; continue; } } m_PolyData->GetPointData()->AddArray(colors); this->Modified(); } diff --git a/Modules/DiffusionCore/Rendering/mitkOdfNormalizationMethodProperty.h b/Modules/DiffusionCore/Rendering/mitkOdfNormalizationMethodProperty.h index ac503e7..53c8979 100644 --- a/Modules/DiffusionCore/Rendering/mitkOdfNormalizationMethodProperty.h +++ b/Modules/DiffusionCore/Rendering/mitkOdfNormalizationMethodProperty.h @@ -1,137 +1,140 @@ /*=================================================================== 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_OdfNormalizationMethodProperty__H_ #define _MITK_OdfNormalizationMethodProperty__H_ #include "mitkEnumerationProperty.h" #include "mitkEnumerationPropertySerializer.h" #include #include "mitkSerializerMacros.h" +#include namespace mitk { enum OdfNormalizationMethod { ODFN_MINMAX, ODFN_MAX, ODFN_NONE, ODFN_GLOBAL_MAX }; /** * Encapsulates the enumeration for ODF normalization. Valid values are * ODFN_MINMAX, ODFN_MAX, ODFN_NONE ODFN_GLOBAL_MAX * Default is ODFN_MINMAX */ class MITKDIFFUSIONCORE_EXPORT OdfNormalizationMethodProperty : public EnumerationProperty { public: mitkClassMacro( OdfNormalizationMethodProperty, EnumerationProperty ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(OdfNormalizationMethodProperty, const IdType&); mitkNewMacro1Param(OdfNormalizationMethodProperty, const std::string&); /** * Returns the current interpolation value as defined by VTK constants. */ virtual int GetNormalization(); /** * Sets the interpolation type to ODFN_MINMAX. */ virtual void SetNormalizationToMinMax(); /** * Sets the interpolation type to ODFN_MAX. */ virtual void SetNormalizationToMax(); /** * Sets the interpolation type to ODFN_NONE. */ virtual void SetNormalizationToNone(); /** * Sets the interpolation type to ODFN_GLOBAL_MAX. */ virtual void SetNormalizationToGlobalMax(); protected: /** Sets reslice interpolation mode to default (VTK_RESLICE_NEAREST). */ OdfNormalizationMethodProperty( ); /** * Constructor. Sets reslice interpolation to the given value. */ OdfNormalizationMethodProperty( const IdType& value ); /** * Constructor. Sets reslice interpolation to the given value. */ OdfNormalizationMethodProperty( const std::string& value ); /** * this function is overridden as protected, so that the user may not add * additional invalid interpolation types. */ bool AddEnum( const std::string& name, const IdType& id ) override; /** * Adds the enumeration types as defined by vtk to the list of known * enumeration values. */ virtual void AddInterpolationTypes(); }; class MITKDIFFUSIONCORE_EXPORT OdfNormalizationMethodPropertySerializer : public EnumerationPropertySerializer { public: mitkClassMacro( OdfNormalizationMethodPropertySerializer, EnumerationPropertySerializer ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) - BaseProperty::Pointer Deserialize(TiXmlElement* element) override + BaseProperty::Pointer Deserialize(const tinyxml2::XMLElement* element) override { if (!element) return nullptr; - const char* sa( element->Attribute("value") ); - std::string s(sa?sa:""); + + std::string s = ""; + s = std::string(element->Attribute("value")); + OdfNormalizationMethodProperty::Pointer property = OdfNormalizationMethodProperty::New(); property->SetValue( s ); return property.GetPointer(); } protected: OdfNormalizationMethodPropertySerializer () {} ~OdfNormalizationMethodPropertySerializer () override {} }; } // end of namespace mitk #endif diff --git a/Modules/DiffusionCore/Rendering/mitkOdfScaleByProperty.h b/Modules/DiffusionCore/Rendering/mitkOdfScaleByProperty.h index 5ab3e6e..49a4d88 100644 --- a/Modules/DiffusionCore/Rendering/mitkOdfScaleByProperty.h +++ b/Modules/DiffusionCore/Rendering/mitkOdfScaleByProperty.h @@ -1,117 +1,120 @@ /*=================================================================== 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_OdfScaleByProperty__H_ #define _MITK_OdfScaleByProperty__H_ #include #include "mitkEnumerationProperty.h" #include "mitkEnumerationPropertySerializer.h" +#include namespace mitk { enum OdfScaleBy { ODFSB_NONE, ODFSB_GFA, ODFSB_PC }; /** * Encapsulates the enumeration for ODF normalization. Valid values are * ODFSB_NONE, ODFSB_GFA, ODFSB_PC * Default is ODFSB_NONE */ class MITKDIFFUSIONCORE_EXPORT OdfScaleByProperty : public EnumerationProperty { public: mitkClassMacro( OdfScaleByProperty, EnumerationProperty ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(OdfScaleByProperty, const IdType&); mitkNewMacro1Param(OdfScaleByProperty, const std::string&); /** * Returns the current interpolation value as defined by VTK constants. */ virtual int GetScaleBy(); virtual void SetScaleByNothing(); virtual void SetScaleByGFA(); virtual void SetScaleByPrincipalCurvature(); protected: /** Sets reslice interpolation mode to default (ODFSB_NONE). */ OdfScaleByProperty( ); /** * Constructor. Sets reslice interpolation to the given value. */ OdfScaleByProperty( const IdType& value ); /** * Constructor. Sets reslice interpolation to the given value. */ OdfScaleByProperty( const std::string& value ); /** * this function is overridden as protected, so that the user may not add * additional invalid interpolation types. */ bool AddEnum( const std::string& name, const IdType& id ) override; /** * Adds the enumeration types as defined by vtk to the list of known * enumeration values. */ virtual void AddInterpolationTypes(); }; class MITKDIFFUSIONCORE_EXPORT OdfScaleByPropertySerializer : public EnumerationPropertySerializer { public: mitkClassMacro( OdfScaleByPropertySerializer, EnumerationPropertySerializer ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) - BaseProperty::Pointer Deserialize(TiXmlElement* element) override + BaseProperty::Pointer Deserialize(const tinyxml2::XMLElement* element) override { if (!element) return nullptr; - const char* sa( element->Attribute("value") ); - std::string s(sa?sa:""); + + std::string s = ""; + s = std::string(element->Attribute("value")); + OdfScaleByProperty::Pointer property = OdfScaleByProperty::New(); property->SetValue( s ); return property.GetPointer(); } protected: OdfScaleByPropertySerializer () {} ~OdfScaleByPropertySerializer () override {} }; } // end of namespace mitk #endif diff --git a/Modules/DiffusionIO/CMakeLists.txt b/Modules/DiffusionIO/CMakeLists.txt index 4fcc91b..9eb3000 100644 --- a/Modules/DiffusionIO/CMakeLists.txt +++ b/Modules/DiffusionIO/CMakeLists.txt @@ -1,10 +1,10 @@ MITK_CREATE_MODULE( SUBPROJECTS MITK-Diffusion INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR} ReaderWriter ObjectFactories DEPENDS MitkConnectomics MitkDiffusionCore MitkSceneSerialization AUTOLOAD_WITH MitkCore PACKAGE_DEPENDS PUBLIC DCMTK - PRIVATE tinyxml + PRIVATE tinyxml2 ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/QmitkTbssRoiAnalysisWidget.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/QmitkTbssRoiAnalysisWidget.cpp index f6f2383..d83dbb9 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/QmitkTbssRoiAnalysisWidget.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/QmitkTbssRoiAnalysisWidget.cpp @@ -1,967 +1,967 @@ /*=================================================================== 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 "QmitkTbssRoiAnalysisWidget.h" #include "mitkImagePixelReadAccessor.h" #include "mitkPixelTypeMultiplex.h" #include #include #include #include QmitkTbssRoiAnalysisWidget::QmitkTbssRoiAnalysisWidget( QWidget * parent ) : QmitkPlotWidget(parent) { m_PlotPicker = new QwtPlotPicker(m_Plot->canvas()); m_PlotPicker->setStateMachine(new QwtPickerDragPointMachine()); m_PlotPicker->setTrackerMode(QwtPicker::ActiveOnly); m_PlottingFiberBundle = false; } void QmitkTbssRoiAnalysisWidget::DoPlotFiberBundles(mitk::FiberBundle *fib, mitk::Image* img, mitk::DataNode* startRoi, mitk::DataNode* endRoi, bool avg, int number) { TractContainerType tracts = CreateTracts(fib, startRoi, endRoi); TractContainerType resampledTracts = ParameterizeTracts(tracts, number); // Now we have the resampled tracts. Next we should use these points to read out the values mitkPixelTypeMultiplex3(PlotFiberBundles,img->GetImageDescriptor()->GetChannelTypeById(0),resampledTracts, img, avg); m_CurrentTracts = resampledTracts; } TractContainerType QmitkTbssRoiAnalysisWidget::CreateTracts(mitk::FiberBundle *fib, mitk::DataNode *startRoi, mitk::DataNode *endRoi) { mitk::PlaneGeometry* startGeometry2D = const_cast(dynamic_cast(startRoi->GetData())->GetPlaneGeometry()); mitk::PlaneGeometry* endGeometry2D = const_cast(dynamic_cast(endRoi->GetData())->GetPlaneGeometry()); mitk::Point3D startCenter = dynamic_cast(startRoi->GetData())->GetWorldControlPoint(0); //center Point of start roi mitk::Point3D endCenter = dynamic_cast(startRoi->GetData())->GetWorldControlPoint(0); //center Point of end roi mitk::FiberBundle::Pointer inStart = fib->ExtractFiberSubset(startRoi, nullptr); mitk::FiberBundle::Pointer inBoth = inStart->ExtractFiberSubset(endRoi, nullptr); int num = inBoth->GetNumFibers(); TractContainerType tracts; vtkSmartPointer fiberPolyData = inBoth->GetFiberPolyData(); vtkCellArray* lines = fiberPolyData->GetLines(); lines->InitTraversal(); // Now find out for each fiber which ROI is encountered first. If this is the startRoi, the direction is ok // Otherwise the plot should be in the reverse direction for( int fiberID( 0 ); fiberID < num; fiberID++ ) { vtkIdType numPointsInCell(0); - vtkIdType* pointsInCell(nullptr); + vtkIdType const* pointsInCell(nullptr); lines->GetNextCell ( numPointsInCell, pointsInCell ); int startId = 0; int endId = 0; mitk::ScalarType minDistStart = std::numeric_limits::max(); mitk::ScalarType minDistEnd = std::numeric_limits::max(); for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) { mitk::ScalarType *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); mitk::Point3D point; point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; mitk::ScalarType distanceToStart = point.EuclideanDistanceTo(startCenter); mitk::ScalarType distanceToEnd = point.EuclideanDistanceTo(endCenter); if(distanceToStart < minDistStart) { minDistStart = distanceToStart; startId = pointInCellID; } if(distanceToEnd < minDistEnd) { minDistEnd = distanceToEnd; endId = pointInCellID; } } /* We found the start and end points of of the part that should be plottet for the current fiber. now we need to plot them. If the endId is smaller than the startId the plot order must be reversed*/ TractType singleTract; PointType point; if(startId < endId) { // Calculate the intersection of the ROI with the startRoi and decide if the startId is part of the roi or must be cut of mitk::ScalarType *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); mitk::Vector3D p0; p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ startId+1 ] ); mitk::Vector3D p1; p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; // Check if p and p2 are both on the same side of the plane mitk::Vector3D normal = startGeometry2D->GetNormal(); mitk::Point3D pStart; pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); bool secondOnPositive = startGeometry2D->IsAbove(pSecond); mitk::Vector3D onPlane; onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId+1 lie on the same side of the plane, so we need need startId-1 to calculate the intersection with the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; } mitk::ScalarType d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); point[0] = d*newPoint[0] + p0[0]; point[1] = d*newPoint[1] + p0[1]; point[2] = d*newPoint[2] + p0[2]; singleTract.push_back(point); if(! (secondOnPositive ^ startOnPositive) ) { /* StartId and startId+1 lie on the same side of the plane so startId is also part of the ROI*/ mitk::ScalarType *start = fiberPolyData->GetPoint( pointsInCell[startId] ); point[0] = start[0]; point[1] = start[1]; point[2] = start[2]; singleTract.push_back(point); } for( int pointInCellID( startId+1 ); pointInCellID < endId ; pointInCellID++) { // push back point mitk::ScalarType *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; singleTract.push_back( point ); } /* endId must be included if endId and endId-1 lie on the same side of the plane defined by endRoi*/ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; mitk::Point3D pLast; pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; mitk::Point3D pBeforeLast; pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; normal = endGeometry2D->GetNormal(); bool lastOnPositive = endGeometry2D->IsAbove(pLast); bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); onPlane[0] = endCenter[0]; onPlane[1] = endCenter[1]; onPlane[2] = endCenter[2]; if(! (lastOnPositive ^ secondLastOnPositive) ) { /* endId and endId-1 lie on the same side of the plane, so we need need endId+1 to calculate the intersection with the planar figure. this should exist since we know that the fiber crosses the planar figure endId is also part of the tract and can be inserted here */ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; singleTract.push_back( point ); p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); point[0] = d*newPoint[0] + p0[0]; point[1] = d*newPoint[1] + p0[1]; point[2] = d*newPoint[2] + p0[2]; singleTract.push_back(point); } else{ // Calculate the intersection of the ROI with the startRoi and decide if the startId is part of the roi or must be cut of mitk::ScalarType *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); mitk::Vector3D p0; p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); mitk::Vector3D p1; p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; // Check if p and p2 are both on the same side of the plane mitk::Vector3D normal = startGeometry2D->GetNormal(); mitk::Point3D pStart; pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); bool secondOnPositive = startGeometry2D->IsAbove(pSecond); mitk::Vector3D onPlane; onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId+1 lie on the same side of the plane, so we need need startId-1 to calculate the intersection with the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; } mitk::ScalarType d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); point[0] = d*newPoint[0] + p0[0]; point[1] = d*newPoint[1] + p0[1]; point[2] = d*newPoint[2] + p0[2]; singleTract.push_back(point); if(! (secondOnPositive ^ startOnPositive) ) { /* StartId and startId+1 lie on the same side of the plane so startId is also part of the ROI*/ mitk::ScalarType *start = fiberPolyData->GetPoint( pointsInCell[startId] ); point[0] = start[0]; point[1] = start[1]; point[2] = start[2]; singleTract.push_back(point); } for( int pointInCellID( startId-1 ); pointInCellID > endId ; pointInCellID--) { // push back point mitk::ScalarType *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; singleTract.push_back( point ); } /* endId must be included if endId and endI+1 lie on the same side of the plane defined by endRoi*/ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; mitk::Point3D pLast; pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; mitk::Point3D pBeforeLast; pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; bool lastOnPositive = endGeometry2D->IsAbove(pLast); bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); normal = endGeometry2D->GetNormal(); onPlane[0] = endCenter[0]; onPlane[1] = endCenter[1]; onPlane[2] = endCenter[2]; if(! (lastOnPositive ^ secondLastOnPositive) ) { /* endId and endId+1 lie on the same side of the plane, so we need need endId-1 to calculate the intersection with the planar figure. this should exist since we know that the fiber crosses the planar figure endId is also part of the tract and can be inserted here */ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; singleTract.push_back( point ); p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); point[0] = d*newPoint[0] + p0[0]; point[1] = d*newPoint[1] + p0[1]; point[2] = d*newPoint[2] + p0[2]; singleTract.push_back(point); } tracts.push_back(singleTract); } return tracts; } void QmitkTbssRoiAnalysisWidget::PlotFiberBetweenRois(mitk::FiberBundle *fib, mitk::Image* img, mitk::DataNode* startRoi, mitk::DataNode* endRoi, bool avg, int number) { if(fib == nullptr || img == nullptr || startRoi == nullptr || endRoi == nullptr) return; m_Fib = fib; m_CurrentImage = img; m_CurrentStartRoi = startRoi; m_CurrentEndRoi = endRoi; DoPlotFiberBundles(fib, img, startRoi, endRoi, avg, number); } void QmitkTbssRoiAnalysisWidget::ModifyPlot(int number, bool avg) { if(m_Fib == nullptr || m_CurrentTbssImage == nullptr || m_CurrentStartRoi == nullptr || m_CurrentEndRoi == nullptr) return; if(m_PlottingFiberBundle) { DoPlotFiberBundles(m_Fib, m_CurrentImage, m_CurrentStartRoi, m_CurrentEndRoi, avg, number); } else { PlotFiber4D(m_CurrentTbssImage, m_Fib, m_CurrentStartRoi, m_CurrentEndRoi, number); } } TractContainerType QmitkTbssRoiAnalysisWidget::ParameterizeTracts(TractContainerType tracts, int number) { TractContainerType resampledTracts; for(auto it = tracts.begin(); it != tracts.end(); ++it) { TractType resampledTract; TractType tract = *it; // Calculate the total length mitk::ScalarType totalLength = 0; if(tract.size() < 2) continue; PointType p0 = tract.at(0); for(unsigned int i = 1; i(resampledTract.size()) <= (number+1); position+=stepSize) { /* In case we walked to far we need to find the next segment we are on and on what relative position on that tract we are on. Small correction for rounding errors */ while(locationBetween > distance+0.001) { if(tractCounter == tract.size()) std::cout << "problem"; // Determine by what distance we are no on the next segment locationBetween = locationBetween - distance; p0 = p1; p1 = tract.at(tractCounter); tractCounter++; distance = p0.EuclideanDistanceTo(p1); } // Direction PointType::VectorType direction = p1-p0; direction.Normalize(); PointType newSample = p0 + direction*locationBetween; resampledTract.push_back(newSample); locationBetween += stepSize; } resampledTracts.push_back(resampledTract); } return resampledTracts; } mitk::Point3D QmitkTbssRoiAnalysisWidget::GetPositionInWorld(int index) { mitk::ScalarType xSum = 0.0; mitk::ScalarType ySum = 0.0; mitk::ScalarType zSum = 0.0; for(auto it = m_CurrentTracts.begin(); it!=m_CurrentTracts.end(); ++it) { TractType tract = *it; PointType p = tract.at(index); xSum += p[0]; ySum += p[1]; zSum += p[2]; } int number = m_CurrentTracts.size(); mitk::ScalarType xPos = xSum / number; mitk::ScalarType yPos = ySum / number; mitk::ScalarType zPos = zSum / number; mitk::Point3D pos; pos[0] = xPos; pos[1] = yPos; pos[2] = zPos; return pos; } std::vector< std::vector > QmitkTbssRoiAnalysisWidget::CalculateGroupProfiles() { MITK_INFO << "make profiles!"; std::vector< std::vector > profiles; int size = m_Projections->GetVectorLength(); for(int s=0; s profile; RoiType::iterator it; it = m_Roi.begin(); while(it != m_Roi.end()) { itk::Index<3> ix = *it; profile.push_back(m_Projections->GetPixel(ix).GetElement(s)); it++; } profiles.push_back(profile); } m_IndividualProfiles = profiles; // Calculate the averages // Here a check could be build in to check whether all profiles have // the same length, but this should normally be the case if the input // data were corrected with the TBSS Module. std::vector< std::vector > groupProfiles; std::vector< std::pair >::iterator it; it = m_Groups.begin(); int c = 0; //the current profile number while(it != m_Groups.end() && profiles.size() > 0) { std::pair p = *it; int size = p.second; //initialize a vector of the right length with zeroes std::vector averageProfile; for(unsigned int i=0; i > groupProfiles = CalculateGroupProfiles(); Plot(groupProfiles); } void QmitkTbssRoiAnalysisWidget::Plot(std::vector > groupProfiles) { this->Clear(); m_Vals.clear(); std::vector v1; std::vector xAxis; for(unsigned int i=0; iSetPlotTitle( title.c_str() ); QPen pen( Qt::SolidLine ); pen.setWidth(2); std::vector< std::pair >::iterator it; it = m_Groups.begin(); int c = 0; //the current profile number QColor colors[4] = {Qt::green, Qt::blue, Qt::yellow, Qt::red}; while(it != m_Groups.end() && groupProfiles.size() > 0) { std::pair< std::string, int > group = *it; pen.setColor(colors[c]); int curveId = this->InsertCurve( group.first.c_str() ); this->SetCurveData( curveId, xAxis, groupProfiles.at(c) ); this->SetCurvePen( curveId, pen ); c++; it++; } auto legend = new QwtLegend; this->SetLegend(legend, QwtPlot::RightLegend, 0.5); std::cout << m_Measure << std::endl; this->m_Plot->setAxisTitle(0, m_Measure.c_str()); this->m_Plot->setAxisTitle(3, "Position"); this->Replot(); } std::vector< std::vector > QmitkTbssRoiAnalysisWidget::CalculateGroupProfilesFibers(mitk::TbssImage::Pointer tbssImage, mitk::FiberBundle *fib, mitk::DataNode* startRoi, mitk::DataNode* endRoi, int number) { TractContainerType tracts = CreateTracts(fib, startRoi, endRoi); TractContainerType resampledTracts = ParameterizeTracts(tracts, number); int nTracts = resampledTracts.size(); this->Clear(); // For every group we have m fibers * n subjects of profiles to fill std::vector< std::vector > profiles; // calculate individual profiles by going through all n subjects int size = m_Projections->GetVectorLength(); for(int s=0; s profile; auto it = resampledTracts[t].begin(); while(it != resampledTracts[t].end()) { PointType p = *it; PointType index; tbssImage->GetGeometry()->WorldToIndex(p, index); itk::Index<3> ix; ix[0] = index[0]; ix[1] = index[1]; ix[2] = index[2]; // Get value from image profile.push_back(m_Projections->GetPixel(ix).GetElement(s)); it++; } profiles.push_back(profile); } } m_IndividualProfiles = profiles; // Now create the group averages (every group contains m fibers * n_i group members std::vector< std::pair >::iterator it; it = m_Groups.begin(); int c = 0; //the current profile number // Calculate the group averages std::vector< std::vector > groupProfiles; while(it != m_Groups.end() && profiles.size() > 0) { std::pair p = *it; int size = p.second; //initialize a vector of the right length with zeroes std::vector averageProfile; for(unsigned int i=0; i > groupProfiles = CalculateGroupProfilesFibers(tbssImage, fib, startRoi, endRoi, number); Plot(groupProfiles); } template void QmitkTbssRoiAnalysisWidget::PlotFiberBundles(const mitk::PixelType, TractContainerType tracts, mitk::Image *img, bool avg) { m_PlottingFiberBundle = true; this->Clear(); auto it = tracts.begin(); std::vector< std::vector > profiles; mitk::ImagePixelReadAccessor imAccess(img,img->GetVolumeData(0)); it = tracts.begin(); while(it != tracts.end()) { TractType tract = *it; auto tractIt = tract.begin(); std::vector profile; while(tractIt != tract.end()) { PointType p = *tractIt; // Get value from image profile.push_back( (mitk::ScalarType) imAccess.GetPixelByWorldCoordinates(p) ); ++tractIt; } profiles.push_back(profile); std::cout << std::endl; ++it; } if(profiles.size() == 0) return; m_IndividualProfiles = profiles; std::string title = "Fiber bundle plot"; this->SetPlotTitle( title.c_str() ); // initialize average profile std::vector averageProfile; std::vector profile = profiles.at(0); // can do this because we checked the size of profiles before for(unsigned int i=0; i profile = *profit; std::vector xAxis; for(unsigned int i=0; iInsertCurve( "" ); this->SetCurveData( curveId, xAxis, profile ); ++profit; id++; } m_Average = averageProfile; if(avg) { // Draw the average profile std::vector xAxis; for(unsigned int i=0; iInsertCurve( "" ); this->SetCurveData( curveId, xAxis, averageProfile ); QPen pen( Qt::SolidLine ); pen.setWidth(3); pen.setColor(Qt::red); this->SetCurvePen( curveId, pen ); id++; } this->Replot(); } void QmitkTbssRoiAnalysisWidget::drawBar(int x) { m_Plot->detachItems(QwtPlotItem::Rtti_PlotMarker, true); QwtPlotMarker *mX = new QwtPlotMarker(); //mX->setLabel(QString::fromLatin1("selected point")); mX->setLabelAlignment(Qt::AlignLeft | Qt::AlignBottom); mX->setLabelOrientation(Qt::Vertical); mX->setLineStyle(QwtPlotMarker::VLine); mX->setLinePen(QPen(Qt::black, 0, Qt::SolidLine)); mX->setXValue(x); mX->attach(m_Plot); this->Replot(); } QmitkTbssRoiAnalysisWidget::~QmitkTbssRoiAnalysisWidget() { delete m_PlotPicker; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTractbasedSpatialStatisticsView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTractbasedSpatialStatisticsView.cpp index 47e0da0..7946ec5 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTractbasedSpatialStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTractbasedSpatialStatisticsView.cpp @@ -1,1142 +1,1142 @@ /*=================================================================== 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. ===================================================================*/ // Qmitk #include "QmitkTractbasedSpatialStatisticsView.h" #include "mitkDataNodeObject.h" #include // Qt #include #include #include #include #include #include #include #include #include "vtkPoints.h" #include #include #include const std::string QmitkTractbasedSpatialStatisticsView::VIEW_ID = "org.mitk.views.tractbasedspatialstatistics"; using namespace berry; QmitkTractbasedSpatialStatisticsView::QmitkTractbasedSpatialStatisticsView() : QmitkAbstractView() , m_Controls( 0 ) , m_Activated(false) { } QmitkTractbasedSpatialStatisticsView::~QmitkTractbasedSpatialStatisticsView() { } void QmitkTractbasedSpatialStatisticsView::PerformChange() { m_Controls->m_RoiPlotWidget->ModifyPlot(m_Controls->m_Segments->value(), m_Controls->m_Average->isChecked()); } void QmitkTractbasedSpatialStatisticsView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { //datamanager selection changed if (!this->IsActivated()) return; // Check which datatypes are selected in the datamanager and enable/disable widgets accordingly bool foundTbssRoi = false; bool foundTbss = false; bool found3dImage = false; bool found4dImage = false; bool foundFiberBundle = false; bool foundStartRoi = false; bool foundEndRoi = false; mitk::TbssRoiImage* roiImage = nullptr; mitk::TbssImage* image = nullptr; mitk::Image* img = nullptr; mitk::FiberBundle* fib = nullptr; mitk::DataNode* start = nullptr; mitk::DataNode* end = nullptr; m_CurrentStartRoi = nullptr; m_CurrentEndRoi = nullptr; for (mitk::DataNode::Pointer node: nodes) { // only look at interesting types // check for valid data mitk::BaseData* nodeData = node->GetData(); if( nodeData ) { if(QString("TbssRoiImage").compare(nodeData->GetNameOfClass())==0) { foundTbssRoi = true; roiImage = static_cast(nodeData); } else if (QString("TbssImage").compare(nodeData->GetNameOfClass())==0) { foundTbss = true; image = static_cast(nodeData); } else if(QString("Image").compare(nodeData->GetNameOfClass())==0) { img = static_cast(nodeData); if(img->GetDimension() == 3) { found3dImage = true; } else if(img->GetDimension() == 4) { found4dImage = true; } } else if (QString("FiberBundle").compare(nodeData->GetNameOfClass())==0) { foundFiberBundle = true; fib = static_cast(nodeData); this->m_CurrentFiberNode = node; } if(QString("PlanarCircle").compare(nodeData->GetNameOfClass())==0) { if(!foundStartRoi) { start = node; this->m_CurrentStartRoi = node; foundStartRoi = true; } else { end = node; this->m_CurrentEndRoi = node; foundEndRoi = true; } } } } this->m_Controls->m_CreateRoi->setEnabled(found3dImage); this->m_Controls->m_ImportFsl->setEnabled(found4dImage); if(foundTbss && foundTbssRoi) { this->Plot(image, roiImage); } else if(found3dImage && foundFiberBundle && foundStartRoi && foundEndRoi) { this->PlotFiberBundle(fib, img, start, end); } else if(found3dImage && foundFiberBundle) { this->PlotFiberBundle(fib, img); } else if(foundTbss && foundStartRoi && foundEndRoi && foundFiberBundle) { this->PlotFiber4D(image, fib, start, end); } if(found3dImage) { this->InitPointsets(); } this->m_Controls->m_Cut->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi); this->m_Controls->m_SegmentLabel->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && (found3dImage || foundTbss)); this->m_Controls->m_Segments->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && (found3dImage || foundTbss)); this->m_Controls->m_Average->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); } void QmitkTractbasedSpatialStatisticsView::InitPointsets() { // Check if PointSetStart exsits, if not create it. m_P1 = this->GetDataStorage()->GetNamedNode("PointSetNode"); if (m_PointSetNode) { //m_PointSetNode = dynamic_cast(m_P1->GetData()); return; } if ((!m_P1) || (!m_PointSetNode)) { // create new ones m_PointSetNode = mitk::PointSet::New(); m_P1 = mitk::DataNode::New(); m_P1->SetData( m_PointSetNode ); m_P1->SetProperty( "name", mitk::StringProperty::New( "PointSet" ) ); m_P1->SetProperty( "opacity", mitk::FloatProperty::New( 1 ) ); m_P1->SetProperty( "helper object", mitk::BoolProperty::New(true) ); // CHANGE if wanted m_P1->SetProperty( "pointsize", mitk::FloatProperty::New( 0.1 ) ); m_P1->SetColor( 1.0, 0.0, 0.0 ); this->GetDataStorage()->Add(m_P1); m_Controls->m_PointWidget->SetPointSetNode(m_P1); auto renderWindowPart = this->GetRenderWindowPart(mitk::WorkbenchUtil::IRenderWindowPartStrategy::OPEN); auto axialSnc = renderWindowPart->GetQmitkRenderWindow("axial")->GetSliceNavigationController(); auto sagittalSnc = renderWindowPart->GetQmitkRenderWindow("sagittal")->GetSliceNavigationController(); auto coronalSnc = renderWindowPart->GetQmitkRenderWindow("coronal")->GetSliceNavigationController(); m_Controls->m_PointWidget->AddSliceNavigationController(axialSnc); m_Controls->m_PointWidget->AddSliceNavigationController(sagittalSnc); m_Controls->m_PointWidget->AddSliceNavigationController(coronalSnc); } } void QmitkTractbasedSpatialStatisticsView::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::QmitkTractbasedSpatialStatisticsViewControls; m_Controls->setupUi( parent ); this->CreateConnections(); } // Table for the FSL TBSS import m_GroupModel = new QmitkTbssTableModel(); m_Controls->m_GroupInfo->setModel(m_GroupModel); } void QmitkTractbasedSpatialStatisticsView::SetFocus() { m_Controls->m_AddGroup->setFocus(); } void QmitkTractbasedSpatialStatisticsView::Activated() { m_Activated = true; } void QmitkTractbasedSpatialStatisticsView::Deactivated() { m_Activated = false; } bool QmitkTractbasedSpatialStatisticsView::IsActivated() const { return m_Activated; } void QmitkTractbasedSpatialStatisticsView::Visible() { } void QmitkTractbasedSpatialStatisticsView::Hidden() { } void QmitkTractbasedSpatialStatisticsView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_CreateRoi), SIGNAL(clicked()), this, SLOT(CreateRoi()) ); connect( (QObject*)(m_Controls->m_ImportFsl), SIGNAL(clicked()), this, SLOT(TbssImport()) ); connect( (QObject*)(m_Controls->m_AddGroup), SIGNAL(clicked()), this, SLOT(AddGroup()) ); connect( (QObject*)(m_Controls->m_RemoveGroup), SIGNAL(clicked()), this, SLOT(RemoveGroup()) ); connect( (QObject*)(m_Controls->m_Clipboard), SIGNAL(clicked()), this, SLOT(CopyToClipboard()) ); connect( m_Controls->m_RoiPlotWidget->m_PlotPicker, SIGNAL(selected(const QPointF&)), SLOT(Clicked(const QPointF&) ) ); connect( m_Controls->m_RoiPlotWidget->m_PlotPicker, SIGNAL(moved(const QPointF&)), SLOT(Clicked(const QPointF&) ) ); connect( (QObject*)(m_Controls->m_Cut), SIGNAL(clicked()), this, SLOT(Cut()) ); connect( (QObject*)(m_Controls->m_Average), SIGNAL(stateChanged(int)), this, SLOT(PerformChange()) ); connect( (QObject*)(m_Controls->m_Segments), SIGNAL(valueChanged(int)), this, SLOT(PerformChange()) ); } } void QmitkTractbasedSpatialStatisticsView::CopyToClipboard() { if(m_Controls->m_RoiPlotWidget->IsPlottingFiber()) { // Working with fiber bundles std::vector > profiles = m_Controls->m_RoiPlotWidget->GetIndividualProfiles(); QString clipboardText; for (std::vector >::iterator it = profiles.begin(); it != profiles.end(); ++it) { for (std::vector::iterator it2 = (*it).begin(); it2 != (*it).end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); } clipboardText.append(QString("\n")); } if(m_Controls->m_Average->isChecked()) { std::vector averages = m_Controls->m_RoiPlotWidget->GetAverageProfile(); clipboardText.append(QString("\nAverage\n")); for (std::vector::iterator it2 = averages.begin(); it2 != averages.end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); } } QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); } else{ // Working with TBSS Data if(m_Controls->m_Average->isChecked()) { std::vector > vals = m_Controls->m_RoiPlotWidget->GetVals(); QString clipboardText; for (std::vector >::iterator it = vals.begin(); it != vals.end(); ++it) { for (std::vector::iterator it2 = (*it).begin(); it2 != (*it).end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); double d = *it2; std::cout << d <setText(clipboardText, QClipboard::Clipboard); } else { std::vector > vals = m_Controls->m_RoiPlotWidget->GetIndividualProfiles(); QString clipboardText; for (std::vector >::iterator it = vals.begin(); it != vals.end(); ++it) { for (std::vector::iterator it2 = (*it).begin(); it2 != (*it).end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); double d = *it2; std::cout << d <setText(clipboardText, QClipboard::Clipboard); } } } void QmitkTractbasedSpatialStatisticsView::RemoveGroup() { QTableView *temp = static_cast(m_Controls->m_GroupInfo); QItemSelectionModel *selectionModel = temp->selectionModel(); QModelIndexList indices = selectionModel->selectedRows(); QModelIndex index; foreach(index, indices) { int row = index.row(); m_GroupModel->removeRows(row, 1, QModelIndex()); } } void QmitkTractbasedSpatialStatisticsView::AddGroup() { QString group("Group"); int number = 0; QPair pair(group, number); QList< QPair >list = m_GroupModel->getList(); if(!list.contains(pair)) { m_GroupModel->insertRows(0, 1, QModelIndex()); QModelIndex index = m_GroupModel->index(0, 0, QModelIndex()); m_GroupModel->setData(index, group, Qt::EditRole); index = m_GroupModel->index(0, 1, QModelIndex()); m_GroupModel->setData(index, number, Qt::EditRole); } } void QmitkTractbasedSpatialStatisticsView::TbssImport() { // Read groups from the interface mitk::TbssImporter::Pointer importer = mitk::TbssImporter::New(); QList< QPair >list = m_GroupModel->getList(); if(list.size() == 0) { QMessageBox msgBox; msgBox.setText("No study group information has been set yet."); msgBox.exec(); return; } std::vector < std::pair > groups; for(int i=0; i pair = list.at(i); std::string s = pair.first.toStdString(); int n = pair.second; std::pair p; p.first = s; p.second = n; groups.push_back(p); } importer->SetGroupInfo(groups); std::string minfo = m_Controls->m_MeasurementInfo->text().toStdString(); importer->SetMeasurementInfo(minfo); std::string name = ""; QList nodes = this->GetDataManagerSelection(); for (auto node: nodes) { if(QString("Image").compare(node->GetData()->GetNameOfClass())==0) { mitk::Image* img = static_cast(node->GetData()); if(img->GetDimension() == 4) { importer->SetImportVolume(img); name = node->GetName(); } } } mitk::TbssImage::Pointer tbssImage; tbssImage = importer->Import(); name += "_tbss"; AddTbssToDataStorage(tbssImage, name); } void QmitkTractbasedSpatialStatisticsView::AddTbssToDataStorage(mitk::Image* image, std::string name) { mitk::LevelWindow levelwindow; levelwindow.SetAuto( image ); mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); levWinProp->SetLevelWindow( levelwindow ); mitk::DataNode::Pointer result = mitk::DataNode::New(); result->SetProperty( "name", mitk::StringProperty::New(name) ); result->SetData( image ); result->SetProperty( "levelwindow", levWinProp ); // add new image to data storage and set as active to ease further processing GetDataStorage()->Add( result ); // show the results mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkTractbasedSpatialStatisticsView::Clicked(const QPointF& pos) { int index = (int)pos.x(); if(m_Roi.size() > 0 && m_CurrentGeometry != nullptr && !m_Controls->m_RoiPlotWidget->IsPlottingFiber() ) { index = std::min( (int)m_Roi.size()-1, std::max(0, index) ); itk::Index<3> ix = m_Roi.at(index); mitk::Vector3D i; i[0] = ix[0]; i[1] = ix[1]; i[2] = ix[2]; mitk::Vector3D w; m_CurrentGeometry->IndexToWorld(i, w); mitk::Point3D origin = m_CurrentGeometry->GetOrigin(); mitk::Point3D p; p[0] = w[0] + origin[0]; p[1] = w[1] + origin[1]; p[2] = w[2] + origin[2]; this->GetRenderWindowPart()->SetSelectedPosition(p); m_Controls->m_RoiPlotWidget->drawBar(index); } else if(m_Controls->m_RoiPlotWidget->IsPlottingFiber() ) { mitk::Point3D point = m_Controls->m_RoiPlotWidget->GetPositionInWorld(index); this->GetRenderWindowPart()->SetSelectedPosition(point); } } void QmitkTractbasedSpatialStatisticsView::Cut() { mitk::BaseData* fibData = m_CurrentFiberNode->GetData(); mitk::FiberBundle* fib = static_cast(fibData); mitk::PlaneGeometry* startGeometry2D = const_cast(dynamic_cast(m_CurrentStartRoi->GetData())->GetPlaneGeometry()); mitk::PlaneGeometry* endGeometry2D = const_cast(dynamic_cast(m_CurrentEndRoi->GetData())->GetPlaneGeometry()); mitk::Point3D startCenter = dynamic_cast(m_CurrentStartRoi->GetData())->GetWorldControlPoint(0); //center Point of start roi mitk::Point3D endCenter = dynamic_cast(m_CurrentEndRoi->GetData())->GetWorldControlPoint(0); //center Point of end roi mitk::FiberBundle::Pointer inStart = fib->ExtractFiberSubset(m_CurrentStartRoi, nullptr); mitk::FiberBundle::Pointer inBoth = inStart->ExtractFiberSubset(m_CurrentEndRoi, nullptr); int num = inBoth->GetNumFibers(); vtkSmartPointer fiberPolyData = inBoth->GetFiberPolyData(); vtkCellArray* lines = fiberPolyData->GetLines(); lines->InitTraversal(); // initialize new vtk polydata vtkSmartPointer points = vtkSmartPointer::New(); vtkSmartPointer polyData = vtkSmartPointer::New(); vtkSmartPointer cells = vtkSmartPointer::New(); int pointIndex=0; // find start and endpoint for( int fiberID( 0 ); fiberID < num; fiberID++ ) { vtkIdType numPointsInCell(0); - vtkIdType* pointsInCell(nullptr); + vtkIdType const* pointsInCell(nullptr); lines->GetNextCell ( numPointsInCell, pointsInCell ); int startId = 0; int endId = 0; float minDistStart = std::numeric_limits::max(); float minDistEnd = std::numeric_limits::max(); vtkSmartPointer polyLine = vtkSmartPointer::New(); int lineIndex=0; for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) { double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); mitk::Point3D point; point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; float distanceToStart = point.EuclideanDistanceTo(startCenter); float distanceToEnd = point.EuclideanDistanceTo(endCenter); if(distanceToStart < minDistStart) { minDistStart = distanceToStart; startId = pointInCellID; } if(distanceToEnd < minDistEnd) { minDistEnd = distanceToEnd; endId = pointInCellID; } } /* We found the start and end points of of the part that should be plottet for the current fiber. now we need to plot them. If the endId is smaller than the startId the plot order must be reversed*/ if(startId < endId) { double *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); mitk::Vector3D p0; p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ startId+1 ] ); mitk::Vector3D p1; p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; // Check if p and p2 are both on the same side of the plane mitk::Vector3D normal = startGeometry2D->GetNormal(); mitk::Point3D pStart; pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool secondOnPositive = startGeometry2D->IsAbove(pSecond); // Calculate intersection with the plane mitk::Vector3D onPlane; onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId+1 lie on the same side of the plane, so we need need startId-1 to calculate the intersection with the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; } double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; double insertPoint[3]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; // First insert the intersection with the start roi points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; if(! (secondOnPositive ^ startOnPositive) ) { /* StartId and startId+1 lie on the same side of the plane so startId is also part of the ROI*/ double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); points->InsertNextPoint(start); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Insert the rest up and to including endId-1 for( int pointInCellID( startId+1 ); pointInCellID < endId ; pointInCellID++) { // create new polyline for new polydata double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); points->InsertNextPoint(p); // add point to line polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } /* endId must be included if endId and endId-1 lie on the same side of the plane defined by endRoi*/ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; mitk::Point3D pLast; pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; mitk::Point3D pBeforeLast; pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; bool lastOnPositive = endGeometry2D->IsAbove(pLast); bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); normal = endGeometry2D->GetNormal(); onPlane[0] = endCenter[0]; onPlane[1] = endCenter[1]; onPlane[2] = endCenter[2]; if(! (lastOnPositive ^ secondLastOnPositive) ) { /* endId and endId-1 lie on the same side of the plane, so we need need endId+1 to calculate the intersection with the planar figure. this should exist since we know that the fiber crosses the planar figure endId is part of the roi so can also be included here*/ p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; double *end = fiberPolyData->GetPoint( pointsInCell[endId] ); points->InsertNextPoint(end); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; //Insert the Last Point (intersection with the end roi) points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Need to reverse walking order else{ double *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); mitk::Vector3D p0; p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); mitk::Vector3D p1; p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; // Check if p and p2 are both on the same side of the plane mitk::Vector3D normal = startGeometry2D->GetNormal(); mitk::Point3D pStart; pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool secondOnPositive = startGeometry2D->IsAbove(pSecond); // Calculate intersection with the plane mitk::Vector3D onPlane; onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId-1 lie on the same side of the plane, so we need need startId+1 to calculate the intersection with the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ startId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; } double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; double insertPoint[3]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; // First insert the intersection with the start roi points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId-1 lie on the same side of the plane so endId is also part of the ROI*/ double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); points->InsertNextPoint(start); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Insert the rest up and to including endId-1 for( int pointInCellID( startId-1 ); pointInCellID > endId ; pointInCellID--) { // create new polyline for new polydata double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); points->InsertNextPoint(p); // add point to line polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } /* startId must be included if startId and startId+ lie on the same side of the plane defined by endRoi*/ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; mitk::Point3D pLast; pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; bool lastOnPositive = endGeometry2D->IsAbove(pLast); mitk::Point3D pBeforeLast; pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); onPlane[0] = endCenter[0]; onPlane[1] = endCenter[1]; onPlane[2] = endCenter[2]; if(! (lastOnPositive ^ secondLastOnPositive) ) { /* endId and endId+1 lie on the same side of the plane, so we need need endId-1 to calculate the intersection with the planar figure. this should exist since we know that the fiber crosses the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; /* endId and endId+1 lie on the same side of the plane so startId is also part of the ROI*/ double *end = fiberPolyData->GetPoint( pointsInCell[endId] ); points->InsertNextPoint(end); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; //Insert the Last Point (intersection with the end roi) points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // add polyline to vtk cell array cells->InsertNextCell(polyLine); } // Add the points to the dataset polyData->SetPoints(points); // Add the lines to the dataset polyData->SetLines(cells); mitk::FiberBundle::Pointer cutBundle = mitk::FiberBundle::New(polyData); mitk::DataNode::Pointer cutNode = mitk::DataNode::New(); cutNode->SetData(cutBundle); std::string name = "fiberCut"; cutNode->SetName(name); GetDataStorage()->Add(cutNode); } void QmitkTractbasedSpatialStatisticsView::CreateRoi() { bool ok; double threshold = QInputDialog::getDouble(m_Controls->m_CreateRoi, tr("Set an FA threshold"), tr("Threshold:"), 0.2, 0.0, 1.0, 2, &ok); if(!ok) return; mitk::Image::Pointer image; QList nodes = this->GetDataManagerSelection(); for (auto node: nodes) { if(QString("Image").compare(node->GetData()->GetNameOfClass())==0) { mitk::Image* img = static_cast(node->GetData()); if(img->GetDimension() == 3) { image = img; } } } if(image.IsNull()) { return; } mitk::TractAnalyzer analyzer; analyzer.SetInputImage(image); analyzer.SetThreshold(threshold); m_PointSetNode = this->m_Controls->m_PointWidget->GetPointSet(); // Set Pointset to analyzer analyzer.SetPointSet(m_PointSetNode); // Run Analyzer try { analyzer.MakeRoi(); } catch (const mitk::Exception& e) { QMessageBox msgBox; msgBox.setText(QString::fromStdString(e.what())); msgBox.exec(); } // Obtain tbss roi image from analyzer mitk::TbssRoiImage::Pointer tbssRoi = analyzer.GetRoiImage(); tbssRoi->SetStructure(m_Controls->m_Structure->text().toStdString()); // get path description and set to interface std::string pathDescription = analyzer.GetPathDescription(); m_Controls->m_PathTextEdit->setPlainText(QString(pathDescription.c_str())); // Add roi image to datastorage AddTbssToDataStorage(tbssRoi, m_Controls->m_RoiName->text().toStdString()); } void QmitkTractbasedSpatialStatisticsView::PlotFiber4D(mitk::TbssImage* image, mitk::FiberBundle* fib, mitk::DataNode* startRoi, mitk::DataNode* endRoi) { if(m_Controls->m_TabWidget->currentWidget() == m_Controls->m_MeasureTAB) { m_CurrentGeometry = image->GetGeometry(); m_Controls->m_RoiPlotWidget->SetGroups(image->GetGroupInfo()); m_Controls->m_RoiPlotWidget->SetProjections(image->GetImage()); m_Controls->m_RoiPlotWidget->SetMeasure( image->GetMeasurementInfo() ); m_Controls->m_RoiPlotWidget->PlotFiber4D(image, fib, startRoi, endRoi, m_Controls->m_Segments->value()); } } void QmitkTractbasedSpatialStatisticsView:: PlotFiberBundle(mitk::FiberBundle *fib, mitk::Image* img, mitk::DataNode* startRoi, mitk::DataNode* endRoi) { bool avg = m_Controls->m_Average->isChecked(); int segments = m_Controls->m_Segments->value(); m_Controls->m_RoiPlotWidget->PlotFiberBetweenRois(fib, img, startRoi ,endRoi, avg, segments); m_Controls->m_RoiPlotWidget->SetPlottingFiber(true); mitk::RenderingManager::GetInstance()->ForceImmediateUpdateAll(); } void QmitkTractbasedSpatialStatisticsView::Plot(mitk::TbssImage* image, mitk::TbssRoiImage* roiImage) { if(m_Controls->m_TabWidget->currentWidget() == m_Controls->m_MeasureTAB) { std::vector< itk::Index<3> > roi = roiImage->GetRoi(); m_Roi = roi; m_CurrentGeometry = image->GetGeometry(); std::string structure = roiImage->GetStructure(); m_Controls->m_RoiPlotWidget->SetGroups(image->GetGroupInfo()); m_Controls->m_RoiPlotWidget->SetProjections(image->GetImage()); m_Controls->m_RoiPlotWidget->SetRoi(roi); m_Controls->m_RoiPlotWidget->SetStructure(structure); m_Controls->m_RoiPlotWidget->SetMeasure( image->GetMeasurementInfo() ); m_Controls->m_RoiPlotWidget->DrawProfiles(); } m_Controls->m_RoiPlotWidget->SetPlottingFiber(false); }