LinearInterpolation.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2020.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/Types.h>
 
#include <cmath> // for modf() (which is an overloaded function in C++)
#include <vector>
 
namespace OpenMS
{
 
  namespace Math
  {
 
    template <typename Key = double, typename Value = Key>
    class LinearInterpolation
    {
 
public:
 
 
      typedef Value value_type;
 
      typedef Key key_type;
      typedef std::vector<value_type> container_type;
 
      typedef value_type      ValueType;
      typedef key_type        KeyType;
      typedef container_type  ContainerType;
 
public:
 
      LinearInterpolation(KeyType scale = 1., KeyType offset = 0.) :
        scale_(scale),
        offset_(offset),
        inside_(),
        outside_(),
        data_()
      {}
 
      LinearInterpolation(LinearInterpolation const & arg) :
        scale_(arg.scale_),
        offset_(arg.offset_),
        inside_(arg.inside_),
        outside_(arg.outside_),
        data_(arg.data_)
      {}
 
      LinearInterpolation & operator=(LinearInterpolation const & arg)
      {
        if (&arg == this)
          return *this;
 
        scale_   = arg.scale_;
        offset_  = arg.offset_;
        inside_  = arg.inside_;
        outside_ = arg.outside_;
        data_    = arg.data_;
 
        return *this;
      }
 
      ~LinearInterpolation() {}
 
      // ----------------------------------------------------------------------
 
 
 
      ValueType value(KeyType arg_pos) const
      {
 
        typedef typename container_type::difference_type DiffType;
 
        // apply the key transformation
        KeyType left_key;
        KeyType pos = key2index(arg_pos);
        KeyType frac = std::modf(pos, &left_key);
        DiffType const left = DiffType(left_key);
 
        // At left margin?
        if (pos < 0)
        {
          if (left /* <= -1 */)
          {
            return 0;
          }
          else        // left == 0
          {
            return data_[0] * (1 + frac);
          }
        }
        else         // pos >= 0
        {
          // At right margin?
          DiffType const back = data_.size() - 1;
          if (left >= back)
          {
            if (left != back)
            {
              return 0;
            }
            else
            {
              return data_[left] * (1 - frac);
            }
          }
          else
          {
            // In between!
            return data_[left + 1] * frac + data_[left] * (1 - frac);
          }
        }
      }
 
      void addValue(KeyType arg_pos, ValueType arg_value)
      {
 
        typedef typename container_type::difference_type DiffType;
 
        // apply the key transformation
        KeyType left_key;
        KeyType const pos = key2index(arg_pos);
        KeyType const frac = std::modf(pos, &left_key);
        DiffType const left = DiffType(left_key);
 
        // At left margin?
        if (pos < 0)
        {
          if (left /* <= -1 */)
          {
            return;
          }
          else        // left == 0
          {
            data_[0] += (1 + frac) * arg_value;
            return;
          }
        }
        else         // pos >= 0
        {
          // At right margin?
          DiffType const back = data_.size() - 1;
          if (left >= back)
          {
            if (left != back)
            {
              return;
            }
            else             // left == back
            {
              data_[left] += (1 - frac) * arg_value;
              return;
            }
          }
          else
          {
            // In between!
            data_[left + 1] += frac * arg_value;
            data_[left] += (1 - frac) * arg_value;
            return;
          }
        }
      }
 
      ValueType derivative(KeyType arg_pos) const
      {
 
        // apply the key transformation
        KeyType const pos = key2index(arg_pos);
 
        SignedSize const size_ = data_.size();
        SignedSize const left = int(pos + 0.5);       // rounds towards zero
 
        if (left < 0)           // quite small
        {
          return 0;
        }
        else
        {
          if (left == 0)             // at the border
          {
            if (pos >= -0.5)               // that is: -0.5 <= pos < +0.5
            {
              return (data_[1] - data_[0]) * (pos + 0.5) + (data_[0]) * (0.5 - pos);
            }
            else             // that is: -1.5 <= pos < -0.5
            {
              return (data_[0]) * (pos + 1.5);
            }
          }
        }
        // "else" case: to the right of the left margin
 
 
        KeyType factor = KeyType(left) - pos + KeyType(0.5);
 
        if (left > size_)           // quite large
        {
          return 0;
        }
        else
        {
          if (left < size_ - 1)             // to the left of the right margin
          {
            // weighted average of derivatives for adjacent intervals
            return (data_[left] - data_[left - 1]) * factor + (data_[left + 1] - data_[left]) * (1. - factor);
          }
          else           // somewhat at the border
          {
            // at the border, first case
            if (left == size_ - 1)
            {
              return (data_[left] - data_[left - 1]) * factor + (-data_[left]) * (1. - factor);
            }
          }
        }
        // else // that is: left == size_
 
        // We pull the last remaining case out of the "if" tree to avoid a
        // compiler warning ...
 
        // at the border, second case
        return (-data_[left - 1]) * factor;
      }
 
 
      // ----------------------------------------------------------------------
 
 
 
      ContainerType & getData()
      {
        return data_;
      }
 
      ContainerType const & getData() const
      {
        return data_;
      }
 
      template <typename SourceContainer>
      void setData(SourceContainer const & data)
      {
        data_ = data;
      }
 
      bool empty() const
      {
        return data_.empty();
      }
 
 
      // ----------------------------------------------------------------------
 
 
 
      KeyType key2index(KeyType pos) const
      {
        if (scale_)
        {
          pos -= offset_;
          pos /= scale_;
          return pos;
        }
        else
        {
          return 0;
        }
      }
 
      KeyType index2key(KeyType pos) const
      {
        pos *= scale_;
        pos += offset_;
        return pos;
      }
 
      KeyType const & getScale() const
      {
        return scale_;
      }
 
      void setScale(KeyType const & scale)
      {
        scale_ = scale;
      }
 
      KeyType const & getOffset() const
      {
        return offset_;
      }
 
      void setOffset(KeyType const & offset)
      {
        offset_ = offset;
      }
 
      void setMapping(KeyType const & scale, KeyType const & inside, KeyType const & outside)
      {
        scale_   = scale;
        inside_  = inside;
        outside_ = outside;
        offset_  = outside - scale * inside;
      }
 
      void setMapping(KeyType const & inside_low, KeyType const & outside_low,
                      KeyType const & inside_high, KeyType const & outside_high)
      {
        if (inside_high != inside_low)
        {
          setMapping((outside_high - outside_low) / (inside_high - inside_low),
                     inside_low, outside_low);
        }
        else
        {
          setMapping(0, inside_low, outside_low);
        }
        return;
      }
 
      KeyType const & getInsideReferencePoint() const
      {
        return inside_;
      }
 
      KeyType const & getOutsideReferencePoint() const
      {
        return outside_;
      }
 
      KeyType supportMin() const
      {
        return index2key(KeyType(empty() ? 0 : -1));
      }
 
      KeyType supportMax() const
      {
        return index2key(KeyType(data_.size()));
      }
 
 
protected:
 
      KeyType scale_;
      KeyType offset_;
      KeyType inside_;
      KeyType outside_;
 
      ContainerType data_;
 
    };
 
  }   // namespace Math
 
} // namespace OpenMS