[ 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 and floating 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 available. 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> <vigra/fixedpoint.hxx><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> <vigra/fixedpoint.hxx><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> <vigra/fixedpoint.hxx><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> <vigra/fixedpoint.hxx><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<TargetType>(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> <vigra/fixedpoint.hxx><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 target 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 target 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 target 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 target 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 target 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 target 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 target 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 target 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 target 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 target 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> <vigra/fixedpoint.hxx><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 hypotenuse. 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) |
html generated using doxygen and Python
|