BilinearInterpolation.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-2017.
//
// 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: $
// --------------------------------------------------------------------------

#ifndef OPENMS_MATH_MISC_BILINEARINTERPOLATION_H
#define OPENMS_MATH_MISC_BILINEARINTERPOLATION_H

#include <OpenMS/DATASTRUCTURES/Matrix.h>

namespace OpenMS
{

  namespace Math
  {

    template <typename Key = double, typename Value = Key>
    class BilinearInterpolation
    {

public:


      typedef Value value_type;

      typedef Key key_type;
      typedef Matrix<value_type> container_type;

      typedef value_type ValueType;
      typedef key_type KeyType;
      typedef container_type ContainerType;

public:


      BilinearInterpolation() :
        scale_0_(1),
        offset_0_(0),
        scale_1_(1),
        offset_1_(0),
        inside_0_(0),
        outside_0_(0),
        inside_1_(0),
        outside_1_(0),
        data_()
      {}

      BilinearInterpolation(BilinearInterpolation const & arg) :
        scale_0_(arg.scale_0_),
        offset_0_(arg.offset_0_),
        scale_1_(arg.scale_1_),
        offset_1_(arg.offset_1_),
        inside_0_(arg.inside_0_),
        outside_0_(arg.outside_0_),
        inside_1_(arg.inside_1_),
        outside_1_(arg.outside_1_),
        data_(arg.data_)
      {}

      BilinearInterpolation & operator=(BilinearInterpolation const & arg)
      {
        if (&arg == this)
          return *this;

        scale_0_ = arg.scale_0_;
        offset_0_ = arg.offset_0_;
        scale_1_ = arg.scale_1_;
        offset_1_ = arg.offset_1_;
        inside_0_ = arg.inside_0_;
        outside_1_ = arg.outside_1_;
        inside_1_ = arg.inside_1_;
        outside_0_ = arg.outside_0_;
        data_ = arg.data_;
        return *this;
      }

      ~BilinearInterpolation()
      {}


      // ----------------------------------------------------------------------



