MedicalVisualization

MR Time Series Processing with ITK

TTP Heart | | Graphic Concepts

Concept paper about 4D post-processing with ITK: PDF

Example 4D processing pipeline to generate MIPt and WI parameter maps:

#include "itkNumericTraits.h"

#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageFileWriter.h"

#include "itkMetaDataDictionary.h"

#include "itkNaryElevateImageFilter.h"
#include "itkUnaryRetractImageFilter.h"

#include "itkStreamingImageFilter.h"

#include "functors.cxx"

int main(int argc, char* argv[])
{
  // check arguments
  if (argc!=2)
  {
    std::cerr << "Usage: executable <dicomDirectory>" << std::endl;
    return EXIT_FAILURE;
  }

  // declare images
  typedef signed short PixelType;
  typedef itk::Image<PixelType, 2> Image2DType;
  typedef itk::Image<PixelType, 3> Image3DType;
  typedef itk::Image<PixelType, 4> Image4DType;

  // create elevation filter
  typedef itk::NaryElevateImageFilter<Image3DType, Image4DType> NaryElevateFilterType;
  NaryElevateFilterType::Pointer elevateFilter=NaryElevateFilterType::New();

  // create name generator and attach to reader
  typedef itk::GDCMSeriesFileNames NamesGeneratorType;
  NamesGeneratorType::Pointer nameGenerator=NamesGeneratorType::New();
  nameGenerator->SetUseSeriesDetails(true);
  nameGenerator->AddSeriesRestriction("0020|0012"); // acquisition number
  nameGenerator->SetDirectory(argv[1]);

  // get series IDs
  typedef std::vector<std::string> SeriesIdContainer;
  const SeriesIdContainer &seriesUID=nameGenerator->GetSeriesUIDs();

  // create reader array
  typedef itk::ImageSeriesReader<Image3DType> ReaderType;
  ReaderType::Pointer *reader=new ReaderType::Pointer[seriesUID.size()];

  // declare series iterators
  SeriesIdContainer::const_iterator seriesItr=seriesUID.begin();
  SeriesIdContainer::const_iterator seriesEnd=seriesUID.end();
  int seriesNum=0;

  // connect each series to elevation filter
  while (seriesItr!=seriesEnd)
  {
    reader[seriesNum]=ReaderType::New();

    typedef itk::GDCMImageIO ImageIOType;
    ImageIOType::Pointer dicomIO=ImageIOType::New();
    reader[seriesNum]->SetImageIO(dicomIO);

    typedef std::vector<std::string> FileNamesContainer;
    FileNamesContainer fileNames;

    fileNames=nameGenerator->GetFileNames(seriesItr->c_str());
    reader[seriesNum]->SetFileNames(fileNames);

    elevateFilter->SetInput(seriesNum,reader[seriesNum]->GetOutput());
    reader[seriesNum]->ReleaseDataFlagOn();

    reader[seriesNum]->Update();
    reader[seriesNum]->GetOutput()->SetMetaDataDictionary(dicomIO->GetMetaDataDictionary());

    seriesNum++;

    ++seriesItr;
  }

  // get statistics from elevate filter
  elevateFilter->Update();
  float dataMin=elevateFilter->GetDataMin();
  float dataMax=elevateFilter->GetDataMax();
  float dataMean=elevateFilter->GetDataMean();
  float noiseThres=elevateFilter->GetNoiseThres();

  // pixel type declarations
  typedef Image4DType::PixelType Pixel4DType;
  typedef Image3DType::PixelType Pixel3DType;

  // create retraction filter #1
  typedef itk::Functor::CalcMin< Pixel4DType, Pixel3DType> MinCalculatorType;
  typedef itk::UnaryRetractImageFilter<Image4DType, Image3DType, MinCalculatorType> UnaryRetractFilterType1;
  UnaryRetractFilterType1::Pointer retractFilter1=UnaryRetractFilterType1::New();

  // create retraction filter #2
  typedef itk::Functor::CalcMinDev< Pixel4DType, Pixel3DType> MinDevCalculatorType;
  typedef itk::UnaryRetractImageFilter<Image4DType, Image3DType, MinDevCalculatorType> UnaryRetractFilterType2;
  UnaryRetractFilterType2::Pointer retractFilter2=UnaryRetractFilterType2::New();

  // connect elevation to retract filters
  retractFilter1->SetInput(elevateFilter->GetOutput());
  retractFilter2->SetInput(elevateFilter->GetOutput());

  // write data to multiframe DICOM
  typedef itk::ImageFileWriter<Image3DType> WriterType;
  WriterType::Pointer writer=WriterType::New();
  writer->SetFileName("out-MIPt.dcm");
  writer->SetInput(retractFilter1->GetOutput());
  writer->Update();
  retractFilter1=NULL;
  writer->SetFileName("out-WI.dcm");
  writer->SetInput(retractFilter2->GetOutput());
  writer->Update();
  retractFilter2=NULL;

  return EXIT_SUCCESS;
}

