From: tsteven4 Date: Sat, 18 Aug 2018 18:50:07 +0000 (-0600) Subject: use continued fractions to create TIFF RATIONAL types. X-Git-Tag: archive/raspbian/1.10.0+ds-2+rpi1~1^2~115^2~1^2~92^2~6 X-Git-Url: https://dgit.raspbian.org/?a=commitdiff_plain;h=3c560d32ea5af9fd5fcbf30a6ff247155c18112d;p=gpsbabel.git use continued fractions to create TIFF RATIONAL types. some references were regenerated as the new algorithm can generate different values for the RATIONAL numerator and denominator, in this case 0.0 -> 0/INT32_MAX. --- diff --git a/exif.cc b/exif.cc index 2983e14bf..c36e99b61 100644 --- a/exif.cc +++ b/exif.cc @@ -48,6 +48,7 @@ #include // for qPrintable #include // for sort #include // for assert +#include // for DBL_EPSILON #include // for fabs, floor #include // for snprintf, SEEK_SET, sscanf #include // for int32_t, int16_t, uint16_t, uint32_t, uint8_t @@ -891,69 +892,67 @@ exif_waypt_from_exif_app(ExifApp* app) return wpt; } -static int exif_gcd(int ui, int vi) +template +class Rational { - int u = abs(ui); - int v = abs(vi); - - /* Modern Euclidean algorithum to find greatest commond divisor */ - /* See Knuth, Seminumerical Algorithms, pg. 320 */ - while (v != 0) { - int r = u % v; - u = v; - v = r; - } - return u; -} +public: + T num; + T den; -// TODO: This algorithm could be improved upon using continued fractions, -// i.e. the domain could be expanded and the accuracy improved. -// We could also achieve an increased domain and accuracy for RATIONAL -// types if we handled them separately. -static void -exif_dec2frac(double val, int32_t* num, int32_t* den) + Rational() = default; + Rational(T n, T d) : num{n}, den{d} {} +}; + +// TODO: we could achieve an increased domain and accuracy for TIFF RATIONAL +// types if we handled them separately (int32_t -> uint32_t, INT32_MAX -> UINT32_MAX). +Rational exif_dec2frac(double val, double tolerance = DBL_EPSILON) { - char sval[16], snum[16]; - char dot = 0; - int den1 = 1; - - assert(val >= 0.0); - if (val < 0.000000001) { - val = 0.0; - } else if (val > 999999999.0) { + constexpr double upper_limit = INT32_MAX; + constexpr double lower_limit = 1.0/upper_limit; + const double pval = fabs(val); + const double tol = fmax(tolerance, DBL_EPSILON); + + if (pval < lower_limit) { + return Rational(0, upper_limit); + } else if (pval > upper_limit) { fatal(MYNAME ": Value (%f) to big for a rational representation!\n", val); + return Rational(copysign(upper_limit, val), 1); } - int num1 = 0; - double vx = fabs(val); - while (vx > 1) { - num1++; - vx = vx / 10; - } - - snprintf(sval, sizeof(sval), "%*.*f", 9, 9 - num1, fabs(val)); - snum[0] = '\0'; + double b; + double remainder = modf(pval, &b); + Rational prev_prev(1.0, 0.0); + Rational prev(b, 1.0); + Rational curr = prev; + + // phi = (1.0+sqrt(5.0))/2.0 is badly approximable and the slowest to converge. + // This is a good test case to see the maximum number of iterations required. + for (int idx = 0; idx < 64; ++idx) { + // Calculate the next simple continued fraction coefficient (b), and remainder (remainder). + if (remainder < lower_limit) { + break; // remainder is nearly zero, use current estimate. + } + remainder = modf(1.0/remainder, &b); - char* cx = sval; - while (*cx) { - if (dot) { - den1 *= 10; + // Convert the truncated simple continued fraction, a.k.a. a convergent, to an ordinary fraction (curr.num/curr.den). + Rational candidate((b * prev.num) + prev_prev.num, (b * prev.den) + prev_prev.den); + if (candidate.num > upper_limit) { + break; // numerator too big, use current estimate. } - if (*cx == '.') { - dot = 1; - } else { - strncat(snum, cx, 1); + if (candidate.den > upper_limit) { + break; // denominator too big, use current estimate. } - cx++; - } + curr = candidate; - num1 = atoi(snum); + if (fabs(pval- (curr.num/curr.den)) < (pval * tol)) { + break; // close enough, use current estimate. + } - int gcd = exif_gcd(num1, den1); - assert(gcd != 0); // Note gcd(0, 0) = 0, but we shouldn't generate num1 = den1 = 0. + prev_prev = prev; + prev = curr; + } - *num = num1 / gcd; - *den = den1 / gcd; + return Rational(round(copysign(curr.num, val)), round(curr.den)); } static ExifTag* @@ -1008,24 +1007,15 @@ exif_put_value(const int ifd_nr, const uint16_t tag_id, const uint16_t type, con case EXIF_TYPE_RAT: case EXIF_TYPE_SRAT: { double val = *(double*)data; - uint32_t* dest = reinterpret_cast(tag->data.data()); + int32_t* dest = reinterpret_cast(tag->data.data()); if ((val < 0.0) && (type == EXIF_TYPE_RAT)) { fatal(MYNAME ": A negative value cannot be stored as type RATIONAL."); } - if ((int)val == val) { - // For integers this expands the domain compared to the sub-optimal exif_dec2frac implementation. - dest[index * 2] = (int)val; - dest[(index * 2) + 1] = 1; - } else { - int32_t Nom, Den; - exif_dec2frac(fabs(val), &Nom, &Den); - if (val < 0.0) { - Nom *= -1; - } - dest[index * 2] = Nom; - dest[(index * 2) + 1] = Den; - } + + Rational rat = exif_dec2frac(val, 1e-11); + dest[index * 2] = rat.num; + dest[(index * 2) + 1] = rat.den; } break; default: { @@ -1060,7 +1050,7 @@ exif_put_coord(const int ifd_nr, const int tag_id, const double val) double vmin; double fractional_part; - fractional_part = modf(val, &vdeg); + fractional_part = modf(val, &vdeg); fractional_part = modf(60.0 * fractional_part, &vmin); double vsec = 60.0 * fractional_part; diff --git a/reference/20180717_080125.jpg.jpg b/reference/20180717_080125.jpg.jpg index ccd67ca9b..e13af24ee 100644 Binary files a/reference/20180717_080125.jpg.jpg and b/reference/20180717_080125.jpg.jpg differ diff --git a/reference/IMG_2065.JPG.jpg b/reference/IMG_2065.JPG.jpg index 5a1318e59..f6135bd43 100644 Binary files a/reference/IMG_2065.JPG.jpg and b/reference/IMG_2065.JPG.jpg differ diff --git a/reference/ricoh-rdc5300.jpg.jpg b/reference/ricoh-rdc5300.jpg.jpg index daf892350..9aa1027cd 100644 Binary files a/reference/ricoh-rdc5300.jpg.jpg and b/reference/ricoh-rdc5300.jpg.jpg differ