      ValueType value(KeyType arg_pos_0, KeyType arg_pos_1) const
      {
        // apply the key transformations
        KeyType const pos_0 = key2index_0(arg_pos_0);
        KeyType const pos_1 = key2index_1(arg_pos_1);

        // ???? should use modf() here!

        SignedSize const size_0 = data_.rows();
        SignedSize const lower_0 = SignedSize(pos_0);           // this rounds towards zero
        SignedSize const size_1 = data_.cols();
        SignedSize const lower_1 = SignedSize(pos_1);           // this rounds towards zero

        // small pos_0
        if (pos_0 <= 0)
        {
          if (lower_0 != 0)
          {
            return 0;
          }
          else        // that is: -1 < pos_0 <= 0
          {             // small pos_1
            if (pos_1 <= 0)
            {
              if (lower_1 != 0)
              {
                return 0;
              }
              else            // that is: -1 < pos_1 <= 0
              {
                return data_(0, 0) * (1. + pos_0) * (1. + pos_1);
              }
            }

            // big pos_1
            if (lower_1 >= size_1 - 1)
            {
              if (lower_1 != size_1 - 1)
              {
                return 0;
              }
              else
              {
                return data_(0, lower_1) * (1. + pos_0) * (size_1 - pos_1);
              }
            }

            // medium pos_1
            KeyType const factor_1 = pos_1 - KeyType(lower_1);
            KeyType const factor_1_complement = KeyType(1.) - factor_1;
            return (
                     data_(0, lower_1 + 1) * factor_1 +
                     data_(0, lower_1) * factor_1_complement
                     ) * (1. + pos_0);
          }
        }

        // big pos_0
        if (lower_0 >= size_0 - 1)
        {
          if (lower_0 != size_0 - 1)
          {
            return 0;
          }
          else        // that is: size_0 - 1 <= pos_0 < size_0
          {             // small pos_1
            if (pos_1 <= 0)
            {
              if (lower_1 != 0)
              {
                return 0;
              }
              else            // that is: -1 < pos_1 <= 0
              {
                return data_(lower_0, 0) * (size_0 - pos_0) * (1. + pos_1);
              }
            }

            // big pos_1
            if (lower_1 >= size_1 - 1)
            {
              if (lower_1 != size_1 - 1)
              {
                return 0;
              }
              else
              {
                return data_(lower_0, lower_1) * (size_0 - pos_0) * (size_1 - pos_1);
              }
            }

            // medium pos_1
            KeyType const factor_1 = pos_1 - KeyType(lower_1);
            KeyType const factor_1_complement = KeyType(1.) - factor_1;
            return (
                     data_(lower_0, lower_1 + 1) * factor_1 +
                     data_(lower_0, lower_1) * factor_1_complement
                     )
                   * (size_0 - pos_0);
          }
        }

        // medium pos_0
        {
          KeyType const factor_0 = pos_0 - KeyType(lower_0);
          KeyType const factor_0_complement = KeyType(1.) - factor_0;

          // small pos_1
          if (pos_1 <= 0)
          {
            if (lower_1 != 0)
            {
              return 0;
            }
            else          // that is: -1 < pos_1 <= 0
            {
              return (
                       data_(lower_0 + 1, 0) * factor_0
                       +
                       data_(lower_0, 0) * factor_0_complement
                       )
                     * (1. + pos_1);
            }
          }

          // big pos_1
          if (lower_1 >= size_1 - 1)
          {
            if (lower_1 != size_1 - 1)
            {
              return 0;
            }
            else
            {
              return (
                       data_(lower_0 + 1, lower_1) * factor_0
                       +
                       data_(lower_0, lower_1) * factor_0_complement
                       )
                     * (size_1 - pos_1);
            }
          }
          KeyType const factor_1 = pos_1 - KeyType(lower_1);
          KeyType const factor_1_complement = KeyType(1.) - factor_1;

          // medium pos_0 and medium pos_1 --> "within" the matrix
          return (
                   data_(lower_0 + 1, lower_1 + 1) * factor_0
                   +
                   data_(lower_0, lower_1 + 1) * factor_0_complement
                   )
                 * factor_1
                 +
                 (
                   data_(lower_0 + 1, lower_1) * factor_0
                   +
                   data_(lower_0, lower_1) * factor_0_complement
                 )
                 * factor_1_complement;
        }
      }

