D:/Programmierung/C/TRTK/include/TRTK/Diffusion.hpp

Go to the documentation of this file.

/*
    Copyright (C) 2010 - 2014 Christoph Haenisch

    Chair of Medical Engineering (mediTEC)
    RWTH Aachen University
    Pauwelsstr. 20
    52074 Aachen
    Germany

    See license.txt for more information.

    Version 0.1.1 (2011-11-29)
*/

#ifndef DIFFUSION_HPP_6431082743
#define DIFFUSION_HPP_6431082743


#include <assert.h>
#include <cmath>
#include <iostream>
#include <vector>

#include <Eigen/Core>

#include "ErrorObj.hpp"
#include "Tools.hpp"


namespace TRTK
{


namespace Diffusion
{


template<class T>
void print(const std::vector<T> & data)
{
    for (unsigned i = 0; i < data.size(); ++i)
    {
        std::cout.precision(4);
        std::cout.width(8);
        std::cout << std::fixed << data[i];
    }
    std::cout << std::endl;
}


enum BorderInterpolation
{
    CONTINUED,      
    CIRCULAR,       
    INTERPOLATED,   
    REFLECTED,      
    ZERO            
};


enum Error
{
    INVALID_ARGUMENT,   
    UNKNOWN_ERROR       
};

template <class T>
const std::vector<T> firstDerivative(const std::vector<T> & signal, const BorderInterpolation interpolation = CONTINUED)
{
    const unsigned size = signal.size();
    const unsigned n = signal.size() - 1;

    std::vector<T> derivative(size);

    if (size == 0)
    {
        return derivative;
    }
    else if (size == 1)
    {
        derivative.push_back(T(0));
        return derivative;
    }
    else
    {
        // compute boundaries

        // left boundary
        // derivative[0] = 0.5 * (signal[1] - signal[-1])

        switch(interpolation)
        {
            case CONTINUED:
                // signal[-1] := signal[0]
                derivative[0] = 0.5 * (signal[1] - signal[0]);
                break;

            case CIRCULAR:
                // signal[-1] := signal[n]
                derivative[0] = 0.5 * (signal[1] - signal[n]);
                break;

            case INTERPOLATED:
                // signal[-1] := signal[0] + (signal[0] - signal[1]) = 2 * signal[0] - signal[1]
                derivative[0] = signal[1] - signal[0];
                break;

            case REFLECTED:
                // signal[-1] := signal[0]
                derivative[0] = 0.5 * (signal[1] - signal[0]);
                break;

            case ZERO:
                // signal[-1] := 0
                derivative[0] = signal[1] / 2.0;
                break;

            default:
                ErrorObj error("Unknown interpolation method.");
                error.setFunctionName("Differentiate");
                error.setErrorCode(INVALID_ARGUMENT);
                throw error;
        }

        // right boundary
        // derivative[n] = 0.5 * (signal[n+1] - signal[n-1])

        switch(interpolation)
        {
            case CONTINUED:
                // signal[n+1] := signal[n]
                derivative[n] = 0.5 * (signal[n] - signal[n-1]);
                break;

            case CIRCULAR:
                // signal[n+1] := signal[0]
                derivative[n] = 0.5 * (signal[0] - signal[n-1]);
                break;

            case INTERPOLATED:
                // signal[n+1] := signal[n] + (signal[n] - signal[n-1]) = 2 * signal[n] - signal[n-1]
                derivative[n] = signal[n] - signal[n-1];
                break;

            case REFLECTED:
                // signal[n+1] := signal[n]
                derivative[n] = 0.5 * (signal[n] - signal[n-1]);
                break;

            case ZERO:
                // signal[n+1] := 0
                derivative[n] = -signal[n-1] / 2.0;
                break;

            default:
                ErrorObj error("Unknown interpolation method.");
                error.setFunctionName("Differentiate");
                error.setErrorCode(INVALID_ARGUMENT);
                throw error;
        }

        // compute inner part

        for (unsigned i = 1; i < n; ++i)
        {
            derivative[i] = 0.5 * (signal[i+1] - signal[i-1]);
        }

        return derivative;
    }
}


template <class T>
const std::vector<T> secondDerivative(const std::vector<T> & signal, const BorderInterpolation interpolation = CONTINUED)
{
    const unsigned size = signal.size();
    const unsigned n = signal.size() - 1;

    std::vector<T> derivative(size);

    if (size == 0)
    {
        return derivative;
    }
    else if (size == 1)
    {
        derivative.push_back(T(0));
        return derivative;
    }
    else
    {
        // compute boundaries

        // left boundary
        // derivative[0] = signal[-1] - 2 * signal[0] + signal[1]

        switch(interpolation)
        {
            case CONTINUED:
                // signal[-1] := signal[0]
                derivative[0] = -signal[0] + signal[1];
                break;

            case CIRCULAR:
                // signal[-1] := signal[n]
                derivative[0] = signal[n] - 2 * signal[0] + signal[1];
                break;

            case INTERPOLATED:
                // signal[-1] := signal[0] + (signal[0] - signal[1]) = 2 * signal[0] - signal[1]
                derivative[0] =  0;
                break;

            case REFLECTED:
                // signal[-1] := signal[0]
                derivative[0] = -signal[0] + signal[1];
                break;

            case ZERO:
                // signal[-1] := 0
                derivative[0] = -2 * signal[0] + signal[1];
                break;

            default:
                ErrorObj error("Unknown interpolation method.");
                error.setFunctionName("Differentiate");
                error.setErrorCode(INVALID_ARGUMENT);
                throw error;
        }

        // right boundary
        // derivative[n] = signal[n-1] - 2 * signal[n] + signal[n+1]

        switch(interpolation)
        {
            case CONTINUED:
                // signal[n+1] := signal[n]
                derivative[n] = signal[n-1] -signal[n];
                break;

            case CIRCULAR:
                // signal[n+1] := signal[0]
                derivative[n] = signal[n-1] - 2 * signal[n] + signal[0];
                break;

            case INTERPOLATED:
                // signal[n+1] := signal[n] + (signal[n] - signal[n-1]) = 2 * signal[n] - signal[n-1]
                derivative[n] = 0;
                break;

            case REFLECTED:
                // signal[n+1] := signal[n]
                derivative[n] = signal[n-1] -signal[n];
                break;

            case ZERO:
                // signal[n+1] := 0
                derivative[n] = signal[n-1] - 2 * signal[n];
                break;

            default:
                ErrorObj error("Unknown interpolation method.");
                error.setFunctionName("Differentiate");
                error.setErrorCode(INVALID_ARGUMENT);
                throw error;
        }

        // compute inner part

        for (unsigned i = 1; i < n; ++i)
        {
            derivative[i] = signal[i-1] - 2 * signal[i] + signal[i+1];
        }

        return derivative;
    }
}


template <class T>
const std::vector<T> linearDiffusion(const std::vector<T> & signal,
                                     T time,
                                     T step_size = 0.25,
                                     BorderInterpolation interpolation_method = CONTINUED)
{
    /* The general formulation of a nonlinear diffusion is
     *
     *      d/dt f(x, t) = div grad f(x, t) = laplace f(x, t)
     *
     * In the 1-dimensional case, the divergence and the gradient simplify to
     * an ordinary differentiation. Thus we get
     *
     *      d/dt f(x, t) = d^2/dx^2 f(x, t)
     *
     * The diffusion equation can be easily solved by using the following
     * approximation
     *
     *      d/dt f(x, t) = [f(x, t+h) - f(x, t) ] / h
     *
     * where h is the step size. With this, we get
     *
     *      d/dt f(x, t+h) = f(x, t) + h *  d^2/dx^2 f(x, t)
     *
     * Iterating until t + n * h = time, yields the solution.
     */


    // If time == 0 return the unchanged signal.

    if (Tools::isZero(time))
    {
        return signal;
    }

    // Modify step_size such, that after n iterations, time is really reached.

    step_size = time / std::ceil(time / step_size);

    // Compute the diffusion.

    const unsigned N = signal.size();
    std::vector<T> B(N);

    std::vector<T> diffused_signal = signal;

    for (T t = step_size; t <= time; t += step_size)
    {
        const std::vector<T> & derivative = secondDerivative(diffused_signal, interpolation_method);

        for (unsigned i = 0; i < signal.size(); ++i)
        {
            diffused_signal[i] = diffused_signal[i] + step_size * derivative[i];
        }
    }

    return diffused_signal;
}


template <class T>
const std::vector<T> nonlinearIsotropicDiffusion(const std::vector<T> & signal,
                                                 T time,
                                                 T step_size = 0.25,
                                                 T alpha = 0.0,
                                                 BorderInterpolation interpolation_method = CONTINUED)
{
    /* The general formulation of a nonlinear diffusion is
     *
     *      d/dt f(x, t) = div [ D(grad f(x, t)) * grad f(x, t) ]
     *
     * where D is the so-called diffusivity which we define to be
     *
     *      D(x) = 1 / abs(x)^alpha
     *
     * In the 1-dimensional case, the divergence and the gradient simplify to
     * an ordinary differentiation. Thus we get
     *
     *      d/dt f(x, t) = d/dx [ D(d/dx f(x, t)) * d/dx f(x, t) ]
     *
     * The diffusion equation can be easily solved by using the following
     * approximation
     *
     *      d/dt f(x, t) = [f(x, t+h) - f(x, t) ] / h
     *
     * where h is the step size. With this, we get (A, B, and C are used within
     * the implementation)
     *
     *      d/dt f(x, t+h) = f(x, t) + h *  { d/dx [ D(d/dx f(x, t)) * d/dx f(x, t) ] }
     *                     = f(x, t) + h *  { d/dx [ D(A)            * A            ] }
     *                     = f(x, t) + h *  { d/dx [ B                              ] }
     *                     = f(x, t) + h *  { C                                       }
     *
     * Iterating until t + n * h = time, yields the solution.
     */


    // If time == 0 return the unchanged signal.

    if (Tools::isZero(time))
    {
        return signal;
    }

    // Modify step_size such, that after n iterations, time is really reached.

    step_size = time / std::ceil(time / step_size);

    // Compute the diffusion.

    const unsigned N = signal.size();
    std::vector<T> B(N);

    std::vector<T> diffused_signal = signal;

    for (T t = step_size; t <= time; t += step_size)
    {
        const std::vector<T> & A = firstDerivative(diffused_signal, interpolation_method);

        for (unsigned i = 0; i < N; ++i)
        {
            const double epsilon = 1e-14;
            B[i] = 1 / std::pow(std::abs(A[i] + epsilon), alpha) * A[i];
        }

        const std::vector<T> & C = firstDerivative(B, interpolation_method);

        for (unsigned i = 0; i < signal.size(); ++i)
        {
            diffused_signal[i] = diffused_signal[i] + step_size * C[i];
        }
    }

    return diffused_signal;
}


template <class T>
const std::vector<T> nonlinearIsotropicDiffusion(const std::vector<T> & signal,
                                                 T (* diffusivity)(T),
                                                 T time,
                                                 T step_size = 0.25,
                                                 BorderInterpolation interpolation_method = CONTINUED)
{
    /* The general formulation of a nonlinear diffusion is
     *
     *      d/dt f(x, t) = div [ D(grad f(x, t)) * grad f(x, t) ]
     *
     * where D is the so-called diffusivity. In the 1-dimensional case, the
     * divergence and the gradient simplify to an ordinary differentiation.
     * Thus we get
     *
     *      d/dt f(x, t) = d/dx [ D(d/dx f(x, t)) * d/dx f(x, t) ]
     *
     * The diffusion equation can be easily solved by using the following
     * approximation
     *
     *      d/dt f(x, t) = [f(x, t+h) - f(x, t) ] / h
     *
     * where h is the step size. With this, we get (A, B, and C are used within
     * the implementation)
     *
     *      d/dt f(x, t+h) = f(x, t) + h *  { d/dx [ D(d/dx f(x, t)) * d/dx f(x, t) ] }
     *                     = f(x, t) + h *  { d/dx [ D(A)            * A            ] }
     *                     = f(x, t) + h *  { d/dx [ B                              ] }
     *                     = f(x, t) + h *  { C                                       }
     *
     * Iterating until t + n * h = time, yields the solution.
     */


    // If time == 0 return the unchanged signal.

    if (Tools::isZero(time))
    {
        return signal;
    }

    // Modify step_size such, that after n iterations, time is really reached.

    step_size = time / std::ceil(time / step_size);

    // Compute the diffusion.

    const unsigned N = signal.size();
    std::vector<T> B(N);

    std::vector<T> diffused_signal = signal;

    for (T t = step_size; t <= time; t += step_size)
    {
        const std::vector<T> & A = firstDerivative(diffused_signal, interpolation_method);

        for (unsigned i = 0; i < N; ++i)
        {
            B[i] = diffusivity(A[i]) * A[i];
        }

        const std::vector<T> & C = firstDerivative(B, interpolation_method);

        for (unsigned i = 0; i < signal.size(); ++i)
        {
            diffused_signal[i] = diffused_signal[i] + step_size * C[i];
        }
    }

    return diffused_signal;
}


} // namespace Diffusion


} // namespace TRTK


#endif // DIFFUSION_HPP_6431082743