[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

fixedpoint.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2004-2005 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_FIXEDPOINT_HXX
37 #define VIGRA_FIXEDPOINT_HXX
38 
39 #include "mathutil.hxx"
40 #include "static_assert.hxx"
41 #include "error.hxx"
42 #include "numerictraits.hxx"
43 
44 namespace vigra {
45 
46 template <unsigned IntBits, unsigned FractionalBits>
47 class FixedPoint;
48 
49 struct Error_FixedPointTraits_not_specialized_for_this_case;
50 
51 template <class T1, class T2>
52 class FixedPointTraits
53 {
54 public:
55  typedef Error_FixedPointTraits_not_specialized_for_this_case PlusType;
56  typedef Error_FixedPointTraits_not_specialized_for_this_case MinusType;
57  typedef Error_FixedPointTraits_not_specialized_for_this_case MultipliesType;
58 // typedef Error_FixedPointTraits_not_specialized_for_this_case DividesType;
59 };
60 
61 // return type policy:
62 // * try to allocate enough bits to represent the biggest possible result
63 // * in case of add/subtract: if all bits of the internal int are used up,
64 // keep the representation
65 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
66 class FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >
67 {
68  enum { MaxIntBits = (IntBits1 < IntBits2) ? IntBits2 : IntBits1,
69  MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1,
70  PlusMinusIntBits = (MaxIntBits + 1 + MaxFracBits < 32) ?
71  MaxIntBits + 1 : MaxIntBits,
72  MultipliesFracBits = (IntBits1 + IntBits2 < 31)
73  ? (FracBits1 + FracBits2) > (31 - IntBits1 - IntBits2)
74  ? 31 - IntBits1 - IntBits2
75  : FracBits1 + FracBits2
76  : 0
77  };
78 public:
79  typedef FixedPoint<PlusMinusIntBits, MaxFracBits> PlusType;
80  typedef FixedPoint<PlusMinusIntBits, MaxFracBits> MinusType;
81  typedef FixedPoint<IntBits1 + IntBits2, MultipliesFracBits> MultipliesType;
82 // typedef FixedPoint<IntBits1 + FracBits2, FracBits1 + IntBits2> DividesType;
83 };
84 
85 template <unsigned IntBits, unsigned FracBits>
86 struct SquareRootTraits<FixedPoint<IntBits, FracBits> >
87 {
88  enum { SRTotalBits = (IntBits + FracBits + 1) / 2,
89  SRIntBits = (IntBits + 1) / 2,
90  SRFracBits = SRTotalBits - SRIntBits
91  };
92 public:
93  typedef FixedPoint<IntBits, FracBits> Type;
94  typedef FixedPoint<SRIntBits, SRFracBits> SquareRootResult;
95  typedef Type SquareRootArgument;
96 };
97 
98 
99 #ifndef DOXYGEN
100 
101 template <int N>
102 struct FixedPoint_overflow_error__More_than_31_bits_requested
103 : staticAssert::AssertBool<(N < 32)>
104 {};
105 
106 #endif /* DOXYGEN */
107 
108 
109 
110 template <bool Predicate>
111 struct FixedPoint_assignment_error__Target_object_has_too_few_integer_bits
112 : staticAssert::AssertBool<Predicate>
113 {};
114 
115 enum FixedPointNoShift { FPNoShift };
116 
117 namespace detail {
118 
119 template <bool MustRound>
120 struct FPAssignWithRound;
121 
122 template <>
123 struct FPAssignWithRound<false>
124 {
125  template <int N>
126  static inline int exec(int v) { return v << (-N); }
127 };
128 
129 template <>
130 struct FPAssignWithRound<true>
131 {
132  template <int N>
133  static inline int exec(int const v)
134  {
135  return (v + (1 << (N - 1))) >> (N);
136  }
137 };
138 
139 template <bool MustRound>
140 struct FPMulImplementation;
141 
142 template <>
143 struct FPMulImplementation<false>
144 {
145  template <int N>
146  static inline int exec(int l, int r) { return (l * r) << (-N); }
147 };
148 
149 template <>
150 struct FPMulImplementation<true>
151 {
152  template <int N>
153  static inline int exec(int l, int r)
154  {
155  // there is not enough space in the result
156  // => perform calculations that preserve as much accuracy as possible
157  enum { diffl = N / 2, diffr = N - diffl, maskl = (1 << diffl) - 1, maskr = (1 << diffr) - 1 };
158  int shiftl = l >> diffl;
159  int shiftr = r >> diffr;
160 
161  return shiftl * shiftr + (((l & maskl) * shiftr) >> diffl) +
162  (((r & maskr) * shiftl) >> diffr);
163  }
164 };
165 
166 } // namespace detail
167 
168 /********************************************************/
169 /* */
170 /* FixedPoint */
171 /* */
172 /********************************************************/
173 
174 /** Template for fixed point arithmetic.
175 
176  Fixed point arithmetic is used when computations with fractional accuracy
177  must be made at the highest speed possible (e.g. in the inner loop
178  of a volume rendering routine). The speed-up relative to floating
179  point arithmetic can be dramatic, especially when one can avoid
180  conversions between integer and floating point numbers (these are
181  very expensive because integer and floating point arithmetic
182  resides in different pipelines).
183 
184  The template wraps an <tt>int</tt> and uses <tt>IntBits</tt> to
185  represent the integral part of a number, and <tt>FractionalBits</tt>
186  for the fractional part, where <tt>IntBits + FractionalBits < 32</tt>.
187  (The 32rd bit is reserved because FixedPoint is a signed type).
188  These numbers will be automatically allocated in an intelligent way
189  in the result of an arithmetic operation. For example, when two
190  fixed point numbers are multiplied, the required number of integer
191  bits in the result is the sum of the number of integer bits of the
192  arguments, but only when so many bits are available. This is figured out
193  by means of FixedPointTraits, and a compile-time error is raised
194  when no suitable representation can be found. The idea is that the right
195  thing happens automatically as often as possible.
196 
197  <tt>FixedPoint</tt> implements the required interface of an
198  \ref AlgebraicRing and the required numeric and
199  promotion traits. In addition, it supports functions <tt>add</tt>,
200  <tt>sub</tt>, and <tt>mul</tt>, where a particular layout of the result can
201  be enforced.
202 
203  <tt>unsigned char, signed char, unsigned short, signed short, int</tt> can be
204  transformed into a FixedPoint with appropriate layout by means of the factory
205  function <tt>fixedPoint()</tt>.
206 
207  <b>See also:</b>
208  <ul>
209  <li> \ref FixedPointOperations
210  <li> \ref FixedPointTraits
211  </ul>
212 
213  <b>\#include</b> <vigra/fixedpoint.hxx><br>
214  Namespace: vigra
215 */
216 template <unsigned IntBits, unsigned FractionalBits>
217 class FixedPoint
218 {
219 public:
220  enum {
221  INT_BITS = IntBits,
222  FRACTIONAL_BITS = FractionalBits,
223  TOTAL_BITS = IntBits + FractionalBits,
224  MAX = (int)(((unsigned)1 << TOTAL_BITS) - 1),
225  ONE = 1 << FractionalBits,
226  ONE_HALF = ONE >> 1,
227  FRACTIONAL_MASK = ONE - 1,
228  INT_MASK = MAX ^ FRACTIONAL_MASK
229  };
230 
231  Int32 value;
232 
233  FixedPoint()
234  {
235  VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
236  }
237 
238  /** Construct from an int (fractional part will become zero).
239  */
240  explicit FixedPoint(int v)
241  : value(v << FractionalBits)
242  {
243  VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
244  }
245 
246  /** Construct from an int by a bitwise copy. This is normally only used internally.
247  */
248  FixedPoint(int v, FixedPointNoShift)
249  : value(v)
250  {
251  VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
252  }
253 
254  /** Construct from an double and round the fractional part to
255  <tt>FractionalBits</tt> accuracy. A PreconditionViolation exception is raised when
256  the integer part is too small to represent the number.
257  */
258  explicit FixedPoint(double rhs)
259  : value((int)round(rhs * ONE))
260  {
261  VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
262  vigra_precondition(abs(rhs * ONE) <= (double)MAX,
263  "FixedPoint(double rhs): Too few integer bits to convert rhs.");
264  }
265 
266  /** Copy constructor.
267  */
268  FixedPoint(const FixedPoint &other)
269  : value(other.value)
270  {}
271 
272  /** Construct from a FixedPoint with different layout. It rounds as appropriate and raises
273  a compile-time error when the target type has too few integer bits.
274  */
275  template <unsigned Int2, unsigned Frac2>
276  FixedPoint(const FixedPoint<Int2, Frac2> &other)
277  : value(detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<int(Frac2) - int(FractionalBits)>(other.value))
278  {
279  VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
280  VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
281  }
282 
283  /** Assignment from int. The fractional part will become zero.
284  A PreconditionViolation exception is raised when
285  the integer part is too small to represent the number.
286  */
287  FixedPoint &operator=(int rhs)
288  {
289  vigra_precondition(abs(rhs) < (1 << IntBits),
290  "FixedPoint::operator=(int rhs): Too few integer bits to represent rhs.");
291  value = rhs << FractionalBits;
292  return *this;
293  }
294 
295  /** Assignment form double. The fractional part is rounded, and a
296  PreconditionViolation exception is raised when
297  the integer part is too small to represent the number.
298  */
299  FixedPoint &operator=(double rhs)
300  {
301  vigra_precondition(abs(rhs) <= ((1 << IntBits) - 1),
302  "FixedPoint::operator=(double rhs): Too few integer bits to convert rhs.");
303  value = (int)round(rhs * ONE);
304  return *this;
305  }
306 
307  /** Copy assignment.
308  */
309  FixedPoint & operator=(const FixedPoint &other)
310  {
311  value = other.value;
312  return *this;
313  }
314 
315  /** Assignment from a FixedPoint with different layout. It rounds as appropriate and raises
316  a compile-time error when the target type has too few integer bits.
317  */
318  template <unsigned Int2, unsigned Frac2>
319  FixedPoint & operator=(const FixedPoint<Int2, Frac2> &other)
320  {
321  VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
322  value = detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<int(Frac2) - int(FractionalBits)>(other.value);
323  return *this;
324  }
325 
326  /** Negation.
327  */
328  FixedPoint operator-() const
329  {
330  return FixedPoint(-value, FPNoShift);
331  }
332 
333  /** Pre-increment.
334  */
335  FixedPoint & operator++()
336  {
337  value += ONE;
338  return *this;
339  }
340 
341  /** Post-increment.
342  */
343  FixedPoint operator++(int)
344  {
345  FixedPoint old(*this);
346  value += ONE;
347  return old;
348  }
349 
350  /** Pre-decrement.
351  */
352  FixedPoint & operator--()
353  {
354  value -= ONE;
355  return *this;
356  }
357 
358  /** Post-decrement.
359  */
360  FixedPoint operator--(int)
361  {
362  FixedPoint old(*this);
363  value -= ONE;
364  return old;
365  }
366 
367  /** Add-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
368  a compile-time error when the target type has too few integer bits.
369  */
370  template <unsigned Int2, unsigned Frac2>
371  FixedPoint & operator+=(const FixedPoint<Int2, Frac2> &other)
372  {
373  VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
374  value += detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
375  return *this;
376  }
377 
378  /** Subtract-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
379  a compile-time error when the target type has too few integer bits.
380  */
381  template <unsigned Int2, unsigned Frac2>
382  FixedPoint & operator-=(const FixedPoint<Int2, Frac2> &other)
383  {
384  VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
385  value -= detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
386  return *this;
387  }
388 
389  /** Multiply-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
390  a compile-time error when the target type has too few integer bits.
391  */
392  template <unsigned Int2, unsigned Frac2>
393  FixedPoint & operator*=(const FixedPoint<Int2, Frac2> &other)
394  {
395  VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
396  value = detail::FPMulImplementation<(Frac2 > 0)>::template exec<Frac2>(value, other.value);
397  return *this;
398  }
399 };
400 
401 #define VIGRA_FIXED_POINT_FACTORY(T, INTBITS) \
402  inline FixedPoint<INTBITS, 0> fixedPoint(T t) \
403  { \
404  return FixedPoint<INTBITS, 0>(t, FPNoShift); \
405  }
406 
407 VIGRA_FIXED_POINT_FACTORY(unsigned char, 8)
408 VIGRA_FIXED_POINT_FACTORY(signed char, 7)
409 VIGRA_FIXED_POINT_FACTORY(unsigned short, 16)
410 VIGRA_FIXED_POINT_FACTORY(signed short, 15)
411 VIGRA_FIXED_POINT_FACTORY(int, 31)
412 
413 #undef VIGRA_FIXED_POINT_FACTORY
414 
415 template <class T>
416 struct FixedPointCast;
417 
418 #define VIGRA_FIXED_POINT_CAST(type) \
419 template <> \
420 struct FixedPointCast<type> \
421 { \
422  template <unsigned IntBits, unsigned FracBits> \
423  static type cast(FixedPoint<IntBits, FracBits> v) \
424  { \
425  return round(v); \
426  } \
427 };
428 
429 VIGRA_FIXED_POINT_CAST(Int8)
430 VIGRA_FIXED_POINT_CAST(UInt8)
431 VIGRA_FIXED_POINT_CAST(Int16)
432 VIGRA_FIXED_POINT_CAST(UInt16)
433 VIGRA_FIXED_POINT_CAST(Int32)
434 VIGRA_FIXED_POINT_CAST(UInt32)
435 
436 #undef VIGRA_FIXED_POINT_CAST
437 
438 template <>
439 struct FixedPointCast<float>
440 {
441  template <unsigned IntBits, unsigned FracBits>
442  static float cast(FixedPoint<IntBits, FracBits> v)
443  {
444  return (float)v.value / FixedPoint<IntBits, FracBits>::ONE;
445  }
446 };
447 
448 template <>
449 struct FixedPointCast<double>
450 {
451  template <unsigned IntBits, unsigned FracBits>
452  static double cast(FixedPoint<IntBits, FracBits> v)
453  {
454  return (double)v.value / FixedPoint<IntBits, FracBits>::ONE;
455  }
456 };
457 
458 /********************************************************/
459 /* */
460 /* FixedPointOperations */
461 /* */
462 /********************************************************/
463 
464 /** \addtogroup FixedPointOperations Functions for FixedPoint
465 
466  \brief <b>\#include</b> <vigra/fixedpoint.hxx><br>
467 
468  These functions fulfill the requirements of an \ref AlgebraicRing.
469 
470  Namespace: vigra
471  <p>
472 
473  */
474 //@{
475 
476  /** Convert a FixedPoint to a built-in type.
477  If the target is integral, the value is rounded.<br>
478  Usage:
479  \code
480  FixedPoint<16,15> fp(...);
481 
482  double d = fixed_point_cast<double>(fp);
483  \endcode
484  */
485 template <class TARGET, unsigned IntBits, unsigned FracBits>
486 TARGET fixed_point_cast(FixedPoint<IntBits, FracBits> v)
487 {
488  return FixedPointCast<TARGET>::cast(v);
489 }
490 
491  /// equal
492 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
493 inline
494 bool operator==(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
495 {
496  enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
497  return (l.value << (MaxFracBits - FracBits1)) == (r.value << (MaxFracBits - FracBits2));
498 }
499 
500  /// not equal
501 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
502 inline
503 bool operator!=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
504 {
505  enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
506  return (l.value << (MaxFracBits - FracBits1)) != (r.value << (MaxFracBits - FracBits2));
507 }
508 
509  /// less than
510 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
511 inline
512 bool operator<(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
513 {
514  enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
515  return (l.value << (MaxFracBits - FracBits1)) < (r.value << (MaxFracBits - FracBits2));
516 }
517 
518  /// less or equal
519 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
520 inline
521 bool operator<=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
522 {
523  enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
524  return (l.value << (MaxFracBits - FracBits1)) <= (r.value << (MaxFracBits - FracBits2));
525 }
526 
527  /// greater
528 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
529 inline
530 bool operator>(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
531 {
532  enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
533  return (l.value << (MaxFracBits - FracBits1)) > (r.value << (MaxFracBits - FracBits2));
534 }
535 
536  /// greater or equal
537 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
538 inline
539 bool operator>=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
540 {
541  enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
542  return (l.value << (MaxFracBits - FracBits1)) >= (r.value << (MaxFracBits - FracBits2));
543 }
544 
545  /// addition with automatic determination of the appropriate result type.
546 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
547 inline
548 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType
549 operator+(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
550 {
551  enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
552  return typename
553  FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
554  PlusType((l.value << (MaxFracBits - FracBits1)) + (r.value << (MaxFracBits - FracBits2)), FPNoShift);
555 }
556 
557  /// addition with enforced result type.
558 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
559  unsigned IntBits3, unsigned FracBits3>
560 inline void
561 add(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
562  FixedPoint<IntBits3, FracBits3> & result)
563 {
564  result = l + r;
565 }
566 
567  /// subtraction with automatic determination of the appropriate result type.
568 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
569 inline
570 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::MinusType
571 operator-(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
572 {
573  enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
574  return typename
575  FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
576  MinusType((l.value << (MaxFracBits - FracBits1)) - (r.value << (MaxFracBits - FracBits2)), FPNoShift);
577 }
578 
579  /// subtraction with enforced result type.
580 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
581  unsigned IntBits3, unsigned FracBits3>
582 inline void
583 sub(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
584  FixedPoint<IntBits3, FracBits3> & result)
585 {
586  result = l - r;
587 }
588 
589  /// multiplication with automatic determination of the appropriate result type.
590 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
591 inline
592 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::MultipliesType
593 operator*(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
594 {
595  typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
596  MultipliesType res;
597  mul(l, r, res);
598  return res;
599 }
600 
601  /// multiplication with enforced result type.
602 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
603  unsigned IntBits3, unsigned FracBits3>
604 inline void
605 mul(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
606  FixedPoint<IntBits3, FracBits3> & result)
607 {
608  VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits1 + IntBits2 <= IntBits3)>));
609  enum { diff = FracBits1 + FracBits2 - FracBits3 };
610  result.value = detail::FPMulImplementation<(diff > 0)>::template exec<diff>(l.value, r.value);
611 }
612 
613  /// square root.
614 template <unsigned IntBits, unsigned FracBits>
615 inline typename SquareRootTraits<FixedPoint<IntBits, FracBits> >::SquareRootResult
616 sqrt(FixedPoint<IntBits, FracBits> v)
617 {
618  return typename SquareRootTraits<FixedPoint<IntBits, FracBits> >::SquareRootResult(sqrti(v.value), FPNoShift);
619 }
620 
621  /// absolute value.
622 template <unsigned IntBits, unsigned FracBits>
623 inline FixedPoint<IntBits, FracBits>
624 abs(FixedPoint<IntBits, FracBits> v)
625 {
626  return FixedPoint<IntBits, FracBits>(abs(v.value), FPNoShift);
627 }
628 
629  /// squared norm (same as v*v).
630 template <unsigned IntBits, unsigned FracBits>
631 inline
632 typename FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
633 squaredNorm(FixedPoint<IntBits, FracBits> v)
634 {
635  return v*v;
636 }
637 
638  /// norm (same as abs).
639 template <unsigned IntBits, unsigned FracBits>
640 inline
641 FixedPoint<IntBits, FracBits>
642 norm(FixedPoint<IntBits, FracBits> const & v)
643 {
644  return abs(v);
645 }
646 
647  /// fractional part.
648 template <unsigned IntBits, unsigned FracBits>
649 inline FixedPoint<0, FracBits>
650 frac(FixedPoint<IntBits, FracBits> v)
651 {
652  return FixedPoint<0, FracBits>(v.value & FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK, FPNoShift);
653 }
654 
655  /// dual fractional part: <tt>1 - frac(v)</tt>.
656 template <unsigned IntBits, unsigned FracBits>
657 inline FixedPoint<0, FracBits>
658 dual_frac(FixedPoint<IntBits, FracBits> v)
659 {
660  return FixedPoint<0, FracBits>(FixedPoint<0, FracBits>::ONE -
661  (v.value & FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK), FPNoShift);
662 }
663 
664  /// rounding down.
665 template <unsigned IntBits, unsigned FracBits>
666 inline int
667 floor(FixedPoint<IntBits, FracBits> v)
668 {
669  return(v.value >> FracBits);
670 }
671 
672  /// rounding up.
673 template <unsigned IntBits, unsigned FracBits>
674 inline int
675 ceil(FixedPoint<IntBits, FracBits> v)
676 {
677  return((v.value + FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK) >> FracBits);
678 }
679 
680  /// rounding to the nearest integer.
681 template <unsigned IntBits, unsigned FracBits>
682 inline int
683 round(FixedPoint<IntBits, FracBits> v)
684 {
685  return((v.value + FixedPoint<IntBits, FracBits>::ONE_HALF) >> FracBits);
686 }
687 
688 //@}
689 
690 /********************************************************/
691 /* */
692 /* FixedPoint-Traits */
693 /* */
694 /********************************************************/
695 
696 /** \page FixedPointTraits Numeric and Promote Traits of FixedPoint
697 
698  The numeric and promote traits for FixedPoint follow
699  the general specifications for \ref NumericPromotionTraits and
700  \ref AlgebraicRing. They are implemented in terms of the traits of the basic types by
701  partial template specialization:
702 
703  \code
704 
705  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
706  class FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >
707  {
708  typedef FixedPoint<PlusMinusIntBits, MaxFracBits> PlusType;
709  typedef FixedPoint<PlusMinusIntBits, MaxFracBits> MinusType;
710  typedef FixedPoint<IntBits1 + IntBits2, FracBits1 + FracBits2> MultipliesType;
711  };
712 
713  template <unsigned IntBits, unsigned FracBits>
714  struct NumericTraits<FixedPoint<IntBits, FracBits> >
715  {
716  typedef FixedPoint<IntBits, FracBits> Type;
717  // Promote undefined because it depends on the layout, use FixedPointTraits
718  // RealPromote in AlgebraicRing -- multiplication with double is not supported.
719  // ComplexPromote in AlgebraicRing -- multiplication with double is not supported.
720  typedef Type ValueType;
721 
722  typedef VigraFalseType isIntegral;
723  typedef VigraTrueType isScalar;
724  typedef VigraTrueType isSigned;
725  typedef VigraTrueType isOrdered;
726  typedef VigraFalseType isComplex;
727 
728  ... // etc.
729  };
730 
731  template <unsigned IntBits, unsigned FracBits>
732  struct SquareRootTraits<FixedPoint<IntBits, FracBits> >
733  {
734  typedef FixedPoint<IntBits, FracBits> Type;
735  typedef FixedPoint<SRIntBits, SRFracBits> SquareRootResult;
736  typedef Type SquareRootArgument;
737  };
738 
739  template <unsigned IntBits, unsigned FracBits>
740  struct NormTraits<FixedPoint<IntBits, FracBits> >
741  {
742  typedef FixedPoint<IntBits, FracBits> Type;
743  typedef typename
744  FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
745  SquaredNormType;
746  typedef Type NormType;
747  };
748 
749  template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
750  struct PromoteTraits<FixedPoint<IntBits1, FracBits1>,
751  FixedPoint<IntBits2, FracBits2> >
752  {
753  typedef typename
754  FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType
755  Promote;
756  };
757  \endcode
758 
759  <b>\#include</b> <vigra/fixedpoint.hxx><br>
760  Namespace: vigra
761 
762 */
763 template <unsigned IntBits, unsigned FracBits>
764 struct NumericTraits<FixedPoint<IntBits, FracBits> >
765 {
766  typedef FixedPoint<IntBits, FracBits> Type;
767  //typedef FixedPoint<IntBits, FracBits> Promote;
768  //typedef FixedPoint<IntBits, FracBits> RealPromote;
769  //typedef std::complex<RealPromote> ComplexPromote;
770  typedef Type ValueType;
771 
772  typedef VigraFalseType isIntegral;
773  typedef VigraTrueType isScalar;
774  typedef VigraTrueType isSigned;
775  typedef VigraTrueType isOrdered;
776  typedef VigraFalseType isComplex;
777 
778  static Type zero() { return Type(0, FPNoShift); }
779  static Type one() { return Type(Type::ONE, FPNoShift); }
780  static Type nonZero() { return one(); }
781  static Type epsilon() { return Type(1, FPNoShift); }
782  static Type smallestPositive() { return Type(1, FPNoShift); }
783  static Type max() { return Type( Type::MAX, FPNoShift); }
784  static Type min() { return -max(); }
785 };
786 
787 template <unsigned IntBits, unsigned FracBits>
788 struct NormTraits<FixedPoint<IntBits, FracBits> >
789 {
790  typedef FixedPoint<IntBits, FracBits> Type;
791  typedef typename
792  FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
793  SquaredNormType;
794  typedef Type NormType;
795 };
796 
797 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
798 struct PromoteTraits<FixedPoint<IntBits1, FracBits1>,
799  FixedPoint<IntBits2, FracBits2> >
800 {
801  typedef typename
802  FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType
803  Promote;
804 };
805 
806 /***********************************************************************************/
807 
808 enum FPOverflowHandling { FPOverflowIgnore, FPOverflowSaturate, FPOverflowError };
809 
810 template <int IntBits, FPOverflowHandling OverflowHandling = FPOverflowIgnore>
811 class FixedPoint16;
812 
813 /********************************************************/
814 /* */
815 /* FixedPoint16-Traits */
816 /* */
817 /********************************************************/
818 
819 /** \page FixedPoint16Traits Numeric and Promote Traits of FixedPoint16
820 
821  The numeric and promote traits for FixedPoint16 follow
822  the general specifications for \ref NumericPromotionTraits and
823  \ref AlgebraicRing. They are implemented in terms of the traits of the basic types by
824  partial template specialization:
825 
826  \code
827  template <int IntBits, FPOverflowHandling OverflowHandling>
828  struct NumericTraits<FixedPoint16<IntBits, OverflowHandling> >
829  {
830  typedef FixedPoint16<IntBits, OverflowHandling> Type;
831  typedef Type Promote;
832  // RealPromote undefined -- multiplication with double is not supported.
833  // ComplexPromote undefined -- multiplication with double is not supported.
834  typedef Type ValueType;
835 
836  typedef VigraFalseType isIntegral;
837  typedef VigraTrueType isScalar;
838  typedef VigraTrueType isSigned;
839  typedef VigraTrueType isOrdered;
840  typedef VigraFalseType isComplex;
841 
842  ... // etc.
843  };
844 
845  template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
846  struct PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>,
847  FixedPoint16<IntBits2, OverflowHandling> >
848  {
849  typedef FixedPoint16<MetaMax<IntBits1, IntBits2>::value, OverflowHandling> Promote;
850  ... // etc.
851  };
852 
853  template <int IntBits, FPOverflowHandling OverflowHandling>
854  struct NormTraits<FixedPoint16<IntBits, OverflowHandling> >
855  {
856  typedef FixedPoint16<IntBits, OverflowHandling> Type;
857  typedef typename PromoteTraits<Type, Type>::Promote SquaredNormType;
858  typedef Type NormType;
859  };
860 
861  template <int IntBits, FPOverflowHandling OverflowHandling>
862  struct SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >
863  {
864  typedef FixedPoint16<IntBits, OverflowHandling> Type;
865  typedef FixedPoint16<(IntBits + 1) / 2, OverflowHandling> SquareRootResult;
866  typedef Type SquareRootArgument;
867  };
868  \endcode
869 
870  <b>\#include</b> <vigra/fixedpoint.hxx><br>
871  Namespace: vigra
872 
873 */
874 template <int IntBits, FPOverflowHandling OverflowHandling>
875 struct NumericTraits<FixedPoint16<IntBits, OverflowHandling> >
876 {
877  typedef FixedPoint16<IntBits, OverflowHandling> Type;
878  typedef Type Promote;
879  // RealPromote undefined -- multiplication with double is not supported.
880  // ComplexPromote undefined -- multiplication with double is not supported.
881  typedef Type ValueType;
882 
883  typedef VigraFalseType isIntegral;
884  typedef VigraTrueType isScalar;
885  typedef VigraTrueType isSigned;
886  typedef VigraTrueType isOrdered;
887  typedef VigraFalseType isComplex;
888 
889  static Type zero() { return Type(0, FPNoShift); }
890  static Type one() { return Type(Type::ONE, FPNoShift); }
891  static Type nonZero() { return one(); }
892  static Type epsilon() { return Type(1, FPNoShift); }
893  static Type smallestPositive() { return Type(1, FPNoShift); }
894  static Type max() { return Type( Type::MAX, FPNoShift); }
895  static Type min() { return Type( Type::MIN, FPNoShift); }
896 
897  static Promote toPromote(Type v) { return v; }
898  static Type fromPromote(Promote v) { return v; };
899 };
900 
901 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
902 struct PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>,
903  FixedPoint16<IntBits2, OverflowHandling> >
904 {
905  typedef FixedPoint16<MetaMax<IntBits1, IntBits2>::value, OverflowHandling> Promote;
906  static Promote toPromote(FixedPoint16<IntBits1, OverflowHandling> v) { return Promote(v); }
907  static Promote toPromote(FixedPoint16<IntBits2, OverflowHandling> v) { return Promote(v); }
908 };
909 
910 template <int IntBits, FPOverflowHandling OverflowHandling>
911 struct PromoteTraits<FixedPoint16<IntBits, OverflowHandling>,
912  FixedPoint16<IntBits, OverflowHandling> >
913 {
914  typedef FixedPoint16<IntBits, OverflowHandling> Promote;
915  static Promote toPromote(FixedPoint16<IntBits, OverflowHandling> v) { return v; }
916 };
917 
918 template <int IntBits, FPOverflowHandling OverflowHandling>
919 struct NormTraits<FixedPoint16<IntBits, OverflowHandling> >
920 {
921  typedef FixedPoint16<IntBits, OverflowHandling> Type;
922  typedef typename PromoteTraits<Type, Type>::Promote SquaredNormType;
923  typedef Type NormType;
924 };
925 
926 template <int IntBits, FPOverflowHandling OverflowHandling>
927 struct SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >
928 {
929  typedef FixedPoint16<IntBits, OverflowHandling> Type;
930  typedef FixedPoint16<(IntBits + 1) / 2, OverflowHandling> SquareRootResult;
931  typedef Type SquareRootArgument;
932 };
933 
934 #ifndef DOXYGEN
935 
936 template <bool Compatible>
937 struct FixedPoint_error__Right_shift_operator_has_unsupported_semantics
938 : staticAssert::AssertBool<Compatible>
939 {};
940 
941 #endif /* DOXYGEN */
942 
943 template <bool Predicate>
944 struct FixedPoint16_assignment_error__Target_object_has_too_few_integer_bits
945 : staticAssert::AssertBool<Predicate>
946 {};
947 
948 namespace detail {
949 
950 template<int BeforeIntBits, int AfterIntBits,
951  bool Round = false,
952  bool RightShift = (AfterIntBits >= BeforeIntBits)>
953 struct FP16Align;
954 
955 template<int BeforeIntBits>
956 struct FP16Align<BeforeIntBits, BeforeIntBits, true, true>
957 {
958  static inline Int32 exec(Int32 v)
959  {
960  return v;
961  }
962 };
963 
964 template<int BeforeIntBits>
965 struct FP16Align<BeforeIntBits, BeforeIntBits, false, true>
966 {
967  static inline Int32 exec(Int32 v)
968  {
969  return v;
970  }
971 };
972 
973 template<int BeforeIntBits, int AfterIntBits>
974 struct FP16Align<BeforeIntBits, AfterIntBits, false, true>
975 {
976  static inline Int32 exec(Int32 v)
977  {
978  VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
979  return v >> (AfterIntBits - BeforeIntBits);
980  }
981 };
982 
983 template<int BeforeIntBits, int AfterIntBits>
984 struct FP16Align<BeforeIntBits, AfterIntBits, true, true>
985 {
986  enum { ONE_HALF = 1 << (AfterIntBits - BeforeIntBits - 1) };
987  static inline Int32 exec(Int32 v)
988  {
989  VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
990  return (v + ONE_HALF) >> (AfterIntBits - BeforeIntBits);
991  }
992 };
993 
994 template<int BeforeIntBits, int AfterIntBits, bool Round>
995 struct FP16Align<BeforeIntBits, AfterIntBits, Round, false>
996 {
997  static inline Int32 exec(Int32 v)
998  {
999  return v << (BeforeIntBits - AfterIntBits);
1000  }
1001 };
1002 
1003 template <FPOverflowHandling OverflowHandling = FPOverflowIgnore>
1004 struct FP16OverflowHandling
1005 {
1006  static inline Int32 exec(Int32 v)
1007  {
1008  return v;
1009  }
1010 
1011  static inline Int32 exec(UInt32 v)
1012  {
1013  return v;
1014  }
1015 };
1016 
1017 template <>
1018 struct FP16OverflowHandling<FPOverflowSaturate>
1019 {
1020  static inline Int32 exec(Int32 v)
1021  {
1022  if(v >= 1 << 15)
1023  return (1 << 15) - 1;
1024  if(v < -(1 << 15))
1025  return -(1 << 15);
1026  return v;
1027  }
1028  static inline Int32 exec(UInt32 v)
1029  {
1030  if(v >= 1 << 15)
1031  return (1 << 15) - 1;
1032  return v;
1033  }
1034 };
1035 
1036 template <>
1037 struct FP16OverflowHandling<FPOverflowError>
1038 {
1039  static inline Int32 exec(Int32 v)
1040  {
1041  vigra_precondition(v < (1 << 15) && v >= -(1 << 15),
1042  "FixedPoint16: Operation overflows.");
1043  return v;
1044  }
1045  static inline Int32 exec(UInt32 v)
1046  {
1047  vigra_precondition(v < (1 << 15),
1048  "FixedPoint16: Operation overflows.");
1049  return v;
1050  }
1051 };
1052 
1053 
1054 template <int IntBits1, int IntBits2, int IntBitsOut,
1055  FPOverflowHandling OverflowHandling >
1056 struct FP16AddImpl
1057 {
1058  enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
1059  static inline Int32 exec(Int32 t1, Int32 t2)
1060  {
1061  return FP16OverflowHandling<OverflowHandling>::exec(
1062  FP16Align<MinIntBits, IntBitsOut, /*Round*/ true>::exec(
1063  FP16Align<IntBits1, MinIntBits, /*Round*/ false>::exec(t1) +
1064  FP16Align<IntBits2, MinIntBits, /*Round*/ false>::exec(t2)));
1065  }
1066 };
1067 
1068 template <int IntBits1, int IntBits2, int IntBitsOut,
1069  FPOverflowHandling OverflowHandling >
1070 struct FP16SubImpl
1071 {
1072  enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
1073  static inline Int32 exec(Int32 t1, Int32 t2)
1074  {
1075  return FP16OverflowHandling<OverflowHandling>::exec(
1076  FP16Align<MinIntBits, IntBitsOut, /*Round*/ true>::exec(
1077  FP16Align<IntBits1, MinIntBits, /*Round*/ false>::exec(t1) -
1078  FP16Align<IntBits2, MinIntBits, /*Round*/ false>::exec(t2)));
1079  }
1080 };
1081 
1082 template <int IntBits1, int IntBits2, int IntBitsOut,
1083  FPOverflowHandling OverflowHandling >
1084 struct FP16MulImpl
1085 {
1086  static inline Int32 exec(Int32 t1, Int32 t2)
1087  {
1088  return FP16OverflowHandling<OverflowHandling>::exec(
1089  FP16Align<IntBits1+IntBits2, IntBitsOut+15, /*Round*/ true>::exec(t1*t2));
1090  }
1091 };
1092 
1093 template <int IntBits1, int IntBits2, int IntBitsOut,
1094  FPOverflowHandling OverflowHandling >
1095 struct FP16DivImpl
1096 {
1097  static inline Int32 exec(Int32 t1, Int32 t2)
1098  {
1099  if(t2 == 0)
1100  return (t1 >= 0)
1101  ? (1 << 15) - 1
1102  : -(1 << 15);
1103  return FP16OverflowHandling<OverflowHandling>::exec(
1104  FP16Align<IntBits1-IntBits2, IntBitsOut+1, /*Round*/ true>::exec((t1<<16)/t2));
1105  }
1106 };
1107 
1108 } // namespace detail
1109 
1110 /********************************************************/
1111 /* */
1112 /* FixedPoint16 */
1113 /* */
1114 /********************************************************/
1115 
1116 template <class TARGET, int IntBits, FPOverflowHandling OverflowHandling>
1117 TARGET fixed_point_cast(FixedPoint16<IntBits, OverflowHandling> v);
1118 
1119 /** Template for 16-bit signed fixed point arithmetic.
1120 
1121  Fixed point arithmetic is used when computations with fractional accuracy
1122  must be made at the highest speed possible (e.g. in the inner loop
1123  of a volume rendering routine). The speed-up relative to floating
1124  point arithmetic can be dramatic, especially when one can avoid
1125  conversions between integer and floating point numbers (these are
1126  very expensive because integer and floating point arithmetic
1127  resides in different pipelines).
1128 
1129  The template wraps an <tt>Int16</tt> and uses <tt>IntBits</tt> to
1130  represent the integral part of a number, and <tt>15 - IntBits</tt>
1131  for the fractional part. The 16th bit is reserved because FixedPoint16
1132  is a signed type. Results of expressions with mixed types will preserve
1133  larger number of <tt>IntBits</tt> of the results, in order to minimize
1134  the possibility for overflow. Nonetheless, overflow can occur, and the
1135  template parameter <tt>OverflowHandling</tt> determines how this will be
1136  handled:
1137 
1138  <DL>
1139  <DT>FPOverflowIgnore<DD> (default) Ignore overflow, i.e. use the usual modulo behavior of the
1140  built-in integer types.
1141 
1142  <DT>FPOverflowSaturate<DD> Use the largest or smallest representable number (depending on sign)
1143  in case of overflow.
1144 
1145  <DT>FPOverflowError<DD> Throw <tt>PreconditionViolation</tt> upon overflow. This is useful for
1146  debugging.
1147  </DL>
1148 
1149  The implementation relies on Int32-arithmetic and requires that the right-shift operator
1150  preserves signedness. Although not enforced by the C++ standard, this is implemented
1151  by most of today's processors. This property is checked by a
1152  VIGRA_STATIC_ASSERT(FixedPoint_error__Right_shift_operator_has_unsupported_semantics).
1153 
1154  <tt>FixedPoint16</tt> implements the required interface of an
1155  \ref AlgebraicRing and the required numeric and
1156  promotion traits. In addition, it supports functions <tt>add</tt>,
1157  <tt>sub</tt>, <tt>mul</tt>, and <tt>div</tt>, where a particular layout
1158  of the result can be enforced.
1159 
1160  Built-in numeric types can be converted into <tt>FixedPoint16</tt> by the
1161  appropriate constructors, and from <tt>FixedPoint16</tt> by means of
1162  <tt>fixed_point_cast&lt;TargetType&gt;(fixedPoint)</tt>.
1163 
1164  <b>See also:</b>
1165  <ul>
1166  <li> \ref FixedPoint16Operations
1167  <li> \ref FixedPoint16Traits
1168  </ul>
1169 
1170  <b>\#include</b> <vigra/fixedpoint.hxx><br>
1171  Namespace: vigra
1172 */
1173 template <int IntBits, FPOverflowHandling OverflowHandling>
1174 class FixedPoint16
1175 {
1176 public:
1177  static const Int32 TOTAL_BITS = 15; // bit 16 is sign
1178  static const Int32 INT_BITS = IntBits;
1179  static const Int32 FRACTIONAL_BITS = TOTAL_BITS - INT_BITS;
1180  static const Int32 MAX = (Int32)((1u << TOTAL_BITS) - 1);
1181  static const Int32 MIN = -(Int32)(1u << TOTAL_BITS);
1182  static const Int32 ONE = 1 << FRACTIONAL_BITS;
1183  static const Int32 ONE_HALF = ONE >> 1;
1184  static const Int32 FRACTIONAL_MASK = (1u << FRACTIONAL_BITS) - 1;
1185  static const Int32 INT_MASK = 0xffffffffu ^ FRACTIONAL_MASK;
1186 
1187  static const FixedPoint16 zero, pi, pi_2, mpi_2;
1188 
1189  Int16 value;
1190 
1191  FixedPoint16()
1192  : value(0)
1193  {
1194  VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
1195  }
1196 
1197  /** Construct from an int (fractional part will become zero).
1198  Possible overflow is handled according to the target type's <tt>OverflowHandling</tt>.
1199  */
1200  explicit FixedPoint16(Int32 v)
1201  : value(detail::FP16OverflowHandling<OverflowHandling>::exec(v << FRACTIONAL_BITS))
1202  {
1203  VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
1204  }
1205 
1206  /** Construct from an int by a bitwise copy. This is normally only used internally.
1207  */
1208  FixedPoint16(Int32 v, FixedPointNoShift)
1209  : value(detail::FP16OverflowHandling<OverflowHandling>::exec(v))
1210  {
1211  VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
1212  }
1213 
1214  /** Construct from a double and round the fractional part to
1215  <tt>FRACTIONAL_BITS</tt> accuracy. Possible overflow is handled according
1216  to the target type's <tt>OverflowHandling</tt>.
1217  */
1218  explicit FixedPoint16(double rhs)
1219  : value(detail::FP16OverflowHandling<OverflowHandling>::exec(roundi(rhs * ONE)))
1220  {
1221  VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
1222  }
1223 
1224  /** Copy constructor.
1225  */
1226  FixedPoint16(const FixedPoint16 &other)
1227  : value(other.value)
1228  {
1229  VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
1230  }
1231 
1232  /** Construct from a FixedPoint16 with different layout. It rounds as appropriate and
1233  handles possible overflow according to the target type's <tt>OverflowHandling</tt>.
1234  */
1235  template <int IntBits2, FPOverflowHandling OverflowHandling2>
1236  FixedPoint16(const FixedPoint16<IntBits2, OverflowHandling2> &other)
1237  : value(detail::FP16OverflowHandling<OverflowHandling>::exec(
1238  detail::FP16Align<IntBits2, IntBits, /*Round*/true>::exec(other.value)))
1239  {
1240  VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
1241  }
1242 
1243  /** Assignment from int. The fractional part will become zero.
1244  Possible overflow is handled according to the target type's <tt>OverflowHandling</tt>.
1245  */
1246  FixedPoint16 &operator=(Int32 rhs)
1247  {
1248  value = detail::FP16OverflowHandling<OverflowHandling>::exec(rhs << FRACTIONAL_BITS);
1249  return *this;
1250  }
1251 
1252  /** Assignment form double. The fractional part is rounded, and possible overflow is
1253  handled according to the target type's <tt>OverflowHandling</tt>.
1254  */
1255  FixedPoint16 &operator=(double rhs)
1256  {
1257  value = detail::FP16OverflowHandling<OverflowHandling>::exec(roundi(rhs * ONE));
1258  return *this;
1259  }
1260 
1261  /** Copy assignment.
1262  */
1263  FixedPoint16 & operator=(const FixedPoint16 &other)
1264  {
1265  value = other.value;
1266  return *this;
1267  }
1268 
1269  /** Assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
1270  handled according to the target type's <tt>OverflowHandling</tt>.
1271  */
1272  template <int IntBits2>
1273  FixedPoint16 & operator=(const FixedPoint16<IntBits2, OverflowHandling> &other)
1274  {
1275  value = detail::FP16OverflowHandling<OverflowHandling>::exec(
1276  detail::FP16Align<IntBits2, IntBits, /*Round*/true>::exec(other.value));
1277  return *this;
1278  }
1279 
1280  /** Conversion to float
1281  */
1282  operator float() const
1283  {
1284  return fixed_point_cast<float>(*this);
1285  }
1286 
1287  /** Conversion to double
1288  */
1289  operator double() const
1290  {
1291  return fixed_point_cast<double>(*this);
1292  }
1293 
1294  /** Unary plus.
1295  */
1296  FixedPoint16 operator+() const
1297  {
1298  return *this;
1299  }
1300 
1301  /** Negation.
1302  */
1303  FixedPoint16 operator-() const
1304  {
1305  return FixedPoint16(-value, FPNoShift);
1306  }
1307 
1308  /** Pre-increment.
1309  */
1310  FixedPoint16 & operator++()
1311  {
1312  value = detail::FP16OverflowHandling<OverflowHandling>::exec(value+ONE);
1313  return *this;
1314  }
1315 
1316  /** Post-increment.
1317  */
1318  FixedPoint16 operator++(int)
1319  {
1320  FixedPoint16 old(*this);
1321  ++(*this);
1322  return old;
1323  }
1324 
1325  /** Pre-decrement.
1326  */
1327  FixedPoint16 & operator--()
1328  {
1329  value = detail::FP16OverflowHandling<OverflowHandling>::exec(value-ONE);
1330  return *this;
1331  }
1332 
1333  /** Post-decrement.
1334  */
1335  FixedPoint16 operator--(int)
1336  {
1337  FixedPoint16 old(*this);
1338  --(*this);
1339  return old;
1340  }
1341 
1342  /** Add-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
1343  handled according to the target type's <tt>OverflowHandling</tt>.
1344  */
1345  template <int IntBits2>
1346  FixedPoint16 & operator+=(const FixedPoint16<IntBits2, OverflowHandling> &other)
1347  {
1348  value = detail::FP16AddImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
1349  return *this;
1350  }
1351 
1352  /** Subtract-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
1353  handled according to the target type's <tt>OverflowHandling</tt>.
1354  */
1355  template <int IntBits2>
1356  FixedPoint16 & operator-=(const FixedPoint16<IntBits2, OverflowHandling> &other)
1357  {
1358  value = detail::FP16SubImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
1359  return *this;
1360  }
1361 
1362  /** Multiply-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
1363  handled according to the target type's <tt>OverflowHandling</tt>.
1364  */
1365  template <int IntBits2>
1366  FixedPoint16 & operator*=(const FixedPoint16<IntBits2, OverflowHandling> &other)
1367  {
1368  value = detail::FP16MulImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
1369  return *this;
1370  }
1371 
1372  /** Divide-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is
1373  handled according to the target type's <tt>OverflowHandling</tt>.
1374  */
1375  template <int IntBits2>
1376  FixedPoint16 & operator/=(const FixedPoint16<IntBits2, OverflowHandling> &other)
1377  {
1378  value = detail::FP16DivImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
1379  return *this;
1380  }
1381 };
1382 
1383 template <int IntBits, FPOverflowHandling OverflowHandling>
1384 const FixedPoint16<IntBits, OverflowHandling> FixedPoint16<IntBits, OverflowHandling>::zero(0);
1385 
1386 template <int IntBits, FPOverflowHandling OverflowHandling>
1387 const FixedPoint16<IntBits, OverflowHandling> FixedPoint16<IntBits, OverflowHandling>::pi(M_PI);
1388 
1389 template <int IntBits, FPOverflowHandling OverflowHandling>
1390 const FixedPoint16<IntBits, OverflowHandling> FixedPoint16<IntBits, OverflowHandling>::pi_2(0.5 * M_PI);
1391 
1392 template <int IntBits, FPOverflowHandling OverflowHandling>
1393 const FixedPoint16<IntBits, OverflowHandling> FixedPoint16<IntBits, OverflowHandling>::mpi_2(-0.5 * M_PI);
1394 
1395 namespace detail {
1396 
1397 template <class T>
1398 struct FixedPoint16Cast;
1399 
1400 #define VIGRA_FIXED_POINT_CAST(type) \
1401 template <> \
1402 struct FixedPoint16Cast<type> \
1403 { \
1404  template <int IntBits, FPOverflowHandling OverflowHandling> \
1405  static type cast(FixedPoint16<IntBits, OverflowHandling> v) \
1406  { \
1407  return round(v); \
1408  } \
1409 };
1410 
1411 VIGRA_FIXED_POINT_CAST(Int8)
1412 VIGRA_FIXED_POINT_CAST(UInt8)
1413 VIGRA_FIXED_POINT_CAST(Int16)
1414 VIGRA_FIXED_POINT_CAST(UInt16)
1415 VIGRA_FIXED_POINT_CAST(Int32)
1416 VIGRA_FIXED_POINT_CAST(UInt32)
1417 VIGRA_FIXED_POINT_CAST(Int64)
1418 VIGRA_FIXED_POINT_CAST(UInt64)
1419 
1420 #undef VIGRA_FIXED_POINT_CAST
1421 
1422 template <>
1423 struct FixedPoint16Cast<float>
1424 {
1425  template <int IntBits, FPOverflowHandling OverflowHandling>
1426  static float cast(FixedPoint16<IntBits, OverflowHandling> v)
1427  {
1428  return (float)v.value / FixedPoint16<IntBits, OverflowHandling>::ONE;
1429  }
1430 };
1431 
1432 template <>
1433 struct FixedPoint16Cast<double>
1434 {
1435  template <int IntBits, FPOverflowHandling OverflowHandling>
1436  static double cast(FixedPoint16<IntBits, OverflowHandling> v)
1437  {
1438  return (double)v.value / FixedPoint16<IntBits, OverflowHandling>::ONE;
1439  }
1440 };
1441 
1442 } // namespace detail
1443 
1444 /********************************************************/
1445 /* */
1446 /* FixedPoint16Operations */
1447 /* */
1448 /********************************************************/
1449 
1450 /** \addtogroup FixedPoint16Operations Functions for FixedPoint16
1451 
1452  \brief <b>\#include</b> <vigra/fixedpoint.hxx><br>
1453 
1454  These functions fulfill the requirements of an \ref AlgebraicRing.
1455 
1456  Namespace: vigra
1457  <p>
1458 
1459  */
1460 //@{
1461 
1462  /** Convert a FixedPoint16 to a built-in type.
1463  If the target is integral, the value is rounded.<br>
1464  Usage:
1465  \code
1466  FixedPoint16<16,15> fp(...);
1467 
1468  double d = fixed_point_cast<double>(fp);
1469  \endcode
1470  */
1471 template <class TARGET, int IntBits, FPOverflowHandling OverflowHandling>
1472 TARGET fixed_point_cast(FixedPoint16<IntBits, OverflowHandling> v)
1473 {
1474  return detail::FixedPoint16Cast<TARGET>::cast(v);
1475 }
1476 
1477 
1478  /// equal
1479 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1480 inline
1481 bool operator==(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1482 {
1483  enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
1484  return (l.value << (IntBits1 - MinIntBits)) == (r.value << (IntBits2 - MinIntBits));
1485 }
1486 
1487  /// not equal
1488 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1489 inline
1490 bool operator!=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1491 {
1492  enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
1493  return (l.value << (IntBits1 - MinIntBits)) != (r.value << (IntBits2 - MinIntBits));
1494 }
1495 
1496  /// less than
1497 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1498 inline
1499 bool operator<(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1500 {
1501  enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
1502  return (l.value << (IntBits1 - MinIntBits)) < (r.value << (IntBits2 - MinIntBits));
1503 }
1504 
1505  /// less or equal
1506 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1507 inline
1508 bool operator<=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1509 {
1510  enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
1511  return (l.value << (IntBits1 - MinIntBits)) <= (r.value << (IntBits2 - MinIntBits));
1512 }
1513 
1514  /// greater
1515 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1516 inline
1517 bool operator>(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1518 {
1519  enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
1520  return (l.value << (IntBits1 - MinIntBits)) > (r.value << (IntBits2 - MinIntBits));
1521 }
1522 
1523  /// greater or equal
1524 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1525 inline
1526 bool operator>=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1527 {
1528  enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
1529  return (l.value << (IntBits1 - MinIntBits)) >= (r.value << (IntBits2 - MinIntBits));
1530 }
1531 
1532  /// addition with automatic determination of the appropriate result type.
1533 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1534 inline
1535 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
1536 operator+(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1537 {
1538  typedef typename
1539  PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
1540  Result;
1541  return Result(detail::FP16AddImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
1542 }
1543 
1544  /// addition with enforced result type.
1545 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
1546 inline
1547 FixedPoint16<IntBits3, OverflowHandling> &
1548 add(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
1549  FixedPoint16<IntBits3, OverflowHandling> & result)
1550 {
1551  result.value = detail::FP16AddImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
1552  return result;
1553 }
1554 
1555  /// subtraction with automatic determination of the appropriate result type.
1556 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1557 inline
1558 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
1559 operator-(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1560 {
1561  typedef typename
1562  PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
1563  Result;
1564  return Result(detail::FP16SubImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
1565 }
1566 
1567  /// subtraction with enforced result type.
1568 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
1569 inline FixedPoint16<IntBits3, OverflowHandling> &
1570 sub(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
1571  FixedPoint16<IntBits3, OverflowHandling> & result)
1572 {
1573  result.value = detail::FP16SubImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
1574  return result;
1575 }
1576 
1577  /// multiplication with automatic determination of the appropriate result type.
1578 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1579 inline
1580 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
1581 operator*(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1582 {
1583  typedef typename
1584  PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
1585  Result;
1586  return Result(detail::FP16MulImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
1587 }
1588 
1589  /// multiplication with enforced result type.
1590 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
1591 inline
1592 FixedPoint16<IntBits3, OverflowHandling> &
1593 mul(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
1594  FixedPoint16<IntBits3, OverflowHandling> & result)
1595 {
1596  result.value = detail::FP16MulImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
1597  return result;
1598 }
1599 
1600  /// division with automatic determination of the appropriate result type.
1601 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
1602 inline
1603 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
1604 operator/(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
1605 {
1606  typedef typename
1607  PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
1608  Result;
1609  return Result(detail::FP16DivImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
1610 }
1611 
1612  /// division with enforced result type.
1613 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
1614 inline
1615 FixedPoint16<IntBits3, OverflowHandling> &
1616 div(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
1617  FixedPoint16<IntBits3, OverflowHandling> & result)
1618 {
1619  result.value = detail::FP16DivImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
1620  return result;
1621 }
1622 
1623  /// square root.
1624 template <int IntBits, FPOverflowHandling OverflowHandling>
1625 inline typename SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >::SquareRootResult
1626 sqrt(FixedPoint16<IntBits, OverflowHandling> v)
1627 {
1628  typedef typename SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >::SquareRootResult Result;
1629  enum { Shift = 15 + IntBits - 2*Result::INT_BITS };
1630  return Result(sqrti(v.value << Shift), FPNoShift);
1631 }
1632 
1633 #ifndef VIGRA_NO_HYPOT
1634  using ::hypot;
1635 #endif
1636 
1637  /// Length of hypotenuse.
1638 template <int IntBits, FPOverflowHandling OverflowHandling>
1639 inline FixedPoint16<IntBits, OverflowHandling>
1640 hypot(FixedPoint16<IntBits, OverflowHandling> v1, FixedPoint16<IntBits, OverflowHandling> v2)
1641 {
1642  UInt32 l = abs(v1.value), r = abs(v2.value);
1643  // sq(l) + sq(r) <= 2**31, so overflow handling after sqrti is sufficient
1644  return FixedPoint16<IntBits, OverflowHandling>(
1645  detail::FP16OverflowHandling<OverflowHandling>::exec(sqrti(sq(l) + sq(r))),
1646  FPNoShift);
1647 }
1648 
1649 using std::atan2;
1650 
1651  /// Arctangent. Accuracy better than 1/3 degree (9 significant bits).
1652 template <int IntBits, FPOverflowHandling OverflowHandling>
1653 FixedPoint16<2, OverflowHandling>
1654 atan2(FixedPoint16<IntBits, OverflowHandling> y, FixedPoint16<IntBits, OverflowHandling> x)
1655 {
1656  enum { ResIntBits = 2 };
1657  typedef FixedPoint16<ResIntBits, OverflowHandling> FP;
1658  static const Int32 Pi_4 = 25736, // = roundi(0.25 * M_PI * (1 << 15)), at 15 frac bits
1659  Pi3_4 = 77208, // = roundi(0.75 * M_PI * (1 << 15)),
1660  c1 = 6497, // = roundi(0.19826763260224867 * (1 << 15)),
1661  c2 = -1047730238; // = roundi(-0.9757748231899761 * (1 << 30));
1662 
1663  // coefficients c1 and c2 minimize
1664  //
1665  // NIntegrate[(c1 r^3 + c2 r + Pi/4 - a)^4 /. r -> (Cos[a] - Sin[a])/(Cos[a] + Sin[a]), {a, 0, Pi/2}]
1666  //
1667  // Thanks to Jim Shima, http://www.dspguru.com/comp.dsp/tricks/alg/fxdatan2.htm
1668 
1669  if(x.value == 0)
1670  return (y.value > 0)
1671  ? FP::pi_2
1672  : (y.value < 0)
1673  ? FP::mpi_2
1674  : FP::zero;
1675 
1676  Int32 abs_y = abs(y.value);
1677  Int32 r, angle;
1678  if(x.value > 0)
1679  {
1680  if(y.value == 0)
1681  return FP::zero;
1682  r = ((x.value - abs_y) << 15) / (x.value + abs_y); // 15 frac bits
1683  angle = Pi_4;
1684  }
1685  else
1686  {
1687  if(y.value == 0)
1688  return FP::pi;
1689  r = ((x.value + abs_y) << 15) / (abs_y - x.value); // 15 frac bits
1690  angle = Pi3_4;
1691  }
1692 
1693  angle += r*((c2 + c1 * (sq(r) >> 15)) >> 15) >> 15;
1694 
1695  return (y.value > 0)
1696  ? FP(detail::FP16Align<0, ResIntBits, true>::exec( angle), FPNoShift)
1697  : FP(detail::FP16Align<0, ResIntBits, true>::exec(-angle), FPNoShift);
1698 }
1699 
1700  /// absolute value.
1701 template <int IntBits, FPOverflowHandling OverflowHandling>
1702 inline FixedPoint16<IntBits, OverflowHandling>
1703 abs(FixedPoint16<IntBits, OverflowHandling> v)
1704 {
1705  return FixedPoint16<IntBits, OverflowHandling>(abs(v.value), FPNoShift);
1706 }
1707 
1708  /// squared norm (same as v*v).
1709 template <int IntBits, FPOverflowHandling OverflowHandling>
1710 inline
1711 typename NormTraits<FixedPoint16<IntBits, OverflowHandling> >::SquaredNormType
1712 squaredNorm(FixedPoint16<IntBits, OverflowHandling> v)
1713 {
1714  return v*v;
1715 }
1716 
1717  /// norm (same as abs).
1718 template <int IntBits, FPOverflowHandling OverflowHandling>
1719 inline
1720 typename NormTraits<FixedPoint16<IntBits, OverflowHandling> >::NormType
1721 norm(FixedPoint16<IntBits, OverflowHandling> const & v)
1722 {
1723  return abs(v);
1724 }
1725 
1726  /// fractional part. (difference between v and its floor)
1727 template <int IntBits, FPOverflowHandling OverflowHandling>
1728 inline FixedPoint16<IntBits, OverflowHandling>
1729 frac(FixedPoint16<IntBits, OverflowHandling> v)
1730 {
1731  return FixedPoint16<IntBits, OverflowHandling>(
1732  v.value - (v.value & FixedPoint16<IntBits, OverflowHandling>::INT_MASK),
1733  FPNoShift);
1734 }
1735 
1736  /// dual fractional part. (1 - frac(v))
1737 template <int IntBits, FPOverflowHandling OverflowHandling>
1738 inline FixedPoint16<IntBits, OverflowHandling>
1739 dual_frac(FixedPoint16<IntBits, OverflowHandling> v)
1740 {
1741  return FixedPoint16<IntBits, OverflowHandling>(
1742  FixedPoint16<IntBits, OverflowHandling>::ONE - v.value + (v.value & FixedPoint16<IntBits, OverflowHandling>::INT_MASK),
1743  FPNoShift);
1744 }
1745 
1746  /// rounding down.
1747 template <int IntBits, FPOverflowHandling OverflowHandling>
1748 inline Int32
1749 floor(FixedPoint16<IntBits, OverflowHandling> v)
1750 {
1751  return(v.value >> FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
1752 }
1753 
1754  /// rounding up.
1755 template <int IntBits, FPOverflowHandling OverflowHandling>
1756 inline Int32
1757 ceil(FixedPoint16<IntBits, OverflowHandling> v)
1758 {
1759  return((v.value + FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_MASK) >>
1760  FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
1761 }
1762 
1763  /// rounding to the nearest integer.
1764 template <int IntBits, FPOverflowHandling OverflowHandling>
1765 inline Int32
1766 round(FixedPoint16<IntBits, OverflowHandling> v)
1767 {
1768  return((v.value + FixedPoint16<IntBits, OverflowHandling>::ONE_HALF) >>
1769  FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
1770 }
1771 
1772  /// rounding to the nearest integer.
1773 template <int IntBits, FPOverflowHandling OverflowHandling>
1774 inline Int32
1775 roundi(FixedPoint16<IntBits, OverflowHandling> v)
1776 {
1777  return round(v);
1778 }
1779 
1780 //@}
1781 
1782 } // namespace vigra
1783 
1784 namespace std {
1785 
1786 template <int IntBits, vigra::FPOverflowHandling OverflowHandling>
1787 ostream & operator<<(ostream & s, vigra::FixedPoint16<IntBits, OverflowHandling> v)
1788 {
1789  s << vigra::fixed_point_cast<float>(v);
1790  return s;
1791 }
1792 
1793 } // namespace std
1794 
1795 #endif // VIGRA_FIXEDPOINT_HXX
vigra
Definition: accessor.hxx:43
vigra::FixedPoint
Definition: fixedpoint.hxx:47

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.10.0