itk::NaryElevateImageFilter.h:

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkNaryElevateImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2007/02/05 10:22:00 $
  Version:   $Revision: 1.0 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#ifndef __itkNaryElevateImageFilter_h
#define __itkNaryElevateImageFilter_h

#include "itkImageToImageFilter.h"
#include "itkArray.h"

namespace itk
{

namespace Functor
{

// this functor is used if the corresponding template function is left undefined
template <class TInput, class TOutput>
class DefaultElevateFunction
{
public:

  DefaultElevateFunction() {}
  ~DefaultElevateFunction() {}

  inline void operator()(const Array<TInput> &A,const Array<double> &T,Array<TOutput> &R)
  {
    int size=A.size();
    for (int i=0; i<size; i++) R[i]=static_cast<TOutput>(A[i]); // just copy without change
  }
};

}

/** \class NaryElevateImageFilter
 * \brief Implements generic elevation of a time series of several nD inputs to a single n+1D output.
 *
 * Contributed by Stefan Roettger, Siemens Medical MR, Mar. 2007
 *
 * This class takes a time series of inputs with the same dimension and type and
 * generates an output image of the next higher dimension which contains the entire time series.
 * For example, a set of 3D inputs will be composed into a 4D output with
 * the highest dimension corresponding to the time axis.
 * Additionally, a functor style template parameter can be used for
 * per-pixel-wise manipulation of the time curves of the input data.
 * The functor can be ommited, then the data is just copied.
 * The filter also provides the min, the max, the mean and
 * an approximation of the noise level of the generated output.
 *
 * A typical times series produced in clinical practice is a MR breast angiography.
 * Gadolinium contrast agent is infused and T1-weighted 3D MR-images are taken every minute.
 * The procedure usually generates 5-10 time points which show the perfusion of the breast tissue.
 * The itkNaryElevateImageFilter takes the 3D breast images and constructs a single 4D image
 * containing the information of the entire scan.
 * After this, the itkUnaryRetractFilter can be used to calculate certain characteristics of
 * the time curves which give insight into the malignancy of the breast tissue.
 * Characteristic operations are MIPt (Maximum Intensity [Projection] over time) and
 * WI (Wash-In or temporal derivative), for example.
 * A tumor usually shows higher WI values than normal tissue, so that this parameter
 * (amongst a variety of others) can be used to diagnose breast cancer.
 *
 * \ingroup IntensityImageFilters Multithreaded
 */


template <class TInputImage, class TOutputImage, class TFunction =
  typename Functor::DefaultElevateFunction<
    typename TInputImage::PixelType,
    typename TOutputImage::PixelType> >
 class ITK_EXPORT NaryElevateImageFilter:
  public ImageToImageFilter<TInputImage, TOutputImage>
{
public:

  /** Standard class typedefs. */
  typedef NaryElevateImageFilter Self;
  typedef ImageToImageFilter<TInputImage, TOutputImage> Superclass;
  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self> ConstPointer;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(NaryElevateImageFilter, ImageToImageFilter);

  /** Some typedefs. */
  typedef TFunction FunctorType;
  typedef TInputImage InputImageType;
  typedef typename InputImageType::Pointer     InputImagePointer;
  typedef typename InputImageType::RegionType  InputImageRegionType;
  typedef typename InputImageType::PixelType   InputImagePixelType;
  typedef typename InputImageType::SizeType    InputImageSizeType;
  typedef typename InputImageType::IndexType   InputImageIndexType;
  typedef typename InputImageType::PointType   InputImageOriginType;
  typedef typename InputImageType::SpacingType InputImageSpacingType;
  typedef TOutputImage OutputImageType;
  typedef typename OutputImageType::Pointer     OutputImagePointer;
  typedef typename OutputImageType::RegionType  OutputImageRegionType;
  typedef typename OutputImageType::PixelType   OutputImagePixelType;
  typedef typename OutputImageType::SizeType    OutputImageSizeType;
  typedef typename OutputImageType::IndexType   OutputImageIndexType;
  typedef typename OutputImageType::PointType   OutputImageOriginType;
  typedef typename OutputImageType::SpacingType OutputImageSpacingType;
  typedef Array<InputImagePixelType> NaryInputType;
  typedef Array<OutputImagePixelType> NaryOutputType;

  /** ImageDimension constants. */
  itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
  itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);

