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

vigra/fixedpoint.hxx
00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 2004-2005 by Ullrich Koethe                  */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 #ifndef VIGRA_FIXEDPOINT_HXX
00037 #define VIGRA_FIXEDPOINT_HXX
00038 
00039 #include "mathutil.hxx"
00040 #include "static_assert.hxx"
00041 #include "error.hxx"
00042 #include "numerictraits.hxx"
00043 
00044 namespace vigra {
00045 
00046 template <unsigned IntBits, unsigned FractionalBits>
00047 class FixedPoint;
00048 
00049 struct Error_FixedPointTraits_not_specialized_for_this_case;
00050 
00051 template <class T1, class T2>
00052 class FixedPointTraits
00053 {
00054 public:
00055     typedef Error_FixedPointTraits_not_specialized_for_this_case PlusType;
00056     typedef Error_FixedPointTraits_not_specialized_for_this_case MinusType;
00057     typedef Error_FixedPointTraits_not_specialized_for_this_case MultipliesType;
00058 //    typedef Error_FixedPointTraits_not_specialized_for_this_case DividesType;
00059 };
00060 
00061 // return type policy: 
00062 //     * try to allocate enough bits to represent the biggest possible result
00063 //     * in case of add/subtract: if all bits of the internal int are used up, 
00064 //                                keep the representation
00065 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00066 class FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >
00067 {
00068     enum { MaxIntBits  = (IntBits1 < IntBits2) ? IntBits2 : IntBits1,
00069            MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1,
00070            PlusMinusIntBits = (MaxIntBits + 1 + MaxFracBits < 32) ?
00071                                MaxIntBits + 1 : MaxIntBits,
00072            MultipliesFracBits = (IntBits1 + IntBits2 < 31) 
00073                                     ? (FracBits1 + FracBits2) > (31 - IntBits1 - IntBits2) 
00074                                             ? 31 - IntBits1 - IntBits2
00075                                             : FracBits1 + FracBits2
00076                                     : 0
00077          };
00078 public:
00079     typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               PlusType;
00080     typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               MinusType;
00081     typedef FixedPoint<IntBits1 + IntBits2, MultipliesFracBits>  MultipliesType;
00082 //    typedef FixedPoint<IntBits1 + FracBits2, FracBits1 + IntBits2>  DividesType;
00083 };
00084 
00085 template <unsigned IntBits, unsigned FracBits>
00086 struct SquareRootTraits<FixedPoint<IntBits, FracBits> >
00087 {
00088     enum { SRTotalBits = (IntBits + FracBits + 1) / 2,
00089            SRIntBits   = (IntBits + 1) / 2,
00090            SRFracBits  = SRTotalBits - SRIntBits
00091          };
00092 public:
00093     typedef FixedPoint<IntBits, FracBits>      Type;
00094     typedef FixedPoint<SRIntBits, SRFracBits>  SquareRootResult;
00095     typedef Type                               SquareRootArgument;
00096 };
00097 
00098 
00099 #ifndef DOXYGEN
00100 
00101 template <int N>
00102 struct FixedPoint_overflow_error__More_than_31_bits_requested
00103 : staticAssert::AssertBool<(N < 32)>
00104 {};
00105 
00106 #endif /* DOXYGEN */
00107 
00108 
00109 
00110 template <bool Predicate>
00111 struct FixedPoint_assignment_error__Target_object_has_too_few_integer_bits
00112 : staticAssert::AssertBool<Predicate>
00113 {};
00114 
00115 enum FixedPointNoShift { FPNoShift };
00116 
00117 namespace detail {
00118 
00119 template <bool MustRound>
00120 struct FPAssignWithRound;
00121 
00122 template <>
00123 struct FPAssignWithRound<false>
00124 {
00125     template <int N>
00126     static inline int exec(int v) { return v << (-N); }
00127 };
00128 
00129 template <>
00130 struct FPAssignWithRound<true>
00131 {
00132     template <int N>
00133     static inline int exec(int const v)
00134     {
00135         return (v + (1 << (N - 1))) >> (N);
00136     }
00137 };
00138 
00139 template <bool MustRound>
00140 struct FPMulImplementation;
00141 
00142 template <>
00143 struct FPMulImplementation<false>
00144 {
00145     template <int N>
00146     static inline int exec(int l, int r) { return (l * r) << (-N); }
00147 };
00148 
00149 template <>
00150 struct FPMulImplementation<true>
00151 {
00152     template <int N>
00153     static inline int exec(int l, int r)
00154     {
00155         // there is not enough space in the result
00156         // => perform calculations that preserve as much accuracy as possible
00157         enum { diffl = N / 2, diffr = N  - diffl, maskl = (1 << diffl) - 1, maskr = (1 << diffr) - 1 };
00158         int shiftl = l >> diffl;
00159         int shiftr = r >> diffr;
00160 
00161         return shiftl * shiftr + (((l & maskl) * shiftr) >> diffl) +
00162                                  (((r & maskr) * shiftl) >> diffr);
00163     }
00164 };
00165 
00166 } // namespace detail
00167 
00168 /********************************************************/
00169 /*                                                      */
00170 /*                      FixedPoint                      */
00171 /*                                                      */
00172 /********************************************************/
00173 
00174 /** Template for fixed point arithmetic.
00175 
00176     Fixed point arithmetic is used when computations with fractional accuracy
00177     must be made at the highest speed possible (e.g. in the inner loop
00178     of a volume rendering routine). The speed-up relative to floating
00179     point arithmetic can be dramatic, especially when one can avoid
00180     conversions between integer anfloating point numbers (these are 
00181     very expensive because integer and floating point arithmetic
00182     resides in different pipelines). 
00183     
00184     The template wraps an <tt>int</tt> and uses <tt>IntBits</tt> to
00185     represent the integral part of a number, and <tt>FractionalBits</tt>
00186     for the fractional part, where <tt>IntBits + FractionalBits < 32</tt>.
00187     (The 32rd bit is reserved because FixedPoint is a signed type).
00188     These numbers will be automatically allocated in an intelligent way
00189     in the result of an arithmetic operation. For example, when two 
00190     fixed point numbers are multiplied, the required number of integer
00191     bits in the result is the sum of the number of integer bits of the
00192     arguments, but only when so many bits are avaiable. This is figured out
00193     by means of FixedPointTraits, and a compile-time error is raised
00194     when no suitable representation can be found. The idea is that the right
00195     thing happens automatically as often as possible.
00196 
00197     <tt>FixedPoint</tt> implements the required interface of an
00198     \ref AlgebraicRing and the required numeric and
00199     promotion traits. In addition, it supports functions <tt>add</tt>, 
00200     <tt>sub</tt>, and <tt>mul</tt>, where a particular layout of the result can
00201     be enforced. 
00202     
00203     <tt>unsigned char, signed char, unsigned short, signed short, int</tt> can be
00204     transformed into a FixedPoint with appropriate layout by means of the factory
00205     function <tt>fixedPoint()</tt>.
00206 
00207     <b>See also:</b>
00208     <ul>
00209     <li> \ref FixedPointOperations
00210     <li> \ref FixedPointTraits
00211     </ul>
00212 
00213     <b>\#include</b> <<a href="fixedpoint_8hxx-source.html">vigra/fixedpoint.hxx</a>><br>
00214     Namespace: vigra
00215 */
00216 template <unsigned IntBits, unsigned FractionalBits>
00217 class FixedPoint
00218 {
00219 public:
00220     enum {
00221         INT_BITS        = IntBits,
00222         FRACTIONAL_BITS = FractionalBits,
00223         TOTAL_BITS      = IntBits + FractionalBits,
00224         MAX             = (int)(((unsigned)1 << TOTAL_BITS) - 1),
00225         ONE             = 1 << FractionalBits,
00226         ONE_HALF        = ONE >> 1,
00227         FRACTIONAL_MASK = ONE - 1,
00228         INT_MASK        = MAX ^ FRACTIONAL_MASK
00229     };
00230 
00231     Int32 value;
00232 
00233     FixedPoint()
00234     {
00235         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
00236     }
00237 
00238         /** Construct from an int (fractional part will become zero).
00239         */
00240     explicit FixedPoint(int v)
00241     : value(v << FractionalBits)
00242     {
00243         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
00244     }
00245 
00246         /** Construct from an int by a bitwise copy. This is normally only used internally.
00247         */
00248     FixedPoint(int v, FixedPointNoShift)
00249     : value(v)
00250     {
00251         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
00252     }
00253 
00254         /** Construct from an double and round the fractional part to 
00255             <tt>FractionalBits</tt> accuracy. A PreconditionViolation exception is raised when
00256             the integer part is too small to represent the number.
00257         */
00258     explicit FixedPoint(double rhs)
00259     : value((int)round(rhs * ONE))
00260     {
00261         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
00262         vigra_precondition(abs(rhs * ONE) <= (double)MAX,
00263             "FixedPoint(double rhs): Too few integer bits to convert rhs.");
00264     }
00265 
00266         /** Copy constructor.
00267         */
00268     FixedPoint(const FixedPoint &other)
00269     : value(other.value)
00270     {}
00271 
00272         /** Construct from a FixedPoint with different layout. It rounds as appropriate and raises
00273             a compile-time error when the target type has too few integer bits.
00274         */
00275     template <unsigned Int2, unsigned Frac2>
00276     FixedPoint(const FixedPoint<Int2, Frac2> &other)
00277     : value(detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value))
00278     {
00279         VIGRA_STATIC_ASSERT((FixedPoint_overflow_error__More_than_31_bits_requested<(IntBits + FractionalBits)>));
00280         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
00281     }
00282 
00283         /** Assignment from int. The fractional part will become zero.  
00284             A PreconditionViolation exception is raised when
00285             the integer part is too small to represent the number.
00286         */
00287     FixedPoint &operator=(int rhs)
00288     {
00289         vigra_precondition(abs(rhs) < (1 << IntBits),
00290             "FixedPoint::operator=(int rhs): Too few integer bits to represent rhs.");
00291         value = rhs << FractionalBits;
00292         return *this;
00293     }
00294 
00295         /** Assignment form double. The fractional part is rounded, and a 
00296             PreconditionViolation exception is raised when
00297             the integer part is too small to represent the number.
00298         */
00299     FixedPoint &operator=(double rhs)
00300     {
00301         vigra_precondition(abs(rhs) <= ((1 << IntBits) - 1),
00302             "FixedPoint::operator=(double rhs): Too few integer bits to convert rhs.");
00303         value = (int)round(rhs * ONE);
00304         return *this;
00305     }
00306 
00307         /** Copy assignment.
00308         */
00309     FixedPoint & operator=(const FixedPoint &other)
00310     {
00311         value = other.value;
00312         return *this;
00313     }
00314 
00315         /** Assignment from a FixedPoint with different layout. It rounds as appropriate and raises
00316             a compile-time error when the target type has too few integer bits.
00317         */
00318     template <unsigned Int2, unsigned Frac2>
00319     FixedPoint & operator=(const FixedPoint<Int2, Frac2> &other)
00320     {
00321         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
00322         value = detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
00323         return *this;
00324     }
00325 
00326         /** Negation.
00327         */
00328     FixedPoint operator-() const
00329     {
00330         return FixedPoint(-value, FPNoShift);
00331     }
00332 
00333         /** Pre-increment.
00334         */
00335     FixedPoint & operator++()
00336     {
00337         value += ONE;
00338         return *this;
00339     }
00340 
00341         /** Post-increment.
00342         */
00343     FixedPoint operator++(int)
00344     {
00345         FixedPoint old(*this);
00346         value += ONE;
00347         return old;
00348     }
00349 
00350         /** Pre-decrement.
00351         */
00352     FixedPoint & operator--()
00353     {
00354         value -= ONE;
00355         return *this;
00356     }
00357 
00358         /** Post-decrement.
00359         */
00360     FixedPoint operator--(int)
00361     {
00362         FixedPoint old(*this);
00363         value -= ONE;
00364         return old;
00365     }
00366 
00367         /** Add-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
00368             a compile-time error when the target type has too few integer bits.
00369         */
00370     template <unsigned Int2, unsigned Frac2>
00371     FixedPoint & operator+=(const FixedPoint<Int2, Frac2> &other)
00372     {
00373         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
00374         value += detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
00375         return *this;
00376     }
00377 
00378         /** Subtract-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
00379             a compile-time error when the target type has too few integer bits.
00380         */
00381     template <unsigned Int2, unsigned Frac2>
00382     FixedPoint & operator-=(const FixedPoint<Int2, Frac2> &other)
00383     {
00384         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
00385         value -= detail::FPAssignWithRound<(Frac2 > FractionalBits)>::template exec<Frac2 - FractionalBits>(other.value);
00386         return *this;
00387     }
00388     
00389         /** Multiply-assignment from a FixedPoint with different layout. It rounds as appropriate and raises
00390             a compile-time error when the target type has too few integer bits.
00391         */
00392     template <unsigned Int2, unsigned Frac2>
00393     FixedPoint & operator*=(const FixedPoint<Int2, Frac2> &other)
00394     {
00395         VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits >= Int2)>));
00396         value = detail::FPMulImplementation<(Frac2 > 0)>::template exec<Frac2>(value, other.value);
00397         return *this;
00398     }
00399 };
00400 
00401 #define VIGRA_FIXED_POINT_FACTORY(T, INTBITS) \
00402     inline FixedPoint<INTBITS, 0> fixedPoint(T t) \
00403     { \
00404         return FixedPoint<INTBITS, 0>(t, FPNoShift); \
00405     }
00406 
00407 VIGRA_FIXED_POINT_FACTORY(unsigned char, 8)
00408 VIGRA_FIXED_POINT_FACTORY(signed char, 7)
00409 VIGRA_FIXED_POINT_FACTORY(unsigned short, 16)
00410 VIGRA_FIXED_POINT_FACTORY(signed short, 15)
00411 VIGRA_FIXED_POINT_FACTORY(int, 31)
00412 
00413 #undef VIGRA_FIXED_POINT_FACTORY
00414 
00415 template <class T>
00416 struct FixedPointCast;
00417 
00418 #define VIGRA_FIXED_POINT_CAST(type) \
00419 template <> \
00420 struct FixedPointCast<type> \
00421 { \
00422     template <unsigned IntBits, unsigned FracBits> \
00423     static type cast(FixedPoint<IntBits, FracBits> v) \
00424     { \
00425         return round(v); \
00426     } \
00427 };
00428 
00429 VIGRA_FIXED_POINT_CAST(Int8)
00430 VIGRA_FIXED_POINT_CAST(UInt8)
00431 VIGRA_FIXED_POINT_CAST(Int16)
00432 VIGRA_FIXED_POINT_CAST(UInt16)
00433 VIGRA_FIXED_POINT_CAST(Int32)
00434 VIGRA_FIXED_POINT_CAST(UInt32)
00435 
00436 #undef VIGRA_FIXED_POINT_CAST
00437 
00438 template <>
00439 struct FixedPointCast<float>
00440 {
00441     template <unsigned IntBits, unsigned FracBits>
00442     static float cast(FixedPoint<IntBits, FracBits> v)
00443     {
00444         return (float)v.value / FixedPoint<IntBits, FracBits>::ONE;
00445     }
00446 };
00447 
00448 template <>
00449 struct FixedPointCast<double>
00450 {
00451     template <unsigned IntBits, unsigned FracBits>
00452     static double cast(FixedPoint<IntBits, FracBits> v)
00453     {
00454         return (double)v.value / FixedPoint<IntBits, FracBits>::ONE;
00455     }
00456 };
00457 
00458 /********************************************************/
00459 /*                                                      */
00460 /*                 FixedPointOperations                 */
00461 /*                                                      */
00462 /********************************************************/
00463 
00464 /** \addtogroup FixedPointOperations Functions for FixedPoint
00465 
00466     \brief     <b>\#include</b> <<a href="fixedpoint_8hxx-source.html">vigra/fixedpoint.hxx</a>><br>
00467 
00468     These functions fulfill the requirements of an \ref AlgebraicRing.
00469 
00470     Namespace: vigra
00471     <p>
00472 
00473  */
00474 //@{
00475 
00476     /** Convert a FixedPoint to a built-in type.
00477         If the target is integral, the value is rounded.<br>
00478         Usage:
00479         \code
00480         FixedPoint<16,15> fp(...);
00481         
00482         double d = fixed_point_cast<double>(fp);
00483         \endcode
00484     */
00485 template <class TARGET, unsigned IntBits, unsigned FracBits>
00486 TARGET fixed_point_cast(FixedPoint<IntBits, FracBits> v)
00487 {
00488     return FixedPointCast<TARGET>::cast(v);
00489 }
00490 
00491     /// equal
00492 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00493 inline
00494 bool operator==(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00495 {
00496     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
00497     return (l.value << (MaxFracBits - FracBits1)) == (r.value << (MaxFracBits - FracBits2));
00498 }
00499 
00500     /// not equal
00501 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00502 inline
00503 bool operator!=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00504 {
00505     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
00506     return (l.value << (MaxFracBits - FracBits1)) != (r.value << (MaxFracBits - FracBits2));
00507 }
00508 
00509     /// less than
00510 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00511 inline
00512 bool operator<(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00513 {
00514     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
00515     return (l.value << (MaxFracBits - FracBits1)) < (r.value << (MaxFracBits - FracBits2));
00516 }
00517 
00518     /// less or equal
00519 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00520 inline
00521 bool operator<=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00522 {
00523     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
00524     return (l.value << (MaxFracBits - FracBits1)) <= (r.value << (MaxFracBits - FracBits2));
00525 }
00526 
00527     /// greater
00528 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00529 inline
00530 bool operator>(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00531 {
00532     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
00533     return (l.value << (MaxFracBits - FracBits1)) > (r.value << (MaxFracBits - FracBits2));
00534 }
00535 
00536     /// greater or equal
00537 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00538 inline
00539 bool operator>=(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00540 {
00541     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
00542     return (l.value << (MaxFracBits - FracBits1)) >= (r.value << (MaxFracBits - FracBits2));
00543 }
00544 
00545     /// addition with automatic determination of the appropriate result type.
00546 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00547 inline
00548 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType
00549 operator+(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00550 {
00551     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
00552     return typename
00553         FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
00554         PlusType((l.value << (MaxFracBits - FracBits1)) + (r.value << (MaxFracBits - FracBits2)), FPNoShift);
00555 }
00556 
00557     /// addition with enforced result type.
00558 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
00559           unsigned IntBits3, unsigned FracBits3>
00560 inline void
00561 add(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
00562     FixedPoint<IntBits3, FracBits3> & result)
00563 {
00564     result = l + r;
00565 }
00566 
00567     /// subtraction with automatic determination of the appropriate result type.
00568 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00569 inline
00570 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::MinusType
00571 operator-(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00572 {
00573     enum { MaxFracBits = (FracBits1 < FracBits2) ? FracBits2 : FracBits1 };
00574     return typename
00575         FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
00576         MinusType((l.value << (MaxFracBits - FracBits1)) - (r.value << (MaxFracBits - FracBits2)), FPNoShift);
00577 }
00578 
00579     /// subtraction with enforced result type.
00580 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
00581           unsigned IntBits3, unsigned FracBits3>
00582 inline void
00583 sub(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
00584     FixedPoint<IntBits3, FracBits3> & result)
00585 {
00586     result = l - r;
00587 }
00588 
00589     /// multiplication with automatic determination of the appropriate result type.
00590 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00591 inline
00592 typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::MultipliesType
00593 operator*(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r)
00594 {
00595     typename FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::
00596         MultipliesType res;
00597     mul(l, r, res);
00598     return res;
00599 }
00600 
00601     /// multiplication with enforced result type.
00602 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2,
00603           unsigned IntBits3, unsigned FracBits3>
00604 inline void
00605 mul(FixedPoint<IntBits1, FracBits1> l, FixedPoint<IntBits2, FracBits2> r,
00606     FixedPoint<IntBits3, FracBits3> & result)
00607 {
00608     VIGRA_STATIC_ASSERT((FixedPoint_assignment_error__Target_object_has_too_few_integer_bits<(IntBits1 + IntBits2 <= IntBits3)>));
00609     enum { diff = FracBits1 + FracBits2 - FracBits3 };
00610     result.value = detail::FPMulImplementation<(diff > 0)>::template exec<diff>(l.value, r.value);
00611 }
00612 
00613     /// square root.
00614 template <unsigned IntBits, unsigned FracBits>
00615 inline typename SquareRootTraits<FixedPoint<IntBits, FracBits> >::SquareRootResult
00616 sqrt(FixedPoint<IntBits, FracBits> v)
00617 {
00618     return typename SquareRootTraits<FixedPoint<IntBits, FracBits> >::SquareRootResult(sqrti(v.value), FPNoShift);
00619 }
00620 
00621     /// absolute value.
00622 template <unsigned IntBits, unsigned FracBits>
00623 inline FixedPoint<IntBits, FracBits>
00624 abs(FixedPoint<IntBits, FracBits> v)
00625 {
00626     return FixedPoint<IntBits, FracBits>(abs(v.value), FPNoShift);
00627 }
00628 
00629     /// squared norm (same as v*v).
00630 template <unsigned IntBits, unsigned FracBits>
00631 inline
00632 typename FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
00633 squaredNorm(FixedPoint<IntBits, FracBits> v)
00634 {
00635     return v*v;
00636 }
00637 
00638     /// norm (same as abs).
00639 template <unsigned IntBits, unsigned FracBits>
00640 inline
00641 FixedPoint<IntBits, FracBits>
00642 norm(FixedPoint<IntBits, FracBits> const & v)
00643 {
00644     return abs(v);
00645 }
00646 
00647     /// fractional part.
00648 template <unsigned IntBits, unsigned FracBits>
00649 inline FixedPoint<0, FracBits>
00650 frac(FixedPoint<IntBits, FracBits> v)
00651 {
00652     return FixedPoint<0, FracBits>(v.value & FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK, FPNoShift);
00653 }
00654 
00655     /// dual fractional part: <tt>1 - frac(v)</tt>.
00656 template <unsigned IntBits, unsigned FracBits>
00657 inline FixedPoint<0, FracBits>
00658 dual_frac(FixedPoint<IntBits, FracBits> v)
00659 {
00660     return FixedPoint<0, FracBits>(FixedPoint<0, FracBits>::ONE - 
00661                                    (v.value & FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK), FPNoShift);
00662 }
00663 
00664     /// rounding down.
00665 template <unsigned IntBits, unsigned FracBits>
00666 inline int
00667 floor(FixedPoint<IntBits, FracBits> v)
00668 {
00669     return(v.value >> FracBits);
00670 }
00671 
00672     /// rounding up.
00673 template <unsigned IntBits, unsigned FracBits>
00674 inline int
00675 ceil(FixedPoint<IntBits, FracBits> v)
00676 {
00677     return((v.value + FixedPoint<IntBits, FracBits>::FRACTIONAL_MASK) >> FracBits);
00678 }
00679 
00680     /// rounding to the nearest integer.
00681 template <unsigned IntBits, unsigned FracBits>
00682 inline int
00683 round(FixedPoint<IntBits, FracBits> v)
00684 {
00685     return((v.value + FixedPoint<IntBits, FracBits>::ONE_HALF) >> FracBits);
00686 }
00687 
00688 //@}
00689 
00690 /********************************************************/
00691 /*                                                      */
00692 /*                     FixedPoint-Traits                */
00693 /*                                                      */
00694 /********************************************************/
00695 
00696 /** \page FixedPointTraits Numeric and Promote Traits of FixedPoint
00697 
00698     The numeric and promote traits for FixedPoint follow
00699     the general specifications for \ref NumericPromotionTraits and
00700     \ref AlgebraicRing. They are implemented in terms of the traits of the basic types by
00701     partial template specialization:
00702 
00703     \code
00704 
00705     template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00706     class FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >
00707     {
00708         typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               PlusType;
00709         typedef FixedPoint<PlusMinusIntBits, MaxFracBits>               MinusType;
00710         typedef FixedPoint<IntBits1 + IntBits2, FracBits1 + FracBits2>  MultipliesType;
00711     };
00712 
00713     template <unsigned IntBits, unsigned FracBits>
00714     struct NumericTraits<FixedPoint<IntBits, FracBits> >
00715     {
00716         typedef FixedPoint<IntBits, FracBits> Type;
00717             // Promote undefined because it depends on the layout, use FixedPointTraits
00718             // RealPromote in AlgebraicRing -- multiplication with double is not supported.
00719             // ComplexPromote in AlgebraicRing -- multiplication with double is not supported.
00720         typedef Type ValueType;
00721 
00722         typedef VigraFalseType isIntegral;
00723         typedef VigraTrueType  isScalar;
00724         typedef VigraTrueType  isSigned;
00725         typedef VigraTrueType  isOrdered;
00726         typedef VigraFalseType isComplex;
00727 
00728         ... // etc.
00729     };
00730 
00731     template <unsigned IntBits, unsigned FracBits>
00732     struct SquareRootTraits<FixedPoint<IntBits, FracBits> >
00733     {
00734         typedef FixedPoint<IntBits, FracBits>      Type;
00735         typedef FixedPoint<SRIntBits, SRFracBits>  SquareRootResult;
00736         typedef Type                               SquareRootArgument;
00737     };
00738     
00739     template <unsigned IntBits, unsigned FracBits>
00740     struct NormTraits<FixedPoint<IntBits, FracBits> >
00741     {
00742         typedef FixedPoint<IntBits, FracBits>         Type;
00743         typedef typename 
00744             FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
00745                                                       SquaredNormType;
00746         typedef Type                                  NormType;
00747     };
00748 
00749     template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00750     struct PromoteTraits<FixedPoint<IntBits1, FracBits1>,
00751                          FixedPoint<IntBits2, FracBits2> >
00752     {
00753         typedef typename 
00754             FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType 
00755             Promote;
00756     };
00757     \endcode
00758 
00759     <b>\#include</b> <<a href="fixedpoint_8hxx-source.html">vigra/fixedpoint.hxx</a>><br>
00760     Namespace: vigra
00761 
00762 */
00763 template <unsigned IntBits, unsigned FracBits>
00764 struct NumericTraits<FixedPoint<IntBits, FracBits> >
00765 {
00766     typedef FixedPoint<IntBits, FracBits> Type;
00767         //typedef FixedPoint<IntBits, FracBits> Promote;
00768         //typedef FixedPoint<IntBits, FracBits> RealPromote;
00769         //typedef std::complex<RealPromote> ComplexPromote;
00770     typedef Type ValueType;
00771 
00772     typedef VigraFalseType isIntegral;
00773     typedef VigraTrueType  isScalar;
00774     typedef VigraTrueType  isSigned;
00775     typedef VigraTrueType  isOrdered;
00776     typedef VigraFalseType isComplex;
00777 
00778     static Type zero() { return Type(0, FPNoShift); }
00779     static Type one() { return Type(Type::ONE, FPNoShift); }
00780     static Type nonZero() { return one(); }
00781     static Type epsilon() { return Type(1, FPNoShift); }
00782     static Type smallestPositive() { return Type(1, FPNoShift); }
00783     static Type max() { return Type( Type::MAX, FPNoShift); }
00784     static Type min() { return -max(); }
00785 };
00786 
00787 template <unsigned IntBits, unsigned FracBits>
00788 struct NormTraits<FixedPoint<IntBits, FracBits> >
00789 {
00790     typedef FixedPoint<IntBits, FracBits>         Type;
00791     typedef typename 
00792         FixedPointTraits<FixedPoint<IntBits, FracBits>, FixedPoint<IntBits, FracBits> >::MultipliesType
00793                                                   SquaredNormType;
00794     typedef Type                                  NormType;
00795 };
00796 
00797 template <unsigned IntBits1, unsigned FracBits1, unsigned IntBits2, unsigned FracBits2>
00798 struct PromoteTraits<FixedPoint<IntBits1, FracBits1>,
00799                      FixedPoint<IntBits2, FracBits2> >
00800 {
00801     typedef typename 
00802         FixedPointTraits<FixedPoint<IntBits1, FracBits1>, FixedPoint<IntBits2, FracBits2> >::PlusType 
00803         Promote;
00804 };
00805 
00806 /***********************************************************************************/
00807 
00808 enum FPOverflowHandling { FPOverflowIgnore, FPOverflowSaturate, FPOverflowError };
00809 
00810 template <int IntBits, FPOverflowHandling OverflowHandling = FPOverflowIgnore>
00811 class FixedPoint16;
00812 
00813 /********************************************************/
00814 /*                                                      */
00815 /*                     FixedPoint16-Traits              */
00816 /*                                                      */
00817 /********************************************************/
00818 
00819 /** \page FixedPoint16Traits Numeric and Promote Traits of FixedPoint16
00820 
00821     The numeric and promote traits for FixedPoint16 follow
00822     the general specifications for \ref NumericPromotionTraits and
00823     \ref AlgebraicRing. They are implemented in terms of the traits of the basic types by
00824     partial template specialization:
00825 
00826     \code
00827     template <int IntBits, FPOverflowHandling OverflowHandling>
00828     struct NumericTraits<FixedPoint16<IntBits, OverflowHandling> >
00829     {
00830         typedef FixedPoint16<IntBits, OverflowHandling> Type;
00831         typedef Type                                    Promote;
00832             // RealPromote undefined -- multiplication with double is not supported.
00833             // ComplexPromote undefined -- multiplication with double is not supported.
00834         typedef Type ValueType;
00835 
00836         typedef VigraFalseType isIntegral;
00837         typedef VigraTrueType  isScalar;
00838         typedef VigraTrueType  isSigned;
00839         typedef VigraTrueType  isOrdered;
00840         typedef VigraFalseType isComplex;
00841 
00842         ... // etc.
00843     };
00844 
00845     template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
00846     struct PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>,
00847                          FixedPoint16<IntBits2, OverflowHandling> >
00848     {
00849         typedef FixedPoint16<MetaMax<IntBits1, IntBits2>::value, OverflowHandling> Promote;
00850         ... // etc.
00851     };
00852 
00853     template <int IntBits, FPOverflowHandling OverflowHandling>
00854     struct NormTraits<FixedPoint16<IntBits, OverflowHandling> >
00855     {
00856         typedef FixedPoint16<IntBits, OverflowHandling>     Type;
00857         typedef typename PromoteTraits<Type, Type>::Promote SquaredNormType;
00858         typedef Type                                        NormType;
00859     };
00860 
00861     template <int IntBits, FPOverflowHandling OverflowHandling>
00862     struct SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >
00863     {
00864         typedef FixedPoint16<IntBits, OverflowHandling>            Type;
00865         typedef FixedPoint16<(IntBits + 1) / 2, OverflowHandling>  SquareRootResult;
00866         typedef Type                                               SquareRootArgument;
00867     };
00868     \endcode
00869 
00870     <b>\#include</b> <<a href="fixedpoint_8hxx-source.html">vigra/fixedpoint.hxx</a>><br>
00871     Namespace: vigra
00872 
00873 */
00874 template <int IntBits, FPOverflowHandling OverflowHandling>
00875 struct NumericTraits<FixedPoint16<IntBits, OverflowHandling> >
00876 {
00877     typedef FixedPoint16<IntBits, OverflowHandling> Type;
00878     typedef Type                                    Promote;
00879         // RealPromote undefined -- multiplication with double is not supported.
00880         // ComplexPromote undefined -- multiplication with double is not supported.
00881     typedef Type ValueType;
00882 
00883     typedef VigraFalseType isIntegral;
00884     typedef VigraTrueType  isScalar;
00885     typedef VigraTrueType  isSigned;
00886     typedef VigraTrueType  isOrdered;
00887     typedef VigraFalseType isComplex;
00888 
00889     static Type zero() { return Type(0, FPNoShift); }
00890     static Type one() { return Type(Type::ONE, FPNoShift); }
00891     static Type nonZero() { return one(); }
00892     static Type epsilon() { return Type(1, FPNoShift); }
00893     static Type smallestPositive() { return Type(1, FPNoShift); }
00894     static Type max() { return Type( Type::MAX, FPNoShift); }
00895     static Type min() { return Type( Type::MIN, FPNoShift); }
00896 
00897     static Promote toPromote(Type v) { return v; }
00898     static Type fromPromote(Promote v) { return v; }; 
00899 };
00900 
00901 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
00902 struct PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>,
00903                      FixedPoint16<IntBits2, OverflowHandling> >
00904 {
00905     typedef FixedPoint16<MetaMax<IntBits1, IntBits2>::value, OverflowHandling> Promote;
00906     static Promote toPromote(FixedPoint16<IntBits1, OverflowHandling> v) { return Promote(v); }
00907     static Promote toPromote(FixedPoint16<IntBits2, OverflowHandling> v) { return Promote(v); }
00908 };
00909 
00910 template <int IntBits, FPOverflowHandling OverflowHandling>
00911 struct PromoteTraits<FixedPoint16<IntBits, OverflowHandling>,
00912                      FixedPoint16<IntBits, OverflowHandling> >
00913 {
00914     typedef FixedPoint16<IntBits, OverflowHandling> Promote;
00915     static Promote toPromote(FixedPoint16<IntBits, OverflowHandling> v) { return v; }
00916 };
00917 
00918 template <int IntBits, FPOverflowHandling OverflowHandling>
00919 struct NormTraits<FixedPoint16<IntBits, OverflowHandling> >
00920 {
00921     typedef FixedPoint16<IntBits, OverflowHandling>     Type;
00922     typedef typename PromoteTraits<Type, Type>::Promote SquaredNormType;
00923     typedef Type                                        NormType;
00924 };
00925 
00926 template <int IntBits, FPOverflowHandling OverflowHandling>
00927 struct SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >
00928 {
00929     typedef FixedPoint16<IntBits, OverflowHandling>            Type;
00930     typedef FixedPoint16<(IntBits + 1) / 2, OverflowHandling>  SquareRootResult;
00931     typedef Type                                               SquareRootArgument;
00932 };
00933 
00934 #ifndef DOXYGEN
00935 
00936 template <bool Compatible>
00937 struct FixedPoint_error__Right_shift_operator_has_unsupported_semantics
00938 : staticAssert::AssertBool<Compatible>
00939 {};
00940 
00941 #endif /* DOXYGEN */
00942 
00943 template <bool Predicate>
00944 struct FixedPoint16_assignment_error__Target_object_has_too_few_integer_bits
00945 : staticAssert::AssertBool<Predicate>
00946 {};
00947 
00948 namespace detail {
00949 
00950 template<int BeforeIntBits, int AfterIntBits, 
00951          bool Round = false,
00952          bool RightShift = (AfterIntBits >= BeforeIntBits)>
00953 struct FP16Align;
00954 
00955 template<int BeforeIntBits>
00956 struct FP16Align<BeforeIntBits, BeforeIntBits, true, true>
00957 {
00958     static inline Int32 exec(Int32 v)
00959     {
00960         return v;
00961     }
00962 };
00963 
00964 template<int BeforeIntBits>
00965 struct FP16Align<BeforeIntBits, BeforeIntBits, false, true>
00966 {
00967     static inline Int32 exec(Int32 v)
00968     {
00969         return v;
00970     }
00971 };
00972 
00973 template<int BeforeIntBits, int AfterIntBits>
00974 struct FP16Align<BeforeIntBits, AfterIntBits, false, true>
00975 {
00976     static inline Int32 exec(Int32 v)
00977     {
00978         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
00979         return v >> (AfterIntBits - BeforeIntBits);
00980     }
00981 };
00982 
00983 template<int BeforeIntBits, int AfterIntBits>
00984 struct FP16Align<BeforeIntBits, AfterIntBits, true, true>
00985 {
00986     enum { ONE_HALF = 1 << (AfterIntBits - BeforeIntBits - 1) };
00987     static inline Int32 exec(Int32 v)
00988     {
00989         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
00990         return (v + ONE_HALF) >> (AfterIntBits - BeforeIntBits);
00991     }
00992 };
00993 
00994 template<int BeforeIntBits, int AfterIntBits, bool Round>
00995 struct FP16Align<BeforeIntBits, AfterIntBits, Round, false>
00996 {
00997     static inline Int32 exec(Int32 v)
00998     {
00999         return v << (BeforeIntBits - AfterIntBits);
01000     }
01001 };
01002 
01003 template <FPOverflowHandling OverflowHandling = FPOverflowIgnore>
01004 struct FP16OverflowHandling
01005 {
01006     static inline Int32 exec(Int32 v)
01007     {
01008         return v;
01009     }
01010     
01011     static inline Int32 exec(UInt32 v)
01012     {
01013         return v;
01014     }
01015 };
01016 
01017 template <>
01018 struct FP16OverflowHandling<FPOverflowSaturate>
01019 {
01020     static inline Int32 exec(Int32 v)
01021     {
01022         if(v >= 1 << 15)
01023             return (1 << 15) - 1;
01024         if(v < -(1 << 15))
01025             return -(1 << 15);
01026         return v;
01027     }
01028     static inline Int32 exec(UInt32 v)
01029     {
01030         if(v >= 1 << 15)
01031             return (1 << 15) - 1;
01032         return v;
01033     }
01034 };
01035 
01036 template <>
01037 struct FP16OverflowHandling<FPOverflowError>
01038 {
01039     static inline Int32 exec(Int32 v)
01040     {
01041         vigra_precondition(v < (1 << 15) && v >= -(1 << 15),
01042              "FixedPoint16: Operation overflows.");
01043         return v;
01044     }
01045     static inline Int32 exec(UInt32 v)
01046     {
01047         vigra_precondition(v < (1 << 15),
01048              "FixedPoint16: Operation overflows.");
01049         return v;
01050     }
01051 };
01052 
01053 
01054 template <int IntBits1, int IntBits2, int IntBitsOut, 
01055           FPOverflowHandling OverflowHandling >
01056 struct FP16AddImpl
01057 {
01058     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
01059     static inline Int32 exec(Int32 t1, Int32 t2)
01060     {
01061        return FP16OverflowHandling<OverflowHandling>::exec(
01062                   FP16Align<MinIntBits, IntBitsOut, /*Round*/ true>::exec(
01063                       FP16Align<IntBits1, MinIntBits, /*Round*/ false>::exec(t1) + 
01064                       FP16Align<IntBits2, MinIntBits, /*Round*/ false>::exec(t2)));
01065     }
01066 };
01067 
01068 template <int IntBits1, int IntBits2, int IntBitsOut, 
01069           FPOverflowHandling OverflowHandling >
01070 struct FP16SubImpl
01071 {
01072     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
01073     static inline Int32 exec(Int32 t1, Int32 t2)
01074     {
01075        return FP16OverflowHandling<OverflowHandling>::exec(
01076                   FP16Align<MinIntBits, IntBitsOut, /*Round*/ true>::exec(
01077                       FP16Align<IntBits1, MinIntBits, /*Round*/ false>::exec(t1) - 
01078                       FP16Align<IntBits2, MinIntBits, /*Round*/ false>::exec(t2)));
01079     }
01080 };
01081 
01082 template <int IntBits1, int IntBits2, int IntBitsOut, 
01083           FPOverflowHandling OverflowHandling >
01084 struct FP16MulImpl
01085 {
01086     static inline Int32 exec(Int32 t1, Int32 t2)
01087     {
01088         return FP16OverflowHandling<OverflowHandling>::exec(
01089                    FP16Align<IntBits1+IntBits2, IntBitsOut+15, /*Round*/ true>::exec(t1*t2));
01090     }
01091 };
01092 
01093 template <int IntBits1, int IntBits2, int IntBitsOut, 
01094           FPOverflowHandling OverflowHandling >
01095 struct FP16DivImpl
01096 {
01097     static inline Int32 exec(Int32 t1, Int32 t2)
01098     {
01099         if(t2 == 0)
01100             return (t1 >= 0)
01101                        ?  (1 << 15) - 1
01102                        : -(1 << 15);
01103         return FP16OverflowHandling<OverflowHandling>::exec(
01104                    FP16Align<IntBits1-IntBits2, IntBitsOut+1, /*Round*/ true>::exec((t1<<16)/t2));
01105     }
01106 };
01107 
01108 } // namespace detail
01109 
01110 /********************************************************/
01111 /*                                                      */
01112 /*                      FixedPoint16                    */
01113 /*                                                      */
01114 /********************************************************/
01115 
01116 template <class TARGET, int IntBits, FPOverflowHandling OverflowHandling>
01117 TARGET fixed_point_cast(FixedPoint16<IntBits, OverflowHandling> v);
01118 
01119 /** Template for 16-bit signed fixed point arithmetic.
01120 
01121     Fixed point arithmetic is used when computations with fractional accuracy
01122     must be made at the highest speed possible (e.g. in the inner loop
01123     of a volume rendering routine). The speed-up relative to floating
01124     point arithmetic can be dramatic, especially when one can avoid
01125     conversions between integer and floating point numbers (these are 
01126     very expensive because integer and floating point arithmetic
01127     resides in different pipelines). 
01128     
01129     The template wraps an <tt>Int16</tt> and uses <tt>IntBits</tt> to
01130     represent the integral part of a number, and <tt>15 - IntBits</tt>
01131     for the fractional part. The 16th bit is reserved because FixedPoint16 
01132     is a signed type. Results of expressions with mixed types will preserve
01133     larger number of <tt>IntBits</tt> of the results, in order to minimize
01134     the possibility for overflow. Nonetheless, overflow can occur, and the 
01135     template parameter <tt>OverflowHandling</tt> determines how this will be
01136     handled:
01137     
01138     <DL>
01139     <DT>FPOverflowIgnore<DD> (default) Ignore overflow, i.e. use the usual modulo behavior of the
01140                              built-in integer types.
01141                        
01142     <DT>FPOverflowSaturate<DD> Use the largest or smallest representable number (depending on sign)
01143                                in case of overflow.
01144 
01145     <DT>FPOverflowError<DD> Throw <tt>PreconditionViolation</tt> upon overflow. This is useful for 
01146                             debugging.
01147     </DL>
01148     
01149     The implementation relies on Int32-arithmetic and requires that the right-shift operator
01150     preserves signedness. Although not enforced by the C++ standard, this is implemented
01151     by most of today's processors. This property is checked by a 
01152     VIGRA_STATIC_ASSERT(FixedPoint_error__Right_shift_operator_has_unsupported_semantics).
01153 
01154     <tt>FixedPoint16</tt> implements the required interface of an
01155     \ref AlgebraicRing and the required numeric and
01156     promotion traits. In addition, it supports functions <tt>add</tt>, 
01157     <tt>sub</tt>, <tt>mul</tt>, and <tt>div</tt>, where a particular layout 
01158     of the result can be enforced. 
01159     
01160     Built-in numeric types can be converted into <tt>FixedPoint16</tt> by the 
01161     appropriate constructors, and from <tt>FixedPoint16</tt> by means of
01162     <tt>fixed_point_cast&lt;TargetType&gt;(fixedPoint)</tt>.
01163 
01164     <b>See also:</b>
01165     <ul>
01166     <li> \ref FixedPoint16Operations
01167     <li> \ref FixedPoint16Traits
01168     </ul>
01169 
01170     <b>\#include</b> <<a href="fixedpoint_8hxx-source.html">vigra/fixedpoint.hxx</a>><br>
01171     Namespace: vigra
01172 */
01173 template <int IntBits, FPOverflowHandling OverflowHandling>
01174 class FixedPoint16
01175 {
01176 public:
01177     static const Int32 TOTAL_BITS      = 15; // bit 16 is sign
01178     static const Int32 INT_BITS        = IntBits;
01179     static const Int32 FRACTIONAL_BITS = TOTAL_BITS - INT_BITS;
01180     static const Int32 MAX             = (Int32)((1u << TOTAL_BITS) - 1);
01181     static const Int32 MIN             = -(Int32)(1u << TOTAL_BITS);
01182     static const Int32 ONE             = 1 << FRACTIONAL_BITS;
01183     static const Int32 ONE_HALF        = ONE >> 1;
01184     static const Int32 FRACTIONAL_MASK = (1u << FRACTIONAL_BITS) - 1;
01185     static const Int32 INT_MASK        = 0xffffffffu ^ FRACTIONAL_MASK;
01186 
01187     Int16 value;
01188 
01189     FixedPoint16()
01190     : value(0)
01191     {
01192         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
01193     }
01194 
01195         /** Construct from an int (fractional part will become zero).
01196             Possible overflow is handled according to the taget type's <tt>OverflowHandling</tt>.
01197         */
01198     explicit FixedPoint16(Int32 v)
01199     : value(detail::FP16OverflowHandling<OverflowHandling>::exec(v << FRACTIONAL_BITS))
01200     {
01201         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
01202     }
01203 
01204         /** Construct from an int by a bitwise copy. This is normally only used internally.
01205         */
01206     FixedPoint16(Int32 v, FixedPointNoShift)
01207     : value(detail::FP16OverflowHandling<OverflowHandling>::exec(v))
01208     {
01209         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
01210     }
01211 
01212         /** Construct from a double and round the fractional part to 
01213             <tt>FRACTIONAL_BITS</tt> accuracy. Possible overflow is handled according 
01214             to the taget type's <tt>OverflowHandling</tt>.
01215         */
01216     explicit FixedPoint16(double rhs)
01217     : value(detail::FP16OverflowHandling<OverflowHandling>::exec(roundi(rhs * ONE)))
01218     {
01219         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
01220     }
01221 
01222         /** Copy constructor.
01223         */
01224     FixedPoint16(const FixedPoint16 &other)
01225     : value(other.value)
01226     {
01227         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
01228     }
01229 
01230         /** Construct from a FixedPoint16 with different layout. It rounds as appropriate and 
01231             handles possible overflow according to the taget type's <tt>OverflowHandling</tt>.
01232         */
01233     template <int IntBits2, FPOverflowHandling OverflowHandling2>
01234     FixedPoint16(const FixedPoint16<IntBits2, OverflowHandling2> &other)
01235     : value(detail::FP16OverflowHandling<OverflowHandling>::exec(
01236             detail::FP16Align<IntBits2, IntBits, /*Round*/true>::exec(other.value)))
01237     {
01238         VIGRA_STATIC_ASSERT((FixedPoint_error__Right_shift_operator_has_unsupported_semantics<((-1 >> 8) == -1)>));
01239     }
01240 
01241         /** Assignment from int. The fractional part will become zero.  
01242             Possible overflow is handled according to the taget type's <tt>OverflowHandling</tt>.
01243         */
01244     FixedPoint16 &operator=(Int32 rhs)
01245     {
01246         value = detail::FP16OverflowHandling<OverflowHandling>::exec(rhs << FRACTIONAL_BITS);
01247         return *this;
01248     }
01249 
01250         /** Assignment form double. The fractional part is rounded, and possible overflow is 
01251             handled according to the taget type's <tt>OverflowHandling</tt>.
01252         */
01253     FixedPoint16 &operator=(double rhs)
01254     {
01255         value = detail::FP16OverflowHandling<OverflowHandling>::exec(roundi(rhs * ONE));
01256         return *this;
01257     }
01258 
01259         /** Copy assignment.
01260         */
01261     FixedPoint16 & operator=(const FixedPoint16 &other)
01262     {
01263         value = other.value;
01264         return *this;
01265     }
01266 
01267         /** Assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
01268             handled according to the taget type's <tt>OverflowHandling</tt>.
01269         */
01270     template <int IntBits2>
01271     FixedPoint16 & operator=(const FixedPoint16<IntBits2, OverflowHandling> &other)
01272     {
01273         value = detail::FP16OverflowHandling<OverflowHandling>::exec(
01274                 detail::FP16Align<IntBits2, IntBits, /*Round*/true>::exec(other.value));
01275         return *this;
01276     }
01277 
01278         /** Conversion to float
01279         */
01280     operator float() const
01281     {
01282         return fixed_point_cast<float>(*this);
01283     }
01284 
01285         /** Conversion to double
01286         */
01287     operator double() const
01288     {
01289         return fixed_point_cast<double>(*this);
01290     }
01291 
01292         /** Unary plus.
01293         */
01294     FixedPoint16 operator+() const
01295     {
01296         return *this;
01297     }
01298 
01299         /** Negation.
01300         */
01301     FixedPoint16 operator-() const
01302     {
01303         return FixedPoint16(-value, FPNoShift);
01304     }
01305 
01306         /** Pre-increment.
01307         */
01308     FixedPoint16 & operator++()
01309     {
01310         value = detail::FP16OverflowHandling<OverflowHandling>::exec(value+ONE);
01311         return *this;
01312     }
01313 
01314         /** Post-increment.
01315         */
01316     FixedPoint16 operator++(int)
01317     {
01318         FixedPoint16 old(*this);
01319         ++(*this);
01320         return old;
01321     }
01322 
01323         /** Pre-decrement.
01324         */
01325     FixedPoint16 & operator--()
01326     {
01327         value = detail::FP16OverflowHandling<OverflowHandling>::exec(value-ONE);
01328         return *this;
01329     }
01330 
01331         /** Post-decrement.
01332         */
01333     FixedPoint16 operator--(int)
01334     {
01335         FixedPoint16 old(*this);
01336         --(*this);
01337         return old;
01338     }
01339 
01340         /** Add-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
01341             handled according to the taget type's <tt>OverflowHandling</tt>.
01342         */
01343     template <int IntBits2>
01344     FixedPoint16 & operator+=(const FixedPoint16<IntBits2, OverflowHandling> &other)
01345     {
01346         value = detail::FP16AddImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
01347         return *this;
01348     }
01349 
01350         /** Subtract-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
01351             handled according to the taget type's <tt>OverflowHandling</tt>.
01352         */
01353     template <int IntBits2>
01354     FixedPoint16 & operator-=(const FixedPoint16<IntBits2, OverflowHandling> &other)
01355     {
01356         value = detail::FP16SubImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
01357         return *this;
01358     }
01359     
01360         /** Multiply-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
01361             handled according to the taget type's <tt>OverflowHandling</tt>.
01362         */
01363     template <int IntBits2>
01364     FixedPoint16 & operator*=(const FixedPoint16<IntBits2, OverflowHandling> &other)
01365     {
01366         value = detail::FP16MulImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
01367         return *this;
01368     }
01369     
01370         /** Divide-assignment from a FixedPoint16 with different layout. It rounds as appropriate, and possible overflow is 
01371             handled according to the taget type's <tt>OverflowHandling</tt>.
01372         */
01373     template <int IntBits2>
01374     FixedPoint16 & operator/=(const FixedPoint16<IntBits2, OverflowHandling> &other)
01375     {
01376         value = detail::FP16DivImpl<IntBits, IntBits2, IntBits, OverflowHandling>::exec(value, other.value);
01377         return *this;
01378     }
01379 };
01380 
01381 namespace detail {
01382 
01383 template <class T>
01384 struct FixedPoint16Cast;
01385 
01386 #define VIGRA_FIXED_POINT_CAST(type) \
01387 template <> \
01388 struct FixedPoint16Cast<type> \
01389 { \
01390     template <int IntBits, FPOverflowHandling OverflowHandling> \
01391     static type cast(FixedPoint16<IntBits, OverflowHandling> v) \
01392     { \
01393         return round(v); \
01394     } \
01395 };
01396 
01397 VIGRA_FIXED_POINT_CAST(Int8)
01398 VIGRA_FIXED_POINT_CAST(UInt8)
01399 VIGRA_FIXED_POINT_CAST(Int16)
01400 VIGRA_FIXED_POINT_CAST(UInt16)
01401 VIGRA_FIXED_POINT_CAST(Int32)
01402 VIGRA_FIXED_POINT_CAST(UInt32)
01403 VIGRA_FIXED_POINT_CAST(Int64)
01404 VIGRA_FIXED_POINT_CAST(UInt64)
01405 
01406 #undef VIGRA_FIXED_POINT_CAST
01407 
01408 template <>
01409 struct FixedPoint16Cast<float>
01410 {
01411     template <int IntBits, FPOverflowHandling OverflowHandling>
01412     static float cast(FixedPoint16<IntBits, OverflowHandling> v)
01413     {
01414         return (float)v.value / FixedPoint16<IntBits, OverflowHandling>::ONE;
01415     }
01416 };
01417 
01418 template <>
01419 struct FixedPoint16Cast<double>
01420 {
01421     template <int IntBits, FPOverflowHandling OverflowHandling>
01422     static double cast(FixedPoint16<IntBits, OverflowHandling> v)
01423     {
01424         return (double)v.value / FixedPoint16<IntBits, OverflowHandling>::ONE;
01425     }
01426 };
01427 
01428 } // namespace detail
01429 
01430 /********************************************************/
01431 /*                                                      */
01432 /*                 FixedPoint16Operations               */
01433 /*                                                      */
01434 /********************************************************/
01435 
01436 /** \addtogroup FixedPoint16Operations Functions for FixedPoint16
01437 
01438     \brief     <b>\#include</b> <<a href="fixedpoint_8hxx-source.html">vigra/fixedpoint.hxx</a>><br>
01439 
01440     These functions fulfill the requirements of an \ref AlgebraicRing.
01441 
01442     Namespace: vigra
01443     <p>
01444 
01445  */
01446 //@{
01447 
01448     /** Convert a FixedPoint16 to a built-in type.
01449         If the target is integral, the value is rounded.<br>
01450         Usage:
01451         \code
01452         FixedPoint16<16,15> fp(...);
01453         
01454         double d = fixed_point_cast<double>(fp);
01455         \endcode
01456     */
01457 template <class TARGET, int IntBits, FPOverflowHandling OverflowHandling>
01458 TARGET fixed_point_cast(FixedPoint16<IntBits, OverflowHandling> v)
01459 {
01460     return detail::FixedPoint16Cast<TARGET>::cast(v);
01461 }
01462 
01463 
01464     /// equal
01465 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01466 inline
01467 bool operator==(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01468 {
01469     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
01470     return (l.value << (IntBits1 - MinIntBits)) == (r.value << (IntBits2 - MinIntBits));
01471 }
01472 
01473     /// not equal
01474 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01475 inline
01476 bool operator!=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01477 {
01478     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
01479     return (l.value << (IntBits1 - MinIntBits)) != (r.value << (IntBits2 - MinIntBits));
01480 }
01481 
01482     /// less than
01483 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01484 inline
01485 bool operator<(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01486 {
01487     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
01488     return (l.value << (IntBits1 - MinIntBits)) < (r.value << (IntBits2 - MinIntBits));
01489 }
01490 
01491     /// less or equal
01492 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01493 inline
01494 bool operator<=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01495 {
01496     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
01497     return (l.value << (IntBits1 - MinIntBits)) <= (r.value << (IntBits2 - MinIntBits));
01498 }
01499 
01500     /// greater
01501 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01502 inline
01503 bool operator>(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01504 {
01505     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
01506     return (l.value << (IntBits1 - MinIntBits)) > (r.value << (IntBits2 - MinIntBits));
01507 }
01508 
01509     /// greater or equal
01510 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01511 inline
01512 bool operator>=(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01513 {
01514     enum { MinIntBits = MetaMin<IntBits1, IntBits2>::value };
01515     return (l.value << (IntBits1 - MinIntBits)) >= (r.value << (IntBits2 - MinIntBits));
01516 }
01517 
01518     /// addition with automatic determination of the appropriate result type.
01519 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01520 inline
01521 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
01522 operator+(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01523 {
01524     typedef typename
01525         PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
01526         Result;
01527     return Result(detail::FP16AddImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
01528 }
01529 
01530     /// addition with enforced result type.
01531 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
01532 inline 
01533 FixedPoint16<IntBits3, OverflowHandling> &
01534 add(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
01535     FixedPoint16<IntBits3, OverflowHandling> & result)
01536 {
01537     result.value = detail::FP16AddImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
01538     return result;
01539 }
01540 
01541     /// subtraction with automatic determination of the appropriate result type.
01542 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01543 inline
01544 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
01545 operator-(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01546 {
01547     typedef typename
01548         PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
01549         Result;
01550     return Result(detail::FP16SubImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
01551 }
01552 
01553     /// subtraction with enforced result type.
01554 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
01555 inline FixedPoint16<IntBits3, OverflowHandling> &
01556 sub(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
01557     FixedPoint16<IntBits3, OverflowHandling> & result)
01558 {
01559     result.value = detail::FP16SubImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
01560     return result;
01561 }
01562 
01563     /// multiplication with automatic determination of the appropriate result type.
01564 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01565 inline
01566 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
01567 operator*(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01568 {
01569     typedef typename
01570         PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
01571         Result;
01572     return Result(detail::FP16MulImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
01573 }
01574 
01575     /// multiplication with enforced result type.
01576 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
01577 inline 
01578 FixedPoint16<IntBits3, OverflowHandling> &
01579 mul(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
01580     FixedPoint16<IntBits3, OverflowHandling> & result)
01581 {
01582     result.value = detail::FP16MulImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
01583     return result;
01584 }
01585 
01586     /// division with automatic determination of the appropriate result type.
01587 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2>
01588 inline
01589 typename PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
01590 operator/(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r)
01591 {
01592     typedef typename
01593         PromoteTraits<FixedPoint16<IntBits1, OverflowHandling>, FixedPoint16<IntBits2, OverflowHandling> >::Promote
01594         Result;
01595     return Result(detail::FP16DivImpl<IntBits1, IntBits2, Result::INT_BITS, OverflowHandling>::exec(l.value, r.value), FPNoShift);
01596 }
01597 
01598     /// division with enforced result type.
01599 template <int IntBits1, FPOverflowHandling OverflowHandling, int IntBits2, int IntBits3>
01600 inline 
01601 FixedPoint16<IntBits3, OverflowHandling> &
01602 div(FixedPoint16<IntBits1, OverflowHandling> l, FixedPoint16<IntBits2, OverflowHandling> r,
01603     FixedPoint16<IntBits3, OverflowHandling> & result)
01604 {
01605     result.value = detail::FP16DivImpl<IntBits1, IntBits2, IntBits3, OverflowHandling>::exec(l.value, r.value);
01606     return result;
01607 }
01608 
01609     /// square root.
01610 template <int IntBits, FPOverflowHandling OverflowHandling>
01611 inline typename SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >::SquareRootResult
01612 sqrt(FixedPoint16<IntBits, OverflowHandling> v)
01613 {
01614     typedef typename SquareRootTraits<FixedPoint16<IntBits, OverflowHandling> >::SquareRootResult Result;
01615     enum { Shift = 15 + IntBits - 2*Result::INT_BITS };
01616     return Result(sqrti(v.value << Shift), FPNoShift);
01617 }
01618 
01619 #ifndef VIGRA_NO_HYPOT
01620     using ::hypot;
01621 #endif
01622 
01623     /// Length of hypothenuse. 
01624 template <int IntBits, FPOverflowHandling OverflowHandling>
01625 inline FixedPoint16<IntBits, OverflowHandling>
01626 hypot(FixedPoint16<IntBits, OverflowHandling> v1, FixedPoint16<IntBits, OverflowHandling> v2)
01627 {
01628     UInt32 l = abs(v1.value), r = abs(v2.value);   
01629     // sq(l) + sq(r) <= 2**31, so overflow handling after sqrti is sufficient
01630     return FixedPoint16<IntBits, OverflowHandling>(
01631                    detail::FP16OverflowHandling<OverflowHandling>::exec(sqrti(sq(l) + sq(r))), 
01632                    FPNoShift);
01633 }
01634 
01635 using std::atan2;
01636 
01637     /// Arctangent. Accuracy better than 1/3 degree (9 significant bits).
01638 template <int IntBits, FPOverflowHandling OverflowHandling>
01639 FixedPoint16<2, OverflowHandling>
01640 atan2(FixedPoint16<IntBits, OverflowHandling> y, FixedPoint16<IntBits, OverflowHandling> x)
01641 {
01642     enum { ResIntBits = 2 };
01643     typedef FixedPoint16<ResIntBits, OverflowHandling> FP;
01644     static const FP zero(0), pi(M_PI), pi_2(0.5 * M_PI), mpi_2(-0.5 * M_PI); 
01645     static const Int32 Pi_4  = roundi(0.25 * M_PI * (1 << 15)), // 15 frac bits
01646                        Pi3_4 = roundi(0.75 * M_PI * (1 << 15)),
01647                        c1    = roundi(0.19826763260224867 * (1 << 15)), 
01648                        c2    = roundi(-0.9757748231899761 * (1 << 30)); 
01649     // coefficients c1 and c2 minimize
01650     //
01651     // NIntegrate[(c1 r^3 + c2 r + Pi/4 - a)^4 /. r -> (Cos[a] - Sin[a])/(Cos[a] + Sin[a]), {a, 0, Pi/2}]
01652     //
01653     // Thanks to Jim Shima, http://www.dspguru.com/comp.dsp/tricks/alg/fxdatan2.htm
01654 
01655     if(x.value == 0)
01656         return (y.value > 0)
01657                    ? pi_2
01658                    : (y.value < 0)
01659                          ? mpi_2
01660                          : zero;
01661                          
01662     Int32 abs_y = abs(y.value);
01663     Int32 r, angle;
01664     if(x.value > 0)
01665     {
01666         if(y.value == 0)
01667             return zero;
01668         r = ((x.value - abs_y) << 15) / (x.value + abs_y); // 15 frac bits
01669         angle = Pi_4;
01670     }
01671     else
01672     {
01673         if(y.value == 0)
01674             return pi;
01675         r = ((x.value + abs_y) << 15) / (abs_y - x.value); // 15 frac bits
01676         angle = Pi3_4;
01677     }
01678     
01679     angle += r*((c2 + c1 * (sq(r) >> 15)) >> 15) >> 15;
01680 
01681     return (y.value > 0)
01682                ? FP(detail::FP16Align<0, ResIntBits, true>::exec( angle), FPNoShift)
01683                : FP(detail::FP16Align<0, ResIntBits, true>::exec(-angle), FPNoShift);
01684 }
01685 
01686     /// absolute value.
01687 template <int IntBits, FPOverflowHandling OverflowHandling>
01688 inline FixedPoint16<IntBits, OverflowHandling>
01689 abs(FixedPoint16<IntBits, OverflowHandling> v)
01690 {
01691     return FixedPoint16<IntBits, OverflowHandling>(abs(v.value), FPNoShift);
01692 }
01693 
01694     /// squared norm (same as v*v).
01695 template <int IntBits, FPOverflowHandling OverflowHandling>
01696 inline
01697 typename NormTraits<FixedPoint16<IntBits, OverflowHandling> >::SquaredNormType
01698 squaredNorm(FixedPoint16<IntBits, OverflowHandling> v)
01699 {
01700     return v*v;
01701 }
01702 
01703     /// norm (same as abs).
01704 template <int IntBits, FPOverflowHandling OverflowHandling>
01705 inline 
01706 typename NormTraits<FixedPoint16<IntBits, OverflowHandling> >::NormType
01707 norm(FixedPoint16<IntBits, OverflowHandling> const & v)
01708 {
01709     return abs(v);
01710 }
01711 
01712     /// fractional part. (difference between v and its floor)
01713 template <int IntBits, FPOverflowHandling OverflowHandling>
01714 inline FixedPoint16<IntBits, OverflowHandling>
01715 frac(FixedPoint16<IntBits, OverflowHandling> v)
01716 {
01717     return FixedPoint16<IntBits, OverflowHandling>(
01718            v.value - (v.value & FixedPoint16<IntBits, OverflowHandling>::INT_MASK), 
01719            FPNoShift);
01720 }
01721 
01722     /// dual fractional part. (1 - frac(v))
01723 template <int IntBits, FPOverflowHandling OverflowHandling>
01724 inline FixedPoint16<IntBits, OverflowHandling>
01725 dual_frac(FixedPoint16<IntBits, OverflowHandling> v)
01726 {
01727     return FixedPoint16<IntBits, OverflowHandling>(
01728            FixedPoint16<IntBits, OverflowHandling>::ONE - v.value + (v.value & FixedPoint16<IntBits, OverflowHandling>::INT_MASK), 
01729            FPNoShift);
01730 }
01731 
01732     /// rounding down.
01733 template <int IntBits, FPOverflowHandling OverflowHandling>
01734 inline Int32
01735 floor(FixedPoint16<IntBits, OverflowHandling> v)
01736 {
01737     return(v.value >> FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
01738 }
01739 
01740     /// rounding up.
01741 template <int IntBits, FPOverflowHandling OverflowHandling>
01742 inline Int32
01743 ceil(FixedPoint16<IntBits, OverflowHandling> v)
01744 {
01745     return((v.value + FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_MASK) >> 
01746                       FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
01747 }
01748 
01749     /// rounding to the nearest integer.
01750 template <int IntBits, FPOverflowHandling OverflowHandling>
01751 inline Int32
01752 round(FixedPoint16<IntBits, OverflowHandling> v)
01753 {
01754     return((v.value + FixedPoint16<IntBits, OverflowHandling>::ONE_HALF) >> 
01755                       FixedPoint16<IntBits, OverflowHandling>::FRACTIONAL_BITS);
01756 }
01757 
01758     /// rounding to the nearest integer.
01759 template <int IntBits, FPOverflowHandling OverflowHandling>
01760 inline Int32
01761 roundi(FixedPoint16<IntBits, OverflowHandling> v)
01762 {
01763     return round(v);
01764 }
01765 
01766 //@}
01767 
01768 } // namespace vigra
01769 
01770 namespace std {
01771 
01772 template <int IntBits, vigra::FPOverflowHandling OverflowHandling>
01773 ostream & operator<<(ostream & s, vigra::FixedPoint16<IntBits, OverflowHandling> v)
01774 {
01775     s << vigra::fixed_point_cast<float>(v);
01776     return s;
01777 }
01778 
01779 } // namespace std
01780 
01781 #endif // VIGRA_FIXEDPOINT_HXX

© 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.7.0 (Thu Aug 25 2011)