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

Go to the documentation of this file.

/*
    Three-dimensional interpolation for irregularly-spaced data.

    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.2 (2013-08-20)
*/

#ifndef INTERPOLATION_3D_HPP_4326730820
#define INTERPOLATION_3D_HPP_4326730820


#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>

#include <flann/flann.hpp>

#include "Coordinate.hpp"
#include "Tools.hpp"

/*

All data items are stored in a kd-tree, internally. If an algorithm needs to
access them, the tree must be traversed. Since the tree behaves like normal
container, this can be easily done with iterators .

The algorithms are called via a function pointer which is set by the constructor
or a setter-method.

*/

namespace TRTK
{


template<class ValueType = double,
         class PointType = Coordinate<double>,
         class DataType  = Coordinate<double> >
class Interpolation3D
{
public:

    typedef ValueType value_type;
    typedef PointType point_type;
    typedef DataType  data_type;

    typedef PointType Point;            

    enum InterpolationMethod
    {
        INVERSE_DISTANCE_WEIGHTING,     
        INVERSE_DISTANCE_WEIGHTING_KNN, 
        NEAREST_NEIGHBOR                
    };

    Interpolation3D(const DataType & default_value = DataType(),
                    InterpolationMethod interpolation_method = INVERSE_DISTANCE_WEIGHTING_KNN);

    Interpolation3D(const std::vector<Point> & data_points,
                    const std::vector<DataType> & data,
                    const DataType & default_value = DataType(),
                    InterpolationMethod interpolation_method = INVERSE_DISTANCE_WEIGHTING_KNN);

    ~Interpolation3D();

    DataType operator()(const Point & point) const;
    DataType operator()(ValueType x, ValueType y, ValueType z) const;

    DataType getData(const Point & point) const;
    DataType getData(ValueType x, ValueType y, ValueType z) const;

    void setData(const std::vector<Point> & data_points,
                 const std::vector<DataType> & data);

    void setExponent(ValueType value = 2.0);

    void setInterpolationMethod(InterpolationMethod interpolation_method);

    void setNumberOfNearestNeighbors(unsigned number = 10);

    void setRadius(ValueType radius);

protected:

    DataType (Interpolation3D::*interpolation_function)(ValueType, ValueType, ValueType) const;

    DataType getDataInverseDistanceWeighting(ValueType x, ValueType y, ValueType z) const;
    DataType getDataInverseDistanceWeightingKNN(ValueType x, ValueType y, ValueType z) const;
    DataType getDataNearestNeighbor(ValueType x, ValueType y, ValueType z) const;

    ValueType exponent;
    ValueType radius;

    unsigned max_number_of_nearest_neighbors;

    flann::Index<flann::L2<ValueType> > * tree;

    const DataType data_type_default_value;

    std::vector<DataType> data;

