mitkConvolutionHelper.h

Go to the documentation of this file.

/*============================================================================
 
The Medical Imaging Interaction Toolkit (MITK)
 
Copyright (c) German Cancer Research Center (DKFZ)
All rights reserved.
 
Use of this source code is governed by a 3-clause BSD license that can be
found in the LICENSE file.
 
============================================================================*/
#ifndef mitkConvolutionHelper_h
#define mitkConvolutionHelper_h
 
#include "itkArray.h"
#include "mitkAIFBasedModelBase.h"
#include <iostream>
#include "MitkPharmacokineticsExports.h"
 
namespace  mitk {
  namespace convolution {
 
inline    itk::Array<double> wrap1d(itk::Array<double> kernel)
        {
            int dim = kernel.GetNumberOfElements();
            itk::Array<double> wrappedKernel(dim);
            wrappedKernel.fill(0.);
                    for(int i=0; i< dim; ++i)
                    {
                        wrappedKernel.SetElement(i, kernel.GetElement((i+(dim/2))%dim));
                    }
 
            return wrappedKernel;
        }
 
   inline  itk::Array<double> zeropadding1d(itk::Array<double> unpaddedSpectrum, int paddedDimension)
        {
 
            int initialDimension = unpaddedSpectrum.GetNumberOfElements();
 
            itk::Array<double> paddedSpectrum(paddedDimension);
            paddedSpectrum.fill(0.);
 
            if(paddedDimension > initialDimension)
            {
                unsigned int padding = paddedDimension - initialDimension;
 
                for(int i=0; i<initialDimension ;++i)
                {
                    paddedSpectrum.SetElement(i+padding/2, unpaddedSpectrum.GetElement(i));
                }
            }
            return paddedSpectrum;
        }
 
  inline  itk::Array<double> unpadAndScale(itk::Array<double> convolutionResult, int initialDimension)
        {
            int transformationDimension = convolutionResult.size();
            unsigned int padding = transformationDimension - initialDimension;
 
            itk::Array<double>  scaledResult(initialDimension);
            scaledResult.fill(0.0);
 
            for(int i = 0; i<initialDimension; ++i)
            {
                double value = convolutionResult(i+padding/2) / transformationDimension;
                scaledResult.SetElement(i,value);
            }
            return scaledResult;
        }
  inline  void prepareConvolution(const itk::Array<double>& kernel, const itk::Array<double>& spectrum, itk::Array<double>& preparedKernel, itk::Array<double>& preparedSpectrum ){
        int convolutionDimensions = kernel.GetSize() + spectrum.GetSize();
 
//        itk::Array<double> paddedKernel = zeropadding1d(kernel,convolutionDimensions);
        preparedKernel=zeropadding1d(kernel,convolutionDimensions);
 
        preparedSpectrum = zeropadding1d(spectrum,convolutionDimensions);
//        preparedKernel = wrap1d(paddedKernel);
    }
 
    }
 
  inline itk::Array<double> convoluteAIFWithExponential(mitk::ModelBase::TimeGridType timeGrid, mitk::AIFBasedModelBase::AterialInputFunctionType aif, double lambda)
  {
      typedef itk::Array<double> ConvolutionResultType;
      ConvolutionResultType convolution(timeGrid.GetSize());
      convolution.fill(0.0);
 
      convolution(0) = 0;
      for(unsigned int i = 0; i< (timeGrid.GetSize()-1); ++i)
      {
          double dt = timeGrid(i+1) - timeGrid(i);
          double m = (aif(i+1) - aif(i))/dt;
          double edt = exp(-lambda *dt);
 
          convolution(i+1) =edt * convolution(i)
                           + (aif(i) - m*timeGrid(i))/lambda * (1 - edt )
                           + m/(lambda * lambda) * ((lambda * timeGrid(i+1) - 1) - edt*(lambda*timeGrid(i) -1));
 
      }
      return convolution;
  }
 
 
  inline itk::Array<double> convoluteAIFWithConstant(mitk::ModelBase::TimeGridType timeGrid, mitk::AIFBasedModelBase::AterialInputFunctionType aif, double constant)
  {
      typedef itk::Array<double> ConvolutionResultType;
      ConvolutionResultType convolution(timeGrid.GetSize());
      convolution.fill(0.0);
 
      convolution(0) = 0;
      for(unsigned int i = 0; i< (timeGrid.GetSize()-1); ++i)
      {
          double dt = timeGrid(i+1) - timeGrid(i);
          double m = (aif(i+1) - aif(i))/dt;
 
          convolution(i+1) = convolution(i) + constant * (aif(i)*dt + m*timeGrid(i)*dt + m/2*(timeGrid(i+1)*timeGrid(i+1) - timeGrid(i)*timeGrid(i)));
 
      }
      return convolution;
  }
 
}
 
#endif