  /** Get functions. */
  itkGetMacro(DataMin,float);
  itkGetMacro(DataMax,float);
  itkGetMacro(DataMean,float);
  itkGetMacro(NoiseThres,float);

  // get the actual functor
  FunctorType& GetFunctor() {return m_Functor;}

  // set the actual functor
  void SetFunctor(FunctorType& functor)
  {
    m_Functor = functor;
    this->Modified();
  }

protected:

  NaryElevateImageFilter();
  virtual ~NaryElevateImageFilter() {};

  void GenerateOutputInformation();
  void GenerateData();

private:

  NaryElevateImageFilter(const Self&); // intentionally not implemented
  void operator=(const Self&);         // intentionally not implemented

  float m_DataMin;
  float m_DataMax;
  float m_DataMean;
  float m_NoiseThres;

  FunctorType m_Functor;
};

}

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkNaryElevateImageFilter.txx"
#endif

#endif

itk::NaryElevateImageFilter.txx:

#ifndef _itkNaryElevateImageFilter_txx
#define _itkNaryElevateImageFilter_txx

#include "stdio.h"

#include "itkNaryElevateImageFilter.h"

#include "itkMetaDataDictionary.h"

#include "itkImageRegionIteratorWithIndex.h"

namespace itk
{

// helper function to get the value of a specific tag
double gettagvalue(std::string &entryId, MetaDataDictionary &metaDict)
{
  ExceptionObject ex("invalid tag value");

  MetaDataDictionary::ConstIterator tagItr=metaDict.Find(entryId);

  if (tagItr!=metaDict.End())
  {
    typedef MetaDataObject<std::string> MetaDataStringType;
    MetaDataObjectBase::Pointer entry=tagItr->second;
    MetaDataStringType::Pointer entryvalue=
      dynamic_cast<MetaDataStringType *>(entry.GetPointer());

    if (entryvalue)
    {
      std::string tagvalue=entryvalue->GetMetaDataObjectValue();

      float value;
      if (sscanf(tagvalue.c_str(),"%g",&value)!=1) throw ex;
      else return(value);
    }
  }
  else throw ex;

  return(0.0);
}

// constructor
template <class TInputImage, class TOutputImage, class TFunction>
NaryElevateImageFilter<TInputImage, TOutputImage, TFunction>::
NaryElevateImageFilter()
{
  m_DataMin=0.0f;
  m_DataMax=0.0f;
  m_DataMean=0.0f;
  m_NoiseThres=0.0f;
}

// check inputs
template <class TInputImage, class TOutputImage, class TFunction>
void
NaryElevateImageFilter<TInputImage, TOutputImage, TFunction>::
GenerateOutputInformation()
{
  const unsigned int dimI=InputImageDimension;
  const unsigned int dimO=OutputImageDimension;

  if (dimO<dimI+1) itkExceptionMacro(<< "inadequate output image dimension");

  OutputImagePointer outputPtr = this->GetOutput(0);

  const unsigned int numberOfInputImages =
    static_cast<int>(this->GetNumberOfInputs());

  for (int i=0; i<numberOfInputImages; i++)
  {
    InputImagePointer inputPtr =
      dynamic_cast<InputImageType *>(ProcessObject::GetInput(i));

    InputImageRegionType regionI=inputPtr->GetRequestedRegion();

    InputImageSizeType    sizeI    = regionI.GetSize();
    InputImageIndexType   indexI   = regionI.GetIndex();
    InputImageOriginType  originI  = inputPtr->GetOrigin();
    InputImageSpacingType spacingI = inputPtr->GetSpacing();

    OutputImageSizeType    sizeO;
    OutputImageIndexType   indexO;
    OutputImageOriginType  originO;
    OutputImageSpacingType spacingO;

    if (i==0)
    {
      for (int j=0; j<dimI; j++)
      {
        sizeO[j]    = sizeI[j];
        indexO[j]   = indexI[j];
        originO[j]  = originI[j];
        spacingO[j] = spacingI[j];
      }

      sizeO[dimI]=numberOfInputImages;
      indexO[dimI]=0;
      originO[dimI]=0;
      spacingO[dimI]=1;

      for (int j=dimI+1; j<dimO; j++)
      {
        sizeO[j]    = 1;
        indexO[j]   = 0;
        originO[j]  = 0;
        spacingO[j] = 1;
      }

      OutputImageRegionType regionO;
      regionO.SetSize(sizeO);
      regionO.SetIndex(indexO);

      outputPtr->SetOrigin(originO);
      outputPtr->SetSpacing(spacingO);

      outputPtr->SetRegions(regionO);
    }
    else
      for (int j=0; j<dimI; j++)
        if (sizeO[j]!=sizeI[j] || indexO[j]!=indexI[j])
          itkExceptionMacro(<< "mismatch of nary input images");
  }

  // pass all tags from first time point:

  InputImagePointer inputPtr0 =
    dynamic_cast<InputImageType *>(ProcessObject::GetInput(0));

  outputPtr->SetMetaDataDictionary(inputPtr0->GetMetaDataDictionary());

  // retrieve all time points as list:

  typedef MetaDataDictionary DictionaryType;

  std::list<double> list;

  double firstdate;
  double firsttime;

  int firsthour;
  int firstminutes;
  double firstseconds;

  bool dontusedict=false;

  for (int i=0; i<numberOfInputImages; i++)
  {
    InputImagePointer inputPtr =
      dynamic_cast<InputImageType *>(ProcessObject::GetInput(i));

    DictionaryType &metaDict=inputPtr->GetMetaDataDictionary();

    if (metaDict.Begin()==metaDict.End() || dontusedict)
    {
      ExceptionObject ex("invalid time series");
      if (i>0 && !dontusedict) throw ex;

      list.push_back(i); // assume that time points are given in minutes
      dontusedict=true;
      continue;
    }

    std::string entryIdAD="0008|0022"; // acquisition date
    std::string entryIdAT="0008|0032"; // acquisition time
    std::string entryIdAN="0020|0012"; // aquisition number

    double entryAD=gettagvalue(entryIdAD,metaDict);
    double entryAT=gettagvalue(entryIdAT,metaDict);

    if (i==0)
    {
      firstdate=entryAD;
      firsttime=entryAT;

      firsthour=(int)(firsttime/10000)%100;
      firstminutes=(int)(firsttime/100)%100;
      firstseconds=firsttime-100*((int)firsttime/100);

      list.push_back(0.0);
    }
    else
    {
      double date=entryAD;
      double time=entryAT;

      int hour=(int)(time/10000)%100;
      int minutes=(int)(time/100)%100;
      double seconds=time-100*((int)time/100);

      // assume that if date differs we have just passed midnight
      if (date!=firstdate)
      {
        hour+=24;
        firstdate=date;
      }

      double difference=(seconds-firstseconds)+
                        60.0*(minutes-firstminutes)+
                        3600.0*(hour-firsthour);

      list.push_back(difference/60.0); // minutes
    }
  }

  // encapsulate list into MetaDataObject and append to MetaDataDictionary:

  typedef MetaDataObject< std::list<double> > MetaDataListType;
  MetaDataListType::Pointer value=MetaDataListType::New();
  value->SetMetaDataObjectValue(list);

  std::string internalId="___internal-4D-filter-time-point-list"; // internal tag
  DictionaryType &metaDict=outputPtr->GetMetaDataDictionary();
  metaDict[internalId]=value;
}

// elevate inputs
template <class TInputImage, class TOutputImage, class TFunction>
void
NaryElevateImageFilter<TInputImage, TOutputImage, TFunction>::
GenerateData()
{
  const unsigned int dimI=InputImageDimension;
  const unsigned int dimO=OutputImageDimension;

  OutputImagePointer outputPtr = this->GetOutput(0);

  const unsigned int numberOfInputImages =
    static_cast<int>(this->GetNumberOfInputs());

  typedef ImageRegionConstIteratorWithIndex<InputImageType> ImageRegionConstIteratorWithIndexType;
  std::vector<ImageRegionConstIteratorWithIndexType *> inputItrVector(numberOfInputImages);

  for (int i=0; i<numberOfInputImages; i++)
  {
    InputImagePointer inputPtr =
      dynamic_cast<InputImageType *>(ProcessObject::GetInput(i));

    inputPtr->Update();

    ImageRegionConstIteratorWithIndexType *inputItr =
      new ImageRegionConstIteratorWithIndexType(inputPtr, inputPtr->GetBufferedRegion());

    inputItrVector[i] = reinterpret_cast<ImageRegionConstIteratorWithIndexType *>(inputItr);
    inputItrVector[i]->GoToBegin();
  }

  // get time point list from meta dictionary:

  typedef MetaDataDictionary DictionaryType;
  DictionaryType &metaDict=outputPtr->GetMetaDataDictionary();

  std::string internalId="___internal-4D-filter-time-point-list"; // internal tag
  MetaDataDictionary::ConstIterator tagItr=metaDict.Find(internalId);

  ExceptionObject ex("invalid time point list");

  Array<double> timePoints;
  timePoints.SetSize(numberOfInputImages);

  if (tagItr!=metaDict.End())
  {
    typedef MetaDataObject< std::list<double> > MetaDataListType;

    MetaDataObjectBase::Pointer entry = tagItr->second;
    MetaDataListType::Pointer entryvalue = dynamic_cast<MetaDataListType *>( entry.GetPointer() );

    // check whether or not the type of the entry value is correct
    if (entryvalue)
    {
      std::list<double> list;
      list=entryvalue->GetMetaDataObjectValue();

      if (list.size()!=numberOfInputImages) throw ex;

      std::list<double>::const_iterator listItr=list.begin();

      for (int k=0; k<numberOfInputImages; k++)
      {
        timePoints[k]=*listItr;
        listItr++;
      }
    }
    else throw ex;
  }
  else throw ex;

  // process inputs by retrieving one time series after another:

  outputPtr->Allocate();

  m_DataMin=NumericTraits<float>::max();
  m_DataMax=NumericTraits<float>::min();

  m_DataMean=0.0f;
  unsigned int meancount=0;

  m_NoiseThres=0.0f;
  unsigned int noisecount=0;

  InputImageIndexType indexI;
  OutputImageIndexType indexO;

  for (int i=dimI+1; i<dimO; i++) indexO[i]=0;

  while(!inputItrVector[0]->IsAtEnd())
  {
    indexI=inputItrVector[0]->GetIndex();

    for (int j=0; j<dimI; j++) indexO[j]=indexI[j];

    InputImagePixelType base;

    float mean_fvalue=0.0f;

    float mingrad=NumericTraits<float>::max();
    float rndgrad=0.0f;

    float last_fvalue=0.0f,lastlast_fvalue=0.0f;

    NaryInputType values;
    values.set_size(numberOfInputImages);

    // get actual time series
    for (int inputNumber=0; inputNumber<numberOfInputImages; inputNumber++)
    {
      InputImagePixelType value;
      value=inputItrVector[inputNumber]->Get();

      values[inputNumber]=value;

      ++(*inputItrVector[inputNumber]);
    }

    // push actual time series through elevate functor
    NaryOutputType result;
    result.SetSize(numberOfInputImages);
    m_Functor(values,timePoints,result);

    // check actual time series
    for (int inputNumber=0; inputNumber<numberOfInputImages; inputNumber++)
    {
      OutputImagePixelType value;
      value=result[inputNumber];

      if (inputNumber==0) base=value;

      float fvalue;
      fvalue=static_cast<float>(value);

      if (fvalue<m_DataMin) m_DataMin=fvalue;
      if (fvalue>m_DataMax) m_DataMax=fvalue;

      mean_fvalue+=fvalue;

      if (inputNumber>1)
      {
        float grad=last_fvalue-lastlast_fvalue;
        if (grad<0.0f) grad=-grad;

        if (grad<mingrad)
        {
          mingrad=grad;
          rndgrad=fvalue-last_fvalue;
          if (rndgrad<0.0f) rndgrad=-rndgrad;
        }
      }

      lastlast_fvalue=last_fvalue;
      last_fvalue=fvalue;
    }

    // store resulting time series
    for (int inputNumber=0; inputNumber<numberOfInputImages; inputNumber++)
    {
      indexO[dimI]=inputNumber;
      outputPtr->SetPixel(indexO, result[inputNumber]);
    }

    mean_fvalue/=numberOfInputImages;

    if (++meancount>1) m_DataMean*=float(meancount-1)/meancount;
    m_DataMean+=mean_fvalue/meancount;

    m_NoiseThres+=rndgrad;
    noisecount++;
  }

  // clean up:

  for (int i=0; i<numberOfInputImages; i++) delete inputItrVector[i];

  m_NoiseThres/=noisecount;
}

}

