From a40ac65bf79a7911aa211907af99392793cb9262 Mon Sep 17 00:00:00 2001 From: David Tschumperle Date: Wed, 15 Nov 2023 15:37:02 +0100 Subject: [PATCH 1/4] Reimplement 'cimg::_hypoth()' with high precision. --- CImg.h | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/CImg.h b/CImg.h index 1928f797..ffb8d9c9 100644 --- a/CImg.h +++ b/CImg.h @@ -6844,16 +6844,28 @@ namespace cimg_library { return std::sqrt(x*x + y*y); } + //! Return sqrt(x^2 + y^2 + z^2). template inline T hypot(const T x, const T y, const T z) { return std::sqrt(x*x + y*y + z*z); } + //! Return sqrt(x^2 + y^2) (better precision). template - inline T _hypot(const T x, const T y) { // Slower but more precise version - T nx = cimg::abs(x), ny = cimg::abs(y), t; - if (nx0) { t/=nx; return nx*std::sqrt(1 + t*t); } + inline T _hypot(const T x, const T y) { + T nx = cimg::abs(x), ny = cimg::abs(y); + if (nx>ny) cimg::swap(nx,ny); + if (nx>0) return nx*std::sqrt(1 + cimg::sqr(ny/nx)); + return 0; + } + + //! Return sqrt(x^2 + y^2 + z^2) (better precision). + template + inline T _hypot(const T x, const T y, const T z) { + T nx = cimg::abs(x), ny = cimg::abs(y), nz = cimg::abs(z); + if (nx>ny) cimg::swap(nx,ny); + if (nx>nz) cimg::swap(nx,nz); + if (nx>0) return nx*std::sqrt(1 + cimg::sqr(ny/nx) + cimg::sqr(nz/nx)); return 0; } @@ -32298,7 +32310,7 @@ namespace cimg_library { (T)(2*X*W + 2*Y*Z),(T)(X*X - Y*Y + Z*Z - W*W),(T)(2*Z*W - 2*X*Y), (T)(2*Y*W - 2*X*Z),(T)(2*X*Y + 2*Z*W),(T)(X*X - Y*Y - Z*Z + W*W)); } - N = cimg::hypot((double)x,(double)y,(double)z); + N = cimg::_hypot((double)x,(double)y,(double)z); if (N>0) { X = x/N; Y = y/N; Z = z/N; } else { X = Y = 0; Z = 1; } const double ang = w*cimg::PI/180, c = std::cos(ang), omc = 1 - c, s = std::sin(ang); @@ -44532,7 +44544,7 @@ namespace cimg_library { const float patch_penalization, const bool allow_identity, const float max_score) { // 2D version - if (!allow_identity && cimg::hypot((float)x1-x2,(float)y1-y2)::inf(); const T *p1 = img1.data(x1*psizec,y1), *p2 = img2.data(x2*psizec,y2); const unsigned int psizewc = psizew*psizec; From cf8a2452871a4afd2aead106fe6ed2ead76432db Mon Sep 17 00:00:00 2001 From: David Tschumperle Date: Wed, 15 Nov 2023 16:17:08 +0100 Subject: [PATCH 2/4] _cimg_math_parser(): Add function 'hypot()'. --- CImg.h | 80 +++++++++++++++++++++++++++++++++------------------------- 1 file changed, 46 insertions(+), 34 deletions(-) diff --git a/CImg.h b/CImg.h index ffb8d9c9..25bd3007 100644 --- a/CImg.h +++ b/CImg.h @@ -21513,40 +21513,6 @@ namespace cimg_library { p1 = _cimg_mp_size(arg1); _cimg_mp_scalar3(mp_vector_normp,arg1,p1,arg2); } - - if (!std::strncmp(ss,"norm",4) && ss4::vector(0,0,0,arg1).move_to(l_opcode); - for (++s; s::sequence(_cimg_mp_size(arg2),arg2 + 1,arg2 + (ulongT)_cimg_mp_size(arg2)). - move_to(l_opcode); - else CImg::vector(arg2).move_to(l_opcode); - is_sth&=_cimg_mp_is_const_scalar(arg2); - s = ns; - } - (l_opcode>'y').move_to(opcode); - op = val==2?_mp_vector_norm2:val==1?_mp_vector_norm1:!val?_mp_vector_norm0: - cimg::type::is_inf(val)?_mp_vector_norminf:_mp_vector_normp; - opcode[0] = (ulongT)op; - opcode[2] = opcode._height; - if (is_sth) _cimg_mp_const_scalar(op(*this)); - if (opcode._height==5) { // Single argument - if (arg1) { _cimg_mp_scalar1(mp_abs,opcode[4]); } - else { _cimg_mp_scalar2(mp_neq,opcode[4],0); } - } - opcode[1] = pos = scalar(); - opcode.move_to(code); - return_new_comp = true; - _cimg_mp_return(pos); - } break; case 'o' : @@ -22767,6 +22733,43 @@ namespace cimg_library { break; } + if ((!std::strncmp(ss,"norm",4) && ss4::vector(0,0,0,arg1).move_to(l_opcode); + for (++s; s::sequence(_cimg_mp_size(arg2),arg2 + 1,arg2 + (ulongT)_cimg_mp_size(arg2)). + move_to(l_opcode); + else CImg::vector(arg2).move_to(l_opcode); + is_sth&=_cimg_mp_is_const_scalar(arg2); + s = ns; + } + (l_opcode>'y').move_to(opcode); + op = val==2?(is_hypot && opcode._height<8?_mp_vector_hypot:_mp_vector_norm2): + val==1?_mp_vector_norm1:!val?_mp_vector_norm0: + cimg::type::is_inf(val)?_mp_vector_norminf:_mp_vector_normp; + opcode[0] = (ulongT)op; + opcode[2] = opcode._height; + if (is_sth) _cimg_mp_const_scalar(op(*this)); + if (opcode._height==5) { // Single argument + if (arg1) { _cimg_mp_scalar1(mp_abs,opcode[4]); } + else { _cimg_mp_scalar2(mp_neq,opcode[4],0); } + } + opcode[1] = pos = scalar(); + opcode.move_to(code); + return_new_comp = true; + _cimg_mp_return(pos); + } + if (!std::strncmp(ss,"max(",4) || !std::strncmp(ss,"min(",4) || !std::strncmp(ss,"maxabs(",7) || !std::strncmp(ss,"minabs(",7) || !std::strncmp(ss,"med(",4) || !std::strncmp(ss,"kth(",4) || @@ -28282,6 +28285,15 @@ namespace cimg_library { return p?cimg::abs(val):(val!=0); } + static double _mp_vector_hypot(_cimg_math_parser& mp) { + switch ((unsigned int)mp.opcode[2]) { + case 5 : return cimg::abs(_mp_arg(4)); + case 6 : return cimg::_hypot(_mp_arg(4),_mp_arg(5)); + case 7 : return cimg::_hypot(_mp_arg(4),_mp_arg(5),_mp_arg(6)); + }; + return _mp_vector_norm2(mp); + } + static double mp_vector_print(_cimg_math_parser& mp) { const bool print_string = (bool)mp.opcode[4]; cimg_pragma_openmp(critical(mp_vector_print)) From 9e30cea109c970b7b6bbb3530eafc48fc801cbfc Mon Sep 17 00:00:00 2001 From: David Tschumperle Date: Fri, 17 Nov 2023 14:13:47 +0100 Subject: [PATCH 3/4] Fix 'cimg::_hypot()'. --- CImg.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CImg.h b/CImg.h index 25bd3007..3c5b2a60 100644 --- a/CImg.h +++ b/CImg.h @@ -6854,7 +6854,7 @@ namespace cimg_library { template inline T _hypot(const T x, const T y) { T nx = cimg::abs(x), ny = cimg::abs(y); - if (nx>ny) cimg::swap(nx,ny); + if (nx0) return nx*std::sqrt(1 + cimg::sqr(ny/nx)); return 0; } @@ -6863,8 +6863,8 @@ namespace cimg_library { template inline T _hypot(const T x, const T y, const T z) { T nx = cimg::abs(x), ny = cimg::abs(y), nz = cimg::abs(z); - if (nx>ny) cimg::swap(nx,ny); - if (nx>nz) cimg::swap(nx,nz); + if (nx0) return nx*std::sqrt(1 + cimg::sqr(ny/nx) + cimg::sqr(nz/nx)); return 0; } From bde802732c5872765b0f8408c806fc8918d44179 Mon Sep 17 00:00:00 2001 From: David Tschumperle Date: Fri, 17 Nov 2023 14:21:53 +0100 Subject: [PATCH 4/4] . --- CImg.h | 41 +++++++++++++---------------------------- 1 file changed, 13 insertions(+), 28 deletions(-) diff --git a/CImg.h b/CImg.h index 3c5b2a60..7e7ac873 100644 --- a/CImg.h +++ b/CImg.h @@ -6841,7 +6841,11 @@ namespace cimg_library { //! Return sqrt(x^2 + y^2). template inline T hypot(const T x, const T y) { +#if cimg_use_cpp11==1 && !defined(_MSC_VER) + return std::hypot(x,y); +#else return std::sqrt(x*x + y*y); +#endif } //! Return sqrt(x^2 + y^2 + z^2). @@ -6850,25 +6854,6 @@ namespace cimg_library { return std::sqrt(x*x + y*y + z*z); } - //! Return sqrt(x^2 + y^2) (better precision). - template - inline T _hypot(const T x, const T y) { - T nx = cimg::abs(x), ny = cimg::abs(y); - if (nx0) return nx*std::sqrt(1 + cimg::sqr(ny/nx)); - return 0; - } - - //! Return sqrt(x^2 + y^2 + z^2) (better precision). - template - inline T _hypot(const T x, const T y, const T z) { - T nx = cimg::abs(x), ny = cimg::abs(y), nz = cimg::abs(z); - if (nx0) return nx*std::sqrt(1 + cimg::sqr(ny/nx) + cimg::sqr(nz/nx)); - return 0; - } - //! Return the factorial of n inline double factorial(const int n) { if (n<0) return cimg::type::nan(); @@ -24273,7 +24258,7 @@ namespace cimg_library { } static double mp_complex_abs(_cimg_math_parser& mp) { - return cimg::_hypot(_mp_arg(2),_mp_arg(3)); + return cimg::hypot(_mp_arg(2),_mp_arg(3)); } static double mp_complex_conj(_cimg_math_parser& mp) { @@ -24413,7 +24398,7 @@ namespace cimg_library { static double mp_complex_sqrt(_cimg_math_parser& mp) { const double real = _mp_arg(2), imag = _mp_arg(3), - r = std::sqrt(cimg::_hypot(real,imag)), + r = std::sqrt(cimg::hypot(real,imag)), theta = std::atan2(imag,real)/2; double *ptrd = &_mp_arg(1) + 1; ptrd[0] = r*std::cos(theta); @@ -28288,8 +28273,8 @@ namespace cimg_library { static double _mp_vector_hypot(_cimg_math_parser& mp) { switch ((unsigned int)mp.opcode[2]) { case 5 : return cimg::abs(_mp_arg(4)); - case 6 : return cimg::_hypot(_mp_arg(4),_mp_arg(5)); - case 7 : return cimg::_hypot(_mp_arg(4),_mp_arg(5),_mp_arg(6)); + case 6 : return cimg::hypot(_mp_arg(4),_mp_arg(5)); + case 7 : return cimg::hypot(_mp_arg(4),_mp_arg(5),_mp_arg(6)); }; return _mp_vector_norm2(mp); } @@ -31379,7 +31364,7 @@ namespace cimg_library { rv1[i] = c*rv1[i]; if ((cimg::abs(f) + anorm)==anorm) break; g = S[i]; - h = cimg::_hypot(f,g); + h = cimg::hypot(f,g); S[i] = h; h = 1/h; c = g*h; @@ -31399,7 +31384,7 @@ namespace cimg_library { g = rv1[nm]; h = rv1[k]; f = ((y - z)*(y + z) + (g - h)*(g + h))/std::max(epsilon,(Ttfloat)2*h*y); - g = cimg::_hypot(f,(Ttfloat)1); + g = cimg::hypot(f,(Ttfloat)1); f = ((x - z)*(x + z) + h*((y/(f + (f>=0?g:-g))) - h))/std::max(epsilon,(Ttfloat)x); c = s = 1; for (int j = l; j<=nm; ++j) { @@ -31407,7 +31392,7 @@ namespace cimg_library { g = rv1[i]; h = s*g; g = c*g; - t y1 = S[i], z1 = cimg::_hypot(f,h); + t y1 = S[i], z1 = cimg::hypot(f,h); rv1[j] = z1; c = f/std::max(epsilon,(Ttfloat)z1); s = h/std::max(epsilon,(Ttfloat)z1); @@ -31420,7 +31405,7 @@ namespace cimg_library { V(j,jj) = x2*c + z2*s; V(i,jj) = z2*c - x2*s; } - z1 = cimg::_hypot(f,h); + z1 = cimg::hypot(f,h); S[j] = z1; if (z1) { z1 = 1/std::max(epsilon,(Ttfloat)z1); @@ -32322,7 +32307,7 @@ namespace cimg_library { (T)(2*X*W + 2*Y*Z),(T)(X*X - Y*Y + Z*Z - W*W),(T)(2*Z*W - 2*X*Y), (T)(2*Y*W - 2*X*Z),(T)(2*X*Y + 2*Z*W),(T)(X*X - Y*Y - Z*Z + W*W)); } - N = cimg::_hypot((double)x,(double)y,(double)z); + N = cimg::hypot((double)x,(double)y,(double)z); if (N>0) { X = x/N; Y = y/N; Z = z/N; } else { X = Y = 0; Z = 1; } const double ang = w*cimg::PI/180, c = std::cos(ang), omc = 1 - c, s = std::sin(ang);