diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkBrainMaskExtractionImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkBrainMaskExtractionImageFilter.txx index 812ac8ffdf..8a3061d02d 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::GetPrimaryOutput()); 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 diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkGroupDiffusionHeadersFilter.h b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkGroupDiffusionHeadersFilter.h index a4f37682bb..4c629cf3f8 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkGroupDiffusionHeadersFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkGroupDiffusionHeadersFilter.h @@ -1,105 +1,109 @@ /*=================================================================== 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 _mitkGroupDiffusionHeadersFilter_h #define _mitkGroupDiffusionHeadersFilter_h #include "DiffusionCoreExports.h" #include "mitkBaseProcess.h" #include "mitkDiffusionImageHeaderInformation.h" namespace mitk { /** * @brief Groups Headers containing Diffusion Information * @ingroup Process */ class DiffusionCore_EXPORT GroupDiffusionHeadersFilter : public BaseProcess { public: mitkClassMacro( GroupDiffusionHeadersFilter, BaseProcess ); itkNewMacro( Self ); typedef mitk::DiffusionImageHeaderInformation HeaderType; typedef HeaderType::Pointer HeaderPointer; typedef std::vector InputType; typedef std::vector OutputType; /** * Sets the input of this process object * @param input the input */ virtual void SetInput( InputType input ); /** * Sets the input n'th of this process object * @param idx the number associated with the given input */ virtual void SetInput( const unsigned int& idx, InputType input ); /** * Sets the input n'th of this process object * @param idx the number associated with the given input */ virtual void SetNthOutput( const unsigned int& idx, InputType output ); /** * @returns the input tree of the process object */ InputType GetInput( void ); /** * @param idx the index of the input to return * @returns the input object with the given index */ InputType GetInput( const unsigned int& idx ); OutputType GetOutput(); + virtual DataObjectPointer MakeOutput ( DataObjectPointerArraySizeType ){} + + virtual DataObjectPointer MakeOutput(const DataObjectIdentifierType& ) {} + virtual void GenerateOutputInformation(); virtual void Update(); protected: /** * A default constructor */ GroupDiffusionHeadersFilter(); /** * The destructor */ virtual ~GroupDiffusionHeadersFilter(); OutputType m_Output; InputType m_Input; private: void operator=( const Self& ); //purposely not implemented } ; } //end of namespace mitk #endif