From fc5ac5bce7f4aa71d96b1696bf1632e32007ac0e Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 12:27:43 +0100
Subject: [PATCH 01/24] math_parser(): Add new function 'warp()' which allows
 to get a warped version a vector seen as an image.

---
 CImg.h | 80 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 80 insertions(+)

diff --git a/CImg.h b/CImg.h
index 7aaba370..079fc4fa 100644
--- a/CImg.h
+++ b/CImg.h
@@ -22693,6 +22693,63 @@ namespace cimg_library {
               _cimg_mp_scalar1(mp_image_whds,p1);
             }
 
+            if (!std::strncmp(ss,"warp(",5)) { // Image warping
+              _cimg_mp_op("Function 'warp()'");
+              pos = 1;
+              for (s = ss5; s<se; ++s, ++pos) {
+                ns = s; while (ns<se && (*ns!=',' || level[ns - expr._data]!=clevel1) &&
+                               (*ns!=')' || level[ns - expr._data]!=clevel)) ++ns;
+                arg2 = compile(s,ns,depth1,0,block_flags);
+                if (pos==1 || pos==6) _cimg_mp_check_type(arg2,pos,2,0);
+                else if (pos<11) _cimg_mp_check_const_scalar(arg2,pos,3);
+                else _cimg_mp_check_type(arg2,pos,1,0);
+                CImg<ulongT>::vector(arg2).move_to(l_opcode);
+                s = ns;
+              }
+              (l_opcode>'y').move_to(opcode);
+              // opcode = [ A, wA,hA,dA,sA, W, wW,hW,dW,sW, mode, interp, boundary_cond ]
+              //          [ 0   1  2  3  4  5   6  7  8  9    10      11             12 ]
+
+              if (opcode.height()<10) compile(s,se1,depth1,0,block_flags); // Not enough arguments -> throw exception
+              arg1 = (unsigned int)opcode[0]; // Image to warp
+              arg2 = (unsigned int)opcode[5]; // Warp map
+              p1 = _cimg_mp_size(arg1);
+              p2 = _cimg_mp_size(arg2);
+              p3 = opcode.height();
+              opcode.resize(1,13,1,1,0);
+              if (p3<11) opcode[10] = 0;
+              if (p3<12) opcode[11] = 1;
+              if (p3<13) opcode[12] = 0;
+              arg3 = (unsigned int)mem[opcode[1]]; opcode[1] = arg3;
+              arg4 = (unsigned int)mem[opcode[2]]; opcode[2] = arg4;
+              arg5 = (unsigned int)mem[opcode[3]]; opcode[3] = arg5;
+              arg6 = (unsigned int)mem[opcode[4]]; opcode[4] = arg6;
+              if (arg3*arg4*arg5*arg6!=std::max(1U,p1))
+                throw CImgArgumentException("[" cimg_appname "_math_parser] "
+                                            "CImg<%s>::%s: %s: Input vector size (%lu values) and its specified "
+                                            "geometry (%u,%u,%u,%u) (%lu values) do not match.",
+                                            pixel_type(),_cimg_mp_calling_function,s_op,
+                                            std::max(p1,1U),arg3,arg4,arg5,arg6,(ulongT)arg3*arg4*arg5*arg6);
+              arg3 = (unsigned int)mem[opcode[6]]; opcode[6] = arg3;
+              arg4 = (unsigned int)mem[opcode[7]]; opcode[7] = arg4;
+              arg5 = (unsigned int)mem[opcode[8]]; opcode[8] = arg5;
+              arg6 = (unsigned int)mem[opcode[9]]; opcode[9] = arg6;
+              if (arg3*arg4*arg5*arg6!=std::max(1U,p2))
+                throw CImgArgumentException("[" cimg_appname "_math_parser] "
+                                            "CImg<%s>::%s: %s: Warp vector size (%lu values) and its specified "
+                                            "geometry (%u,%u,%u,%u) (%lu values) do not match.",
+                                            pixel_type(),_cimg_mp_calling_function,s_op,
+                                            std::max(p2,1U),arg3,arg4,arg5,arg6,(ulongT)arg3*arg4*arg5*arg6);
+
+              pos = vector(arg3*arg4*arg5*(unsigned int)opcode[4]);
+              opcode.resize(1,15,1,1,0,0,0,1);
+              opcode[0] = (ulongT)mp_vector_warp;
+              opcode[1] = (ulongT)pos;
+              opcode.move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
             if (!std::strncmp(ss,"while(",6)) { // While...do
               _cimg_mp_op("Function 'while()'");
               s0 = *ss5=='('?ss6:ss8;
@@ -28418,6 +28475,29 @@ namespace cimg_library {
         return val?(_mp_arg(2)?1:val):0;
       }
 
+      static double mp_vector_warp(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int
+          wA = (unsigned int)mp.opcode[3],
+          hA = (unsigned int)mp.opcode[4],
+          dA = (unsigned int)mp.opcode[5],
+          sA = (unsigned int)mp.opcode[6],
+          wW = (unsigned int)mp.opcode[8],
+          hW = (unsigned int)mp.opcode[9],
+          dW = (unsigned int)mp.opcode[10],
+          sW = (unsigned int)mp.opcode[11];
+        const int
+          mode = (int)_mp_arg(12),
+          interpolation = (int)_mp_arg(13),
+          boundary_conditions = (int)_mp_arg(14);
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          *const ptrw = &_mp_arg(7) + 1;
+        CImg<doubleT>(ptrd,wW,hW,dW,sA,true) = CImg<doubleT>(ptrs,wA,hA,dA,sA,true).
+          get_warp(CImg<doubleT>(ptrw,wW,hW,dW,sW,true),mode,interpolation,boundary_conditions);
+        return cimg::type<double>::nan();
+      }
+
 #define _cimg_mp_vfunc(func) \
       const longT sizd = (longT)mp.opcode[2];\
       const unsigned int nbargs = (unsigned int)(mp.opcode[3] - 4)/2; \

From 7d4f269218250c6351979b605ece10ef892ac98b Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 13:22:17 +0100
Subject: [PATCH 02/24] .

---
 CImg.h | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/CImg.h b/CImg.h
index 079fc4fa..31c17f87 100644
--- a/CImg.h
+++ b/CImg.h
@@ -22707,7 +22707,7 @@ namespace cimg_library {
                 s = ns;
               }
               (l_opcode>'y').move_to(opcode);
-              // opcode = [ A, wA,hA,dA,sA, W, wW,hW,dW,sW, mode, interp, boundary_cond ]
+              // opcode = [ A, wA,hA,dA,sA, B, wB,hB,dB,sB, mode, interp, boundary_cond ]
               //          [ 0   1  2  3  4  5   6  7  8  9    10      11             12 ]
 
               if (opcode.height()<10) compile(s,se1,depth1,0,block_flags); // Not enough arguments -> throw exception
@@ -28482,10 +28482,10 @@ namespace cimg_library {
           hA = (unsigned int)mp.opcode[4],
           dA = (unsigned int)mp.opcode[5],
           sA = (unsigned int)mp.opcode[6],
-          wW = (unsigned int)mp.opcode[8],
-          hW = (unsigned int)mp.opcode[9],
-          dW = (unsigned int)mp.opcode[10],
-          sW = (unsigned int)mp.opcode[11];
+          wB = (unsigned int)mp.opcode[8],
+          hB = (unsigned int)mp.opcode[9],
+          dB = (unsigned int)mp.opcode[10],
+          sB = (unsigned int)mp.opcode[11];
         const int
           mode = (int)_mp_arg(12),
           interpolation = (int)_mp_arg(13),
@@ -28493,8 +28493,8 @@ namespace cimg_library {
         const double
           *const ptrs = &_mp_arg(2) + 1,
           *const ptrw = &_mp_arg(7) + 1;
-        CImg<doubleT>(ptrd,wW,hW,dW,sA,true) = CImg<doubleT>(ptrs,wA,hA,dA,sA,true).
-          get_warp(CImg<doubleT>(ptrw,wW,hW,dW,sW,true),mode,interpolation,boundary_conditions);
+        CImg<doubleT>(ptrd,wB,hB,dB,sA,true) = CImg<doubleT>(ptrs,wA,hA,dA,sA,true).
+          get_warp(CImg<doubleT>(ptrw,wB,hB,dB,sB,true),mode,interpolation,boundary_conditions);
         return cimg::type<double>::nan();
       }
 

From 2cd3ad913033a71c8a728481b081f0c0d310e025 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 13:48:57 +0100
Subject: [PATCH 03/24] .

---
 CImg.h | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/CImg.h b/CImg.h
index 31c17f87..017b45c2 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21529,10 +21529,10 @@ namespace cimg_library {
             break;
 
           case 'p' :
-            if (!std::strncmp(ss,"permut(",7)) { // Number of permutations
-              _cimg_mp_op("Function 'permut()'");
-              s1 = ss7; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
-              arg1 = compile(ss7,s1,depth1,0,block_flags);
+            if (!std::strncmp(ss,"perm(",5)) { // Number of permutations
+              _cimg_mp_op("Function 'perm()'");
+              s1 = ss5; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(ss5,s1,depth1,0,block_flags);
               s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
               arg2 = compile(s1,s2,depth1,0,block_flags);
               arg3 = compile(++s2,se1,depth1,0,block_flags);

From dcfd9d896baa136d41ff10847b82c622a14bb38c Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 13:51:46 +0100
Subject: [PATCH 04/24] .

---
 CImg.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CImg.h b/CImg.h
index 017b45c2..81dcbebc 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21535,7 +21535,7 @@ namespace cimg_library {
               arg1 = compile(ss5,s1,depth1,0,block_flags);
               s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
               arg2 = compile(s1,s2,depth1,0,block_flags);
-              arg3 = compile(++s2,se1,depth1,0,block_flags);
+              arg3 = s2<se1?compile(++s2,se1,depth1,0,block_flags):0;
               _cimg_mp_check_type(arg1,1,1,0);
               _cimg_mp_check_type(arg2,2,1,0);
               _cimg_mp_check_type(arg3,3,1,0);

From 35a036a3c9f7b15138793612d9c80a1a0800e131 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 14:15:37 +0100
Subject: [PATCH 05/24] cimg_math_parser(): New function 'permute()'.

---
 CImg.h | 71 +++++++++++++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 63 insertions(+), 8 deletions(-)

diff --git a/CImg.h b/CImg.h
index 81dcbebc..61486625 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21544,6 +21544,43 @@ namespace cimg_library {
               _cimg_mp_scalar3(mp_permutations,arg1,arg2,arg3);
             }
 
+            if (!std::strncmp(ss,"permute(",8)) { // Permute axes
+              _cimg_mp_op("Function 'permute()'");
+              s1 = ss8; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(ss8,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg2 = compile(s1,s2,depth1,0,block_flags);
+              s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg3 = compile(s2,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg4 = compile(s1,s2,depth1,0,block_flags);
+              s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg5 = compile(s2,s1,depth1,0,block_flags);
+              arg6 = compile(++s1,se1,depth1,0,block_flags);
+              _cimg_mp_check_type(arg1,1,2,0);
+              _cimg_mp_check_const_scalar(arg2,2,3);
+              _cimg_mp_check_const_scalar(arg3,3,3);
+              _cimg_mp_check_const_scalar(arg4,4,3);
+              _cimg_mp_check_const_scalar(arg5,5,3);
+              _cimg_mp_check_type(arg6,6,2,0);
+              p1 = _cimg_mp_size(arg1);
+              p2 = _cimg_mp_size(arg6);
+              arg2 = (unsigned int)mem[arg2];
+              arg3 = (unsigned int)mem[arg3];
+              arg4 = (unsigned int)mem[arg4];
+              arg5 = (unsigned int)mem[arg5];
+              if (arg2*arg3*arg4*arg5!=std::max(1U,p1))
+                throw CImgArgumentException("[" cimg_appname "_math_parser] "
+                                            "CImg<%s>::%s: %s: Input vector size (%lu values) and its specified "
+                                            "geometry (%u,%u,%u,%u) (%lu values) do not match.",
+                                            pixel_type(),_cimg_mp_calling_function,s_op,
+                                            std::max(p1,1U),arg2,arg3,arg4,arg5,(ulongT)arg2*arg3*arg4*arg5);
+              pos = vector(arg2*arg3*arg4*arg5);
+              CImg<ulongT>::vector((ulongT)mp_vector_permute,pos,arg1,arg2,arg3,arg4,arg5,arg6,p2).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
             if (!std::strncmp(ss,"polygon(",8)) { // Polygon/line drawing
               if (!is_inside_critical) is_parallelizable = false;
               _cimg_mp_op("Function 'polygon()'");
@@ -22740,7 +22777,6 @@ namespace cimg_library {
                                             "geometry (%u,%u,%u,%u) (%lu values) do not match.",
                                             pixel_type(),_cimg_mp_calling_function,s_op,
                                             std::max(p2,1U),arg3,arg4,arg5,arg6,(ulongT)arg3*arg4*arg5*arg6);
-
               pos = vector(arg3*arg4*arg5*(unsigned int)opcode[4]);
               opcode.resize(1,15,1,1,0,0,0,1);
               opcode[0] = (ulongT)mp_vector_warp;
@@ -28129,6 +28165,15 @@ namespace cimg_library {
         return cimg::type<double>::nan();
       }
 
+      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_init(_cimg_math_parser& mp) {
         unsigned int
           ptrs = 4U,
@@ -28345,13 +28390,23 @@ 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_permute(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int
+          wA = (unsigned int)mp.opcode[3],
+          hA = (unsigned int)mp.opcode[4],
+          dA = (unsigned int)mp.opcode[5],
+          sA = (unsigned int)mp.opcode[6],
+          sizp = (unsigned int)mp.opcode[8];
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          *const ptrp = &_mp_arg(7) + 1;
+        CImg<charT> str(sizp + 1);
+        for (unsigned int p = 0; p<sizp; ++p) str[p] = (char)ptrp[p];
+        str.back() = 0;
+        CImg<doubleT>(ptrd,wA,hA,dA,sA,true) = CImg<doubleT>(ptrs,wA,hA,dA,sA,true).
+          get_permute_axes(str);
+        return cimg::type<double>::nan();
       }
 
       static double mp_vector_print(_cimg_math_parser& mp) {

From c0d265478ba6988b1b148bfb4e1e079d9edf2c1f Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 17:09:19 +0100
Subject: [PATCH 06/24] math_parser(): Add new function 'histogram()'.

---
 CImg.h | 41 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)

diff --git a/CImg.h b/CImg.h
index 61486625..305107bd 100644
--- a/CImg.h
+++ b/CImg.h
@@ -20904,6 +20904,31 @@ namespace cimg_library {
               } else { if (ss2!=se1) break; p1 = ~0U; }
               _cimg_mp_scalar1(mp_image_h,p1);
             }
+
+            if (!std::strncmp(ss,"histogram(",10)) { // Compute histogram
+              _cimg_mp_op("Function 'histogram()'");
+              s0 = ss + 10;
+              s1 = s0; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(s0,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg2 = compile(s1,s2,depth1,0,block_flags);
+              arg3 = arg4 = ~0U;
+              if (s2<se1) {
+                s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+                arg3 = compile(s2,s1,depth1,0,block_flags);
+                arg4 = s1<se1?compile(++s1,se1,depth1,0,block_flags):~0U;
+              }
+              _cimg_mp_check_type(arg1,1,2,0);
+              _cimg_mp_check_const_scalar(arg2,2,3);
+              if (arg3!=~0U) _cimg_mp_check_type(arg3,3,1,0);
+              if (arg4!=~0U) _cimg_mp_check_type(arg4,3,1,0);
+              p1 = _cimg_mp_size(arg1);
+              arg2 = (unsigned int)mem[arg2];
+              pos = vector(arg2);
+              CImg<ulongT>::vector((ulongT)mp_vector_histogram,pos,arg1,p1,arg2,arg3,arg4).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
             break;
 
           case 'i' :
@@ -28165,6 +28190,22 @@ namespace cimg_library {
         return cimg::type<double>::nan();
       }
 
+      static double mp_vector_histogram(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int
+          siz = (unsigned int)mp.opcode[3],
+          nb_levels = (unsigned int)mp.opcode[4];
+        const double *const ptrs = &_mp_arg(2) + 1;
+        CImg<doubleT> img(ptrs,siz,1,1,1,true);
+        double min_value, max_value;
+        if ((unsigned int)mp.opcode[5]==~0U || (unsigned int)mp.opcode[6]==~0U)
+          min_value = img.min_max(max_value);
+        if ((unsigned int)mp.opcode[5]!=~0U) min_value = _mp_arg(5);
+        if ((unsigned int)mp.opcode[6]!=~0U) max_value = _mp_arg(6);
+        CImg<doubleT>(ptrd,nb_levels,1,1,1,true) = img.get_histogram(nb_levels,min_value,max_value);
+        return cimg::type<double>::nan();
+      }
+
       static double _mp_vector_hypot(_cimg_math_parser& mp) {
         switch ((unsigned int)mp.opcode[2]) {
           case 5 : return cimg::abs(_mp_arg(4));

From deec4620e1ccae33bdc9a188b323bfd2861a6d3b Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 17:16:26 +0100
Subject: [PATCH 07/24] .

---
 CImg.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CImg.h b/CImg.h
index 305107bd..b469b73e 100644
--- a/CImg.h
+++ b/CImg.h
@@ -28197,7 +28197,7 @@ namespace cimg_library {
           nb_levels = (unsigned int)mp.opcode[4];
         const double *const ptrs = &_mp_arg(2) + 1;
         CImg<doubleT> img(ptrs,siz,1,1,1,true);
-        double min_value, max_value;
+        double min_value = 0, max_value = 0;
         if ((unsigned int)mp.opcode[5]==~0U || (unsigned int)mp.opcode[6]==~0U)
           min_value = img.min_max(max_value);
         if ((unsigned int)mp.opcode[5]!=~0U) min_value = _mp_arg(5);

From 177183d3af15830c761dcd09174e25487a9a31a8 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 17:48:56 +0100
Subject: [PATCH 08/24] math_parser(): Add new function 'equalize()'.

---
 CImg.h | 41 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)

diff --git a/CImg.h b/CImg.h
index b469b73e..04230df6 100644
--- a/CImg.h
+++ b/CImg.h
@@ -20548,6 +20548,31 @@ namespace cimg_library {
               _cimg_mp_return_nan();
             }
 
+            if (!std::strncmp(ss,"equalize(",9)) { // Equalize
+              _cimg_mp_op("Function 'equalize()'");
+              s0 = ss + 9;
+              s1 = s0; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(s0,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg2 = compile(s1,s2,depth1,0,block_flags);
+              arg3 = arg4 = ~0U;
+              if (s2<se1) {
+                s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+                arg3 = compile(s2,s1,depth1,0,block_flags);
+                arg4 = s1<se1?compile(++s1,se1,depth1,0,block_flags):~0U;
+              }
+              _cimg_mp_check_type(arg1,1,2,0);
+              _cimg_mp_check_const_scalar(arg2,2,3);
+              if (arg3!=~0U) _cimg_mp_check_type(arg3,3,1,0);
+              if (arg4!=~0U) _cimg_mp_check_type(arg4,3,1,0);
+              p1 = _cimg_mp_size(arg1);
+              arg2 = (unsigned int)mem[arg2];
+              pos = vector(p1);
+              CImg<ulongT>::vector((ulongT)mp_vector_equalize,pos,arg1,p1,arg2,arg3,arg4).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
 #if cimg_use_cpp11==1
             if (!std::strncmp(ss,"erf(",4)) { // Error function
               _cimg_mp_op("Function 'erf()'");
@@ -28190,6 +28215,22 @@ namespace cimg_library {
         return cimg::type<double>::nan();
       }
 
+      static double mp_vector_equalize(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int
+          siz = (unsigned int)mp.opcode[3],
+          nb_levels = (unsigned int)mp.opcode[4];
+        const double *const ptrs = &_mp_arg(2) + 1;
+        CImg<doubleT> img(ptrs,siz,1,1,1,true);
+        double min_value = 0, max_value = 0;
+        if ((unsigned int)mp.opcode[5]==~0U || (unsigned int)mp.opcode[6]==~0U)
+          min_value = img.min_max(max_value);
+        if ((unsigned int)mp.opcode[5]!=~0U) min_value = _mp_arg(5);
+        if ((unsigned int)mp.opcode[6]!=~0U) max_value = _mp_arg(6);
+        CImg<doubleT>(ptrd,siz,1,1,1,true) = img.get_equalize(nb_levels,min_value,max_value);
+        return cimg::type<double>::nan();
+      }
+
       static double mp_vector_histogram(_cimg_math_parser& mp) {
         double *const ptrd = &_mp_arg(1) + 1;
         const unsigned int

From d2e10163810e6a874489be097f2a4b8cddc8cb62 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 19:21:42 +0100
Subject: [PATCH 09/24] math_parser(): Add new function 'mirror()'.

---
 CImg.h | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)

diff --git a/CImg.h b/CImg.h
index 04230df6..a449c6c0 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21377,6 +21377,42 @@ namespace cimg_library {
               _cimg_mp_return(pos);
             }
 
+            if (!std::strncmp(ss,"mirror(",7)) { // Mirror image
+              _cimg_mp_op("Function 'mirror()'");
+              s1 = ss7; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(ss7,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg2 = compile(s1,s2,depth1,0,block_flags);
+              s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg3 = compile(s2,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg4 = compile(s1,s2,depth1,0,block_flags);
+              s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg5 = compile(s2,s1,depth1,0,block_flags);
+              arg6 = compile(++s1,se1,depth1,0,block_flags);
+              _cimg_mp_check_type(arg1,1,2,0);
+              _cimg_mp_check_const_scalar(arg2,2,3);
+              _cimg_mp_check_const_scalar(arg3,3,3);
+              _cimg_mp_check_const_scalar(arg4,4,3);
+              _cimg_mp_check_const_scalar(arg5,5,3);
+              p1 = _cimg_mp_size(arg1);
+              p2 = _cimg_mp_size(arg6);
+              arg2 = (unsigned int)mem[arg2];
+              arg3 = (unsigned int)mem[arg3];
+              arg4 = (unsigned int)mem[arg4];
+              arg5 = (unsigned int)mem[arg5];
+              if (arg2*arg3*arg4*arg5!=std::max(1U,p1))
+                throw CImgArgumentException("[" cimg_appname "_math_parser] "
+                                            "CImg<%s>::%s: %s: Input vector size (%lu values) and its specified "
+                                            "geometry (%u,%u,%u,%u) (%lu values) do not match.",
+                                            pixel_type(),_cimg_mp_calling_function,s_op,
+                                            std::max(p1,1U),arg2,arg3,arg4,arg5,(ulongT)arg2*arg3*arg4*arg5);
+              pos = vector(arg2*arg3*arg4*arg5);
+              CImg<ulongT>::vector((ulongT)mp_vector_mirror,pos,arg1,arg2,arg3,arg4,arg5,arg6,p2).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
             if (!std::strncmp(ss,"mproj(",6)) { // Project matrix onto dictionary
               _cimg_mp_op("Function 'mproj()'");
               s1 = ss6; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
@@ -28399,6 +28435,26 @@ namespace cimg_library {
         return cimg::type<double>::nan();
       }
 
+      static double mp_vector_mirror(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int
+          wA = (unsigned int)mp.opcode[3],
+          hA = (unsigned int)mp.opcode[4],
+          dA = (unsigned int)mp.opcode[5],
+          sA = (unsigned int)mp.opcode[6],
+          sizp = (unsigned int)mp.opcode[8];
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          *const ptrp = &_mp_arg(7) + 1;
+        CImg<charT> str(std::max(1U,sizp) + 1);
+        if (!sizp) str[0] = _mp_arg(7);
+        else for (unsigned int p = 0; p<sizp; ++p) str[p] = (char)ptrp[p];
+        str.back() = 0;
+        CImg<doubleT>(ptrd,wA,hA,dA,sA,true) = CImg<doubleT>(ptrs,wA,hA,dA,sA,true).
+          get_mirror(str);
+        return cimg::type<double>::nan();
+      }
+
       static double mp_vector_neq(_cimg_math_parser& mp) {
         return !mp_vector_eq(mp);
       }

From a587e2af6b541a0fbe7bff40520563687e37f5b1 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 20:23:42 +0100
Subject: [PATCH 10/24] cimg_math_parser(): Add new functions 'mirror()' and
 'cumulate()'.

---
 CImg.h | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 68 insertions(+)

diff --git a/CImg.h b/CImg.h
index a449c6c0..364063a5 100644
--- a/CImg.h
+++ b/CImg.h
@@ -19848,6 +19848,50 @@ namespace cimg_library {
               _cimg_mp_scalar3(mp_cut,arg1,arg2,arg3);
             }
 
+            if (!std::strncmp(ss,"cumulate(",9)) { // Mirror image
+              _cimg_mp_op("Function 'cumulate()'");
+              s0 = ss + 9;
+              s1 = s0; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(s0,s1,depth1,0,block_flags);
+              _cimg_mp_check_type(arg1,1,2,0);
+              p1 = _cimg_mp_size(arg1);
+
+              arg2 = p1;
+              arg3 = arg4 = arg5 = 1;
+              arg6 = p2 = ~0U;
+              if (s1<se1) {
+                s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+                arg2 = compile(s1,s2,depth1,0,block_flags);
+                s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+                arg3 = compile(s2,s1,depth1,0,block_flags);
+                s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+                arg4 = compile(s1,s2,depth1,0,block_flags);
+                s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+                arg5 = compile(s2,s1,depth1,0,block_flags);
+                arg6 = s1<se1?compile(++s1,se1,depth1,0,block_flags):~0U;
+                _cimg_mp_check_const_scalar(arg2,2,3);
+                _cimg_mp_check_const_scalar(arg3,3,3);
+                _cimg_mp_check_const_scalar(arg4,4,3);
+                _cimg_mp_check_const_scalar(arg5,5,3);
+                arg2 = (unsigned int)mem[arg2];
+                arg3 = (unsigned int)mem[arg3];
+                arg4 = (unsigned int)mem[arg4];
+                arg5 = (unsigned int)mem[arg5];
+                p2 = arg6!=~0U?_cimg_mp_size(arg6):~0U;
+              }
+
+              if (arg2*arg3*arg4*arg5!=std::max(1U,p1))
+                throw CImgArgumentException("[" cimg_appname "_math_parser] "
+                                            "CImg<%s>::%s: %s: Input vector size (%lu values) and its specified "
+                                            "geometry (%u,%u,%u,%u) (%lu values) do not match.",
+                                            pixel_type(),_cimg_mp_calling_function,s_op,
+                                            std::max(p1,1U),arg2,arg3,arg4,arg5,(ulongT)arg2*arg3*arg4*arg5);
+              pos = vector(arg2*arg3*arg4*arg5);
+              CImg<ulongT>::vector((ulongT)mp_vector_cumulate,pos,arg1,arg2,arg3,arg4,arg5,arg6,p2).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
             if (!std::strncmp(ss,"convolve(",9) || !std::strncmp(ss,"correlate(",10)) { // Convolve & Correlate
               is_sth = *ss2=='n'; // is_convolve?
               _cimg_mp_op(is_sth?"Function 'convolve()'":"Function 'correlate()'");
@@ -28199,6 +28243,29 @@ namespace cimg_library {
         return cimg::type<double>::nan();
       }
 
+      static double mp_vector_cumulate(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int
+          wA = (unsigned int)mp.opcode[3],
+          hA = (unsigned int)mp.opcode[4],
+          dA = (unsigned int)mp.opcode[5],
+          sA = (unsigned int)mp.opcode[6],
+          sizp = (unsigned int)mp.opcode[8];
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          *const ptrp = sizp!=~0U?&_mp_arg(7) + 1:0;
+        CImg<charT> str;
+        if (ptrp) {
+          str.assign(std::max(1U,sizp) + 1);
+          if (!sizp) str[0] = _mp_arg(7);
+          else for (unsigned int p = 0; p<sizp; ++p) str[p] = (char)ptrp[p];
+          str.back() = 0;
+        }
+        CImg<doubleT>(ptrd,wA,hA,dA,sA,true) = CImg<doubleT>(ptrs,wA,hA,dA,sA,true).
+          get_cumulate(str);
+        return cimg::type<double>::nan();
+      }
+
       static double mp_vector_draw(_cimg_math_parser& mp) {
         double *const ptrd = &_mp_arg(1) + 1;
         const double *const ptrs = &_mp_arg(7) + 1;
@@ -40812,6 +40879,7 @@ namespace cimg_library {
        \note \c axes may contains multiple characters, e.g. \c "xyz"
     **/
     CImg<T>& cumulate(const char *const axes) {
+      if (!axes) return cumulate();
       for (const char *s = axes; *s; ++s) cumulate(*s);
       return *this;
     }

From 55eaa031fcb625e774c62a770ac68fdeb4a674ab Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 21:15:41 +0100
Subject: [PATCH 11/24] cimg_math_parser(): New function 'index()'.

---
 CImg.h | 87 +++++++++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 74 insertions(+), 13 deletions(-)

diff --git a/CImg.h b/CImg.h
index 364063a5..81ad872a 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21051,6 +21051,48 @@ namespace cimg_library {
               _cimg_mp_return(pos);
             }
 
+            if (!std::strncmp(ss,"index(",6)) { // Index colors
+              _cimg_mp_op("Function 'index()'");
+              s1 = ss6; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(ss6,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg2 = compile(s1,s2,depth1,0,block_flags);
+              s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg3 = compile(s2,s1,depth1,0,block_flags);
+              arg4 = arg5 = 0;
+              if (s1<se1) {
+                s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+                arg4 = compile(s1,s2,depth1,0,block_flags);
+                arg5 = s2<se1?compile(++s2,se1,depth1,0,block_flags):0;
+              }
+              _cimg_mp_check_type(arg1,1,2,0);
+              _cimg_mp_check_type(arg2,2,2,0);
+              _cimg_mp_check_const_scalar(arg3,3,3);
+              _cimg_mp_check_type(arg4,4,1,0);
+              _cimg_mp_check_const_scalar(arg5,5,2);
+              p1 = _cimg_mp_size(arg1);
+              p2 = _cimg_mp_size(arg2);
+              arg3 = (unsigned int)mem[arg3];
+              arg5 = (unsigned int)mem[arg5];
+              p3 = p2/arg3; // Number of color entries
+              if (p2!=p3*arg3)
+                throw CImgArgumentException("[" cimg_appname "_math_parser] "
+                                            "CImg<%s>::%s: %s: Colormap size (%lu values) and specified "
+                                            "dimension of colormap entries (%u) do not match.",
+                                            pixel_type(),_cimg_mp_calling_function,s_op,
+                                            std::max(p2,1U),arg3);
+              if (p1%arg3)
+                throw CImgArgumentException("[" cimg_appname "_math_parser] "
+                                            "CImg<%s>::%s: %s: Input vector size (%lu values) and specified "
+                                            "dimension of colormap entries (%u) do not match.",
+                                            pixel_type(),_cimg_mp_calling_function,s_op,
+                                            std::max(p1,1U),arg3);
+              pos = vector(arg5?p1:p1/arg3);
+              CImg<ulongT>::vector((ulongT)mp_vector_index,pos,arg1,p1,arg2,p2,arg3,arg4,arg5).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
             if (!std::strncmp(ss,"inrange(",8)) { // Check value range
               _cimg_mp_op("Function 'inrange()'");
               s1 = ss8; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
@@ -28359,6 +28401,25 @@ namespace cimg_library {
         return _mp_vector_norm2(mp);
       }
 
+      static double mp_vector_index(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int
+          sizA = (unsigned int)mp.opcode[3],
+          sizP = (unsigned int)mp.opcode[5],
+          dim_colors = (unsigned int)mp.opcode[6],
+          nb_colors = sizP/dim_colors,
+          wA = sizA/dim_colors;
+        const double dithering = _mp_arg(7);
+        const bool map_colors = (bool)mp.opcode[8];
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          *const ptrp = &_mp_arg(4) + 1;
+        CImg<doubleT> colormap(ptrp,nb_colors,1,1,dim_colors,true);
+        CImg<doubleT>(ptrd,wA,1,1,map_colors?dim_colors:1,true) = CImg<doubleT>(ptrs,wA,1,1,dim_colors,true).
+          get_index(colormap,dithering,map_colors);
+        return cimg::type<double>::nan();
+      }
+
       static double mp_vector_init(_cimg_math_parser& mp) {
         unsigned int
           ptrs = 4U,
@@ -34235,7 +34296,7 @@ namespace cimg_library {
     /**
        \param colormap Multi-valued colormap used as the basis for multi-valued pixel indexing.
        \param dithering Level of dithering (0=disable, 1=standard level).
-       \param map_indexes Tell if the values of the resulting image are the colormap indices or the colormap vectors.
+       \param map_colors Tell if the values of the resulting image are the colormap indices or the colormap vectors.
        \note
        - \p img.index(colormap,dithering,1) is equivalent to <tt>img.index(colormap,dithering,0).map(colormap)</tt>.
        \par Example
@@ -34247,14 +34308,14 @@ namespace cimg_library {
        \image html ref_index.jpg
     **/
     template<typename t>
-    CImg<T>& index(const CImg<t>& colormap, const float dithering=1, const bool map_indexes=false) {
-      return get_index(colormap,dithering,map_indexes).move_to(*this);
+    CImg<T>& index(const CImg<t>& colormap, const float dithering=1, const bool map_colors=false) {
+      return get_index(colormap,dithering,map_colors).move_to(*this);
     }
 
     //! Index multi-valued pixels regarding to a specified colormap \newinstance.
     template<typename t>
     CImg<typename CImg<t>::Tuint>
-    get_index(const CImg<t>& colormap, const float dithering=1, const bool map_indexes=true) const {
+    get_index(const CImg<t>& colormap, const float dithering=1, const bool map_colors=true) const {
       if (colormap._spectrum!=_spectrum)
         throw CImgArgumentException(_cimg_instance
                                     "index(): Instance and specified colormap (%u,%u,%u,%u,%p) "
@@ -34267,7 +34328,7 @@ namespace cimg_library {
       const ulongT
         whd = (ulongT)_width*_height*_depth,
         pwhd = (ulongT)colormap._width*colormap._height*colormap._depth;
-      CImg<tuint> res(_width,_height,_depth,map_indexes?_spectrum:1);
+      CImg<tuint> res(_width,_height,_depth,map_colors?_spectrum:1);
       if (dithering>0) { // Dithered versions
         tuint *ptrd = res._data;
         const float ndithering = cimg::cut(dithering,0,1)/16;
@@ -34293,7 +34354,7 @@ namespace cimg_library {
               }
               const Tfloat err0 = ((*(ptrs0++)=val0) - (Tfloat)*ptrmin0)*ndithering;
               *ptrs0+=7*err0; *(ptrsn0 - 1)+=3*err0; *(ptrsn0++)+=5*err0; *ptrsn0+=err0;
-              if (map_indexes) *(ptrd++) = (tuint)*ptrmin0; else *(ptrd++) = (tuint)(ptrmin0 - colormap._data);
+              if (map_colors) *(ptrd++) = (tuint)*ptrmin0; else *(ptrd++) = (tuint)(ptrmin0 - colormap._data);
             }
             cimg::swap(cache_current,cache_next);
           }
@@ -34328,7 +34389,7 @@ namespace cimg_library {
               *(ptrsn0 - 1)+=3*err0; *(ptrsn1 - 1)+=3*err1;
               *(ptrsn0++)+=5*err0; *(ptrsn1++)+=5*err1;
               *ptrsn0+=err0; *ptrsn1+=err1;
-              if (map_indexes) { *(ptrd++) = (tuint)*ptrmin0; *(ptrd1++) = (tuint)*ptrmin1; }
+              if (map_colors) { *(ptrd++) = (tuint)*ptrmin0; *(ptrd1++) = (tuint)*ptrmin1; }
               else *(ptrd++) = (tuint)(ptrmin0 - colormap._data);
             }
             cimg::swap(cache_current,cache_next);
@@ -34373,7 +34434,7 @@ namespace cimg_library {
               *(ptrsn0++)+=5*err0; *(ptrsn1++)+=5*err1; *(ptrsn2++)+=5*err2;
               *ptrsn0+=err0; *ptrsn1+=err1; *ptrsn2+=err2;
 
-              if (map_indexes) {
+              if (map_colors) {
                 *(ptrd++) = (tuint)*ptrmin0; *(ptrd1++) = (tuint)*ptrmin1; *(ptrd2++) = (tuint)*ptrmin2;
               } else *(ptrd++) = (tuint)(ptrmin0 - colormap._data);
             }
@@ -34407,7 +34468,7 @@ namespace cimg_library {
                 *_ptrs+=7*err; *(_ptrsn++)+=3*err; *(_ptrsn++)+=5*err; *_ptrsn+=err;
                 _ptrmin+=pwhd; _ptrs+=cwhd - 1; _ptrsn+=cwhd - 2;
               }
-              if (map_indexes) {
+              if (map_colors) {
                 tuint *_ptrd = ptrd++;
                 cimg_forC(*this,c) { *_ptrd = (tuint)*ptrmin; _ptrd+=whd; ptrmin+=pwhd; }
               }
@@ -34430,7 +34491,7 @@ namespace cimg_library {
                 const Tfloat pval0 = (Tfloat)*(ptrp0++) - val0, dist = pval0*pval0;
                 if (dist<distmin) { ptrmin0 = ptrp0 - 1; distmin = dist; }
               }
-              if (map_indexes) *(ptrd++) = (tuint)*ptrmin0; else *(ptrd++) = (tuint)(ptrmin0 - colormap._data);
+              if (map_colors) *(ptrd++) = (tuint)*ptrmin0; else *(ptrd++) = (tuint)(ptrmin0 - colormap._data);
             }
           }
         } break;
@@ -34448,7 +34509,7 @@ namespace cimg_library {
                   dist = pval0*pval0 + pval1*pval1;
                 if (dist<distmin) { ptrmin0 = ptrp0 - 1; distmin = dist; }
               }
-              if (map_indexes) { *(ptrd++) = (tuint)*ptrmin0; *(ptrd1++) = (tuint)*(ptrmin0 + pwhd); }
+              if (map_colors) { *(ptrd++) = (tuint)*ptrmin0; *(ptrd1++) = (tuint)*(ptrmin0 + pwhd); }
               else *(ptrd++) = (tuint)(ptrmin0 - colormap._data);
             }
           }
@@ -34471,7 +34532,7 @@ namespace cimg_library {
                   dist = pval0*pval0 + pval1*pval1 + pval2*pval2;
                 if (dist<distmin) { ptrmin0 = ptrp0 - 1; distmin = dist; }
               }
-              if (map_indexes) {
+              if (map_colors) {
                 *(ptrd++) = (tuint)*ptrmin0;
                 *(ptrd1++) = (tuint)*(ptrmin0 + pwhd);
                 *(ptrd2++) = (tuint)*(ptrmin0 + 2*pwhd);
@@ -34491,7 +34552,7 @@ namespace cimg_library {
                 cimg_forC(*this,c) { dist+=cimg::sqr((Tfloat)*_ptrs - (Tfloat)*_ptrp); _ptrs+=whd; _ptrp+=pwhd; }
                 if (dist<distmin) { ptrmin = ptrp; distmin = dist; }
               }
-              if (map_indexes) {
+              if (map_colors) {
                 tuint *_ptrd = ptrd++;
                 cimg_forC(*this,c) { *_ptrd = (tuint)*ptrmin; _ptrd+=whd; ptrmin+=pwhd; }
               }

From f72db3024070c2f03a67d144c2d3788bed65b03a Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 21:36:43 +0100
Subject: [PATCH 12/24] cimg_math_parser(): New function 'normalize()'.

---
 CImg.h | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/CImg.h b/CImg.h
index 81ad872a..a85555e7 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21671,6 +21671,25 @@ namespace cimg_library {
               _cimg_mp_const_scalar((double)arg1);
             }
 
+            if (!std::strncmp(ss,"normalize(",10)) { // Normalize
+              _cimg_mp_op("Function 'normalize()'");
+              s0 = ss + 10;
+              s1 = s0; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(s0,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg2 = compile(s1,s2,depth1,0,block_flags);
+              s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg3 = compile(s2,se1,depth1,0,block_flags);
+              _cimg_mp_check_type(arg2,2,1,0);
+              _cimg_mp_check_type(arg3,3,1,0);
+              if (_cimg_mp_is_scalar(arg1)) _cimg_mp_return(arg2);
+              p1 = _cimg_mp_size(arg1);
+              pos = vector(p1);
+              CImg<ulongT>::vector((ulongT)mp_vector_normalize,pos,arg1,p1,arg2,arg3).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
             if (!std::strncmp(ss,"normp(",6)) { // Lp norm, with variable argument p.
               _cimg_mp_op("Function 'normp()'");
               s1 = ss6; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
@@ -28587,6 +28606,17 @@ namespace cimg_library {
         return !mp_vector_eq(mp);
       }
 
+      static double mp_vector_normalize(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int siz = (unsigned int)mp.opcode[3];
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          min_value = _mp_arg(4),
+          max_value = _mp_arg(5);
+        CImg<doubleT>(ptrd,siz,1,1,1,true) = CImg<doubleT>(ptrs,siz,1,1,1,true).get_normalize(min_value,max_value);
+        return cimg::type<double>::nan();
+      }
+
       static double _mp_vector_norm0(_cimg_math_parser& mp) {
         const unsigned int siz = (unsigned int)mp.opcode[2];
         double res = 0;

From 429cb0ffd5c9deeb85b54bc4fdf9723a75ca35c6 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sat, 30 Dec 2023 21:48:41 +0100
Subject: [PATCH 13/24] cimg_math_parser(): New function 'noise()'.

---
 CImg.h | 56 +++++++++++++++++++++++++++++++++++++++++++-------------
 1 file changed, 43 insertions(+), 13 deletions(-)

diff --git a/CImg.h b/CImg.h
index a85555e7..802d4c5a 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21671,6 +21671,23 @@ namespace cimg_library {
               _cimg_mp_const_scalar((double)arg1);
             }
 
+            if (!std::strncmp(ss,"noise(",6)) { // Add noise
+              _cimg_mp_op("Function 'noise()'");
+              s1 = ss6; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(ss6,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg2 = compile(s1,s2,depth1,0,block_flags);
+              arg3 = s2<se1?compile(++s2,se1,depth1,0,block_flags):0;
+              _cimg_mp_check_type(arg1,1,2,0);
+              _cimg_mp_check_type(arg2,2,1,0);
+              _cimg_mp_check_type(arg3,3,1,0);
+              p1 = _cimg_mp_size(arg1);
+              pos = vector(p1);
+              CImg<ulongT>::vector((ulongT)mp_vector_noise,pos,arg1,p1,arg2,arg3).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
             if (!std::strncmp(ss,"normalize(",10)) { // Normalize
               _cimg_mp_op("Function 'normalize()'");
               s0 = ss + 10;
@@ -28606,6 +28623,19 @@ namespace cimg_library {
         return !mp_vector_eq(mp);
       }
 
+      static double mp_vector_noise(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const unsigned int
+          siz = (unsigned int)mp.opcode[3],
+          noise_type = (unsigned int)_mp_arg(5);
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          amplitude = _mp_arg(4);
+        CImg<doubleT>(ptrd,siz,1,1,1,true) = CImg<doubleT>(ptrs,siz,1,1,1,true).get_noise(amplitude,noise_type);
+        return cimg::type<double>::nan();
+      }
+
+
       static double mp_vector_normalize(_cimg_math_parser& mp) {
         double *const ptrd = &_mp_arg(1) + 1;
         const unsigned int siz = (unsigned int)mp.opcode[3];
@@ -33817,13 +33847,13 @@ namespace cimg_library {
        \endcode
        \image html ref_noise.jpg
     **/
-    CImg<T>& noise(const double sigma, const unsigned int noise_type=0) {
+    CImg<T>& noise(const double amplitude, const unsigned int noise_type=0) {
       if (is_empty()) return *this;
       const Tfloat vmin = (Tfloat)cimg::type<T>::min(), vmax = (Tfloat)cimg::type<T>::max();
-      Tfloat nsigma = (Tfloat)sigma, m = 0, M = 0;
-      if (nsigma==0 && noise_type!=3) return *this;
-      if (nsigma<0 || noise_type==2) m = (Tfloat)min_max(M);
-      if (nsigma<0) nsigma = (Tfloat)(-nsigma*(M-m)/100.);
+      Tfloat namplitude = (Tfloat)amplitude, m = 0, M = 0;
+      if (namplitude==0 && noise_type!=3) return *this;
+      if (namplitude<0 || noise_type==2) m = (Tfloat)min_max(M);
+      if (namplitude<0) namplitude = (Tfloat)(-namplitude*(M-m)/100.);
       switch (noise_type) {
       case 0 : { // Gaussian noise
         cimg_pragma_openmp(parallel cimg_openmp_if_size(size(),131072)) {
@@ -33834,7 +33864,7 @@ namespace cimg_library {
 #endif
           cimg_pragma_openmp(for)
             cimg_rofoff(*this,off) {
-            Tfloat val = (Tfloat)(_data[off] + nsigma*cimg::grand(&rng));
+            Tfloat val = (Tfloat)(_data[off] + namplitude*cimg::grand(&rng));
             if (val>vmax) val = vmax;
             if (val<vmin) val = vmin;
             _data[off] = (T)val;
@@ -33851,7 +33881,7 @@ namespace cimg_library {
 #endif
           cimg_pragma_openmp(for)
             cimg_rofoff(*this,off) {
-            Tfloat val = (Tfloat)(_data[off] + nsigma*cimg::rand(-1,1,&rng));
+            Tfloat val = (Tfloat)(_data[off] + namplitude*cimg::rand(-1,1,&rng));
             if (val>vmax) val = vmax;
             if (val<vmin) val = vmin;
             _data[off] = (T)val;
@@ -33860,7 +33890,7 @@ namespace cimg_library {
         }
       } break;
       case 2 : { // Salt & Pepper noise
-        if (nsigma<0) nsigma = -nsigma;
+        if (namplitude<0) namplitude = -namplitude;
         if (M==m) {
           if (cimg::type<T>::is_float()) { --m; ++M; }
           else { m = (Tfloat)cimg::type<T>::min(); M = (Tfloat)cimg::type<T>::max(); }
@@ -33872,7 +33902,7 @@ namespace cimg_library {
           rng+=omp_get_thread_num();
 #endif
           cimg_pragma_openmp(for)
-            cimg_rofoff(*this,off) if (cimg::rand(100,&rng)<nsigma) _data[off] = (T)(cimg::rand(1,&rng)<0.5?M:m);
+            cimg_rofoff(*this,off) if (cimg::rand(100,&rng)<namplitude) _data[off] = (T)(cimg::rand(1,&rng)<0.5?M:m);
           cimg::srand(rng);
           }
       } break;
@@ -33900,8 +33930,8 @@ namespace cimg_library {
             cimg_rofoff(*this,off) {
             const Tfloat
               val0 = (Tfloat)_data[off]/sqrt2,
-              re = (Tfloat)(val0 + nsigma*cimg::grand(&rng)),
-              im = (Tfloat)(val0 + nsigma*cimg::grand(&rng));
+              re = (Tfloat)(val0 + namplitude*cimg::grand(&rng)),
+              im = (Tfloat)(val0 + namplitude*cimg::grand(&rng));
             Tfloat val = cimg::hypot(re,im);
             if (val>vmax) val = vmax;
             if (val<vmin) val = vmin;
@@ -33921,8 +33951,8 @@ namespace cimg_library {
     }
 
     //! Add random noise to pixel values \newinstance.
-    CImg<T> get_noise(const double sigma, const unsigned int noise_type=0) const {
-      return (+*this).noise(sigma,noise_type);
+    CImg<T> get_noise(const double amplitude, const unsigned int noise_type=0) const {
+      return (+*this).noise(amplitude,noise_type);
     }
 
     //! Linearly normalize pixel values.

From cfb9f02cb61aa0edaaaffa4f18210ff395aa68f2 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 10:50:49 +0100
Subject: [PATCH 14/24] Add variation of method 'CImg<T>::rand()' that allow to
 specify the probability density function.

---
 CImg.h | 44 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/CImg.h b/CImg.h
index 802d4c5a..aa552154 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33811,6 +33811,50 @@ namespace cimg_library {
       return (+*this).rand(val_min,val_max);
     }
 
+    //! Fill image with random values following specified distribution and range.
+    /**
+       \param pdf Probability density function.
+       \param val_min Minimal authorized random value.
+       \param val_max Maximal authorized random value.
+       \param prec Precision of generated values. Set to '0' for automatic precision. Negative value means percentage.
+     **/
+    template<typename t>
+    CImg<T>& rand(const CImg<t>& pdf, const T& val_min, const T& val_max, const int prec=-400) {
+      typedef _cimg_tfloat tfloat;
+      const unsigned int siz = (unsigned int)pdf.size();
+      if (siz<2) return fill(val_min);
+      const tfloat
+        delta = (tfloat)val_max - (tfloat)val_min,
+        delta_over_siz1 = delta/(siz - 1);
+      const unsigned int nprec = !prec?siz:prec<0?(unsigned int)(-siz*prec/100):(unsigned int)prec;
+      CImg<tfloat> cdf = pdf.get_max((t)0).cumulate();
+      cdf*=(nprec - 1)/cdf.back();
+
+      // Compute inverse cdf.
+      CImg<tfloat> icdf(nprec);
+      unsigned int k = 0;
+      tfloat p = 0;
+      cimg_forX(icdf,x) {
+        while (k<siz && (!(p = cdf[k]) || p<x)) ++k;
+        icdf[x] = val_min + k*delta_over_siz1;
+      }
+
+      // Generate random numbers.
+      cimg_for(*this,ptr,T) {
+        const unsigned int
+          _ind = (unsigned int)cimg::rand(0,(double)nprec),
+          ind = _ind==nprec?0:_ind;
+        *ptr = (T)icdf[ind];
+      }
+      return *this;
+    }
+
+    //! Fill image with random values following specified distribution and range \newinstance.
+    template<typename t>
+    CImg<T> get_rand(const CImg<t>& pdf, const T& val_min, const T& val_max) const {
+      return (+*this).rand(pdf,val_min,val_max);
+    }
+
     //! Round pixel values.
     /**
        \param y Rounding precision.

From 73995db31cbfa8278d46c856235fbed306c67072 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 10:55:11 +0100
Subject: [PATCH 15/24] .

---
 CImg.h | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/CImg.h b/CImg.h
index aa552154..13d61127 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33819,21 +33819,20 @@ namespace cimg_library {
        \param prec Precision of generated values. Set to '0' for automatic precision. Negative value means percentage.
      **/
     template<typename t>
-    CImg<T>& rand(const CImg<t>& pdf, const T& val_min, const T& val_max, const int prec=-400) {
+    CImg<T>& rand(const CImg<t>& pdf, const T& val_min, const T& val_max, const int prec=65536) {
       typedef _cimg_tfloat tfloat;
       const unsigned int siz = (unsigned int)pdf.size();
       if (siz<2) return fill(val_min);
       const tfloat
         delta = (tfloat)val_max - (tfloat)val_min,
         delta_over_siz1 = delta/(siz - 1);
-      const unsigned int nprec = !prec?siz:prec<0?(unsigned int)(-siz*prec/100):(unsigned int)prec;
-      CImg<tfloat> cdf = pdf.get_max((t)0).cumulate();
-      cdf*=(nprec - 1)/cdf.back();
+      const unsigned int nprec = !prec?4*siz:prec<0?(unsigned int)(-siz*prec/100):(unsigned int)prec;
 
       // Compute inverse cdf.
-      CImg<tfloat> icdf(nprec);
+      CImg<tfloat> cdf = pdf.get_max((t)0).cumulate(), icdf(nprec);
       unsigned int k = 0;
       tfloat p = 0;
+      cdf*=(nprec - 1)/cdf.back();
       cimg_forX(icdf,x) {
         while (k<siz && (!(p = cdf[k]) || p<x)) ++k;
         icdf[x] = val_min + k*delta_over_siz1;

From e1a4b5b6da54f770a1b420ca66b0c8587168bc10 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 10:56:01 +0100
Subject: [PATCH 16/24] .

---
 CImg.h | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/CImg.h b/CImg.h
index 13d61127..e3727212 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33816,7 +33816,8 @@ namespace cimg_library {
        \param pdf Probability density function.
        \param val_min Minimal authorized random value.
        \param val_max Maximal authorized random value.
-       \param prec Precision of generated values. Set to '0' for automatic precision. Negative value means percentage.
+       \param prec Precision of generated values. Set to '0' for automatic precision.
+         Negative value means 'percentage of the pdf size'.
      **/
     template<typename t>
     CImg<T>& rand(const CImg<t>& pdf, const T& val_min, const T& val_max, const int prec=65536) {

From 701f338be55ec9889d68821da0f2372462e0c232 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 10:56:46 +0100
Subject: [PATCH 17/24] .

---
 CImg.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CImg.h b/CImg.h
index e3727212..50b08bf0 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33817,7 +33817,7 @@ namespace cimg_library {
        \param val_min Minimal authorized random value.
        \param val_max Maximal authorized random value.
        \param prec Precision of generated values. Set to '0' for automatic precision.
-         Negative value means 'percentage of the pdf size'.
+         A negative value means 'percentage of the pdf size'.
      **/
     template<typename t>
     CImg<T>& rand(const CImg<t>& pdf, const T& val_min, const T& val_max, const int prec=65536) {

From a80602e4e8ff288276f1bdcc9c132c543150d403 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 11:04:54 +0100
Subject: [PATCH 18/24] .

---
 CImg.h | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/CImg.h b/CImg.h
index 50b08bf0..7266c37c 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33840,11 +33840,20 @@ namespace cimg_library {
       }
 
       // Generate random numbers.
-      cimg_for(*this,ptr,T) {
-        const unsigned int
-          _ind = (unsigned int)cimg::rand(0,(double)nprec),
-          ind = _ind==nprec?0:_ind;
-        *ptr = (T)icdf[ind];
+      cimg_pragma_openmp(parallel cimg_openmp_if_size(size(),524288)) {
+        cimg_uint64 rng = (cimg::_rand(),cimg::rng());
+
+#if cimg_use_openmp!=0
+        rng+=omp_get_thread_num();
+#endif
+        cimg_pragma_openmp(for)
+          cimg_rofoff(*this,off) {
+          const unsigned int
+            _ind = (unsigned int)cimg::rand(0,(double)nprec,&rng),
+            ind = _ind==nprec?0:_ind;
+          _data[off] = (T)icdf[ind];
+        }
+        cimg::srand(rng);
       }
       return *this;
     }

From 91e38349b160c2bd141372b5f1fe20ee3287765a Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 11:48:09 +0100
Subject: [PATCH 19/24] Command 'rand' now accepts a '[pdf]' parameter that
 sets the probability density function followed by the random values.

---
 CImg.h | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/CImg.h b/CImg.h
index 7266c37c..d72197fe 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33813,27 +33813,27 @@ namespace cimg_library {
 
     //! Fill image with random values following specified distribution and range.
     /**
-       \param pdf Probability density function.
        \param val_min Minimal authorized random value.
        \param val_max Maximal authorized random value.
-       \param prec Precision of generated values. Set to '0' for automatic precision.
+       \param pdf Probability density function.
+       \param precision Precision of generated values. Set to '0' for automatic precision.
          A negative value means 'percentage of the pdf size'.
      **/
     template<typename t>
-    CImg<T>& rand(const CImg<t>& pdf, const T& val_min, const T& val_max, const int prec=65536) {
+    CImg<T>& rand(const T& val_min, const T& val_max, const CImg<t>& pdf, const int precision=65536) {
       typedef _cimg_tfloat tfloat;
       const unsigned int siz = (unsigned int)pdf.size();
       if (siz<2) return fill(val_min);
       const tfloat
         delta = (tfloat)val_max - (tfloat)val_min,
         delta_over_siz1 = delta/(siz - 1);
-      const unsigned int nprec = !prec?4*siz:prec<0?(unsigned int)(-siz*prec/100):(unsigned int)prec;
+      const unsigned int prec = !precision?4*siz:precision<0?(unsigned int)(-siz*precision/100):(unsigned int)precision;
 
       // Compute inverse cdf.
-      CImg<tfloat> cdf = pdf.get_max((t)0).cumulate(), icdf(nprec);
+      CImg<tfloat> cdf = pdf.get_max((t)0).cumulate(), icdf(prec);
       unsigned int k = 0;
       tfloat p = 0;
-      cdf*=(nprec - 1)/cdf.back();
+      cdf*=(prec - 1)/cdf.back();
       cimg_forX(icdf,x) {
         while (k<siz && (!(p = cdf[k]) || p<x)) ++k;
         icdf[x] = val_min + k*delta_over_siz1;
@@ -33849,8 +33849,8 @@ namespace cimg_library {
         cimg_pragma_openmp(for)
           cimg_rofoff(*this,off) {
           const unsigned int
-            _ind = (unsigned int)cimg::rand(0,(double)nprec,&rng),
-            ind = _ind==nprec?0:_ind;
+            _ind = (unsigned int)cimg::rand(0,(double)prec,&rng),
+            ind = _ind==prec?0:_ind;
           _data[off] = (T)icdf[ind];
         }
         cimg::srand(rng);
@@ -33860,8 +33860,8 @@ namespace cimg_library {
 
     //! Fill image with random values following specified distribution and range \newinstance.
     template<typename t>
-    CImg<T> get_rand(const CImg<t>& pdf, const T& val_min, const T& val_max) const {
-      return (+*this).rand(pdf,val_min,val_max);
+    CImg<T> get_rand(const T& val_min, const T& val_max, const CImg<t>& pdf, const int precision=65536) const {
+      return (+*this).rand(val_min,val_max,pdf,precision);
     }
 
     //! Round pixel values.

From d1759d8d128c8c62c2d056540d288dfa670fa9cc Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 12:11:50 +0100
Subject: [PATCH 20/24] .

---
 CImg.h | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/CImg.h b/CImg.h
index d72197fe..8a921a45 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33822,12 +33822,13 @@ namespace cimg_library {
     template<typename t>
     CImg<T>& rand(const T& val_min, const T& val_max, const CImg<t>& pdf, const int precision=65536) {
       typedef _cimg_tfloat tfloat;
-      const unsigned int siz = (unsigned int)pdf.size();
-      if (siz<2) return fill(val_min);
+      const unsigned int
+        siz = (unsigned int)pdf.size(),
+        prec = precision<0?(unsigned int)(-siz*precision/100):(unsigned int)precision;
+      if (siz<2 || precision<2) return fill(val_min);
       const tfloat
         delta = (tfloat)val_max - (tfloat)val_min,
         delta_over_siz1 = delta/(siz - 1);
-      const unsigned int prec = !precision?4*siz:precision<0?(unsigned int)(-siz*precision/100):(unsigned int)precision;
 
       // Compute inverse cdf.
       CImg<tfloat> cdf = pdf.get_max((t)0).cumulate(), icdf(prec);

From 9d0a4a3ddf65a23e72c5ff0e01a3341bbecbb4c2 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 12:16:11 +0100
Subject: [PATCH 21/24] .

---
 CImg.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CImg.h b/CImg.h
index 8a921a45..319357a4 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33837,6 +33837,7 @@ namespace cimg_library {
       cdf*=(prec - 1)/cdf.back();
       cimg_forX(icdf,x) {
         while (k<siz && (!(p = cdf[k]) || p<x)) ++k;
+        if (k==siz) --k;
         icdf[x] = val_min + k*delta_over_siz1;
       }
 

From 5b0484ed093aa613fcd2d32428ae9fa00d35ff85 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Wed, 3 Jan 2024 14:18:46 +0100
Subject: [PATCH 22/24] .

---
 CImg.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CImg.h b/CImg.h
index 319357a4..f42f5b3c 100644
--- a/CImg.h
+++ b/CImg.h
@@ -33837,7 +33837,7 @@ namespace cimg_library {
       cdf*=(prec - 1)/cdf.back();
       cimg_forX(icdf,x) {
         while (k<siz && (!(p = cdf[k]) || p<x)) ++k;
-        if (k==siz) --k;
+        if (k>=siz) { while (x<icdf.width()) icdf[x++] = val_max; break; }
         icdf[x] = val_min + k*delta_over_siz1;
       }
 

From 6198cf5a4dceb23660ab6b1cb6f494a8e0497742 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Sun, 7 Jan 2024 20:50:03 +0100
Subject: [PATCH 23/24] .

---
 CImg.h | 75 ++++++++++++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 68 insertions(+), 7 deletions(-)

diff --git a/CImg.h b/CImg.h
index f42f5b3c..35ee00fc 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21878,6 +21878,48 @@ namespace cimg_library {
               _cimg_mp_scalar1(mp_rad2deg,arg1);
             }
 
+            if ((cimg_sscanf(ss,"rand%u%c",&(arg1=~0U),&sep)==2 && sep=='(' && arg1>0) ||
+                !std::strncmp(ss,"rand(#",6) ||
+                (!std::strncmp(ss,"rand",4) && ss4<se1 && *ss4!='(' && (s=std::strchr(ss4,'('))!=0)) { // Random vector
+              _cimg_mp_op("Function 'rand()'");
+              p3 = 0;
+              if (arg1==~0U && *ss4!='(') {
+                arg1 = compile(ss4,s++,depth1,0,block_flags);
+                _cimg_mp_check_const_scalar(arg1,0,3);
+                arg1 = (unsigned int)mem[arg1];
+              } else s = std::strchr(ss4,'(') + 1;
+
+              if (arg1==~0U && *s=='#') { // Number of elements specified as first argument with '#'
+                s0 = ++s; while (s0<se1 && (*s0!=',' || level[s0 - expr._data]!=clevel1)) ++s0;
+                arg1 = compile(s,s0++,depth1,0,block_flags);
+                _cimg_mp_check_const_scalar(arg1,1,3);
+                arg1 = (unsigned int)mem[arg1];
+                p3 = 1;
+                s = s0;
+              }
+
+              s1 = s + 1; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg2 = compile(s,s1,depth1,0,block_flags);
+              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+              arg3 = compile(s1,s2,depth1,0,block_flags);
+              _cimg_mp_check_type(arg2,p3 + 1,1,0);
+              _cimg_mp_check_type(arg3,p3 + 2,1,0);
+              arg4 = arg5 = ~0U; p2 = 0;
+              if (s2<se1) {
+                s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+                arg4 = compile(s2,s1,depth1,0,block_flags);
+                arg5 = s1<se1?compile(++s1,se1,depth1,0,block_flags):~0U;
+                p2 = _cimg_mp_size(arg4);
+                _cimg_mp_check_type(arg4,p3 + 3,2,0);
+                if (arg5!=~0U) _cimg_mp_check_type(arg5,p3 + 4,1,0);
+              }
+
+              pos = vector(arg1);
+              CImg<ulongT>::vector((ulongT)mp_vector_rand,pos,arg1,arg2,arg3,arg4,p2,arg5).move_to(code);
+              return_new_comp = true;
+              _cimg_mp_return(pos);
+            }
+
             if (!std::strncmp(ss,"ref(",4)) { // Variable declaration
               _cimg_mp_op("Function 'ref()'");
               s1 = ss4; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
@@ -28772,20 +28814,39 @@ namespace cimg_library {
         return cimg::type<double>::nan();
       }
 
+      static double mp_vector_rand(_cimg_math_parser& mp) {
+        double *const ptrd = &_mp_arg(1) + 1;
+        const double *const ptr_pdf = (unsigned int)mp.opcode[5]!=~0U?&_mp_arg(5) + 1:0;
+        const unsigned int
+          siz = (unsigned int)mp.opcode[2],
+          siz_pdf = (unsigned int)mp.opcode[6],
+          prec = (unsigned int)mp.opcode[7]!=~0U?(unsigned int)std::abs(_mp_arg(7)):65536;
+        const double
+          val_min = _mp_arg(3),
+          val_max = _mp_arg(4);
+        if (!ptr_pdf)
+          CImg<doubleT>(ptrd,siz,1,1,1,true).rand(val_min,val_max);
+        else {
+          CImg<doubleT> pdf(ptr_pdf,siz_pdf,1,1,1,true);
+          CImg<doubleT>(ptrd,siz,1,1,1,true).rand(val_min,val_max,pdf,prec);
+        }
+        return cimg::type<double>::nan();
+      }
+
       static double mp_vector_resize(_cimg_math_parser& mp) {
         double *const ptrd = &_mp_arg(1) + 1;
-        const unsigned int p1 = (unsigned int)mp.opcode[2], p2 = (unsigned int)mp.opcode[4];
+        const unsigned int siz = (unsigned int)mp.opcode[2], p2 = (unsigned int)mp.opcode[4];
         const int
           interpolation = (int)_mp_arg(5),
           boundary_conditions = (int)_mp_arg(6);
         if (p2) { // Resize vector
           const double *const ptrs = &_mp_arg(3) + 1;
-          CImg<doubleT>(ptrd,p1,1,1,1,true) = CImg<doubleT>(ptrs,p2,1,1,1,true).
-            get_resize(p1,1,1,1,interpolation,boundary_conditions);
+          CImg<doubleT>(ptrd,siz,1,1,1,true) = CImg<doubleT>(ptrs,p2,1,1,1,true).
+            get_resize(siz,1,1,1,interpolation,boundary_conditions);
         } else { // Resize scalar
           const double value = _mp_arg(3);
-          CImg<doubleT>(ptrd,p1,1,1,1,true) = CImg<doubleT>(1,1,1,1,value).resize(p1,1,1,1,interpolation,
-                                                                                  boundary_conditions);
+          CImg<doubleT>(ptrd,siz,1,1,1,true) = CImg<doubleT>(1,1,1,1,value).resize(siz,1,1,1,interpolation,
+                                                                                   boundary_conditions);
         }
         return cimg::type<double>::nan();
       }
@@ -28825,8 +28886,8 @@ namespace cimg_library {
       static double mp_vector_reverse(_cimg_math_parser& mp) {
         double *const ptrd = &_mp_arg(1) + 1;
         const double *const ptrs = &_mp_arg(2) + 1;
-        const unsigned int p1 = (unsigned int)mp.opcode[3];
-        CImg<doubleT>(ptrd,p1,1,1,1,true) = CImg<doubleT>(ptrs,p1,1,1,1,true).get_mirror('x');
+        const unsigned int siz = (unsigned int)mp.opcode[3];
+        CImg<doubleT>(ptrd,siz,1,1,1,true) = CImg<doubleT>(ptrs,siz,1,1,1,true).get_mirror('x');
         return cimg::type<double>::nan();
       }
 

From 1b4660ee66df76f0cc80e45c976a2879d6fcbec8 Mon Sep 17 00:00:00 2001
From: David Tschumperle <David.Tschumperle@unicaen.fr>
Date: Mon, 8 Jan 2024 08:54:52 +0100
Subject: [PATCH 24/24] .

---
 CImg.h | 32 ++++++++++++++++++--------------
 1 file changed, 18 insertions(+), 14 deletions(-)

diff --git a/CImg.h b/CImg.h
index 35ee00fc..b9864d91 100644
--- a/CImg.h
+++ b/CImg.h
@@ -21898,22 +21898,26 @@ namespace cimg_library {
                 s = s0;
               }
 
-              s1 = s + 1; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
-              arg2 = compile(s,s1,depth1,0,block_flags);
-              s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
-              arg3 = compile(s1,s2,depth1,0,block_flags);
+              arg2 = p2 = 0; arg3 = 1; arg4 = arg5 = ~0U;
+              if (s<se1) {
+                s1 = s + 1; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+                arg3 = compile(s,s1,depth1,0,block_flags);
+                if (s1<se1) {
+                  s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+                  arg2 = arg3;
+                  arg3 = compile(s1,s2,depth1,0,block_flags);
+                  if (s2<se1) {
+                    s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+                    arg4 = compile(s2,s1,depth1,0,block_flags);
+                    arg5 = s1<se1?compile(++s1,se1,depth1,0,block_flags):~0U;
+                    p2 = _cimg_mp_size(arg4);
+                  }
+                }
+              }
               _cimg_mp_check_type(arg2,p3 + 1,1,0);
               _cimg_mp_check_type(arg3,p3 + 2,1,0);
-              arg4 = arg5 = ~0U; p2 = 0;
-              if (s2<se1) {
-                s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
-                arg4 = compile(s2,s1,depth1,0,block_flags);
-                arg5 = s1<se1?compile(++s1,se1,depth1,0,block_flags):~0U;
-                p2 = _cimg_mp_size(arg4);
-                _cimg_mp_check_type(arg4,p3 + 3,2,0);
-                if (arg5!=~0U) _cimg_mp_check_type(arg5,p3 + 4,1,0);
-              }
-
+              if (arg4!=~0U) _cimg_mp_check_type(arg4,p3 + 3,2,0);
+              if (arg5!=~0U) _cimg_mp_check_type(arg5,p3 + 4,1,0);
               pos = vector(arg1);
               CImg<ulongT>::vector((ulongT)mp_vector_rand,pos,arg1,arg2,arg3,arg4,p2,arg5).move_to(code);
               return_new_comp = true;