MedicalVisualization
MR Time Series Processing with ITK
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;
}
#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
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
#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
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
#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