      void addValue(KeyType arg_pos_0, KeyType arg_pos_1, ValueType arg_value)
      {

        typedef typename container_type::difference_type DiffType;

        // apply key transformation _0
        KeyType const pos_0 = key2index_0(arg_pos_0);
        KeyType lower_0_key;
        KeyType const frac_0 = std::modf(pos_0, &lower_0_key);
        DiffType const lower_0 = DiffType(lower_0_key);

        // Small pos_0 ?
        if (pos_0 < 0)
        {
          if (lower_0)
          {
            return;
          }
          else        // lower_0 == 0
          {             // apply key transformation _1
            KeyType const pos_1 = key2index_1(arg_pos_1);
            KeyType lower_1_key;
            KeyType const frac_1 = std::modf(pos_1, &lower_1_key);
            DiffType const lower_1 = DiffType(lower_1_key);

            // Small pos_1 ?
            if (pos_1 < 0)
            {
              if (lower_1)
              {
                return;
              }
              else            // lower_1 == 0
              {
                data_(0, 0) += arg_value * (1 + frac_0) * (1 + frac_1);
                return;
              }
            }
            else             // pos_1 >= 0
            {
              DiffType const back_1 = data_.cols() - 1;
              // big pos_1
              if (lower_1 >= back_1)
              {
                if (lower_1 != back_1)
                {
                  return;
                }
                else                 // lower_1 == back_1
                {
                  data_(0, lower_1) += arg_value * (1 + frac_0) * (1 - frac_1);
                  return;
                }
              }
              else
              {
                // medium pos_1
                KeyType const tmp_prod = KeyType(arg_value * (1. + frac_0));
                data_(0, lower_1 + 1) += tmp_prod * frac_1;
                data_(0, lower_1) += tmp_prod * (1. - frac_1);
                return;
              }
            }
          }
        }
        else         // pos_0 >= 0
        {
          DiffType const back_0 = data_.rows() - 1;
          if (lower_0 >= back_0)
          {
            if (lower_0 != back_0)
            {
              return;
            }
            else             // lower_0 == back_0
            {

              KeyType const tmp_prod = KeyType(arg_value * (1. - frac_0));

              // apply key transformation _1
              KeyType const pos_1 = key2index_1(arg_pos_1);
              KeyType lower_1_key;
              KeyType const frac_1 = std::modf(pos_1, &lower_1_key);
              DiffType const lower_1 = DiffType(lower_1_key);

              // Small pos_1 ?
              if (pos_1 < 0)
              {
                if (lower_1)
                {
                  return;
                }
                else              // lower_1 == 0
                {
                  data_(lower_0, 0) += tmp_prod * (1 + frac_1);
                  return;
                }
              }
              else               // pos_1 >= 0
              {
                DiffType const back_1 = data_.cols() - 1;
                // big pos_1
                if (lower_1 >= back_1)
                {
                  if (lower_1 != back_1)
                  {
                    return;
                  }
                  else                   // lower_1 == back_1
                  {
                    data_(lower_0, lower_1) += tmp_prod * (1 - frac_1);
                    return;
                  }
                }
                else
                {
                  // medium pos_1
                  data_(lower_0, lower_1 + 1) += tmp_prod * frac_1;
                  data_(lower_0, lower_1) += tmp_prod * (1 - frac_1);
                  return;
                }
              }
            }
          }
          else           // lower_0 < back_0
          {

            // Medium pos_0 !

            // apply key transformation _1
            KeyType const pos_1 = key2index_1(arg_pos_1);
            KeyType lower_1_key;
            KeyType const frac_1 = std::modf(pos_1, &lower_1_key);
            DiffType const lower_1 = DiffType(lower_1_key);

            // Small pos_1 ?
            if (pos_1 < 0)
            {
              if (lower_1)
              {
                return;
              }
              else            // lower_1 == 0
              {
                KeyType const tmp_prod = KeyType(arg_value * (1 + frac_1));
                data_(lower_0 + 1, 0) += tmp_prod * frac_0;
                data_(lower_0, 0) += tmp_prod * (1 - frac_0);
                return;
              }
            }
            else             // pos_1 >= 0
            {
              DiffType const back_1 = data_.cols() - 1;
              // big pos_1
              if (lower_1 >= back_1)
              {
                if (lower_1 != back_1)
                {
                  return;
                }
                else                 // lower_1 == back_1
                {
                  KeyType const tmp_prod = KeyType(arg_value * (1 - frac_1));
                  data_(lower_0 + 1, lower_1) += tmp_prod * frac_0;
                  data_(lower_0, lower_1) += tmp_prod * (1 - frac_0);
                  return;
                }
              }
              else
              {
                // Medium pos_1 !

                // medium pos_0 and medium pos_1 --> "within" the matrix
                KeyType tmp_prod = KeyType(arg_value * frac_0);
                data_(lower_0 + 1, lower_1 + 1) += tmp_prod * frac_1;
                data_(lower_0 + 1, lower_1) += tmp_prod * (1 - frac_1);
                tmp_prod = KeyType(arg_value * (1 - frac_0));
                data_(lower_0, lower_1 + 1) += tmp_prod * frac_1;
                data_(lower_0, lower_1) += tmp_prod * (1 - frac_1);
                return;
              }
            }
          }
        }
      }


      // ----------------------------------------------------------------------



      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_0(KeyType pos) const
      {
        if (scale_0_)
        {
          pos -= offset_0_;
          pos /= scale_0_;
          return pos;
        }
        else
        {
          return 0;
        }
      }