    flann::Matrix<ValueType> data_points;
    mutable flann::Matrix<ValueType> query_point;
    mutable flann::Matrix<ValueType> distances;
    mutable flann::Matrix<int> indices;
};


template<class ValueType, class PointType, class DataType>
Interpolation3D<ValueType, PointType, DataType>::Interpolation3D(const DataType & default_value, InterpolationMethod interpolation_method) :
    exponent(2),
    radius(1),
    tree(NULL),
    data_type_default_value(default_value)
{
    query_point = flann::Matrix<ValueType>(new ValueType[3], 1, 3);

    setInterpolationMethod(interpolation_method);
    setNumberOfNearestNeighbors(10);
}


template<class ValueType, class PointType, class DataType>
Interpolation3D<ValueType, PointType, DataType>::Interpolation3D(const std::vector<Point> & data_points,
                                                                 const std::vector<DataType> & data,
                                                                 const DataType & default_value,
                                                                 InterpolationMethod interpolation_method) :
    exponent(2),
    radius(1),
    tree(NULL),
    data_type_default_value(default_value)
{
    query_point = flann::Matrix<ValueType>(new ValueType[3], 1, 3);

    setData(data_points, data);
    setInterpolationMethod(interpolation_method);
    setNumberOfNearestNeighbors(10);
}


template<class ValueType, class PointType, class DataType>
Interpolation3D<ValueType, PointType, DataType>::~Interpolation3D()
{
    if (tree != NULL) delete tree;

    data_points.free();
    query_point.free();
    distances.free();
    indices.free();
}


template<class ValueType, class PointType, class DataType>
inline DataType Interpolation3D<ValueType, PointType, DataType>::operator()(const Point & point) const
{
    return (this->*interpolation_function)(point.x(), point.y(), point.z());
}


template<class ValueType, class PointType, class DataType>
inline DataType Interpolation3D<ValueType, PointType, DataType>::operator()(ValueType x, ValueType y, ValueType z) const
{
    return (this->*interpolation_function)(x, y, z);
}


template<class ValueType, class PointType, class DataType>
inline DataType Interpolation3D<ValueType, PointType, DataType>::getData(const Point & point) const
{
    return (this->*interpolation_function)(point.x(), point.y(), point.z());
}


template<class ValueType, class PointType, class DataType>
inline DataType Interpolation3D<ValueType, PointType, DataType>::getData(ValueType x, ValueType y, ValueType z) const
{
    return (this->*interpolation_function)(x, y, z);
}


template<class ValueType, class PointType, class DataType>
DataType Interpolation3D<ValueType, PointType, DataType>::getDataInverseDistanceWeighting(ValueType x, ValueType y, ValueType z) const
{
    if (this->tree->size() == 0)
    {
        throw std::range_error("Interpolation3D::getDataInverseDistanceWeighting(): Tree is empty.");
    }

    query_point[0][0] = x;
    query_point[0][1] = y;
    query_point[0][2] = z;

    // Find the first nearest neighbor.

    unsigned number_of_found_nearest_neighbors = tree->radiusSearch(query_point, indices, distances, radius, flann::SearchParams(128));

    if (number_of_found_nearest_neighbors == 0)
    {
        throw std::runtime_error("Interpolation3D::getDataInverseDistanceWeighting(): Radius too small. No nearest neighbor found.");
    }

    // If the found point is the query point itself, just return the associated data...

    if (Tools::isZero(distances[0][0]))
    {
        return data[indices[0][0]];
    }

    // ... otherwise return an interpolated datum from the first n nearest neighbors.

    unsigned number_of_neighbors = std::min(number_of_found_nearest_neighbors, max_number_of_nearest_neighbors);

    ValueType normalization_coefficient = 0;
    DataType data = data_type_default_value;

    for (unsigned i = 0; i < number_of_neighbors; ++i)
    {
        using std::pow;
        ValueType weight = pow(distances[0][i], -exponent);

        normalization_coefficient += weight;

        data += this->data[indices[0][i]] * weight;
    }

    return data * (1 / normalization_coefficient);
}


template<class ValueType, class PointType, class DataType>
DataType Interpolation3D<ValueType, PointType, DataType>::getDataInverseDistanceWeightingKNN(ValueType x, ValueType y, ValueType z) const
{
    if (this->tree->size() == 0)
    {
        throw std::range_error("Interpolation3D::getDataInverseDistanceWeightingKNN(): Tree is empty.");
    }

    query_point[0][0] = x;
    query_point[0][1] = y;
    query_point[0][2] = z;

    // Find the first nearest neighbor.

    tree->knnSearch(query_point, indices, distances, max_number_of_nearest_neighbors, flann::SearchParams(128));

    // If the found point is the query point itself, just return the associated data...

    if (Tools::isZero(distances[0][0]))
    {
        return data[indices[0][0]];
    }

    // ... otherwise return an interpolated datum from the first n nearest neighbors.

    const unsigned tree_size = tree->size();

    ValueType normalization_coefficient = 0;
    DataType data = data_type_default_value;

    for (unsigned i = 0; i < max_number_of_nearest_neighbors && i < tree_size; ++i)
    {
        using std::pow;
        ValueType weight = pow(distances[0][i], -exponent);

        normalization_coefficient += weight;

        data += this->data[indices[0][i]] * weight;
    }

    return data * (1 / normalization_coefficient);
}


template<class ValueType, class PointType, class DataType>
DataType Interpolation3D<ValueType, PointType, DataType>::getDataNearestNeighbor(ValueType x, ValueType y, ValueType z) const
{
    if (this->tree->size() == 0)
    {
        throw std::range_error("Interpolation3D::getDataNearestNeighbor(): Tree is empty.");
    }

    query_point[0][0] = x;
    query_point[0][1] = y;
    query_point[0][2] = z;

    tree->knnSearch(query_point, indices, distances, 1, flann::SearchParams(128));

    return data[indices[0][0]];
}


template<class ValueType, class PointType, class DataType>
void Interpolation3D<ValueType, PointType, DataType>::setData(const std::vector<Point> & data_points,
                                                              const std::vector<DataType> & data)
{
    if (data_points.size() != data.size())
    {
        throw std::logic_error("Interpolation3D::setData(): 'data_points' and 'data' must have the same size.");
    }

    if (data_points.size() == 0)
    {
        throw std::runtime_error("Interpolation3D::setData(): The input data may not be emtpy.");
    }

    // Copy the data points as well as the associated data.

    const unsigned number_points = data_points.size();

    this->data_points.free();
    this->data_points = flann::Matrix<ValueType>(new ValueType[3 * number_points], number_points, 3);;

    for (unsigned i = 0; i < number_points; ++i)
    {
        this->data_points[i][0] = data_points[i][0];
        this->data_points[i][1] = data_points[i][1];
        this->data_points[i][2] = data_points[i][2];
    }

    this->data = data;

    // Setup tree.

    if (tree != NULL) delete tree;

    tree = new flann::Index<flann::L2<ValueType> >(this->data_points, flann::KDTreeSingleIndexParams());
    tree->buildIndex();
}


template<class ValueType, class PointType, class DataType>
void Interpolation3D<ValueType, PointType, DataType>::setExponent(ValueType value)
{
    this->exponent = value;
}


template<class ValueType, class PointType, class DataType>
void Interpolation3D<ValueType, PointType, DataType>::setInterpolationMethod(InterpolationMethod interpolation_method)
{
    switch (interpolation_method)
    {
        case INVERSE_DISTANCE_WEIGHTING:
        {
            interpolation_function = &Interpolation3D::getDataInverseDistanceWeighting;
            break;
        }

        case INVERSE_DISTANCE_WEIGHTING_KNN:
        {
            interpolation_function = &Interpolation3D::getDataInverseDistanceWeightingKNN;
            break;
        }

        case NEAREST_NEIGHBOR:
        {
            interpolation_function = &Interpolation3D::getDataNearestNeighbor;
            break;
        }

        default:
        {
            throw std::logic_error("Interpolation3D::setInterpolationMethod(): Unknown interpolation method.");
        }
    }
}


template<class ValueType, class PointType, class DataType>
void Interpolation3D<ValueType, PointType, DataType>::setNumberOfNearestNeighbors(unsigned number)
{
    if (number == 0)
    {
        throw std::logic_error("Interpolation3D::setNumberOfNearestNeighbors(): "
                               "Number of nearest neighbors must be greater than zero.");
    }

    this->max_number_of_nearest_neighbors = number;

    indices.free();
    distances.free();

    indices = flann::Matrix<int>(new int[number], 1, number);
    distances = flann::Matrix<ValueType>(new ValueType[number], 1, number);

}


template<class ValueType, class PointType, class DataType>
void Interpolation3D<ValueType, PointType, DataType>::setRadius(ValueType radius)
{
    if (radius <= 0)
    {
        throw std::logic_error("Interpolation3D::setRadius(): Radius must be greater than zero.");
    }

    this->radius = radius;
}


}


#endif // INTERPOLATION_3D_HPP_4326730820