#endif

itk::UnaryRetractImageFilter.h:

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkUnaryRetractImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2007/02/05 11:15:00 $
  Version:   $Revision: 1.0 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#ifndef __itkUnaryRetractImageFilter_h
#define __itkUnaryRetractImageFilter_h

#include "itkImageToImageFilter.h"

namespace itk
{

/** \class UnaryRetractImageFilter
 * \brief Implements generic retraction of a nD time series producing a n-1D output.
 *
 * Contributed by Stefan Roettger, Siemens Medical MR, Mar. 2007
 *
 * This class interpretes the input image as a time series and
 * generates an output image of the next lower dimension.
 * The retraction operation is specified by a functor style template parameter.
 * It receives the time curves of each pixel and produces a single scalar value
 * effectively removing the time dimension of the input image.
 * For example, a 4D times series (generated by the itkNaryElevateImageFilter)
 * will be transformed into a 3D output image by applying the functor to each pixel.
 *
 * A typical times series produced in clinical practice is a MR breast angiography.
 * Gadolinium contrast agent is infused and T1-weighted 3D MR-images are taken every minute.
 * The procedure usually generates 5-10 time points which show the perfusion of the breast tissue.
 * The itkNaryElevateImageFilter can be used to a single 4D image from the original 3D breast images.
 * After this, the itkUnaryRetractFilter can be used to calculate certain characteristics of
 * the time curves which give insight into the malignancy of the breast tissue.
 * Characteristic functor operations are MIPt (Maximum Intensity [Projection] over time) and
 * WI (Wash-In or temporal derivative), for example.
 * A tumor usually shows higher WI values than normal tissue, so that this parameter
 * (amongst a variety of others) can be used to diagnose breast cancer.
 *
 * \ingroup IntensityImageFilters Multithreaded
 */


template <class TInputImage, class TOutputImage, class TFunction>
class ITK_EXPORT UnaryRetractImageFilter:
  public ImageToImageFilter<TInputImage, TOutputImage>
{
public:

  /** Standard class typedefs. */
  typedef UnaryRetractImageFilter Self;
  typedef ImageToImageFilter<TInputImage, TOutputImage> Superclass;
  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self> ConstPointer;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(UnaryRetractImageFilter, ImageToImageFilter);

  /** Some typedefs. */
  typedef TFunction FunctorType;
  typedef TInputImage InputImageType;
  typedef typename InputImageType::Pointer     InputImagePointer;
  typedef typename InputImageType::RegionType  InputImageRegionType;
  typedef typename InputImageType::PixelType   InputImagePixelType;
  typedef typename InputImageType::SizeType    InputImageSizeType;
  typedef typename InputImageType::IndexType   InputImageIndexType;
  typedef typename InputImageType::PointType   InputImageOriginType;
  typedef typename InputImageType::SpacingType InputImageSpacingType;
  typedef TOutputImage OutputImageType;
  typedef typename OutputImageType::Pointer     OutputImagePointer;
  typedef typename OutputImageType::RegionType  OutputImageRegionType;
  typedef typename OutputImageType::PixelType   OutputImagePixelType;
  typedef typename OutputImageType::SizeType    OutputImageSizeType;
  typedef typename OutputImageType::IndexType   OutputImageIndexType;
  typedef typename OutputImageType::PointType   OutputImageOriginType;
  typedef typename OutputImageType::SpacingType OutputImageSpacingType;
  typedef Array<InputImagePixelType> NaryArrayType;

  /** ImageDimension constants. */
  itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
  itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);

