diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkBrainMaskExtractionImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkBrainMaskExtractionImageFilter.txx index 812ac8ffdf..ed41c65f13 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkBrainMaskExtractionImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkBrainMaskExtractionImageFilter.txx @@ -1,443 +1,443 @@ /*=================================================================== 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 __itkBrainMaskExtractionImageFilter_txx #define __itkBrainMaskExtractionImageFilter_txx #include "itkBrainMaskExtractionImageFilter.h" #include #include #include #include #include #include #include #include #include #include namespace itk { template< class TOutputImagePixelType > BrainMaskExtractionImageFilter< TOutputImagePixelType > ::BrainMaskExtractionImageFilter() { // At least 1 inputs is necessary for a vector image. // For images added one at a time we need at least six this->SetNumberOfRequiredInputs( 1 ); } template< class TOutputImagePixelType > void BrainMaskExtractionImageFilter< TOutputImagePixelType > ::GenerateData() { // gaussian filter typename itk::RecursiveGaussianImageFilter::Pointer gaussian = itk::RecursiveGaussianImageFilter::New(); gaussian->SetInput( this->GetInput(0) ); gaussian->SetSigma( 1.0 ); try { gaussian->Update(); } catch( itk::ExceptionObject &e) { std::cerr << e; return; } // threshold the image typename itk::BinaryThresholdImageFilter::Pointer threshold = itk::BinaryThresholdImageFilter::New(); threshold->SetInput( gaussian->GetOutput() ); int seuil = static_cast( ComputeHistogram( gaussian->GetOutput() ) ); threshold->SetLowerThreshold( seuil ); std::cout << "Thresholding..." << std::flush; try { threshold->Update(); } catch( itk::ExceptionObject &e) { std::cerr << e; return; } std::cout << "Done." << std::endl; #ifdef DEBUG_ME { WriterType::Pointer writer = WriterType::New(); writer->SetInput( threshold->GetOutput() ); writer->SetFileName( "AfterThreshold.hdr" ); writer->Update(); } #endif // erode to remove background noise typedef itk::BinaryBallStructuringElement StructuralElementType; StructuralElementType ball; typename itk::BinaryErodeImageFilter::Pointer erode = itk::BinaryErodeImageFilter::New(); ball.SetRadius( 3 ); erode->SetInput( threshold->GetOutput() ); erode->SetKernel( ball ); std::cout << "Eroding..." << std::flush; try { erode->Update(); } catch( itk::ExceptionObject &e) { std::cerr << e; return; } std::cout << "Done." << std::endl; #ifdef DEBUG_ME { WriterType::Pointer writer = WriterType::New(); writer->SetInput( erode->GetOutput() ); writer->SetFileName( "AfterErode.hdr" ); writer->Update(); } #endif typedef BinaryCrossStructuringElement CrossType; typedef BinaryDilateImageFilter DilateFilterType; typedef AndImageFilter AndFilterType; typename OutputImageType::Pointer M0 = threshold->GetOutput(); typename OutputImageType::Pointer Mn = OutputImageType::New(); Mn->SetRegions( M0->GetLargestPossibleRegion() ); Mn->SetSpacing( M0->GetSpacing() ); Mn->SetOrigin( M0->GetOrigin() ); Mn->SetDirection( M0->GetDirection() ); Mn->Allocate(); typename OutputImageType::Pointer Mnplus1 = erode->GetOutput(); CrossType cross; cross.SetRadius( 1 ); //unsigned long rad[3]={3,3,3}; //ball2.SetRadius( rad ); std::cout << "Conditional reconstruction..." << std::flush; int iter = 0; do { std::cout << "Iteration: " << iter++ << std::endl; CopyImage( Mn, Mnplus1); typename DilateFilterType::Pointer dilater = DilateFilterType::New(); dilater->SetInput( Mn ); dilater->SetKernel( cross ); try { dilater->Update(); } catch( itk::ExceptionObject &e) { std::cerr << e; return; } typename AndFilterType::Pointer andfilter = AndFilterType::New(); andfilter->SetInput(0, M0); andfilter->SetInput(1, dilater->GetOutput() ); try { andfilter->Update(); } catch( itk::ExceptionObject &e) { std::cerr << e; return; } Mnplus1 = andfilter->GetOutput(); /* #ifdef DEBUG_ME { WriterType::Pointer writer = WriterType::New(); writer->SetInput( andfilter->GetOutput() ); char filename[512]; sprintf( filename, "CondReconstruction_iter_%d.hdr", iter); writer->SetFileName( filename ); writer->Update(); } #endif*/ } while( !CompareImages( Mn, Mnplus1) ); std::cout << "Done." << std::endl; #ifdef DEBUG_ME { WriterType::Pointer writer = WriterType::New(); writer->SetInput( Mn ); writer->SetFileName( "AfterCondReconstruction.hdr" ); writer->Update(); } #endif // now fill the holes typename itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::Pointer filler = itk::VotingBinaryIterativeHoleFillingImageFilter< OutputImageType >::New(); filler->SetInput( Mn ); filler->SetMaximumNumberOfIterations (1000); std::cout << "Filling the holes..." << std::flush; try { filler->Update(); } catch( itk::ExceptionObject &e) { std::cerr << e; return; } std::cout << "Done." << std::endl; typename OutputImageType::Pointer outputImage = - static_cast< OutputImageType * >(this->ProcessObject::GetOutput()); + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); outputImage->SetSpacing( filler->GetOutput()->GetSpacing() ); // Set the image spacing outputImage->SetOrigin( filler->GetOutput()->GetOrigin() ); // Set the image origin outputImage->SetLargestPossibleRegion( filler->GetOutput()->GetLargestPossibleRegion()); outputImage->SetBufferedRegion( filler->GetOutput()->GetLargestPossibleRegion() ); outputImage->SetDirection( filler->GetOutput()->GetDirection() ); outputImage->Allocate(); itk::ImageRegionIterator itIn( filler->GetOutput(), filler->GetOutput()->GetLargestPossibleRegion() ); itk::ImageRegionIterator itOut( outputImage, outputImage->GetLargestPossibleRegion() ); while( !itIn.IsAtEnd() ) { itOut.Set(itIn.Get()); ++itIn; ++itOut; } } template< class TOutputImagePixelType > void BrainMaskExtractionImageFilter< TOutputImagePixelType > ::CopyImage( typename OutputImageType::Pointer target, typename OutputImageType::Pointer source) { itk::ImageRegionConstIterator itIn( source, source->GetLargestPossibleRegion() ); itk::ImageRegionIterator itOut( target, target->GetLargestPossibleRegion() ); while( !itOut.IsAtEnd() ) { itOut.Set( itIn.Get() ); ++itIn; ++itOut; } } template< class TOutputImagePixelType > bool BrainMaskExtractionImageFilter< TOutputImagePixelType > ::CompareImages( typename OutputImageType::Pointer im1, typename OutputImageType::Pointer im2) { itk::ImageRegionConstIterator itIn( im1, im1->GetLargestPossibleRegion() ); itk::ImageRegionConstIterator itOut( im2, im2->GetLargestPossibleRegion() ); while( !itOut.IsAtEnd() ) { if( itOut.Value() != itIn.Value() ) { return false; } ++itOut; ++itIn; } return true; } template< class TOutputImagePixelType > int BrainMaskExtractionImageFilter< TOutputImagePixelType > ::ComputeHistogram( typename InputImageType::Pointer image) { // IMPORTANT: IMAGE MUST BE UNSIGNED SHORT int N=65535; int* histogram = new int[N]; for( int i=0; i itIn( image, image->GetLargestPossibleRegion() ); long totVoxels = 0; int max = -1; int min = 9999999; while( !itIn.IsAtEnd() ) { histogram[ (int)(itIn.Value()) ]++; if( itIn.Value()>max ) { max = itIn.Value(); } if( itIn.Value()EPS ) for( seuil = min; seuil<=max; seuil++) { //seuil = newseuil; // compute the classes: double mean0 = 0.0; double mean1 = 0.0; //double std0 = 0.0; //double std1 = 0.0; double num0 = 0.0; double num1 = 0.0; for( int i=min; i 65535 || newseuil<0) { std::cerr << "Error: threshold is too low or high, exiting" << std::endl; return -1; }*/ double Vm = num0 * (mean0 - mean)*(mean0 - mean) + num1*(mean1 - mean)*(mean1 - mean); if( Vm > V ) { V = Vm; S = seuil; //std::cout << "New seuil: " << S << std::endl; //getchar(); } } delete [] histogram; std::cout << "Seuil: " << S << std::endl; return S; } } #endif // __itkBrainMaskExtractionImageFilter_txx