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 <cstring>
51 : #include <numeric>
52 : #include <optional>
53 :
54 : /************************************************************************/
55 : /* HalfToFloat() */
56 : /* */
57 : /* 16-bit floating point number to 32-bit one. */
58 : /************************************************************************/
59 :
60 162020 : GUInt32 CPLHalfToFloat(GUInt16 iHalf)
61 : {
62 :
63 162020 : GUInt32 iSign = (iHalf >> 15) & 0x00000001;
64 162020 : int iExponent = (iHalf >> 10) & 0x0000001f;
65 162020 : GUInt32 iMantissa = iHalf & 0x000003ff;
66 :
67 162020 : if (iExponent == 0)
68 : {
69 261 : if (iMantissa == 0)
70 : {
71 : /* --------------------------------------------------------------------
72 : */
73 : /* Plus or minus zero. */
74 : /* --------------------------------------------------------------------
75 : */
76 :
77 256 : return iSign << 31;
78 : }
79 : else
80 : {
81 : /* --------------------------------------------------------------------
82 : */
83 : /* Denormalized number -- renormalize it. */
84 : /* --------------------------------------------------------------------
85 : */
86 :
87 37 : while (!(iMantissa & 0x00000400))
88 : {
89 32 : iMantissa <<= 1;
90 32 : iExponent -= 1;
91 : }
92 :
93 5 : iExponent += 1;
94 5 : iMantissa &= ~0x00000400U;
95 : }
96 : }
97 161759 : else if (iExponent == 31)
98 : {
99 11 : if (iMantissa == 0)
100 : {
101 : /* --------------------------------------------------------------------
102 : */
103 : /* Positive or negative infinity. */
104 : /* --------------------------------------------------------------------
105 : */
106 :
107 5 : return (iSign << 31) | 0x7f800000;
108 : }
109 : else
110 : {
111 : /* --------------------------------------------------------------------
112 : */
113 : /* NaN -- preserve sign and significand bits. */
114 : /* --------------------------------------------------------------------
115 : */
116 :
117 6 : return (iSign << 31) | 0x7f800000 | (iMantissa << 13);
118 : }
119 : }
120 :
121 : /* -------------------------------------------------------------------- */
122 : /* Normalized number. */
123 : /* -------------------------------------------------------------------- */
124 :
125 161753 : iExponent = iExponent + (127 - 15);
126 161753 : iMantissa = iMantissa << 13;
127 :
128 : /* -------------------------------------------------------------------- */
129 : /* Assemble sign, exponent and mantissa. */
130 : /* -------------------------------------------------------------------- */
131 :
132 : /* coverity[overflow_sink] */
133 161753 : return (iSign << 31) | (static_cast<GUInt32>(iExponent) << 23) | iMantissa;
134 : }
135 :
136 : /************************************************************************/
137 : /* TripleToFloat() */
138 : /* */
139 : /* 24-bit floating point number to 32-bit one. */
140 : /************************************************************************/
141 :
142 400 : GUInt32 CPLTripleToFloat(GUInt32 iTriple)
143 : {
144 :
145 400 : GUInt32 iSign = (iTriple >> 23) & 0x00000001;
146 400 : int iExponent = (iTriple >> 16) & 0x0000007f;
147 400 : GUInt32 iMantissa = iTriple & 0x0000ffff;
148 :
149 400 : if (iExponent == 0)
150 : {
151 0 : if (iMantissa == 0)
152 : {
153 : /* --------------------------------------------------------------------
154 : */
155 : /* Plus or minus zero. */
156 : /* --------------------------------------------------------------------
157 : */
158 :
159 0 : return iSign << 31;
160 : }
161 : else
162 : {
163 : /* --------------------------------------------------------------------
164 : */
165 : /* Denormalized number -- renormalize it. */
166 : /* --------------------------------------------------------------------
167 : */
168 :
169 0 : while (!(iMantissa & 0x00010000))
170 : {
171 0 : iMantissa <<= 1;
172 0 : iExponent -= 1;
173 : }
174 :
175 0 : iExponent += 1;
176 0 : iMantissa &= ~0x00010000U;
177 : }
178 : }
179 400 : else if (iExponent == 127)
180 : {
181 0 : if (iMantissa == 0)
182 : {
183 : /* --------------------------------------------------------------------
184 : */
185 : /* Positive or negative infinity. */
186 : /* --------------------------------------------------------------------
187 : */
188 :
189 0 : return (iSign << 31) | 0x7f800000;
190 : }
191 : else
192 : {
193 : /* --------------------------------------------------------------------
194 : */
195 : /* NaN -- preserve sign and significand bits. */
196 : /* --------------------------------------------------------------------
197 : */
198 :
199 0 : return (iSign << 31) | 0x7f800000 | (iMantissa << 7);
200 : }
201 : }
202 :
203 : /* -------------------------------------------------------------------- */
204 : /* Normalized number. */
205 : /* -------------------------------------------------------------------- */
206 :
207 400 : iExponent = iExponent + (127 - 63);
208 400 : iMantissa = iMantissa << 7;
209 :
210 : /* -------------------------------------------------------------------- */
211 : /* Assemble sign, exponent and mantissa. */
212 : /* -------------------------------------------------------------------- */
213 :
214 : /* coverity[overflow_sink] */
215 400 : return (iSign << 31) | (static_cast<GUInt32>(iExponent) << 23) | iMantissa;
216 : }
217 :
218 : /************************************************************************/
219 : /* FloatToHalf() */
220 : /************************************************************************/
221 :
222 240432 : GUInt16 CPLFloatToHalf(GUInt32 iFloat32, bool &bHasWarned)
223 : {
224 240432 : GUInt32 iSign = (iFloat32 >> 31) & 0x00000001;
225 240432 : GUInt32 iExponent = (iFloat32 >> 23) & 0x000000ff;
226 240432 : GUInt32 iMantissa = iFloat32 & 0x007fffff;
227 :
228 240432 : if (iExponent == 255)
229 : {
230 10 : if (iMantissa == 0)
231 : {
232 : /* --------------------------------------------------------------------
233 : */
234 : /* Positive or negative infinity. */
235 : /* --------------------------------------------------------------------
236 : */
237 :
238 4 : return static_cast<GUInt16>((iSign << 15) | 0x7C00);
239 : }
240 : else
241 : {
242 : /* --------------------------------------------------------------------
243 : */
244 : /* NaN -- preserve sign and significand bits. */
245 : /* --------------------------------------------------------------------
246 : */
247 6 : if (iMantissa >> 13)
248 4 : return static_cast<GUInt16>((iSign << 15) | 0x7C00 |
249 4 : (iMantissa >> 13));
250 :
251 2 : return static_cast<GUInt16>((iSign << 15) | 0x7E00);
252 : }
253 : }
254 :
255 240422 : if (iExponent <= 127 - 15)
256 : {
257 : // Zero, float32 denormalized number or float32 too small normalized
258 : // number
259 302 : if (13 + 1 + 127 - 15 - iExponent >= 32)
260 297 : return static_cast<GUInt16>(iSign << 15);
261 :
262 : // Return a denormalized number
263 : return static_cast<GUInt16>(
264 5 : (iSign << 15) |
265 5 : ((iMantissa | 0x00800000) >> (13 + 1 + 127 - 15 - iExponent)));
266 : }
267 240120 : if (iExponent - (127 - 15) >= 31)
268 : {
269 1 : if (!bHasWarned)
270 : {
271 1 : bHasWarned = true;
272 1 : float fVal = 0.0f;
273 1 : memcpy(&fVal, &iFloat32, 4);
274 1 : CPLError(
275 : CE_Failure, CPLE_AppDefined,
276 : "Value %.8g is beyond range of float16. Converted to %sinf",
277 1 : fVal, (fVal > 0) ? "+" : "-");
278 : }
279 1 : return static_cast<GUInt16>((iSign << 15) | 0x7C00); // Infinity
280 : }
281 :
282 : /* -------------------------------------------------------------------- */
283 : /* Normalized number. */
284 : /* -------------------------------------------------------------------- */
285 :
286 240119 : iExponent = iExponent - (127 - 15);
287 240119 : iMantissa = iMantissa >> 13;
288 :
289 : /* -------------------------------------------------------------------- */
290 : /* Assemble sign, exponent and mantissa. */
291 : /* -------------------------------------------------------------------- */
292 :
293 : // coverity[overflow_sink]
294 240119 : return static_cast<GUInt16>((iSign << 15) | (iExponent << 10) | iMantissa);
295 : }
296 :
297 0 : GUInt16 CPLConvertFloatToHalf(float fFloat32)
298 : {
299 : GUInt32 nFloat32;
300 0 : std::memcpy(&nFloat32, &fFloat32, sizeof nFloat32);
301 0 : bool bHasWarned = true;
302 0 : return CPLFloatToHalf(nFloat32, bHasWarned);
303 : }
304 :
305 0 : float CPLConvertHalfToFloat(GUInt16 nHalf)
306 : {
307 0 : GUInt32 nFloat32 = CPLHalfToFloat(nHalf);
308 : float fFloat32;
309 0 : std::memcpy(&fFloat32, &nFloat32, sizeof fFloat32);
310 0 : return fFloat32;
311 : }
312 :
313 : namespace
314 : {
315 :
316 : template <typename T> struct Fraction
317 : {
318 : using value_type = T;
319 :
320 : T num;
321 : T denom;
322 : };
323 :
324 : /** Approximate a floating point number as a fraction, using the method describe
325 : * in Richards, Ian (1981). Continued Fractions Without Tears. Mathematics
326 : * Magazine, Vol. 54, No. 4. https://doi.org/10.2307/2689627
327 : *
328 : * If the fraction cannot be approximated within the specified error tolerance
329 : * in a certain amount of iterations, a warning will be raised and std::nullopt
330 : * will be returned.
331 : *
332 : * @param x the number to approximate as a fraction
333 : * @param err the maximum allowable absolute error in the approximation
334 : *
335 : * @return the approximated value, or std::nullopt
336 : *
337 : */
338 58 : std::optional<Fraction<std::uint64_t>> FloatToFraction(double x, double err)
339 : {
340 : using inttype = std::uint64_t;
341 58 : int nMaxIter = 1000;
342 :
343 58 : double sign = std::signbit(x) ? -1 : 1;
344 :
345 58 : double g(std::abs(x));
346 58 : inttype a(0);
347 58 : inttype b(1);
348 58 : inttype c(1);
349 58 : inttype d(0);
350 :
351 : Fraction<std::uint64_t> ret;
352 :
353 112 : for (int i = 0; i < nMaxIter; i++)
354 : {
355 112 : inttype s = static_cast<inttype>(std::floor(g));
356 112 : ret.num = a + s * c;
357 112 : ret.denom = b + s * d;
358 :
359 112 : a = c;
360 112 : b = d;
361 112 : c = ret.num;
362 112 : d = ret.denom;
363 112 : g = 1.0 / (g - static_cast<double>(s));
364 :
365 112 : double approx = sign * static_cast<double>(ret.num) /
366 112 : static_cast<double>(ret.denom);
367 :
368 112 : if (std::abs(approx - x) < err)
369 : {
370 58 : return ret;
371 : }
372 : }
373 :
374 0 : CPLError(CE_Warning, CPLE_AppDefined,
375 : "Failed to approximate %g as a fraction with error < %g in %d "
376 : "iterations",
377 : x, err, nMaxIter);
378 0 : return std::nullopt;
379 : }
380 : } // namespace
381 :
382 : /** Return the largest value by which two input values can be
383 : * divided, with the result being an integer. If no suitable
384 : * value can be found, zero will be returned.
385 : */
386 35 : double CPLGreatestCommonDivisor(double a, double b)
387 : {
388 35 : if (a == b)
389 : {
390 0 : return a;
391 : }
392 :
393 : // Check if one resolution is an integer factor of the other.
394 : // This is fast and succeeds in some cases where the method below fails.
395 35 : if (a > b && std::abs(std::round(a / b) - a / b) < 1e-8)
396 : {
397 3 : return b;
398 : }
399 32 : if (b > a && std::abs(std::round(b / a) - b / a) < 1e-8)
400 : {
401 3 : return a;
402 : }
403 :
404 29 : auto approx_a = FloatToFraction(a, 1e-10);
405 29 : if (!approx_a.has_value())
406 : {
407 0 : CPLError(CE_Failure, CPLE_AppDefined,
408 : "Could not approximate resolution %.18g as a fraction", a);
409 0 : return 0;
410 : }
411 :
412 29 : auto approx_b = FloatToFraction(b, 1e-10);
413 29 : if (!approx_b.has_value())
414 : {
415 0 : CPLError(CE_Failure, CPLE_AppDefined,
416 : "Could not approximate resolution %.18g as a fraction", b);
417 0 : return 0;
418 : }
419 :
420 29 : double sign = std::signbit(a) ? -1 : 1;
421 :
422 29 : auto &frac_a = approx_a.value();
423 29 : auto &frac_b = approx_b.value();
424 :
425 29 : auto common_denom = std::lcm(frac_a.denom, frac_b.denom);
426 :
427 : auto num_a = static_cast<std::uint64_t>(
428 29 : frac_a.num * std::round(common_denom / frac_a.denom));
429 : auto num_b = static_cast<std::uint64_t>(
430 29 : frac_b.num * std::round(common_denom / frac_b.denom));
431 :
432 29 : auto common_num = std::gcd(num_a, num_b);
433 :
434 29 : auto common = sign * static_cast<double>(common_num) /
435 29 : static_cast<double>(common_denom);
436 :
437 29 : auto disaggregation_factor = std::max(a / common, b / common);
438 29 : if (disaggregation_factor > 10000)
439 : {
440 4 : CPLError(CE_Failure, CPLE_AppDefined,
441 : "Common resolution between %.18g and %.18g calculated at "
442 : "%.18g which "
443 : "would cause excessive disaggregation",
444 : a, b, common);
445 4 : return 0;
446 : }
447 :
448 25 : return common;
449 : }
|