      KeyType index2key_0(KeyType pos) const
      {
        pos *= scale_0_;
        pos += offset_0_;
        return pos;
      }

      KeyType key2index_1(KeyType pos) const
      {
        if (scale_1_)
        {
          pos -= offset_1_;
          pos /= scale_1_;
          return pos;
        }
        else
        {
          return 0;
        }
      }

      KeyType index2key_1(KeyType pos) const
      {
        pos *= scale_1_;
        pos += offset_1_;
        return pos;
      }

      KeyType const & getScale_0() const
      {
        return scale_0_;
      }

      KeyType const & getScale_1() const
      {
        return scale_1_;
      }

      void setScale_0(KeyType const & scale)
      {
        scale_0_ = scale;
      }

      void setScale_1(KeyType const & scale)
      {
        scale_1_ = scale;
      }

      KeyType const & getOffset_0() const
      {
        return offset_0_;
      }

      KeyType const & getOffset_1() const
      {
        return offset_1_;
      }

      void setOffset_0(KeyType const & offset)
      {
        offset_0_ = offset;
      }

      void setOffset_1(KeyType const & offset)
      {
        offset_1_ = offset;
      }

      void setMapping_0(KeyType const & scale, KeyType const & inside_low, KeyType const & outside_low)
      {
        scale_0_ = scale;
        inside_0_ = inside_low;
        outside_0_ = outside_low;
        offset_0_ = outside_low - scale * inside_low;
        return;
      }

      void setMapping_0(KeyType const & inside_low, KeyType const & outside_low,
                        KeyType const & inside_high, KeyType const & outside_high)
      {
        if (inside_high != inside_low)
        {
          setMapping_0((outside_high - outside_low) / (inside_high - inside_low),
                       inside_low, outside_low);
        }
        else
        {
          setMapping_0(0, inside_low, outside_low);
        }
        return;
      }

      void setMapping_1(KeyType const & scale, KeyType const & inside_low, KeyType const & outside_low)
      {
        scale_1_ = scale;
        inside_1_ = inside_low;
        outside_1_ = outside_low;
        offset_1_ = outside_low - scale * inside_low;
        return;
      }

      void setMapping_1(KeyType const & inside_low, KeyType const & outside_low,
                        KeyType const & inside_high, KeyType const & outside_high)
      {
        if (inside_high != inside_low)
        {
          setMapping_1((outside_high - outside_low) / (inside_high - inside_low),
                       inside_low, outside_low);
        }
        else
        {
          setMapping_1(0, inside_low, outside_low);
        }
        return;
      }

      KeyType const & getInsideReferencePoint_0() const
      {
        return inside_0_;
      }

      KeyType const & getInsideReferencePoint_1() const
      {
        return inside_1_;
      }

      KeyType const & getOutsideReferencePoint_0() const
      {
        return outside_0_;
      }

      KeyType const & getOutsideReferencePoint_1() const
      {
        return outside_1_;
      }

      KeyType supportMin_0() const
      {
        return index2key_0(empty() ? KeyType(0.) : KeyType(-1.));
      }

      KeyType supportMin_1() const
      {
        return index2key_1(empty() ? KeyType(0.) : KeyType(-1.));
      }

      KeyType supportMax_0() const
      {
        return index2key_0(KeyType(data_.rows()));
      }

      KeyType supportMax_1() const
      {
        return index2key_1(KeyType(data_.cols()));
      }


protected:

      KeyType scale_0_;
      KeyType offset_0_;
      KeyType scale_1_;
      KeyType offset_1_;
      KeyType inside_0_;
      KeyType outside_0_;
      KeyType inside_1_;
      KeyType outside_1_;
      ContainerType data_;
    };

  }   // namespace Math

} // namespace OpenMS

#endif // OPENMS_MATH_MISC_BILINEARINTERPOLATION_H