  // get the actual functor
  FunctorType& GetFunctor() {return m_Functor;}

  // set the actual functor
  void SetFunctor(FunctorType& functor)
  {
    m_Functor = functor;
    this->Modified();
  }

protected:

  UnaryRetractImageFilter();
  virtual ~UnaryRetractImageFilter() {};

  void GenerateOutputInformation();
  void GenerateData();

private:

  UnaryRetractImageFilter(const Self&); // intentionally not implemented
  void operator=(const Self&);          // intentionally not implemented

  FunctorType m_Functor;
};

}

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkUnaryRetractImageFilter.txx"
#endif

#endif

itk::UnaryRetractImageFilter.txx:

#ifndef _itkUnaryRetractImageFilter_txx
#define _itkUnaryRetractImageFilter_txx

#include "itkUnaryRetractImageFilter.h"

#include "itkImageLinearConstIteratorWithIndex.h"

namespace itk
{

template <class TInputImage, class TOutputImage, class TFunction>
UnaryRetractImageFilter<TInputImage, TOutputImage, TFunction>::
UnaryRetractImageFilter() {}

// check input
template <class TInputImage, class TOutputImage, class TFunction>
void
UnaryRetractImageFilter<TInputImage, TOutputImage, TFunction>::
GenerateOutputInformation()
{
  const unsigned int dimI=InputImageDimension;
  const unsigned int dimO=OutputImageDimension;

  if (dimI<1 || dimO<dimI-1) itkExceptionMacro(<< "inadequate output image dimension");

  OutputImagePointer outputPtr = this->GetOutput(0);

  InputImagePointer inputPtr =
    dynamic_cast<InputImageType *>(ProcessObject::GetInput(0));

  InputImageRegionType regionI=inputPtr->GetRequestedRegion();

  InputImageSizeType    sizeI    = regionI.GetSize();
  InputImageIndexType   indexI   = regionI.GetIndex();
  InputImageOriginType  originI  = inputPtr->GetOrigin();
  InputImageSpacingType spacingI = inputPtr->GetSpacing();

  OutputImageSizeType    sizeO;
  OutputImageIndexType   indexO;
  OutputImageOriginType  originO;
  OutputImageSpacingType spacingO;

  for (int i=0; i<dimI-1; i++)
  {
    sizeO[i]    = sizeI[i];
    indexO[i]   = indexI[i];
    originO[i]  = originI[i];
    spacingO[i] = spacingI[i];
  }

  for (int i=dimI; i<dimO; i++)
  {
    sizeO[i]    = 1;
    indexO[i]   = 0;
    originO[i]  = 0;
    spacingO[i] = 1;
  }

  OutputImageRegionType regionO;
  regionO.SetSize(sizeO);
  regionO.SetIndex(indexO);

  outputPtr->SetOrigin(originO);
  outputPtr->SetSpacing(spacingO);

  outputPtr->SetRegions(regionO);

  // pass tags:

  InputImagePointer inputPtr0 =
    dynamic_cast<InputImageType *>(ProcessObject::GetInput(0));

  outputPtr->SetMetaDataDictionary(inputPtr0->GetMetaDataDictionary());
}

// retract input
template <class TInputImage, class TOutputImage, class TFunction>
void
UnaryRetractImageFilter<TInputImage, TOutputImage, TFunction>::
GenerateData()
{
  const unsigned int dimI=InputImageDimension;
  const unsigned int dimO=OutputImageDimension;

  OutputImagePointer outputPtr = this->GetOutput(0);

  InputImagePointer inputPtr =
    dynamic_cast<InputImageType *>(ProcessObject::GetInput(0));

  inputPtr->Update();

  outputPtr->Allocate();

  const unsigned int runLength=inputPtr->GetBufferedRegion().GetSize()[dimI-1];

  // get time point list from meta dictionary:

  typedef MetaDataDictionary DictionaryType;
  DictionaryType &metaDict=inputPtr->GetMetaDataDictionary();

  std::string internalId="___internal-4D-filter-time-point-list"; // internal tag
  MetaDataDictionary::ConstIterator tagItr=metaDict.Find(internalId);

  ExceptionObject ex("invalid time point list");

  Array<double> timePoints;
  timePoints.SetSize(runLength);

  if (tagItr!=metaDict.End())
  {
    typedef MetaDataObject< std::list<double> > MetaDataListType;

    MetaDataObjectBase::Pointer entry = tagItr->second;
    MetaDataListType::Pointer entryvalue = dynamic_cast<MetaDataListType *>( entry.GetPointer() );

    // check whether or not the type of the entry value is correct
    if (entryvalue)
    {
      std::list<double> list;
      list=entryvalue->GetMetaDataObjectValue();

      if (list.size()!=runLength) throw ex;

      std::list<double>::const_iterator listItr=list.begin();

      for (int k=0; k<runLength; k++)
      {
        timePoints[k]=*listItr;
        listItr++;
      }
    }
    else throw ex;
  }
  else throw ex;

  typedef ImageLinearConstIteratorWithIndex<InputImageType> IteratorType;

  // process input by retrieving one time series after another:

  IteratorType itr(inputPtr, inputPtr->GetBufferedRegion());
  itr.SetDirection(dimI-1); // walk along highest dimension
  itr.GoToBegin();

  InputImageIndexType indexI;
  OutputImageIndexType indexO;

  for (int i=dimI-1; i<dimO; i++) indexO[i]=0;

  NaryArrayType naryInputArray;
  naryInputArray.SetSize(runLength);

  while(!itr.IsAtEnd())
  {
    itr.GoToBeginOfLine();

    unsigned int actualRun=0;

    // get actual time series
    while(!itr.IsAtEndOfLine())
    {
      naryInputArray[actualRun++] = itr.Get();
      ++itr;
    }

    indexI=itr.GetIndex();

    for (int i=0; i<dimI-1; i++) indexO[i]=indexI[i];

    // push actual time series through retract functor and store result
    outputPtr->SetPixel(indexO, static_cast<OutputImagePixelType>(m_Functor(naryInputArray,timePoints)));

    itr.NextLine();
  }
}

}

#endif


TTP Heart | | Graphic Concepts

Options: