diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h index 14ffcc9abf..9bea0276a6 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h @@ -1,102 +1,102 @@ /*=================================================================== /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkDiffusionTensorPrincipalDirectionImageFilter_h_ #define __itkDiffusionTensorPrincipalDirectionImageFilter_h_ #include "DiffusionCoreExports.h" #include #include #include #include #include #include #include #include namespace itk{ /** \brief Extracts principal eigenvectors of the input tensors */ template< class TTensorPixelType, class TPDPixelType=float> class DiffusionTensorPrincipalDirectionImageFilter : public ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< Vector< TPDPixelType, 3 >, 3 > > { public: typedef DiffusionTensorPrincipalDirectionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< Vector< TPDPixelType, 3 >, 3 > > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(DiffusionTensorPrincipalDirectionImageFilter, ImageToImageFilter) typedef TTensorPixelType TensorComponentType; typedef TPDPixelType DirectionPixelType; typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef itk::Image ItkUcharImgType; typedef vnl_vector_fixed< double, 3 > DirectionType; void SetImage( const InputImageType *image ); // input itkSetMacro( MaskImage, ItkUcharImgType::Pointer) itkSetMacro( NormalizeVectors, bool) // output itkGetMacro( OutputFiberBundle, mitk::FiberBundleX::Pointer) itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer) protected: DiffusionTensorPrincipalDirectionImageFilter(); ~DiffusionTensorPrincipalDirectionImageFilter() {} void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); - void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int ); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType ); void AfterThreadedGenerateData(); private: bool m_NormalizeVectors; ///< Normalizes the output vector to length 1 mitk::FiberBundleX::Pointer m_OutputFiberBundle; ///< Vector field representation of the output vectors ItkUcharImgType::Pointer m_NumDirectionsImage; ///< Image containing the number of fiber directions per voxel ItkUcharImgType::Pointer m_MaskImage; ///< Extraction is only performed inside of the binary mask float m_MaxEigenvalue; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionTensorPrincipalDirectionImageFilter.txx" #endif #endif //__itkDiffusionTensorPrincipalDirectionImageFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.txx index 296185dd4f..ee854f5c8d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.txx @@ -1,236 +1,236 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkDiffusionTensorPrincipalDirectionImageFilter_txx #define __itkDiffusionTensorPrincipalDirectionImageFilter_txx #include #include #include #include "itkDiffusionTensorPrincipalDirectionImageFilter.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionIterator.h" #include "itkArray.h" #include "vnl/vnl_vector.h" #include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { //#define QBALL_RECON_PI M_PI template< class TTensorPixelType, class TPDPixelType> DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType, TPDPixelType> ::DiffusionTensorPrincipalDirectionImageFilter() : m_NormalizeVectors(true) , m_MaxEigenvalue(0.0) { this->SetNumberOfRequiredInputs( 1 ); this->SetNumberOfOutputs(1); } template< class TTensorPixelType, class TPDPixelType> void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType, TPDPixelType> ::BeforeThreadedGenerateData() { typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); Vector spacing = inputImagePointer->GetSpacing(); mitk::Point3D origin = inputImagePointer->GetOrigin(); itk::Matrix direction = inputImagePointer->GetDirection(); ImageRegion<3> imageRegion = inputImagePointer->GetLargestPossibleRegion(); if (m_MaskImage.IsNull()) { m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( spacing ); m_MaskImage->SetOrigin( origin ); m_MaskImage->SetDirection( direction ); m_MaskImage->SetRegions( imageRegion ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } m_NumDirectionsImage = ItkUcharImgType::New(); m_NumDirectionsImage->SetSpacing( spacing ); m_NumDirectionsImage->SetOrigin( origin ); m_NumDirectionsImage->SetDirection( direction ); m_NumDirectionsImage->SetRegions( imageRegion ); m_NumDirectionsImage->Allocate(); m_NumDirectionsImage->FillBuffer(0); itk::Vector< TPDPixelType, 3 > nullVec; nullVec.Fill(0.0); typename OutputImageType::Pointer outputImage = OutputImageType::New(); outputImage->SetSpacing( spacing ); outputImage->SetOrigin( origin ); outputImage->SetDirection( direction ); outputImage->SetRegions( imageRegion ); outputImage->Allocate(); outputImage->FillBuffer(nullVec); this->SetNthOutput(0, outputImage); } template< class TTensorPixelType, class TPDPixelType> void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType, TPDPixelType> ::AfterThreadedGenerateData() { vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); typename OutputImageType::Pointer directionImage = static_cast< OutputImageType* >( this->ProcessObject::GetPrimaryOutput() ); ImageRegionConstIterator< OutputImageType > it(directionImage, directionImage->GetLargestPossibleRegion() ); mitk::Vector3D spacing = directionImage->GetSpacing(); double minSpacing = spacing[0]; if (spacing[1]GetPixel(index)==0) { ++it; continue; } itk::Vector< float, 3 > pixel = directionImage->GetPixel(index); DirectionType dir; dir[0] = pixel[0]; dir[1] = pixel[1]; dir[2] = pixel[2]; if (!m_NormalizeVectors && m_MaxEigenvalue>0) dir /= m_MaxEigenvalue; vtkSmartPointer container = vtkSmartPointer::New(); itk::ContinuousIndex center; center[0] = index[0]; center[1] = index[1]; center[2] = index[2]; itk::Point worldCenter; directionImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); itk::Point worldStart; worldStart[0] = worldCenter[0]-dir[0]/2 * minSpacing; worldStart[1] = worldCenter[1]-dir[1]/2 * minSpacing; worldStart[2] = worldCenter[2]-dir[2]/2 * minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+dir[0]/2 * minSpacing; worldEnd[1] = worldCenter[1]+dir[1]/2 * minSpacing; worldEnd[2] = worldCenter[2]+dir[2]/2 * minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); ++it; } vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); directionsPolyData->SetPoints(m_VtkPoints); directionsPolyData->SetLines(m_VtkCellArray); m_OutputFiberBundle = mitk::FiberBundleX::New(directionsPolyData); } template< class TTensorPixelType, class TPDPixelType> void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType, TPDPixelType> -::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int ) +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType ) { typedef itk::DiffusionTensor3D TensorType; typedef ImageRegionConstIterator< InputImageType > InputIteratorType; typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); - typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput()); + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); ImageRegionIterator< OutputImageType > outIt(outputImage, outputRegionForThread); InputIteratorType inIt(inputImagePointer, outputRegionForThread ); ImageRegionIterator< ItkUcharImgType > nIt(m_NumDirectionsImage, outputRegionForThread ); while( !inIt.IsAtEnd() ) { typename InputImageType::IndexType index = inIt.GetIndex(); if (m_MaskImage->GetPixel(index)==0) { ++inIt; ++nIt; ++outIt; continue; } typename InputImageType::PixelType b = inIt.Get(); TensorType tensor = b.GetDataPointer(); typename OutputImageType::PixelType dir; typename TensorType::EigenValuesArrayType eigenvalues; typename TensorType::EigenVectorsMatrixType eigenvectors; if(tensor.GetTrace()!=0) { tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); vnl_vector_fixed vec; vec[0] = eigenvectors(2,0); vec[1] = eigenvectors(2,1); vec[2] = eigenvectors(2,2); if (!m_NormalizeVectors) vec *= eigenvalues[2]; if (eigenvalues[2]>m_MaxEigenvalue) m_MaxEigenvalue = eigenvalues[2]; dir[0] = (TPDPixelType)vec[0]; dir[1] = (TPDPixelType)vec[1]; dir[2] = (TPDPixelType)vec[2]; outIt.Set( dir ); m_NumDirectionsImage->SetPixel(index, 1); } ++outIt; ++inIt; ++nIt; } std::cout << "One Thread finished extraction" << std::endl; } template< class TTensorPixelType, class TPDPixelType> void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType, TPDPixelType> ::PrintSelf(std::ostream& os, Indent indent) const { } } #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp index 6c93d53412..406d1f73a4 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp @@ -1,992 +1,994 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPartialVolumeAnalysisClusteringCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageToItk.h" #include "itkScalarImageToHistogramGenerator.h" #include "itkListSample.h" #define _USE_MATH_DEFINES #include #define PVA_PI M_PI namespace mitk { PartialVolumeAnalysisClusteringCalculator::PartialVolumeAnalysisClusteringCalculator() : m_MaxIt(100), m_StepsNumIntegration(100) { } PartialVolumeAnalysisClusteringCalculator::~PartialVolumeAnalysisClusteringCalculator() { } PartialVolumeAnalysisClusteringCalculator::ParamsType* PartialVolumeAnalysisClusteringCalculator::InitialGuess(HistType h) const { int s = h.xVals.size(); int s30 = 0.3 * s; int s70 = 0.7 * s; ParamsType* params = new ParamsType(); params->means[0] = h.xVals(s30); params->means[1] = h.xVals(s70); params->sigmas[0] = (h.xVals(s-1) - h.xVals(0)) / 5.0; params->sigmas[1] = (h.xVals(s-1) - h.xVals(0)) / 5.0; params->ps[0] = 0.4; params->ps[1] = 0.4; return params; } PartialVolumeAnalysisClusteringCalculator::ParamsType* PartialVolumeAnalysisClusteringCalculator::Cluster(HistType h, ParamsType *initialGuess) const { ParamsType *params = new ParamsType(); params->Initialize(initialGuess); double ll = 9999999999999999999.9; double oll; double sml = (h.xVals[1] - h.xVals[0]) / 1000; int arraysz = h.xVals.size(); for (unsigned int it = 0; it < m_MaxIt; it++) { // wie sehen basisfunktionen zu aktuellen params aus? ClusterResultType curves = CalculateCurves(*params,h.xVals); // summe der basisfunktionen for(int i=0; i2 && oll-ll < arraysz/2*1e-15 ) { break; } for(int j=0; j<2; j++) { VecType array, p(arraysz); array = curves.vals[j]; for(int i=0; ips[j] = 0; params->means[j] = 0; for(int i=0; ips[j] += p(i); params->means[j] += h.xVals(i)*p(i); } params->means[j] /= params->ps[j]; VecType vr = h.xVals; for(int i=0; imeans[j]; } params->sigmas[j] = 0; for(int i=0; isigmas[j] += vr(i)*vr(i)*p(i); } params->sigmas[j] /= params->ps[j]; params->sigmas[j] += sml; } double p3 = 0; for(int i=0; ips[j] = params->ps[j] + 1e-3; sum += params->ps[j]; } sum += p3; for(int j=0; j<2; j++) { params->ps[j] = params->ps[j] / sum; } p3 /= sum; } return params; } void PartialVolumeAnalysisClusteringCalculator::Normalize(ParamsType params, ClusterResultType* curves) const { double sum1=0, sum2=0, sum3=0; int arraysz = curves->vals[0].size(); for(int i=0; ivals[0](i); sum2 += curves->vals[1](i); sum3 += curves->mixedVals[0](i); } sum1 /= params.ps[0]; sum2 /= params.ps[1]; sum3 /= 1.0 -params.ps[0] -params.ps[1]; for(int i=0; ivals[0](i) /= sum1; curves->vals[1](i) /= sum2; curves->mixedVals[0](i) /= sum3; } for(int i=0; icombiVals(i) = curves->mixedVals[0](i) + curves->vals[0](i) + curves->vals[1](i); } } PartialVolumeAnalysisClusteringCalculator::ClusterResultType PartialVolumeAnalysisClusteringCalculator::CalculateCurves(ParamsType params, VecType xVals) const { int arraysz = xVals.size(); ClusterResultType result(arraysz); for( int j=0; j<2; j++) { for(int i=0; i 0.0000000000000001) { itprob.Set( aposteriori ); maxp = aposteriori > maxp ? aposteriori : maxp; } else { itprob.Set(0.0f); } } ++itimage; ++itprob; } itprob = itprob.Begin(); itdisp = itdisp.Begin(); while( !itprob.IsAtEnd() ) { if(itprob.Get()) { typename DisplayImageType::PixelType rgba; rgba.Set(255.0f, 0.0f, 0.0f, 255.0f*(itprob.Get()/maxp)); itdisp.Set( rgba ); } ++itprob; ++itdisp; } outImage1->InitializeByItk(probimage.GetPointer()); outImage1->SetVolume(probimage->GetBufferPointer()); outImage2->InitializeByItk(displayimage.GetPointer()); outImage2->SetVolume(displayimage->GetBufferPointer()); } double* PartialVolumeAnalysisClusteringCalculator::PerformQuantification( mitk::Image::ConstPointer image, mitk::Image::Pointer clusteredImage, mitk::Image::Pointer mask) const { double *retval = new double[2]; AccessFixedDimensionByItk_3( image.GetPointer(), InternalQuantify, 3, clusteredImage.GetPointer(), retval, mask ); return retval; } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisClusteringCalculator::InternalQuantify( const itk::Image< TPixel, VImageDimension > *image, mitk::Image::Pointer clusteredImage, double* retval, mitk::Image::Pointer mask ) const { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > ProbImageType; typedef itk::Image< unsigned char, VImageDimension > MaskImageType; typedef mitk::ImageToItk CastFilterType; typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput( clusteredImage ); castFilter->Update(); typename ProbImageType::Pointer clusterImage = castFilter->GetOutput(); typename MaskImageType::Pointer itkmask = 0; if(mask.IsNotNull()) { typedef mitk::ImageToItk CastFilterType2; typename CastFilterType2::Pointer castFilter2 = CastFilterType2::New(); castFilter2->SetInput( mask ); castFilter2->Update(); itkmask = castFilter2->GetOutput(); } else { itkmask = MaskImageType::New(); itkmask->SetSpacing( clusterImage->GetSpacing() ); // Set the image spacing itkmask->SetOrigin( clusterImage->GetOrigin() ); // Set the image origin itkmask->SetDirection( clusterImage->GetDirection() ); // Set the image direction itkmask->SetRegions( clusterImage->GetLargestPossibleRegion() ); itkmask->Allocate(); itkmask->FillBuffer(1); } itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itprob(clusterImage, clusterImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itmask(itkmask, itkmask->GetLargestPossibleRegion()); itimage = itimage.Begin(); itprob = itprob.Begin(); itmask = itmask.Begin(); double totalProb = 0; double measurement = 0; double error = 0; while( !itimage.IsAtEnd() && !itprob.IsAtEnd() && !itmask.IsAtEnd() ) { double valImag = itimage.Get(); double valProb = itprob.Get(); double valMask = itmask.Get(); typename ProbImageType::PixelType prop = valProb * valMask; totalProb += prop; measurement += valImag * prop; error += valImag * valImag * prop; ++itimage; ++itprob; ++itmask; } measurement = measurement / totalProb; error = error / totalProb; retval[0] = measurement; retval[1] = sqrt( error - measurement*measurement ); } mitk::Image::Pointer PartialVolumeAnalysisClusteringCalculator::CaculateAngularErrorImage( mitk::Image::Pointer comp1, mitk::Image::Pointer comp2, mitk::Image::Pointer probImg) const { // cast input images to itk typedef itk::Image ImageType; typedef mitk::ImageToItk CastType; CastType::Pointer caster = CastType::New(); caster->SetInput(comp1); caster->Update(); ImageType::Pointer comp1Image = caster->GetOutput(); caster = CastType::New(); caster->SetInput(comp2); caster->Update(); ImageType::Pointer comp2Image = caster->GetOutput(); caster = CastType::New(); caster->SetInput(probImg); caster->Update(); ImageType::Pointer probImage = caster->GetOutput(); // figure out maximum probability for fiber class float maxProb = 0; itk::ImageRegionConstIterator itprob(probImage, probImage->GetLargestPossibleRegion()); itprob = itprob.Begin(); while( !itprob.IsAtEnd() ) { maxProb = itprob.Get() > maxProb ? itprob.Get() : maxProb; ++itprob; } // generate a list sample of angles at positions // where the fiber-prob is higher than .2*maxprob typedef float MeasurementType; const unsigned int MeasurementVectorLength = 2; typedef itk::Vector< MeasurementType , MeasurementVectorLength > MeasurementVectorType; typedef itk::Statistics::ListSample< MeasurementVectorType > ListSampleType; ListSampleType::Pointer listSample = ListSampleType::New(); listSample->SetMeasurementVectorSize( MeasurementVectorLength ); itk::ImageRegionIterator it1(comp1Image, comp1Image->GetLargestPossibleRegion()); itk::ImageRegionIterator it2(comp2Image, comp2Image->GetLargestPossibleRegion()); it1 = it1.Begin(); it2 = it2.Begin(); itprob = itprob.Begin(); while( !itprob.IsAtEnd() ) { if(itprob.Get() > 0.2 * maxProb) { MeasurementVectorType mv; mv[0] = ( MeasurementType ) it1.Get(); mv[1] = ( MeasurementType ) it2.Get(); listSample->PushBack(mv); } ++it1; ++it2; ++itprob; } // generate a histogram from the list sample typedef float HistogramMeasurementType; typedef itk::Statistics::Histogram< HistogramMeasurementType, itk::Statistics::DenseFrequencyContainer2 > HistogramType; typedef itk::Statistics::SampleToHistogramFilter< ListSampleType, HistogramType > GeneratorType; GeneratorType::Pointer generator = GeneratorType::New(); - GeneratorType::HistogramType::SizeType size; + GeneratorType::HistogramType::SizeType size(2); size.Fill(30); generator->SetHistogramSize( size ); generator->SetInput( listSample ); generator->SetMarginalScale( 10.0 ); generator->Update(); // look for frequency mode in the histogram GeneratorType::HistogramType::ConstPointer histogram = generator->GetOutput(); GeneratorType::HistogramType::ConstIterator iter = histogram->Begin(); float maxFreq = 0; MeasurementVectorType maxValue; while ( iter != histogram->End() ) { if(iter.GetFrequency() > maxFreq) { maxFreq = iter.GetFrequency(); maxValue[0] = iter.GetMeasurementVector()[0]; maxValue[1] = iter.GetMeasurementVector()[1]; } ++iter; } // generate return image that contains the angular // error of the voxels to the histogram max measurement ImageType::Pointer returnImage = ImageType::New(); returnImage->SetSpacing( comp1Image->GetSpacing() ); // Set the image spacing returnImage->SetOrigin( comp1Image->GetOrigin() ); // Set the image origin returnImage->SetDirection( comp1Image->GetDirection() ); // Set the image direction returnImage->SetRegions( comp1Image->GetLargestPossibleRegion() ); returnImage->Allocate(); itk::ImageRegionConstIterator cit1(comp1Image, comp1Image->GetLargestPossibleRegion()); itk::ImageRegionConstIterator cit2(comp2Image, comp2Image->GetLargestPossibleRegion()); itk::ImageRegionIterator itout(returnImage, returnImage->GetLargestPossibleRegion()); cit1 = cit1.Begin(); cit2 = cit2.Begin(); itout = itout.Begin(); vnl_vector v(3); v[0] = cos( maxValue[0] ) * sin( maxValue[1] ); v[1] = sin( maxValue[0] ) * sin( maxValue[1] ); v[2] = cos( maxValue[1] ); // MITK_INFO << "max vector: " << v; while( !cit1.IsAtEnd() ) { vnl_vector v1(3); v1[0] = cos( cit1.Get() ) * sin( cit2.Get() ); v1[1] = sin( cit1.Get() ) * sin( cit2.Get() ); v1[2] = cos( cit2.Get() ); itout.Set(fabs(angle(v,v1))); // MITK_INFO << "ang_error " << v1 << ": " << fabs(angle(v,v1)); ++cit1; ++cit2; ++itout; } mitk::Image::Pointer retval = mitk::Image::New(); retval->InitializeByItk(returnImage.GetPointer()); retval->SetVolume(returnImage->GetBufferPointer()); return retval; } PartialVolumeAnalysisClusteringCalculator::HelperStructPerformRGBClusteringRetval* PartialVolumeAnalysisClusteringCalculator::PerformRGBQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2) const { HelperStructRGBChannels *rgbChannels = new HelperStructRGBChannels(); HelperStructPerformClusteringRetval *resultr = PerformQuantiles(image, histogram, p2, 999999999.0 ); rgbChannels->r = resultr->clusteredImage; HelperStructPerformClusteringRetval *resultg = PerformQuantiles(image, histogram, -999999999.0, p1 ); rgbChannels->g = resultg->clusteredImage; HelperStructPerformClusteringRetval *resultb = PerformQuantiles(image, histogram, p1, p2 ); rgbChannels->b = resultb->clusteredImage; mitk::Image::Pointer outImage = mitk::Image::New(); switch(rgbChannels->r->GetDimension()) { case 2: InternalGenerateRGB<2>(rgbChannels, outImage); break; case 3: InternalGenerateRGB<3>(rgbChannels, outImage); break; case 4: InternalGenerateRGB<4>(rgbChannels, outImage); break; default: InternalGenerateRGB<3>(rgbChannels, outImage); } HelperStructPerformRGBClusteringRetval *retval = new HelperStructPerformRGBClusteringRetval(); retval->rgbChannels = rgbChannels; retval->rgb = outImage; retval->params = resultr->params; retval->result = resultr->result; retval->hist = resultr->hist; delete resultr; delete resultg; return retval; } PartialVolumeAnalysisClusteringCalculator::HelperStructPerformClusteringRetval* PartialVolumeAnalysisClusteringCalculator::PerformQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2 ) const { HelperStructPerformClusteringRetval *retval = new HelperStructPerformClusteringRetval(); retval->hist = new HistType(); retval->hist->InitByMitkHistogram(histogram); double sum = histogram->GetTotalFrequency(); double* q = new double[2]; q[0] = histogram->Quantile(0, p1); q[1] = histogram->Quantile(0, p2); mitk::Image::Pointer outImage1 = mitk::Image::New(); mitk::Image::Pointer outImage2 = mitk::Image::New(); AccessFixedDimensionByItk_3( image.GetPointer(), InternalGenerateQuantileImage, 3, q, outImage1, outImage2); retval->clusteredImage = outImage1; retval->displayImage = outImage2; delete[] q; return retval; } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisClusteringCalculator::InternalGenerateQuantileImage( const itk::Image< TPixel, VImageDimension > *image, double* q, mitk::Image::Pointer outImage1, mitk::Image::Pointer outImage2 ) const { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< itk::RGBAPixel, VImageDimension > DisplayImageType; typedef itk::Image< float, VImageDimension > ProbImageType; typename ProbImageType::Pointer probimage = ProbImageType::New(); probimage->SetSpacing( image->GetSpacing() ); // Set the image spacing probimage->SetOrigin( image->GetOrigin() ); // Set the image origin probimage->SetDirection( image->GetDirection() ); // Set the image direction probimage->SetRegions( image->GetLargestPossibleRegion() ); probimage->Allocate(); probimage->FillBuffer(0); typename DisplayImageType::Pointer displayimage = DisplayImageType::New(); displayimage->SetSpacing( image->GetSpacing() ); // Set the image spacing displayimage->SetOrigin( image->GetOrigin() ); // Set the image origin displayimage->SetDirection( image->GetDirection() ); // Set the image direction displayimage->SetRegions( image->GetLargestPossibleRegion() ); displayimage->Allocate(); typename DisplayImageType::PixelType rgba; rgba.Set(0.0f, 0.0f, 0.0f, 0.0f); displayimage->FillBuffer(rgba); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itk::ImageRegionIterator itprob(probimage, probimage->GetLargestPossibleRegion()); itk::ImageRegionIterator itdisp(displayimage, displayimage->GetLargestPossibleRegion()); itimage = itimage.Begin(); itprob = itprob.Begin(); while( !itimage.IsAtEnd() ) { if(itimage.Get() > q[0] && itimage.Get() < q[1]) { itprob.Set(1.0f); } ++itimage; ++itprob; } itprob = itprob.Begin(); itdisp = itdisp.Begin(); while( !itprob.IsAtEnd() ) { if(itprob.Get()) { typename DisplayImageType::PixelType rgba; rgba.Set(255.0f, 0.0f, 0.0f, 255.0f); itdisp.Set( rgba ); } ++itprob; ++itdisp; } outImage1->InitializeByItk(probimage.GetPointer()); outImage1->SetVolume(probimage->GetBufferPointer()); outImage2->InitializeByItk(displayimage.GetPointer()); outImage2->SetVolume(displayimage->GetBufferPointer()); } } diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp index ef029319fb..9a3a5925c3 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp @@ -1,1245 +1,1246 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPartialVolumeAnalysisHistogramCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkExtractImageFilter.h" #include #include #include #include #include #include "itkResampleImageFilter.h" #include "itkGaussianInterpolateImageFunction.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkGaussianInterpolateImageFunction.h" #include "itkBSplineInterpolateImageFunction.h" #include "itkNearestNeighborInterpolateImageFunction.h" #include "itkImageMaskSpatialObject.h" #include "itkRegionOfInterestImageFilter.h" #include "itkListSample.h" #include #include namespace mitk { PartialVolumeAnalysisHistogramCalculator::PartialVolumeAnalysisHistogramCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_NumberOfBins(256), m_UpsamplingFactor(1), m_GaussianSigma(0), m_ForceUpdate(false), m_PlanarFigureThickness(0) { m_EmptyHistogram = HistogramType::New(); - HistogramType::SizeType histogramSize; + m_EmptyHistogram->SetMeasurementVectorSize(1); + HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } PartialVolumeAnalysisHistogramCalculator::~PartialVolumeAnalysisHistogramCalculator() { } void PartialVolumeAnalysisHistogramCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); m_ImageStatisticsTimeStamp.Modified(); m_ImageStatisticsCalculationTriggerBool = true; } } void PartialVolumeAnalysisHistogramCalculator::AddAdditionalResamplingImage( const mitk::Image *image ) { m_AdditionalResamplingImages.push_back(image); this->Modified(); m_ImageStatisticsTimeStamp.Modified(); m_ImageStatisticsCalculationTriggerBool = true; } void PartialVolumeAnalysisHistogramCalculator::SetModified( ) { this->Modified(); m_ImageStatisticsTimeStamp.Modified(); m_ImageStatisticsCalculationTriggerBool = true; m_MaskedImageStatisticsTimeStamp.Modified(); m_MaskedImageStatisticsCalculationTriggerBool = true; m_PlanarFigureStatisticsTimeStamp.Modified(); m_PlanarFigureStatisticsCalculationTriggerBool = true; } void PartialVolumeAnalysisHistogramCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) { itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); m_MaskedImageStatisticsTimeStamp.Modified(); m_MaskedImageStatisticsCalculationTriggerBool = true; } } void PartialVolumeAnalysisHistogramCalculator::SetPlanarFigure( const mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); m_PlanarFigureStatisticsTimeStamp.Modified(); m_PlanarFigureStatisticsCalculationTriggerBool = true; } } void PartialVolumeAnalysisHistogramCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } bool PartialVolumeAnalysisHistogramCalculator::ComputeStatistics() { MITK_DEBUG << "ComputeStatistics() start"; if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image not set!" ); } if ( m_Image->GetReferenceCount() == 1 ) { MITK_DEBUG << "No Stats calculated; no one else holds a reference on it"; return false; } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = NULL; } // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStamp.GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStamp.GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStamp.GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerBool; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerBool; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerBool; if ( /*prevent calculation without mask*/!m_ForceUpdate &&( m_MaskingMode == MASKING_MODE_NONE || ( ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (m_ImageMask.IsNotNull() && maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (m_PlanarFigure.IsNotNull() && planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) ) ) { MITK_DEBUG << "Returning, statistics already up to date!"; if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; return true; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( ); Statistics *statistics; HistogramType::ConstPointer *histogram; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: statistics = &m_ImageStatistics; histogram = &m_ImageHistogram; m_ImageStatisticsTimeStamp.Modified(); m_ImageStatisticsCalculationTriggerBool = false; break; case MASKING_MODE_IMAGE: statistics = &m_MaskedImageStatistics; histogram = &m_MaskedImageHistogram; m_MaskedImageStatisticsTimeStamp.Modified(); m_MaskedImageStatisticsCalculationTriggerBool = false; break; case MASKING_MODE_PLANARFIGURE: statistics = &m_PlanarFigureStatistics; histogram = &m_PlanarFigureHistogram; m_PlanarFigureStatisticsTimeStamp.Modified(); m_PlanarFigureStatisticsCalculationTriggerBool = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE ) { // Reset state changed flag AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, *statistics, histogram ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), *statistics, histogram ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, *statistics, histogram ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), *statistics, histogram ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory // m_InternalImage = mitk::Image::Pointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } const PartialVolumeAnalysisHistogramCalculator::HistogramType * PartialVolumeAnalysisHistogramCalculator::GetHistogram( ) const { if ( m_Image.IsNull() ) { return NULL; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: return m_ImageHistogram; case MASKING_MODE_IMAGE: return m_MaskedImageHistogram; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogram; } } const PartialVolumeAnalysisHistogramCalculator::Statistics & PartialVolumeAnalysisHistogramCalculator::GetStatistics( ) const { if ( m_Image.IsNull() ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: return m_ImageStatistics; case MASKING_MODE_IMAGE: return m_MaskedImageStatistics; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatistics; } } void PartialVolumeAnalysisHistogramCalculator::ExtractImageAndMask( ) { MITK_DEBUG << "ExtractImageAndMask( ) start"; if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } mitk::Image *timeSliceImage = const_cast(m_Image.GetPointer());//imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; int s = m_AdditionalResamplingImages.size(); m_InternalAdditionalResamplingImages.resize(s); for(int i=0; i(m_AdditionalResamplingImages[i].GetPointer()); } m_InternalImageMask2D = NULL; m_InternalImageMask3D = NULL; break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( 0 ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); InternalMaskImage(timeSliceMaskedImage); if(m_UpsamplingFactor != 1) { InternalResampleImage(m_InternalImageMask3D); } AccessFixedDimensionByItk_1( timeSliceImage, InternalResampleImageFromMask, 3, -1 ); int s = m_AdditionalResamplingImages.size(); m_InternalAdditionalResamplingImages.resize(s); for(int i=0; iIsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const Geometry3D *imageGeometry = timeSliceImage->GetUpdatedGeometry(); if ( imageGeometry == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); if ( planarFigureGeometry == NULL ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } // unsigned int axis = 2; // unsigned int slice = 0; AccessFixedDimensionByItk_3( timeSliceImage, InternalReorientImagePlane, 3, timeSliceImage->GetGeometry(), m_PlanarFigure->GetGeometry(), -1 ); AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2 ); int s = m_AdditionalResamplingImages.size(); for(int i=0; iGetGeometry(), m_PlanarFigure->GetGeometry(), i ); AccessFixedDimensionByItk_1( m_InternalAdditionalResamplingImages[i], InternalCropAdditionalImage, 3, i ); } } } } bool PartialVolumeAnalysisHistogramCalculator::GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } void PartialVolumeAnalysisHistogramCalculator::InternalMaskImage( mitk::Image *image ) { typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject; typedef itk::Image< unsigned char, 3 > ImageType; typedef itk::ImageRegion<3> RegionType; typedef mitk::ImageToItk CastType; CastType::Pointer caster = CastType::New(); caster->SetInput(image); caster->Update(); ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New(); maskSO->SetImage ( caster->GetOutput() ); m_InternalMask3D = maskSO->GetAxisAlignedBoundingBoxRegion(); MITK_DEBUG << "Bounding Box Region: " << m_InternalMask3D; typedef itk::RegionOfInterestImageFilter< ImageType, MaskImage3DType > ROIFilterType; ROIFilterType::Pointer roi = ROIFilterType::New(); roi->SetRegionOfInterest(m_InternalMask3D); roi->SetInput(caster->GetOutput()); roi->Update(); m_InternalImageMask3D = roi->GetOutput(); MITK_DEBUG << "Created m_InternalImageMask3D"; } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::Geometry3D* imggeo, mitk::Geometry3D* planegeo3D, int additionalIndex ) { MITK_DEBUG << "InternalReorientImagePlane() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > FloatImageType; typedef itk::ResampleImageFilter ResamplerType; typename ResamplerType::Pointer resampler = ResamplerType::New(); mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); float upsamp = m_UpsamplingFactor; float gausssigma = m_GaussianSigma; // Spacing typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); spacing[0] = image->GetSpacing()[0] / upsamp; spacing[1] = image->GetSpacing()[1] / upsamp; spacing[2] = image->GetSpacing()[2]; if(m_PlanarFigureThickness) { spacing[2] = image->GetSpacing()[2] / upsamp; } resampler->SetOutputSpacing( spacing ); // Size typename ResamplerType::SizeType size; size[0] = planegeo->GetParametricExtentInMM(0) / spacing[0]; size[1] = planegeo->GetParametricExtentInMM(1) / spacing[1]; size[2] = 1+2*m_PlanarFigureThickness; // klaus add +2*m_PlanarFigureThickness MITK_DEBUG << "setting size2:="<SetSize( size ); // Origin typename mitk::Point3D orig = planegeo->GetOrigin(); typename mitk::Point3D corrorig; planegeo3D->WorldToIndex(orig,corrorig); corrorig[0] += 0.5/upsamp; corrorig[1] += 0.5/upsamp; if(m_PlanarFigureThickness) { float thickyyy = m_PlanarFigureThickness; thickyyy/=upsamp; corrorig[2] -= thickyyy; // klaus add -= (float)m_PlanarFigureThickness/upsamp statt += 0 } planegeo3D->IndexToWorld(corrorig,corrorig); resampler->SetOutputOrigin(corrorig ); // Direction typename ResamplerType::DirectionType direction; typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); for(int c=0; cSetOutputDirection( direction ); // Gaussian interpolation if(gausssigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) { sigma[d] = gausssigma * image->GetSpacing()[d]; } double alpha = 2.0; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); interpolator->SetInputImage( image ); interpolator->SetParameters( sigma, alpha ); resampler->SetInterpolator( interpolator ); } else { // typedef typename itk::BSplineInterpolateImageFunction // InterpolatorType; typedef typename itk::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } // Other resampling options resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); MITK_DEBUG << "Resampling requested image plane ... "; resampler->Update(); MITK_DEBUG << " ... done"; if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } else { unsigned int myIndex = additionalIndex; this->m_InternalAdditionalResamplingImages.push_back(mitk::Image::New()); this->m_InternalAdditionalResamplingImages[myIndex]->InitializeByItk( resampler->GetOutput() ); this->m_InternalAdditionalResamplingImages[myIndex]->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalResampleImageFromMask( const itk::Image< TPixel, VImageDimension > *image, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typename ImageType::Pointer outImage = ImageType::New(); outImage->SetSpacing( m_InternalImageMask3D->GetSpacing() ); // Set the image spacing outImage->SetOrigin( m_InternalImageMask3D->GetOrigin() ); // Set the image origin outImage->SetDirection( m_InternalImageMask3D->GetDirection() ); // Set the image direction outImage->SetRegions( m_InternalImageMask3D->GetLargestPossibleRegion() ); outImage->Allocate(); outImage->FillBuffer(0); typedef itk::InterpolateImageFunction BaseInterpType; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typedef itk::LinearInterpolateImageFunction LinearInterpolatorType; typename BaseInterpType::Pointer interpolator; if(m_GaussianSigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) { sigma[d] = m_GaussianSigma * image->GetSpacing()[d]; } double alpha = 2.0; interpolator = GaussianInterpolatorType::New(); dynamic_cast(interpolator.GetPointer())->SetParameters( sigma, alpha ); } else { interpolator = LinearInterpolatorType::New(); } interpolator->SetInputImage( image ); itk::ImageRegionConstIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itk::ImageRegionIterator itimage(outImage, outImage->GetLargestPossibleRegion()); itmask = itmask.Begin(); itimage = itimage.Begin(); itk::Point< double, 3 > point; itk::ContinuousIndex< double, 3 > index; while( !itmask.IsAtEnd() ) { if(itmask.Get() != 0) { outImage->TransformIndexToPhysicalPoint (itimage.GetIndex(), point); image->TransformPhysicalPointToContinuousIndex(point, index); itimage.Set(interpolator->EvaluateAtContinuousIndex(index)); } ++itmask; ++itimage; } if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( outImage.GetPointer() ); this->m_InternalImage->SetVolume( outImage->GetBufferPointer() ); } else { this->m_InternalAdditionalResamplingImages[additionalIndex] = mitk::Image::New(); this->m_InternalAdditionalResamplingImages[additionalIndex]->InitializeByItk( outImage.GetPointer() ); this->m_InternalAdditionalResamplingImages[additionalIndex]->SetVolume( outImage->GetBufferPointer() ); } } void PartialVolumeAnalysisHistogramCalculator::InternalResampleImage( const MaskImage3DType *image ) { typedef itk::ResampleImageFilter ResamplerType; ResamplerType::Pointer resampler = ResamplerType::New(); // Size ResamplerType::SizeType size; size[0] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[0]; size[1] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[1]; size[2] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[2];; resampler->SetSize( size ); // Origin mitk::Point3D orig = image->GetOrigin(); resampler->SetOutputOrigin(orig ); // Spacing ResamplerType::SpacingType spacing; spacing[0] = image->GetSpacing()[0] / m_UpsamplingFactor; spacing[1] = image->GetSpacing()[1] / m_UpsamplingFactor; spacing[2] = image->GetSpacing()[2] / m_UpsamplingFactor; resampler->SetOutputSpacing( spacing ); resampler->SetOutputDirection( image->GetDirection() ); typedef itk::NearestNeighborInterpolateImageFunction InterpolatorType; InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); // Other resampling options resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); resampler->Update(); m_InternalImageMask3D = resampler->GetOutput(); } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, Statistics &statistics, typename HistogramType::ConstPointer *histogram ) { MITK_DEBUG << "InternalCalculateStatisticsUnmasked()"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned char, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; // Progress listening... typedef itk::SimpleMemberCommand< PartialVolumeAnalysisHistogramCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &PartialVolumeAnalysisHistogramCalculator::UnmaskedStatisticsProgressUpdate ); // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator // does not (yet?) support progress reporting this->InvokeEvent( itk::StartEvent() ); for ( unsigned int i = 0; i < 100; ++i ) { this->UnmaskedStatisticsProgressUpdate(); } // Calculate statistics (separate filter) typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); statisticsFilter->Update(); statisticsFilter->RemoveObserver( observerTag ); this->InvokeEvent( itk::EndEvent() ); statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); statistics.Min = statisticsFilter->GetMinimum(); statistics.Max = statisticsFilter->GetMaximum(); statistics.Mean = statisticsFilter->GetMean(); statistics.Median = 0.0; statistics.Sigma = statisticsFilter->GetSigma(); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); typename ImageType::Pointer inImage = const_cast(image); // Calculate histogram typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( inImage ); histogramGenerator->SetMarginalScale( 10 ); // Defines y-margin width of histogram histogramGenerator->SetNumberOfBins( m_NumberOfBins ); // CT range [-1024, +2048] --> bin size 4 values histogramGenerator->SetHistogramMin( statistics.Min ); histogramGenerator->SetHistogramMax( statistics.Max ); histogramGenerator->Compute(); *histogram = histogramGenerator->GetOutput(); } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned char, VImageDimension > *maskImage, Statistics &statistics, typename HistogramType::ConstPointer *histogram ) { MITK_DEBUG << "InternalCalculateStatisticsMasked() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned char, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; // generate a list sample of angles at positions // where the fiber-prob is higher than .2*maxprob typedef TPixel MeasurementType; const unsigned int MeasurementVectorLength = 1; typedef itk::Vector< MeasurementType , MeasurementVectorLength > MeasurementVectorType; typedef itk::Statistics::ListSample< MeasurementVectorType > ListSampleType; typename ListSampleType::Pointer listSample = ListSampleType::New(); listSample->SetMeasurementVectorSize( MeasurementVectorLength ); itk::ImageRegionConstIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask = itmask.Begin(); itimage = itimage.Begin(); while( !itmask.IsAtEnd() ) { if(itmask.Get() != 0) { // apend to list MeasurementVectorType mv; mv[0] = ( MeasurementType ) itimage.Get(); listSample->PushBack(mv); } ++itmask; ++itimage; } // generate a histogram from the list sample typedef double HistogramMeasurementType; typedef itk::Statistics::Histogram< HistogramMeasurementType, itk::Statistics::DenseFrequencyContainer2 > HistogramType; typedef itk::Statistics::SampleToHistogramFilter< ListSampleType, HistogramType > GeneratorType; typename GeneratorType::Pointer generator = GeneratorType::New(); - typename GeneratorType::HistogramType::SizeType size; + typename GeneratorType::HistogramType::SizeType size(MeasurementVectorLength); size.Fill(m_NumberOfBins); generator->SetHistogramSize( size ); generator->SetInput( listSample ); generator->SetMarginalScale( 10.0 ); generator->Update(); *histogram = generator->GetOutput(); } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalCropAdditionalImage( itk::Image< TPixel, VImageDimension > *image, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::RegionOfInterestImageFilter< ImageType, ImageType > ROIFilterType; typename ROIFilterType::Pointer roi = ROIFilterType::New(); roi->SetRegionOfInterest(m_CropRegion); roi->SetInput(image); roi->Update(); m_InternalAdditionalResamplingImages[additionalIndex] = mitk::Image::New(); m_InternalAdditionalResamplingImages[additionalIndex]->InitializeByItk(roi->GetOutput()); m_InternalAdditionalResamplingImages[additionalIndex]->SetVolume(roi->GetOutput()->GetBufferPointer()); } template < typename TPixel, unsigned int VImageDimension > void PartialVolumeAnalysisHistogramCalculator::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { MITK_DEBUG << "InternalCalculateMaskFromPlanarFigure() start"; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage3DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". MaskImage3DType::Pointer newMaskImage = MaskImage3DType::New(); newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); newMaskImage->Allocate(); newMaskImage->FillBuffer( 1 ); // Generate VTK polygon from (closed) PlanarFigure polyline // (The polyline points are shifted by -0.5 in z-direction to make sure // that the extrusion filter, which afterwards elevates all points by +0.5 // in z-direction, creates a 3D object which is cut by the the plane z=0) const mitk::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::Geometry3D *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); vtkPolyData *polyline = vtkPolyData::New(); polyline->Allocate( 1, 1 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // Create VTK polydata object of polyline contour bool outOfBounds = false; vtkPoints *points = vtkPoints::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image mitk::Point2D point2D = it->Point; planarFigureGeometry2D->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigureGeometry2D->IndexToWorld(point2D, point2D); planarFigureGeometry2D->Map( point2D, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } imageGeometry3D->WorldToIndex( point3D, point3D ); point3D[i0] += 0.5; point3D[i1] += 0.5; // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); } polyline->SetPoints( points ); points->Delete(); if ( outOfBounds ) { polyline->Delete(); throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } unsigned int numberOfPoints = planarFigurePolyline.size(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( vtkIdType i = 0; i < numberOfPoints; ++i ) { ptIds[i] = i; } polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->SetInput( polyline ); extrudeFilter->SetScaleFactor( 1 ); extrudeFilter->SetExtrusionTypeToNormalExtrusion(); extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); // Make a stencil from the extruded polygon vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); polyDataToImageStencil->SetInput( extrudeFilter->GetOutput() ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage3DType > ImageImportType; typedef itk::VTKImageExport< MaskImage3DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( newMaskImage ); vtkImageImport *vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); vtkImporter->Update(); // Apply the generated image stencil to the input image vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); imageStencilFilter->SetInput( vtkImporter->GetOutput() ); imageStencilFilter->SetStencil( polyDataToImageStencil->GetOutput() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkImageExport *vtkExporter = vtkImageExport::New(); vtkExporter->SetInput( imageStencilFilter->GetOutput() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // calculate cropping bounding box m_InternalImageMask3D = itkImporter->GetOutput(); m_InternalImageMask3D->SetDirection(image->GetDirection()); itk::ImageRegionIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itmask = itmask.Begin(); while( !itmask.IsAtEnd() ) { if(itmask.Get() != 0) { typename ImageType::IndexType index = itmask.GetIndex(); for(int thick=0; thick<2*m_PlanarFigureThickness+1; thick++) { index[axis] = thick; m_InternalImageMask3D->SetPixel(index, itmask.Get()); } } ++itmask; } itmask = itmask.Begin(); itk::ImageRegionIterator itimage(image, image->GetLargestPossibleRegion()); itimage = itimage.Begin(); typename ImageType::SizeType lowersize = {{9999999999.0,9999999999.0,9999999999.0}}; typename ImageType::SizeType uppersize = {{0,0,0}}; while( !itmask.IsAtEnd() ) { if(itmask.Get() == 0) { itimage.Set(0); } else { typename ImageType::IndexType index = itimage.GetIndex(); typename ImageType::SizeType signedindex; signedindex[0] = index[0]; signedindex[1] = index[1]; signedindex[2] = index[2]; lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; } ++itmask; ++itimage; } typename ImageType::IndexType index; index[0] = lowersize[0]; index[1] = lowersize[1]; index[2] = lowersize[2]; typename ImageType::SizeType size; size[0] = uppersize[0] - lowersize[0] + 1; size[1] = uppersize[1] - lowersize[1] + 1; size[2] = uppersize[2] - lowersize[2] + 1; m_CropRegion = itk::ImageRegion<3>(index, size); // crop internal image typedef itk::RegionOfInterestImageFilter< ImageType, ImageType > ROIFilterType; typename ROIFilterType::Pointer roi = ROIFilterType::New(); roi->SetRegionOfInterest(m_CropRegion); roi->SetInput(image); roi->Update(); m_InternalImage = mitk::Image::New(); m_InternalImage->InitializeByItk(roi->GetOutput()); m_InternalImage->SetVolume(roi->GetOutput()->GetBufferPointer()); // crop internal mask typedef itk::RegionOfInterestImageFilter< MaskImage3DType, MaskImage3DType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(m_CropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); // Clean up VTK objects polyline->Delete(); extrudeFilter->Delete(); polyDataToImageStencil->Delete(); vtkImporter->Delete(); imageStencilFilter->Delete(); //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? delete[] ptIds; } void PartialVolumeAnalysisHistogramCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void PartialVolumeAnalysisHistogramCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index bcf0ec3836..31dabfb6ef 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,1217 +1,1217 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkExtractImageFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #if ( ( VTK_MAJOR_VERSION <= 5 ) && ( VTK_MINOR_VERSION<=8) ) #include "mitkvtkLassoStencilSource.h" #else #include "vtkLassoStencilSource.h" #endif #include #include namespace mitk { ImageStatisticsCalculator::ImageStatisticsCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_IgnorePixelValue(0.0), m_DoIgnorePixelValue(false), m_IgnorePixelValueChanged(false), m_PlanarFigureAxis (0), m_PlanarFigureSlice (0), m_PlanarFigureCoordinate0 (0), m_PlanarFigureCoordinate1 (0) { m_EmptyHistogram = HistogramType::New(); - m_EmptyHistogram->SetMeasurementVectorSize(2); - HistogramType::SizeType histogramSize(2); + m_EmptyHistogram->SetMeasurementVectorSize(1); + HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } ImageStatisticsCalculator::~ImageStatisticsCalculator() { } void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); unsigned int numberOfTimeSteps = image->GetTimeSteps(); // Initialize vectors to time-size of this image m_ImageHistogramVector.resize( numberOfTimeSteps ); m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); m_ImageStatisticsVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) { m_ImageStatisticsTimeStampVector[t].Modified(); m_ImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) { itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_MaskedImageStatisticsTimeStampVector[t].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_PlanarFigureStatisticsTimeStampVector[t].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetIgnorePixelValue(double value) { if ( m_IgnorePixelValue != value ) { m_IgnorePixelValue = value; if(m_DoIgnorePixelValue) { m_IgnorePixelValueChanged = true; } this->Modified(); } } double ImageStatisticsCalculator::GetIgnorePixelValue() { return m_IgnorePixelValue; } void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) { if ( m_DoIgnorePixelValue != value ) { m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; this->Modified(); } } bool ImageStatisticsCalculator::GetDoIgnorePixelValue() { return m_DoIgnorePixelValue; } bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) { if (m_Image.IsNull() ) { mitkThrow() << "Image not set!"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if ( m_Image->GetReferenceCount() == 1 ) { // Image no longer valid; we are the only ones to still hold a reference on it return false; } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = NULL; } // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; if ( !m_IgnorePixelValueChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; return true; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; m_IgnorePixelValueChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( timeStep ); StatisticsContainer *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return NULL; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep][label]; return m_ImageHistogramVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep][label]; } } const ImageStatisticsCalculator::HistogramContainer & ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyHistogramContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep]; return m_ImageHistogramVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep][label]; return m_ImageStatisticsVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep][label]; } } const ImageStatisticsCalculator::StatisticsContainer & ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatisticsContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep]; return m_ImageStatisticsVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep]; } } void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( timeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; m_InternalImageMask2D = NULL; m_InternalImageMask3D = NULL; if(m_DoIgnorePixelValue) { if( m_InternalImage->GetDimension() == 3 ) { CastToItkImage( timeSliceImage, m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } if( m_InternalImage->GetDimension() == 2 ) { CastToItkImage( timeSliceImage, m_InternalImageMask2D ); m_InternalImageMask2D->FillBuffer(1); } } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { if ( timeStep < m_ImageMask->GetTimeSteps() ) { ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( timeStep ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); m_InternalImage = timeSliceImage; CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D ); } else { throw std::runtime_error( "Error: image mask has not enough time steps!" ); } } else { throw std::runtime_error( "Error: image mask empty!" ); } break; } case MASKING_MODE_PLANARFIGURE: { m_InternalImageMask2D = NULL; if ( m_PlanarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !m_PlanarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const Geometry3D *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); if ( planarFigureGeometry == NULL ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image MaskImage3DType::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // Extract slice with given position and direction from image ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( timeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); m_InternalImage = imageExtractor->GetOutput(); // Compute mask from PlanarFigure AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 2, axis ); } } if(m_DoIgnorePixelValue) { if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); } } } bool ImageStatisticsCalculator::GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer *statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; statisticsContainer->clear(); histogramContainer->clear(); // Progress listening... typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator // does not (yet?) support progress reporting this->InvokeEvent( itk::StartEvent() ); for ( unsigned int i = 0; i < 100; ++i ) { this->UnmaskedStatisticsProgressUpdate(); } // Calculate statistics (separate filter) typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); statisticsFilter->Update(); statisticsFilter->RemoveObserver( observerTag ); this->InvokeEvent( itk::EndEvent() ); // Calculate minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( image ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); Statistics statistics; statistics.Reset(); statistics.Label = 1; statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); statistics.Min = statisticsFilter->GetMinimum(); statistics.Max = statisticsFilter->GetMaximum(); statistics.Mean = statisticsFilter->GetMean(); statistics.Median = 0.0; statistics.Sigma = statisticsFilter->GetSigma(); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); statistics.MinIndex.set_size(image->GetImageDimension()); statistics.MaxIndex.set_size(image->GetImageDimension()); for (int i=0; iGetIndexOfMaximum()[i]; statistics.MinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } statisticsContainer->push_back( statistics ); // Calculate histogram typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); histogramGenerator->SetNumberOfBins( 768 ); histogramGenerator->SetHistogramMin( statistics.Min ); histogramGenerator->SetHistogramMax( statistics.Max ); histogramGenerator->Compute(); histogramContainer->push_back( histogramGenerator->GetOutput() ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(m_IgnorePixelValue == itimage.Get()) { itmask.Set(0); } ++itmask; ++itimage; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::SpacingType SpacingType; typedef itk::LabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; statisticsContainer->clear(); histogramContainer->clear(); // Make sure that mask is set if ( maskImage == NULL ) { itkExceptionMacro( << "Mask image needs to be set!" ); } // Make sure that spacing of mask and image are the same SpacingType imageSpacing = image->GetSpacing(); SpacingType maskSpacing = maskImage->GetSpacing(); PointType zeroPoint; zeroPoint.Fill( 0.0 ); if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" ); } // Make sure that orientation of mask and image are the same typedef typename ImageType::DirectionType DirectionType; DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); long offset[ImageType::ImageDimension]; typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); if ( fabs( misalignment ) > mitk::eps ) { itkExceptionMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; offset[i] = (int) indexCoordDistance + image->GetBufferedRegion().GetIndex()[i]; } // Adapt the origin and region (index/size) of the mask so that the origin of both are the same typename ChangeInformationFilterType::Pointer adaptMaskFilter; adaptMaskFilter = ChangeInformationFilterType::New(); adaptMaskFilter->ChangeOriginOn(); adaptMaskFilter->ChangeRegionOn(); adaptMaskFilter->SetInput( maskImage ); adaptMaskFilter->SetOutputOrigin( image->GetOrigin() ); adaptMaskFilter->SetOutputOffset( offset ); adaptMaskFilter->Update(); typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); // Make sure that mask region is contained within image region if ( !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) { itkExceptionMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); } // If mask region is smaller than image region, extract the sub-sampled region from the original image typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } typename ImageType::ConstPointer adaptedImage; if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); statisticsFilter->Update(); int numberOfBins = ( m_DoIgnorePixelValue && (m_MaskingMode == MASKING_MODE_NONE) ) ? 768 : 384; typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( numberOfBins, statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum() ); // Add progress listening typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate ); unsigned long observerTag = labelStatisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); // Execute filter this->InvokeEvent( itk::StartEvent() ); // Make sure that only the mask region is considered (otherwise, if the mask region is smaller // than the image region, the Update() would result in an exception). labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); // Execute the filter labelStatisticsFilter->Update(); this->InvokeEvent( itk::EndEvent() ); labelStatisticsFilter->RemoveObserver( observerTag ); // Find all relevant labels of mask (other than 0) std::list< int > relevantLabels; bool maskNonEmpty = false; unsigned int i; for ( i = 1; i < 4096; ++i ) { if ( labelStatisticsFilter->HasLabel( i ) ) { relevantLabels.push_back( i ); maskNonEmpty = true; } } if ( maskNonEmpty ) { std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); Statistics statistics; statistics.Label = (*it); statistics.N = labelStatisticsFilter->GetCount( *it ); statistics.Min = labelStatisticsFilter->GetMinimum( *it ); statistics.Max = labelStatisticsFilter->GetMaximum( *it ); statistics.Mean = labelStatisticsFilter->GetMean( *it ); statistics.Median = labelStatisticsFilter->GetMedian( *it ); statistics.Sigma = labelStatisticsFilter->GetSigma( *it ); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); masker->SetOutsideValue( (statistics.Min+statistics.Max)/2 ); masker->SetInput1(adaptedImage); masker->SetInput2(adaptedMaskImage); masker->Update(); // get index of minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( masker->GetOutput() ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); statistics.MinIndex.set_size(adaptedImage->GetImageDimension()); statistics.MaxIndex.set_size(adaptedImage->GetImageDimension()); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); typename MinMaxFilterType::IndexType tempMinIndex = minMaxFilter->GetIndexOfMinimum(); // FIX BUG 14644 //If a PlanarFigure is used for segmentation the //adaptedImage is a single slice (2D). Adding the // 3. dimension. if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { statistics.MaxIndex.set_size(m_Image->GetDimension()); statistics.MaxIndex[m_PlanarFigureCoordinate0]=tempMaxIndex[0]; statistics.MaxIndex[m_PlanarFigureCoordinate1]=tempMaxIndex[1]; statistics.MaxIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; statistics.MinIndex.set_size(m_Image->GetDimension()); statistics.MinIndex[m_PlanarFigureCoordinate0]=tempMinIndex[0]; statistics.MinIndex[m_PlanarFigureCoordinate1]=tempMinIndex[1]; statistics.MinIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; } else { for (int i = 0; ipush_back( statistics ); } } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); statisticsContainer->push_back( Statistics() );; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput( image ); castFilter->Update(); castFilter->GetOutput()->FillBuffer( 1 ); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::Geometry3D *imageGeometry3D = m_Image->GetGeometry( 0 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } m_PlanarFigureCoordinate0= i0; m_PlanarFigureCoordinate1= i1; // store the polyline contour as vtkPoints object bool outOfBounds = false; vtkSmartPointer points = vtkSmartPointer::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image planarFigureGeometry2D->Map( it->Point, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } imageGeometry3D->WorldToIndex( point3D, point3D ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencil( lassoStencil->GetOutput() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkImageExport::New(); // TODO: this is WRONG, should be vtkSmartPointer::New(), but bug # 14455 vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // Store mask m_InternalImageMask2D = itkImporter->GetOutput(); } void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } }