Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: CPL
4 : * Purpose: Floating point conversion functions. Convert 16- and 24-bit
5 : * floating point numbers into the 32-bit IEEE 754 compliant ones.
6 : * Author: Andrey Kiselev, dron@remotesensing.org
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2005, Andrey Kiselev <dron@remotesensing.org>
10 : *
11 : * This code is based on the code from OpenEXR project with the following
12 : * copyright:
13 : *
14 : * Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
15 : * Digital Ltd. LLC
16 : *
17 : * All rights reserved.
18 : *
19 : * Redistribution and use in source and binary forms, with or without
20 : * modification, are permitted provided that the following conditions are
21 : * met:
22 : * * Redistributions of source code must retain the above copyright
23 : * notice, this list of conditions and the following disclaimer.
24 : * * Redistributions in binary form must reproduce the above
25 : * copyright notice, this list of conditions and the following disclaimer
26 : * in the documentation and/or other materials provided with the
27 : * distribution.
28 : * * Neither the name of Industrial Light & Magic nor the names of
29 : * its contributors may be used to endorse or promote products derived
30 : * from this software without specific prior written permission.
31 : *
32 : * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
33 : * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
34 : * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
35 : * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
36 : * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
37 : * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
38 : * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
39 : * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
40 : * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
41 : * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
42 : * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43 : *
44 : ****************************************************************************/
45 :
46 : #include "cpl_float.h"
47 : #include "cpl_error.h"
48 :
49 : #include <algorithm>
50 : #include <cmath>
51 : #include <cstring>
52 : #include <limits>
53 : #include <numeric>
54 : #include <optional>
55 :
56 : /************************************************************************/
57 : /* HalfToFloat() */
58 : /* */
59 : /* 16-bit floating point number to 32-bit one. */
60 : /************************************************************************/
61 :
62 788 : GUInt32 CPLHalfToFloat(GUInt16 iHalf)
63 : {
64 :
65 788 : GUInt32 iSign = (iHalf >> 15) & 0x00000001;
66 788 : int iExponent = (iHalf >> 10) & 0x0000001f;
67 788 : GUInt32 iMantissa = iHalf & 0x000003ff;
68 :
69 788 : if (iExponent == 0)
70 : {
71 8 : if (iMantissa == 0)
72 : {
73 : /* --------------------------------------------------------------------
74 : */
75 : /* Plus or minus zero. */
76 : /* --------------------------------------------------------------------
77 : */
78 :
79 8 : return iSign << 31;
80 : }
81 : else
82 : {
83 : /* --------------------------------------------------------------------
84 : */
85 : /* Denormalized number -- renormalize it. */
86 : /* --------------------------------------------------------------------
87 : */
88 :
89 0 : while (!(iMantissa & 0x00000400))
90 : {
91 0 : iMantissa <<= 1;
92 0 : iExponent -= 1;
93 : }
94 :
95 0 : iExponent += 1;
96 0 : iMantissa &= ~0x00000400U;
97 : }
98 : }
99 780 : else if (iExponent == 31)
100 : {
101 0 : if (iMantissa == 0)
102 : {
103 : /* --------------------------------------------------------------------
104 : */
105 : /* Positive or negative infinity. */
106 : /* --------------------------------------------------------------------
107 : */
108 :
109 0 : return (iSign << 31) | 0x7f800000;
110 : }
111 : else
112 : {
113 : /* --------------------------------------------------------------------
114 : */
115 : /* NaN -- preserve sign and significand bits. */
116 : /* --------------------------------------------------------------------
117 : */
118 :
119 0 : return (iSign << 31) | 0x7f800000 | (iMantissa << 13);
120 : }
121 : }
122 :
123 : /* -------------------------------------------------------------------- */
124 : /* Normalized number. */
125 : /* -------------------------------------------------------------------- */
126 :
127 780 : iExponent = iExponent + (127 - 15);
128 780 : iMantissa = iMantissa << 13;
129 :
130 : /* -------------------------------------------------------------------- */
131 : /* Assemble sign, exponent and mantissa. */
132 : /* -------------------------------------------------------------------- */
133 :
134 : /* coverity[overflow_sink] */
135 780 : return (iSign << 31) | (static_cast<GUInt32>(iExponent) << 23) | iMantissa;
136 : }
137 :
138 : /************************************************************************/
139 : /* TripleToFloat() */
140 : /* */
141 : /* 24-bit floating point number to 32-bit one. */
142 : /************************************************************************/
143 :
144 400 : GUInt32 CPLTripleToFloat(GUInt32 iTriple)
145 : {
146 :
147 400 : GUInt32 iSign = (iTriple >> 23) & 0x00000001;
148 400 : int iExponent = (iTriple >> 16) & 0x0000007f;
149 400 : GUInt32 iMantissa = iTriple & 0x0000ffff;
150 :
151 400 : if (iExponent == 0)
152 : {
153 0 : if (iMantissa == 0)
154 : {
155 : /* --------------------------------------------------------------------
156 : */
157 : /* Plus or minus zero. */
158 : /* --------------------------------------------------------------------
159 : */
160 :
161 0 : return iSign << 31;
162 : }
163 : else
164 : {
165 : /* --------------------------------------------------------------------
166 : */
167 : /* Denormalized number -- renormalize it. */
168 : /* --------------------------------------------------------------------
169 : */
170 :
171 0 : while (!(iMantissa & 0x00010000))
172 : {
173 0 : iMantissa <<= 1;
174 0 : iExponent -= 1;
175 : }
176 :
177 0 : iExponent += 1;
178 0 : iMantissa &= ~0x00010000U;
179 : }
180 : }
181 400 : else if (iExponent == 127)
182 : {
183 0 : if (iMantissa == 0)
184 : {
185 : /* --------------------------------------------------------------------
186 : */
187 : /* Positive or negative infinity. */
188 : /* --------------------------------------------------------------------
189 : */
190 :
191 0 : return (iSign << 31) | 0x7f800000;
192 : }
193 : else
194 : {
195 : /* --------------------------------------------------------------------
196 : */
197 : /* NaN -- preserve sign and significand bits. */
198 : /* --------------------------------------------------------------------
199 : */
200 :
201 0 : return (iSign << 31) | 0x7f800000 | (iMantissa << 7);
202 : }
203 : }
204 :
205 : /* -------------------------------------------------------------------- */
206 : /* Normalized number. */
207 : /* -------------------------------------------------------------------- */
208 :
209 400 : iExponent = iExponent + (127 - 63);
210 400 : iMantissa = iMantissa << 7;
211 :
212 : /* -------------------------------------------------------------------- */
213 : /* Assemble sign, exponent and mantissa. */
214 : /* -------------------------------------------------------------------- */
215 :
216 : /* coverity[overflow_sink] */
217 400 : return (iSign << 31) | (static_cast<GUInt32>(iExponent) << 23) | iMantissa;
218 : }
219 :
220 : /************************************************************************/
221 : /* FloatToHalf() */
222 : /************************************************************************/
223 :
224 0 : GUInt16 CPLFloatToHalf(GUInt32 iFloat32, bool &bHasWarned)
225 : {
226 0 : GUInt32 iSign = (iFloat32 >> 31) & 0x00000001;
227 0 : GUInt32 iExponent = (iFloat32 >> 23) & 0x000000ff;
228 0 : GUInt32 iMantissa = iFloat32 & 0x007fffff;
229 :
230 0 : if (iExponent == 255)
231 : {
232 0 : if (iMantissa == 0)
233 : {
234 : /* --------------------------------------------------------------------
235 : */
236 : /* Positive or negative infinity. */
237 : /* --------------------------------------------------------------------
238 : */
239 :
240 0 : return static_cast<GUInt16>((iSign << 15) | 0x7C00);
241 : }
242 : else
243 : {
244 : /* --------------------------------------------------------------------
245 : */
246 : /* NaN -- preserve sign and significand bits. */
247 : /* --------------------------------------------------------------------
248 : */
249 0 : if (iMantissa >> 13)
250 0 : return static_cast<GUInt16>((iSign << 15) | 0x7C00 |
251 0 : (iMantissa >> 13));
252 :
253 0 : return static_cast<GUInt16>((iSign << 15) | 0x7E00);
254 : }
255 : }
256 :
257 0 : if (iExponent <= 127 - 15)
258 : {
259 : // Zero, float32 denormalized number or float32 too small normalized
260 : // number
261 0 : if (13 + 1 + 127 - 15 - iExponent >= 32)
262 0 : return static_cast<GUInt16>(iSign << 15);
263 :
264 : // Return a denormalized number
265 : return static_cast<GUInt16>(
266 0 : (iSign << 15) |
267 0 : ((iMantissa | 0x00800000) >> (13 + 1 + 127 - 15 - iExponent)));
268 : }
269 0 : if (iExponent - (127 - 15) >= 31)
270 : {
271 0 : if (!bHasWarned)
272 : {
273 0 : bHasWarned = true;
274 0 : float fVal = 0.0f;
275 0 : memcpy(&fVal, &iFloat32, 4);
276 0 : CPLError(
277 : CE_Failure, CPLE_AppDefined,
278 : "Value %.8g is beyond range of float16. Converted to %sinf",
279 0 : fVal, (fVal > 0) ? "+" : "-");
280 : }
281 0 : return static_cast<GUInt16>((iSign << 15) | 0x7C00); // Infinity
282 : }
283 :
284 : /* -------------------------------------------------------------------- */
285 : /* Normalized number. */
286 : /* -------------------------------------------------------------------- */
287 :
288 0 : iExponent = iExponent - (127 - 15);
289 0 : iMantissa = iMantissa >> 13;
290 :
291 : /* -------------------------------------------------------------------- */
292 : /* Assemble sign, exponent and mantissa. */
293 : /* -------------------------------------------------------------------- */
294 :
295 : // coverity[overflow_sink]
296 0 : return static_cast<GUInt16>((iSign << 15) | (iExponent << 10) | iMantissa);
297 : }
298 :
299 0 : GUInt16 CPLConvertFloatToHalf(float fFloat32)
300 : {
301 : GUInt32 nFloat32;
302 0 : std::memcpy(&nFloat32, &fFloat32, sizeof nFloat32);
303 0 : bool bHasWarned = true;
304 0 : return CPLFloatToHalf(nFloat32, bHasWarned);
305 : }
306 :
307 0 : float CPLConvertHalfToFloat(GUInt16 nHalf)
308 : {
309 0 : GUInt32 nFloat32 = CPLHalfToFloat(nHalf);
310 : float fFloat32;
311 0 : std::memcpy(&fFloat32, &nFloat32, sizeof fFloat32);
312 0 : return fFloat32;
313 : }
314 :
315 : namespace
316 : {
317 :
318 : template <typename T> struct Fraction
319 : {
320 : using value_type = T;
321 :
322 : T num;
323 : T denom;
324 : };
325 :
326 : /** Approximate a floating point number as a fraction, using the method describe
327 : * in Richards, Ian (1981). Continued Fractions Without Tears. Mathematics
328 : * Magazine, Vol. 54, No. 4. https://doi.org/10.2307/2689627
329 : *
330 : * If the fraction cannot be approximated within the specified error tolerance
331 : * in a certain amount of iterations, a warning will be raised and std::nullopt
332 : * will be returned.
333 : *
334 : * @param x the number to approximate as a fraction
335 : * @param err the maximum allowable absolute error in the approximation
336 : *
337 : * @return the approximated value, or std::nullopt
338 : *
339 : */
340 64 : std::optional<Fraction<std::uint64_t>> FloatToFraction(double x, double err)
341 : {
342 : using inttype = std::uint64_t;
343 64 : constexpr int MAX_ITER = 1000;
344 :
345 64 : const double sign = std::signbit(x) ? -1 : 1;
346 :
347 64 : double g(std::abs(x));
348 64 : inttype a(0);
349 64 : inttype b(1);
350 64 : inttype c(1);
351 64 : inttype d(0);
352 :
353 : Fraction<std::uint64_t> ret;
354 :
355 118 : for (int i = 0; i < MAX_ITER; i++)
356 : {
357 236 : if (!(g >= 0 &&
358 118 : g <= static_cast<double>(std::numeric_limits<inttype>::max())))
359 : {
360 1 : break;
361 : }
362 117 : const inttype s = static_cast<inttype>(std::floor(g));
363 117 : ret.num = a + s * c;
364 117 : ret.denom = b + s * d;
365 :
366 117 : a = c;
367 117 : b = d;
368 117 : c = ret.num;
369 117 : d = ret.denom;
370 117 : g = 1.0 / (g - static_cast<double>(s));
371 :
372 117 : const double approx = sign * static_cast<double>(ret.num) /
373 117 : static_cast<double>(ret.denom);
374 :
375 117 : if (std::abs(approx - x) < err)
376 : {
377 63 : return ret;
378 : }
379 : }
380 :
381 1 : CPLError(CE_Warning, CPLE_AppDefined,
382 : "Failed to approximate %g as a fraction with error < %g in %d "
383 : "iterations",
384 : x, err, MAX_ITER);
385 1 : return std::nullopt;
386 : }
387 : } // namespace
388 :
389 : /** Return the largest value by which two input values can be
390 : * divided, with the result being an integer. If no suitable
391 : * value can be found, zero will be returned.
392 : */
393 47 : double CPLGreatestCommonDivisor(double a, double b)
394 : {
395 47 : if (a == 0 || !std::isfinite(a) || b == 0 || !std::isfinite(b))
396 : {
397 8 : CPLError(CE_Failure, CPLE_AppDefined,
398 : "Input values must be finite non-null values");
399 8 : return 0;
400 : }
401 :
402 39 : if (a == b)
403 : {
404 0 : return a;
405 : }
406 :
407 : // Check if one resolution is an integer factor of the other.
408 : // This is fast and succeeds in some cases where the method below fails.
409 39 : if (a > b && std::abs(std::round(a / b) - a / b) < 1e-8)
410 : {
411 3 : return b;
412 : }
413 36 : if (b > a && std::abs(std::round(b / a) - b / a) < 1e-8)
414 : {
415 4 : return a;
416 : }
417 :
418 32 : const auto approx_a = FloatToFraction(a, 1e-10);
419 32 : if (!approx_a.has_value())
420 : {
421 0 : CPLError(CE_Failure, CPLE_AppDefined,
422 : "Could not approximate resolution %.18g as a fraction", a);
423 0 : return 0;
424 : }
425 :
426 32 : const auto approx_b = FloatToFraction(b, 1e-10);
427 32 : if (!approx_b.has_value())
428 : {
429 1 : CPLError(CE_Failure, CPLE_AppDefined,
430 : "Could not approximate resolution %.18g as a fraction", b);
431 1 : return 0;
432 : }
433 :
434 31 : const double sign = std::signbit(a) ? -1 : 1;
435 :
436 31 : const auto &frac_a = approx_a.value();
437 31 : const auto &frac_b = approx_b.value();
438 :
439 31 : const auto common_denom = std::lcm(frac_a.denom, frac_b.denom);
440 :
441 : const auto num_a = static_cast<std::uint64_t>(
442 31 : frac_a.num * std::round(common_denom / frac_a.denom));
443 : const auto num_b = static_cast<std::uint64_t>(
444 31 : frac_b.num * std::round(common_denom / frac_b.denom));
445 :
446 31 : const auto common_num = std::gcd(num_a, num_b);
447 :
448 : // coverity[divide_by_zero]
449 31 : const auto common = sign * static_cast<double>(common_num) /
450 31 : static_cast<double>(common_denom);
451 :
452 31 : const auto disaggregation_factor = std::max(a / common, b / common);
453 31 : if (disaggregation_factor > 10000)
454 : {
455 4 : CPLError(CE_Failure, CPLE_AppDefined,
456 : "Common resolution between %.18g and %.18g calculated at "
457 : "%.18g which "
458 : "would cause excessive disaggregation",
459 : a, b, common);
460 4 : return 0;
461 : }
462 :
463 27 : return common;
464 : }
|
