OpenMS
BilinearInterpolation.h
Go to the documentation of this file.
1 // --------------------------------------------------------------------------
2 // OpenMS -- Open-Source Mass Spectrometry
3 // --------------------------------------------------------------------------
4 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
5 // ETH Zurich, and Freie Universitaet Berlin 2002-2023.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Timo Sachsenberg $
32 // $Authors: $
33 // --------------------------------------------------------------------------
34 
35 #pragma once
36 
38 
39 namespace OpenMS
40 {
41 
42  namespace Math
43  {
44 
71  template <typename Key = double, typename Value = Key>
73  {
74 
75 public:
76 
78 
79  typedef Value value_type;
80 
81  typedef Key key_type;
83 
85  typedef key_type KeyType;
88 
89 public:
90 
94 
97  scale_0_(1),
98  offset_0_(0),
99  scale_1_(1),
100  offset_1_(0),
101  inside_0_(0),
102  outside_0_(0),
103  inside_1_(0),
104  outside_1_(0),
105  data_()
106  {}
107 
110  scale_0_(arg.scale_0_),
111  offset_0_(arg.offset_0_),
112  scale_1_(arg.scale_1_),
113  offset_1_(arg.offset_1_),
114  inside_0_(arg.inside_0_),
115  outside_0_(arg.outside_0_),
116  inside_1_(arg.inside_1_),
117  outside_1_(arg.outside_1_),
118  data_(arg.data_)
119  {}
120 
123  {
124  if (&arg == this)
125  return *this;
126 
127  scale_0_ = arg.scale_0_;
128  offset_0_ = arg.offset_0_;
129  scale_1_ = arg.scale_1_;
130  offset_1_ = arg.offset_1_;
131  inside_0_ = arg.inside_0_;
132  outside_1_ = arg.outside_1_;
133  inside_1_ = arg.inside_1_;
134  outside_0_ = arg.outside_0_;
135  data_ = arg.data_;
136  return *this;
137  }
138 
141  {}
142 
144 
145  // ----------------------------------------------------------------------
146 
148 
149 
151  ValueType value(KeyType arg_pos_0, KeyType arg_pos_1) const
152  {
153  // apply the key transformations
154  KeyType const pos_0 = key2index_0(arg_pos_0);
155  KeyType const pos_1 = key2index_1(arg_pos_1);
156 
157  // ???? should use modf() here!
158 
159  SignedSize const size_0 = data_.rows();
160  SignedSize const lower_0 = SignedSize(pos_0); // this rounds towards zero
161  SignedSize const size_1 = data_.cols();
162  SignedSize const lower_1 = SignedSize(pos_1); // this rounds towards zero
163 
164  // small pos_0
165  if (pos_0 <= 0)
166  {
167  if (lower_0 != 0)
168  {
169  return 0;
170  }
171  else // that is: -1 < pos_0 <= 0
172  { // small pos_1
173  if (pos_1 <= 0)
174  {
175  if (lower_1 != 0)
176  {
177  return 0;
178  }
179  else // that is: -1 < pos_1 <= 0
180  {
181  return data_(0, 0) * (1. + pos_0) * (1. + pos_1);
182  }
183  }
184 
185  // big pos_1
186  if (lower_1 >= size_1 - 1)
187  {
188  if (lower_1 != size_1 - 1)
189  {
190  return 0;
191  }
192  else
193  {
194  return data_(0, lower_1) * (1. + pos_0) * (size_1 - pos_1);
195  }
196  }
197 
198  // medium pos_1
199  KeyType const factor_1 = pos_1 - KeyType(lower_1);
200  KeyType const factor_1_complement = KeyType(1.) - factor_1;
201  return (
202  data_(0, lower_1 + 1) * factor_1 +
203  data_(0, lower_1) * factor_1_complement
204  ) * (1. + pos_0);
205  }
206  }
207 
208  // big pos_0
209  if (lower_0 >= size_0 - 1)
210  {
211  if (lower_0 != size_0 - 1)
212  {
213  return 0;
214  }
215  else // that is: size_0 - 1 <= pos_0 < size_0
216  { // small pos_1
217  if (pos_1 <= 0)
218  {
219  if (lower_1 != 0)
220  {
221  return 0;
222  }
223  else // that is: -1 < pos_1 <= 0
224  {
225  return data_(lower_0, 0) * (size_0 - pos_0) * (1. + pos_1);
226  }
227  }
228 
229  // big pos_1
230  if (lower_1 >= size_1 - 1)
231  {
232  if (lower_1 != size_1 - 1)
233  {
234  return 0;
235  }
236  else
237  {
238  return data_(lower_0, lower_1) * (size_0 - pos_0) * (size_1 - pos_1);
239  }
240  }
241 
242  // medium pos_1
243  KeyType const factor_1 = pos_1 - KeyType(lower_1);
244  KeyType const factor_1_complement = KeyType(1.) - factor_1;
245  return (
246  data_(lower_0, lower_1 + 1) * factor_1 +
247  data_(lower_0, lower_1) * factor_1_complement
248  )
249  * (size_0 - pos_0);
250  }
251  }
252 
253  // medium pos_0
254  {
255  KeyType const factor_0 = pos_0 - KeyType(lower_0);
256  KeyType const factor_0_complement = KeyType(1.) - factor_0;
257 
258  // small pos_1
259  if (pos_1 <= 0)
260  {
261  if (lower_1 != 0)
262  {
263  return 0;
264  }
265  else // that is: -1 < pos_1 <= 0
266  {
267  return (
268  data_(lower_0 + 1, 0) * factor_0
269  +
270  data_(lower_0, 0) * factor_0_complement
271  )
272  * (1. + pos_1);
273  }
274  }
275 
276  // big pos_1
277  if (lower_1 >= size_1 - 1)
278  {
279  if (lower_1 != size_1 - 1)
280  {
281  return 0;
282  }
283  else
284  {
285  return (
286  data_(lower_0 + 1, lower_1) * factor_0
287  +
288  data_(lower_0, lower_1) * factor_0_complement
289  )
290  * (size_1 - pos_1);
291  }
292  }
293  KeyType const factor_1 = pos_1 - KeyType(lower_1);
294  KeyType const factor_1_complement = KeyType(1.) - factor_1;
295 
296  // medium pos_0 and medium pos_1 --> "within" the matrix
297  return (
298  data_(lower_0 + 1, lower_1 + 1) * factor_0
299  +
300  data_(lower_0, lower_1 + 1) * factor_0_complement
301  )
302  * factor_1
303  +
304  (
305  data_(lower_0 + 1, lower_1) * factor_0
306  +
307  data_(lower_0, lower_1) * factor_0_complement
308  )
309  * factor_1_complement;
310  }
311  }
312 
316  void addValue(KeyType arg_pos_0, KeyType arg_pos_1, ValueType arg_value)
317  {
318 
319  typedef typename container_type::difference_type DiffType;
320 
321  // apply key transformation _0
322  KeyType const pos_0 = key2index_0(arg_pos_0);
323  KeyType lower_0_key;
324  KeyType const frac_0 = std::modf(pos_0, &lower_0_key);
325  DiffType const lower_0 = DiffType(lower_0_key);
326 
327  // Small pos_0 ?
328  if (pos_0 < 0)
329  {
330  if (lower_0)
331  {
332  return;
333  }
334  else // lower_0 == 0
335  { // apply key transformation _1
336  KeyType const pos_1 = key2index_1(arg_pos_1);
337  KeyType lower_1_key;
338  KeyType const frac_1 = std::modf(pos_1, &lower_1_key);
339  DiffType const lower_1 = DiffType(lower_1_key);
340 
341  // Small pos_1 ?
342  if (pos_1 < 0)
343  {
344  if (lower_1)
345  {
346  return;
347  }
348  else // lower_1 == 0
349  {
350  data_(0, 0) += arg_value * (1 + frac_0) * (1 + frac_1);
351  return;
352  }
353  }
354  else // pos_1 >= 0
355  {
356  DiffType const back_1 = data_.cols() - 1;
357  // big pos_1
358  if (lower_1 >= back_1)
359  {
360  if (lower_1 != back_1)
361  {
362  return;
363  }
364  else // lower_1 == back_1
365  {
366  data_(0, lower_1) += arg_value * (1 + frac_0) * (1 - frac_1);
367  return;
368  }
369  }
370  else
371  {
372  // medium pos_1
373  KeyType const tmp_prod = KeyType(arg_value * (1. + frac_0));
374  data_(0, lower_1 + 1) += tmp_prod * frac_1;
375  data_(0, lower_1) += tmp_prod * (1. - frac_1);
376  return;
377  }
378  }
379  }
380  }
381  else // pos_0 >= 0
382  {
383  DiffType const back_0 = data_.rows() - 1;
384  if (lower_0 >= back_0)
385  {
386  if (lower_0 != back_0)
387  {
388  return;
389  }
390  else // lower_0 == back_0
391  {
392 
393  KeyType const tmp_prod = KeyType(arg_value * (1. - frac_0));
394 
395  // apply key transformation _1
396  KeyType const pos_1 = key2index_1(arg_pos_1);
397  KeyType lower_1_key;
398  KeyType const frac_1 = std::modf(pos_1, &lower_1_key);
399  DiffType const lower_1 = DiffType(lower_1_key);
400 
401  // Small pos_1 ?
402  if (pos_1 < 0)
403  {
404  if (lower_1)
405  {
406  return;
407  }
408  else // lower_1 == 0
409  {
410  data_(lower_0, 0) += tmp_prod * (1 + frac_1);
411  return;
412  }
413  }
414  else // pos_1 >= 0
415  {
416  DiffType const back_1 = data_.cols() - 1;
417  // big pos_1
418  if (lower_1 >= back_1)
419  {
420  if (lower_1 != back_1)
421  {
422  return;
423  }
424  else // lower_1 == back_1
425  {
426  data_(lower_0, lower_1) += tmp_prod * (1 - frac_1);
427  return;
428  }
429  }
430  else
431  {
432  // medium pos_1
433  data_(lower_0, lower_1 + 1) += tmp_prod * frac_1;
434  data_(lower_0, lower_1) += tmp_prod * (1 - frac_1);
435  return;
436  }
437  }
438  }
439  }
440  else // lower_0 < back_0
441  {
442 
443  // Medium pos_0 !
444 
445  // apply key transformation _1
446  KeyType const pos_1 = key2index_1(arg_pos_1);
447  KeyType lower_1_key;
448  KeyType const frac_1 = std::modf(pos_1, &lower_1_key);
449  DiffType const lower_1 = DiffType(lower_1_key);
450 
451  // Small pos_1 ?
452  if (pos_1 < 0)
453  {
454  if (lower_1)
455  {
456  return;
457  }
458  else // lower_1 == 0
459  {
460  KeyType const tmp_prod = KeyType(arg_value * (1 + frac_1));
461  data_(lower_0 + 1, 0) += tmp_prod * frac_0;
462  data_(lower_0, 0) += tmp_prod * (1 - frac_0);
463  return;
464  }
465  }
466  else // pos_1 >= 0
467  {
468  DiffType const back_1 = data_.cols() - 1;
469  // big pos_1
470  if (lower_1 >= back_1)
471  {
472  if (lower_1 != back_1)
473  {
474  return;
475  }
476  else // lower_1 == back_1
477  {
478  KeyType const tmp_prod = KeyType(arg_value * (1 - frac_1));
479  data_(lower_0 + 1, lower_1) += tmp_prod * frac_0;
480  data_(lower_0, lower_1) += tmp_prod * (1 - frac_0);
481  return;
482  }
483  }
484  else
485  {
486  // Medium pos_1 !
487 
488  // medium pos_0 and medium pos_1 --> "within" the matrix
489  KeyType tmp_prod = KeyType(arg_value * frac_0);
490  data_(lower_0 + 1, lower_1 + 1) += tmp_prod * frac_1;
491  data_(lower_0 + 1, lower_1) += tmp_prod * (1 - frac_1);
492  tmp_prod = KeyType(arg_value * (1 - frac_0));
493  data_(lower_0, lower_1 + 1) += tmp_prod * frac_1;
494  data_(lower_0, lower_1) += tmp_prod * (1 - frac_1);
495  return;
496  }
497  }
498  }
499  }
500  }
501 
503 
504  // ----------------------------------------------------------------------
505 
507 
508 
511  {
512  return data_;
513  }
514 
516  ContainerType const & getData() const
517  {
518  return data_;
519  }
520 
526  template <typename SourceContainer>
527  void setData(SourceContainer const & data)
528  {
529  data_ = data;
530  }
531 
533  bool empty() const
534  {
535  return data_.empty();
536  }
537 
539 
540  // ----------------------------------------------------------------------
541 
543 
544 
547  {
548  if (scale_0_)
549  {
550  pos -= offset_0_;
551  pos /= scale_0_;
552  return pos;
553  }
554  else
555  {
556  return 0;
557  }
558  }
559 
562  {
563  pos *= scale_0_;
564  pos += offset_0_;
565  return pos;
566  }
567 
570  {
571  if (scale_1_)
572  {
573  pos -= offset_1_;
574  pos /= scale_1_;
575  return pos;
576  }
577  else
578  {
579  return 0;
580  }
581  }
582 
585  {
586  pos *= scale_1_;
587  pos += offset_1_;
588  return pos;
589  }
590 
592  KeyType const & getScale_0() const
593  {
594  return scale_0_;
595  }
596 
598  KeyType const & getScale_1() const
599  {
600  return scale_1_;
601  }
602 
608  void setScale_0(KeyType const & scale)
609  {
610  scale_0_ = scale;
611  }
612 
618  void setScale_1(KeyType const & scale)
619  {
620  scale_1_ = scale;
621  }
622 
624  KeyType const & getOffset_0() const
625  {
626  return offset_0_;
627  }
628 
630  KeyType const & getOffset_1() const
631  {
632  return offset_1_;
633  }
634 
641  void setOffset_0(KeyType const & offset)
642  {
643  offset_0_ = offset;
644  }
645 
652  void setOffset_1(KeyType const & offset)
653  {
654  offset_1_ = offset;
655  }
656 
670  void setMapping_0(KeyType const & scale, KeyType const & inside_low, KeyType const & outside_low)
671  {
672  scale_0_ = scale;
673  inside_0_ = inside_low;
674  outside_0_ = outside_low;
675  offset_0_ = outside_low - scale * inside_low;
676  return;
677  }
678 
685  void setMapping_0(KeyType const & inside_low, KeyType const & outside_low,
686  KeyType const & inside_high, KeyType const & outside_high)
687  {
688  if (inside_high != inside_low)
689  {
690  setMapping_0((outside_high - outside_low) / (inside_high - inside_low),
691  inside_low, outside_low);
692  }
693  else
694  {
695  setMapping_0(0, inside_low, outside_low);
696  }
697  return;
698  }
699 
713  void setMapping_1(KeyType const & scale, KeyType const & inside_low, KeyType const & outside_low)
714  {
715  scale_1_ = scale;
716  inside_1_ = inside_low;
717  outside_1_ = outside_low;
718  offset_1_ = outside_low - scale * inside_low;
719  return;
720  }
721 
728  void setMapping_1(KeyType const & inside_low, KeyType const & outside_low,
729  KeyType const & inside_high, KeyType const & outside_high)
730  {
731  if (inside_high != inside_low)
732  {
733  setMapping_1((outside_high - outside_low) / (inside_high - inside_low),
734  inside_low, outside_low);
735  }
736  else
737  {
738  setMapping_1(0, inside_low, outside_low);
739  }
740  return;
741  }
742 
745  {
746  return inside_0_;
747  }
748 
751  {
752  return inside_1_;
753  }
754 
757  {
758  return outside_0_;
759  }
760 
763  {
764  return outside_1_;
765  }
766 
769  {
770  return index2key_0(empty() ? KeyType(0.) : KeyType(-1.));
771  }
772 
775  {
776  return index2key_1(empty() ? KeyType(0.) : KeyType(-1.));
777  }
778 
781  {
782  return index2key_0(KeyType(data_.rows()));
783  }
784 
787  {
788  return index2key_1(KeyType(data_.cols()));
789  }
790 
792 
793 protected:
794 
807  };
808 
809  } // namespace Math
810 
811 } // namespace OpenMS
812 
Provides access to bilinearly interpolated values (and derivatives) from discrete data points....
Definition: BilinearInterpolation.h:73
KeyType const & getOffset_1() const
Accessor. "Offset" is the point (in "outside" units) which corresponds to "Data(0,...
Definition: BilinearInterpolation.h:630
ContainerType const & getData() const
Returns the internal random access container storing the data.
Definition: BilinearInterpolation.h:516
Value value_type
Definition: BilinearInterpolation.h:79
KeyType supportMin_1() const
Lower boundary of the support, in "outside" coordinates.
Definition: BilinearInterpolation.h:774
ValueType value(KeyType arg_pos_0, KeyType arg_pos_1) const
Returns the interpolated value ("backward resampling")
Definition: BilinearInterpolation.h:151
KeyType const & getOffset_0() const
Accessor. "Offset" is the point (in "outside" units) which corresponds to "Data(0,...
Definition: BilinearInterpolation.h:624
Key key_type
Definition: BilinearInterpolation.h:81
Matrix< value_type > container_type
Definition: BilinearInterpolation.h:82
ContainerType data_
Definition: BilinearInterpolation.h:805
container_type ContainerType
Definition: BilinearInterpolation.h:86
void setData(SourceContainer const &data)
Assigns data to the internal random access container storing the data.
Definition: BilinearInterpolation.h:527
KeyType outside_0_
Definition: BilinearInterpolation.h:802
void setMapping_1(KeyType const &inside_low, KeyType const &outside_low, KeyType const &inside_high, KeyType const &outside_high)
Specifies the mapping from "outside" to "inside" coordinates by the following data:
Definition: BilinearInterpolation.h:728
KeyType supportMax_1() const
Upper boundary of the support, in "outside" coordinates.
Definition: BilinearInterpolation.h:786
BilinearInterpolation & operator=(BilinearInterpolation const &arg)
Assignment operator.
Definition: BilinearInterpolation.h:122
KeyType inside_1_
Definition: BilinearInterpolation.h:803
BilinearInterpolation(BilinearInterpolation const &arg)
Copy constructor.
Definition: BilinearInterpolation.h:109
KeyType const & getScale_1() const
Accessor. "Scale" is the difference (in "outside" units) between consecutive entries in "Data".
Definition: BilinearInterpolation.h:598
void setMapping_1(KeyType const &scale, KeyType const &inside_low, KeyType const &outside_low)
Specifies the mapping from "outside" to "inside" coordinates by the following data:
Definition: BilinearInterpolation.h:713
bool empty() const
Returns true if getData() is empty.
Definition: BilinearInterpolation.h:533
KeyType supportMax_0() const
Upper boundary of the support, in "outside" coordinates.
Definition: BilinearInterpolation.h:780
KeyType const & getOutsideReferencePoint_1() const
Accessor. See setMapping().
Definition: BilinearInterpolation.h:762
void setOffset_1(KeyType const &offset)
Accessor. "Offset" is the point (in "outside" units) which corresponds to "Data(0,...
Definition: BilinearInterpolation.h:652
KeyType const & getOutsideReferencePoint_0() const
Accessor. See setMapping().
Definition: BilinearInterpolation.h:756
void setMapping_0(KeyType const &scale, KeyType const &inside_low, KeyType const &outside_low)
Specifies the mapping from "outside" to "inside" coordinates by the following data:
Definition: BilinearInterpolation.h:670
KeyType const & getInsideReferencePoint_1() const
Accessor. See setMapping().
Definition: BilinearInterpolation.h:750
value_type ValueType
Definition: BilinearInterpolation.h:84
KeyType scale_1_
Definition: BilinearInterpolation.h:799
BilinearInterpolation()
Constructors and destructor.
Definition: BilinearInterpolation.h:96
KeyType const & getInsideReferencePoint_0() const
Accessor. See setMapping().
Definition: BilinearInterpolation.h:744
KeyType index2key_0(KeyType pos) const
The transformation from "inside" to "outside" coordinates.
Definition: BilinearInterpolation.h:561
KeyType offset_0_
Definition: BilinearInterpolation.h:798
KeyType scale_0_
Data members.
Definition: BilinearInterpolation.h:797
KeyType inside_0_
Definition: BilinearInterpolation.h:801
ContainerType & getData()
Returns the internal random access container storing the data.
Definition: BilinearInterpolation.h:510
void setMapping_0(KeyType const &inside_low, KeyType const &outside_low, KeyType const &inside_high, KeyType const &outside_high)
Specifies the mapping from "outside" to "inside" coordinates by the following data:
Definition: BilinearInterpolation.h:685
void addValue(KeyType arg_pos_0, KeyType arg_pos_1, ValueType arg_value)
Performs bilinear resampling. The arg_value is split up and added to the data points around arg_pos....
Definition: BilinearInterpolation.h:316
key_type KeyType
Definition: BilinearInterpolation.h:85
KeyType offset_1_
Definition: BilinearInterpolation.h:800
~BilinearInterpolation()
Destructor.
Definition: BilinearInterpolation.h:140
KeyType outside_1_
Definition: BilinearInterpolation.h:804
void setScale_1(KeyType const &scale)
Accessor. "Scale" is the difference (in "outside" units) between consecutive entries in "Data".
Definition: BilinearInterpolation.h:618
KeyType key2index_1(KeyType pos) const
The transformation from "outside" to "inside" coordinates.
Definition: BilinearInterpolation.h:569
void setOffset_0(KeyType const &offset)
Accessor. "Offset" is the point (in "outside" units) which corresponds to "Data(0,...
Definition: BilinearInterpolation.h:641
KeyType index2key_1(KeyType pos) const
The transformation from "inside" to "outside" coordinates.
Definition: BilinearInterpolation.h:584
KeyType const & getScale_0() const
Accessor. "Scale" is the difference (in "outside" units) between consecutive entries in "Data".
Definition: BilinearInterpolation.h:592
KeyType supportMin_0() const
Lower boundary of the support, in "outside" coordinates.
Definition: BilinearInterpolation.h:768
KeyType key2index_0(KeyType pos) const
The transformation from "outside" to "inside" coordinates.
Definition: BilinearInterpolation.h:546
void setScale_0(KeyType const &scale)
Accessor. "Scale" is the difference (in "outside" units) between consecutive entries in "Data".
Definition: BilinearInterpolation.h:608
A two-dimensional matrix. Similar to std::vector, but uses a binary operator(,) for element access.
Definition: Matrix.h:77
SizeType rows() const
Number of rows.
Definition: Matrix.h:258
Base::difference_type difference_type
Definition: Matrix.h:87
SizeType cols() const
Number of columns.
Definition: Matrix.h:264
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:48