diff --git a/.github/workflows/cimg.yml b/.github/workflows/cimg.yml
new file mode 100644
index 0000000..90817df
--- /dev/null
+++ b/.github/workflows/cimg.yml
@@ -0,0 +1,43 @@
+name: Update CImg
+on:
+  schedule:
+    - cron: 0 0 3 * * # The first day after this workflow was merged
+  release:
+    types: [published]
+  workflow_dispatch:
+jobs:
+  check-version:
+    name: Check for updates
+    runs-on: ubuntu-latest
+    outputs:
+      latest-tag: ${{ steps.latest-tag.outputs.result }}
+    steps:
+      - name: Fetch the latest tag (could be buggy)
+        uses: actions/github-script@v7
+        id: latest-tag
+        with:
+          result-encoding: string
+          script: |
+              return (await github.rest.repos.listTags({ owner: 'GreycLab', repo: 'CImg', per_page: 28 })).data[27].name
+  pull-file:
+    name: Update CImg.h
+    runs-on: ubuntu-latest
+    needs: check-version
+    if: ${{ needs.check-version.outputs.latest-tag != 'v.3.3.3' }}
+    steps:
+      - uses: actions/checkout@v4
+        name: Checkout this repository
+        with:
+          token: ${{ secrets.CIMG_UPDATE_TOKEN }}
+      - name: Download new CImg version
+        uses: carlosperate/download-file-action@v2
+        with:
+          file-url: 'https://github.com/GreycLab/CImg/raw/${{ needs.check-version.outputs.latest-tag }}/CImg.h'
+          location: 'src'
+      - name: Commit new CImg version (and update script)
+        run: |
+          sed -i 's/v.3.3.3/${{ needs.check-version.outputs.latest-tag }}/' .github/workflows/cimg.yml
+          git config user.name 'GitHub Actions'
+          git config user.email 'actions@github.com'
+          git commit -am "Update CImg to ${{ needs.check-version.outputs.latest-tag }}"
+          git push
diff --git a/src/CImg.h b/src/CImg.h
index 844a313..63b2126 100644
--- a/src/CImg.h
+++ b/src/CImg.h
@@ -54,7 +54,7 @@
 
 // Set version number of the library.
 #ifndef cimg_version
-#define cimg_version 331
+#define cimg_version 333
 
 /*-----------------------------------------------------------
  #
@@ -5786,6 +5786,7 @@ namespace cimg_library {
     inline int system(const char *const command, const char *const module_name=0, const bool is_verbose=false) {
       cimg::unused(module_name);
 #ifdef cimg_no_system_calls
+      cimg::unused(command,is_verbose);
       return -1;
 #else
 
@@ -6841,22 +6842,19 @@ namespace cimg_library {
     //! Return sqrt(x^2 + y^2).
     template<typename T>
     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).
     template<typename T>
     inline T hypot(const T x, const T y, const T z) {
       return std::sqrt(x*x + y*y + z*z);
     }
 
-    template<typename T>
-    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 (nx<ny) { t = nx; nx = ny; } else t = ny;
-      if (nx>0) { t/=nx; return nx*std::sqrt(1 + t*t); }
-      return 0;
-    }
-
     //! Return the factorial of n
     inline double factorial(const int n) {
       if (n<0) return cimg::type<double>::nan();
@@ -19850,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()'");
@@ -19995,11 +20037,13 @@ namespace cimg_library {
             }
 
             if (!std::strncmp(ss,"da_back(",8) ||
-                !std::strncmp(ss,"da_pop(",7)) { // Get latest element in a dynamic array
+                !std::strncmp(ss,"da_pop(",7) ||
+                !std::strncmp(ss,"da_pop_heap(",12)) { // Get latest element in a dynamic array
               if (!is_inside_critical) is_parallelizable = false;
-              const bool is_pop = *ss3=='p';
-              _cimg_mp_op(is_pop?"Function 'da_pop()'":"Function 'da_back()'");
-              s0 = ss + (is_pop?7:8);
+              const bool is_pop = *ss3=='p', is_pop_heap = *ss6=='_';
+              _cimg_mp_op(is_pop_heap?"Function 'da_pop_heap()":
+                          is_pop?"Function 'da_pop()'":"Function 'da_back()'");
+              s0 = ss + (is_pop_heap?12:is_pop?7:8);
               if (*s0=='#') { // Index specified
                 s1 = ++s0; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
                 p1 = compile(s0,s1++,depth1,0,block_flags);
@@ -20010,17 +20054,19 @@ namespace cimg_library {
               p3 = (unsigned int)cimg::mod((int)mem[p1],imglist.width());
               p2 = imglist[p3]._spectrum;
               if (p2>1) pos = vector(p2); else pos = scalar(); // Return vector or scalar result
-              CImg<ulongT>::vector((ulongT)mp_da_back_or_pop,pos,p2,p1,is_pop).move_to(code);
+              CImg<ulongT>::vector((ulongT)mp_da_back_or_pop,pos,p2,p1,is_pop_heap?2:is_pop?1:0).move_to(code);
               return_new_comp = true;
               _cimg_mp_return(pos);
             }
 
             if (!std::strncmp(ss,"da_insert(",10) ||
-                !std::strncmp(ss,"da_push(",8)) { // Insert element(s) in a dynamic array
+                !std::strncmp(ss,"da_push(",8) ||
+                !std::strncmp(ss,"da_push_heap(",13)) { // Insert element(s) in a dynamic array
               if (!is_inside_critical) is_parallelizable = false;
-              const bool is_push = *ss3=='p';
-              _cimg_mp_op(is_push?"Function 'da_push()'":"Function 'da_insert()'");
-              s0 = ss + (is_push?8:10);
+              const bool is_push = *ss3=='p', is_push_heap = *ss7=='_';
+              _cimg_mp_op(is_push_heap?"Function 'da_push_heap()'":
+                          is_push?"Function 'da_push()'":"Function 'da_insert()'");
+              s0 = ss + (is_push_heap?13:is_push?8:10);
               if (*s0=='#') { // Index specified
                 s1 = ++s0; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
                 p1 = compile(s0,s1++,depth1,0,block_flags);
@@ -20031,7 +20077,8 @@ namespace cimg_library {
               if (!is_push) {
                 s0 = s1; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
                 arg1 = compile(s0,s1++,depth1,0,block_flags); // Position
-              } else arg1 = ~0U;
+              } else if (is_push_heap) arg1 = ~0U - 1;
+              else arg1 = ~0U;
 
               CImg<ulongT>::vector((ulongT)mp_da_insert_or_push,_cimg_mp_slot_nan,p1,arg1,0,0).move_to(l_opcode);
               p3 = p1==~0U?2:3;
@@ -20545,6 +20592,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()'");
@@ -20901,6 +20973,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' :
@@ -20954,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;
@@ -21068,6 +21207,15 @@ namespace cimg_library {
                 _cimg_mp_scalar1(mp_isbool,arg1);
               }
 
+              if (!std::strncmp(ss,"isconst(",8)) { // Is constant?
+                _cimg_mp_op("Function 'isconst()'");
+                if (ss8==se1) _cimg_mp_return(0);
+                try { arg1 = compile(ss8,se1,depth1,0,block_flags); }
+                catch(CImgException&) { _cimg_mp_return(0); }
+                if (_cimg_mp_is_const_scalar(arg1)) _cimg_mp_return(1);
+                _cimg_mp_return(0);
+              }
+
               if (!std::strncmp(ss,"isdir(",6)) { // Is directory?
                 _cimg_mp_op("Function 'isdir()'");
                 arg1 = compile(ss6,se1,depth1,0,block_flags);
@@ -21151,6 +21299,11 @@ namespace cimg_library {
                 _cimg_mp_return(1);
               }
 
+              if (!std::strncmp(ss,"ispercentage(",13)) { // Does argument ends with '%'?
+                _cimg_mp_op("Function 'ispercentage()'");
+                _cimg_mp_return((unsigned int)(*se2=='%'));
+              }
+
               if (!std::strncmp(ss,"isvarname(",10)) { // Is variable name?
                 _cimg_mp_op("Function 'isvarname()'");
                 arg1 = compile(ss + 10,se1,depth1,0,block_flags);
@@ -21310,6 +21463,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;
@@ -21482,6 +21671,51 @@ 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;
+              s1 = s0; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+              arg1 = compile(s0,s1,depth1,0,block_flags);
+
+              arg2 = arg4 = 0; arg3 = 1;
+              if (s1<se1) {
+                s2 = ++s1; while (s2<se1 && (*s2!=',' || level[s2 - expr._data]!=clevel1)) ++s2;
+                arg3 = compile(s1,s2,depth1,0,block_flags);
+                if (s2<se1) {
+                  s1 = ++s2; while (s1<se1 && (*s1!=',' || level[s1 - expr._data]!=clevel1)) ++s1;
+                  arg2 = arg3;
+                  arg3 = compile(s2,s1,depth1,0,block_flags);
+                  arg4 = s1<se1?compile(++s1,se1,depth1,0,block_flags):0;
+                }
+              }
+              _cimg_mp_check_type(arg2,2,1,0);
+              _cimg_mp_check_type(arg3,3,1,0);
+              _cimg_mp_check_type(arg4,4,1,0);
+              if (_cimg_mp_is_scalar(arg1)) _cimg_mp_scalar3(mp_lerp,arg2,arg3,arg4);
+              p1 = _cimg_mp_size(arg1);
+              pos = vector(p1);
+              CImg<ulongT>::vector((ulongT)mp_vector_normalize,pos,arg1,p1,arg2,arg3,arg4).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;
@@ -21491,40 +21725,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<se1 && (s = std::strchr(ss4,'('))!=0) { // Lp norm (constant p)
-              _cimg_mp_op("Function 'norm()'");
-              arg1 = s!=ss4?compile(ss4,s,depth1,0,block_flags):2;
-              _cimg_mp_check_const_scalar(arg1,0,0);
-              val = mem[arg1];
-              is_sth = true; // Tell if all arguments are constant
-              CImg<ulongT>::vector(0,0,0,arg1).move_to(l_opcode);
-              for (++s; s<se; ++s) {
-                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 (_cimg_mp_is_vector(arg2))
-                  CImg<ulongT>::sequence(_cimg_mp_size(arg2),arg2 + 1,arg2 + (ulongT)_cimg_mp_size(arg2)).
-                    move_to(l_opcode);
-                else CImg<ulongT>::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<double>::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' :
@@ -21546,13 +21746,13 @@ 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);
+              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);
@@ -21561,6 +21761,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()'");
@@ -21650,6 +21887,52 @@ 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;
+              }
+
+              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);
+              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;
+              _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;
@@ -22710,6 +22993,62 @@ 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, 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
+              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;
@@ -22745,6 +23084,43 @@ namespace cimg_library {
             break;
           }
 
+          if ((!std::strncmp(ss,"norm",4) && ss4<se1 && (s = std::strchr(ss4,'('))!=0) || // Lp norm (constant p)
+              (!std::strncmp(ss,"hypot(",6) && (s = ss5))) {
+            const bool is_hypot = *ss=='h';
+            _cimg_mp_op(is_hypot?"Function 'hypot()'":"Function 'norm()'");
+            arg1 = !is_hypot && s!=ss4?compile(ss4,s,depth1,0,block_flags):2;
+            _cimg_mp_check_const_scalar(arg1,0,0);
+            val = mem[arg1];
+            is_sth = true; // Tell if all arguments are constant
+            CImg<ulongT>::vector(0,0,0,arg1).move_to(l_opcode);
+            for (++s; s<se; ++s) {
+              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 (_cimg_mp_is_vector(arg2))
+                CImg<ulongT>::sequence(_cimg_mp_size(arg2),arg2 + 1,arg2 + (ulongT)_cimg_mp_size(arg2)).
+                  move_to(l_opcode);
+              else CImg<ulongT>::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<double>::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) ||
@@ -24227,9 +24603,11 @@ namespace cimg_library {
       }
 
       static double mp_c2o(_cimg_math_parser& mp) {
-        mp_check_list(mp,"c2o");
         unsigned int ind = (unsigned int)mp.opcode[2];
-        if (ind!=~0U) ind = (unsigned int)cimg::mod((int)_mp_arg(2),mp.imglist.width());
+        if (ind!=~0U) {
+          mp_check_list(mp,"c2o");
+          ind = (unsigned int)cimg::mod((int)_mp_arg(2),mp.imglist.width());
+        }
         const CImg<T> &img = ind==~0U?mp.imgin:mp.imglist[ind];
         const int
           x = (int)_mp_arg(3),
@@ -24248,7 +24626,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) {
@@ -24388,7 +24766,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);
@@ -24578,15 +24956,15 @@ namespace cimg_library {
       }
 
       static double mp_da_back_or_pop(_cimg_math_parser& mp) {
-        const bool is_pop = (bool)mp.opcode[4];
-        const char *const s_op = is_pop?"da_pop":"da_back";
+        const bool is_pop_heap = mp.opcode[4]==2, is_pop = (bool)mp.opcode[4];
+        const char *const s_op = is_pop_heap?"da_pop_heap":is_pop?"da_pop":"da_back";
         mp_check_list(mp,s_op);
         const unsigned int
           dim = (unsigned int)mp.opcode[2],
           ind = (unsigned int)cimg::mod((int)_mp_arg(3),mp.imglist.width());
         double *const ptrd = &_mp_arg(1) + (dim>1?1:0);
         CImg<T> &img = mp.imglist[ind];
-        int siz = img?(int)img[img._height - 1]:0;
+        int siz = img?(int)cimg::float2uint((float)img[img._height - 1]):0;
         if (img && (img._width!=1 || img._depth!=1 || siz<0 || siz>img.height() - 1))
           throw CImgArgumentException("[" cimg_appname "_math_parser] CImg<%s>: Function '%s()': "
                                       "Specified image #%u of size (%d,%d,%d,%d) cannot be used as dynamic array%s.",
@@ -24597,15 +24975,41 @@ namespace cimg_library {
           throw CImgArgumentException("[" cimg_appname "_math_parser] CImg<%s>: Function '%s()': "
                                       "Specified dynamic array #%u contains no elements.",
                                       mp.imgout.pixel_type(),s_op,ind);
+        const int siz1 = siz - 1;
+        if (is_pop_heap) { // Heapify-down
+          if (dim==1) cimg::swap(img[0],img[siz1]);
+          else {
+            T *ptr0 = img.data(), *ptr1 = img.data(0,siz1);
+            cimg_forC(img,c) { cimg::swap(*ptr0,*ptr1); ptr0+=img._height; ptr1+=img._height; }
+          }
+          int index = 0;
+          while (true) {
+            const int child_left = 2*index + 1, child_right = child_left + 1;
+            int smallest = index;
+            if (child_left<siz1 && img[child_left]<img[smallest]) smallest = child_left;
+            if (child_right<siz1 && img[child_right]<img[smallest]) smallest = child_right;
+            if (smallest!=index) {
+              if (dim==1) cimg::swap(img[index],img[smallest]);
+              else {
+                T *ptr0 = img.data(0,index), *ptr1 = img.data(0,smallest);
+                cimg_forC(img,c) { cimg::swap(*ptr0,*ptr1); ptr0+=img._height; ptr1+=img._height; }
+              }
+              index = smallest;
+            } else break;
+          }
+        }
 
         double ret = cimg::type<double>::nan();
-        if (dim==1) ret = img[siz - 1]; // Scalar element
-        else cimg_forC(img,c) ptrd[c] = img(0,siz - 1,0,c); // Vector element
+        if (dim==1) ret = img[siz1]; // Scalar element
+        else {
+          const T *ptrs = img.data(0,siz1);
+          cimg_forC(img,c) { ptrd[c] = *ptrs; ptrs+=img._height; } // Vector element
+        }
         if (is_pop) { // Remove element from array
           --siz;
-          if (img.height()>32 && siz<2*img.height()/3) // Reduce size of dynamic array
+          if (img.height()>32 && siz<img.height()/8) // Reduce size of dynamic array
             img.resize(1,std::max(2*siz + 1,32),1,-100,0);
-          img[img._height - 1] = (T)siz;
+          img[img._height - 1] = (T)cimg::uint2float(siz);
         }
         return ret;
       }
@@ -24627,7 +25031,8 @@ namespace cimg_library {
       }
 
       static double mp_da_insert_or_push(_cimg_math_parser& mp) {
-        const char *const s_op = mp.opcode[3]==~0U?"da_push":"da_insert";
+        const bool is_push_heap = mp.opcode[3]==~0U - 1, is_push = mp.opcode[3]>=~0U - 1;
+        const char *const s_op = is_push_heap?"da_push_heap":is_push?"da_push":"da_insert";
         mp_check_list(mp,s_op);
         const unsigned int
           dim = (unsigned int)mp.opcode[4],
@@ -24637,7 +25042,7 @@ namespace cimg_library {
         CImg<T> &img = mp.imglist[ind];
         const int
           siz = img?(int)cimg::float2uint((float)img[img._height - 1]):0,
-          pos0 = mp.opcode[3]==~0U?siz:(int)_mp_arg(3),
+          pos0 = is_push?siz:(int)_mp_arg(3),
           pos = pos0<0?pos0 + siz:pos0;
 
         if (img && _dim!=img._spectrum)
@@ -24662,11 +25067,32 @@ namespace cimg_library {
           cimg_forC(img,c) std::memmove(img.data(0,pos + nb_elts,0,c),img.data(0,pos,0,c),(siz - pos)*sizeof(T));
 
         if (!dim) // Scalar or vector1() elements
-          for (unsigned int k = 0; k<nb_elts; ++k) img[pos + k] = (T)_mp_arg(6 + k);
+          for (unsigned int k = 0; k<nb_elts; ++k) {
+            int index = pos + k;
+            img[index] = (T)_mp_arg(6 + k);
+            if (is_push_heap) while (index>0) { // Heapify-up
+                const int index_parent = (index - 1)/2;
+                if (img[index]<img[index_parent]) {
+                  cimg::swap(img[index],img[index_parent]);
+                  index = index_parent; }
+                else break;
+              }
+          }
         else // vectorN() elements, with N>1
           for (unsigned int k = 0; k<nb_elts; ++k) {
-            double *ptrs = &_mp_arg(6 + k) + 1;
-            cimg_forC(img,c) img(0,pos + k,0,c) = ptrs[c];
+            int index = pos + k;
+            const double *const ptrs = &_mp_arg(6 + k) + 1;
+            T *ptrd = img.data(0,index);
+            cimg_forC(img,c) { *ptrd = ptrs[c]; ptrd+=img._height; }
+            if (is_push_heap) while (index>0) { // Heapify-up
+                const int index_parent = (index - 1)/2;
+                if (img[index]<img[index_parent]) {
+                  T *ptr0 = img.data(0,index), *ptr1 = img.data(0,index_parent);
+                  cimg_forC(img,c) { cimg::swap(*ptr0,*ptr1); ptr0+=img._height; ptr1+=img._height; }
+                  index = index_parent;
+                }
+                else break;
+              }
           }
         img[img._height - 1] = (T)cimg::uint2float(siz + nb_elts);
         return cimg::type<double>::nan();
@@ -24691,7 +25117,7 @@ namespace cimg_library {
           start0 = mp.opcode[3]==~0U?siz - 1:_mp_arg(3),
           end0 = mp.opcode[4]==~0U?start0:_mp_arg(4),
           start = start0<0?start0 + siz:start0,
-                         end = end0<0?end0 + siz:end0;
+          end = end0<0?end0 + siz:end0;
         if (start<0 || start>=siz || end<0 || end>=siz || start>end)
           throw CImgArgumentException("[" cimg_appname "_math_parser] CImg<%s>: Function 'da_remove()': "
                                       "Invalid starting (%d) and ending (%d) positions "
@@ -24700,7 +25126,7 @@ namespace cimg_library {
         if (end<siz - 1) // Move remaining data in dynamic array
           cimg_forC(img,c) std::memmove(img.data(0,start,0,c),img.data(0,end + 1,0,c),(siz - 1 - end)*sizeof(T));
         siz-=end - start + 1;
-        if (img.height()>32 && siz<2*img.height()/3) // Reduce size of dynamic array
+        if (img.height()>32 && siz<img.height()/8) // Reduce size of dynamic array
           img.resize(1,std::max(2*siz + 1,32),1,-100,0);
         img[img._height - 1] = (T)cimg::uint2float(siz);
         return cimg::type<double>::nan();
@@ -24938,10 +25364,12 @@ namespace cimg_library {
       }
 
       static double mp_ellipse(_cimg_math_parser& mp) {
-        mp_check_list(mp,"ellipse");
         const unsigned int i_end = (unsigned int)mp.opcode[2];
         unsigned int ind = (unsigned int)mp.opcode[3];
-        if (ind!=~0U) ind = (unsigned int)cimg::mod((int)_mp_arg(3),mp.imglist.width());
+        if (ind!=~0U) {
+          mp_check_list(mp,"ellipse");
+          ind = (unsigned int)cimg::mod((int)_mp_arg(3),mp.imglist.width());
+        }
         CImg<T> &img = ind==~0U?mp.imgout:mp.imglist[ind];
         CImg<T> color(img._spectrum,1,1,1,0);
         bool is_invalid_arguments = false, is_outlined = false;
@@ -25288,8 +25716,10 @@ namespace cimg_library {
       }
 
       static double mp_i(_cimg_math_parser& mp) {
-        return (double)mp.imgin.atXYZC((int)mp.mem[_cimg_mp_slot_x],(int)mp.mem[_cimg_mp_slot_y],
-                                       (int)mp.mem[_cimg_mp_slot_z],(int)mp.mem[_cimg_mp_slot_c],(T)0);
+        if (mp.imgin)
+          return (double)mp.imgin((int)mp.mem[_cimg_mp_slot_x],(int)mp.mem[_cimg_mp_slot_y],
+                                  (int)mp.mem[_cimg_mp_slot_z],(int)mp.mem[_cimg_mp_slot_c]);
+        return 0;
       }
 
       static double mp_if(_cimg_math_parser& mp) {
@@ -27103,9 +27533,11 @@ namespace cimg_library {
       }
 
       static double mp_o2c(_cimg_math_parser& mp) {
-        mp_check_list(mp,"o2c");
         unsigned int ind = (unsigned int)mp.opcode[2];
-        if (ind!=~0U) ind = (unsigned int)cimg::mod((int)_mp_arg(2),mp.imglist.width());
+        if (ind!=~0U) {
+          mp_check_list(mp,"o2c");
+          ind = (unsigned int)cimg::mod((int)_mp_arg(2),mp.imglist.width());
+        }
         const CImg<T> &img = ind==~0U?mp.imgin:mp.imglist[ind];
         longT offset = (longT)_mp_arg(3);
         double *ptrd = &_mp_arg(1) + 1;
@@ -27944,6 +28376,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;
@@ -27996,6 +28451,66 @@ 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
+          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,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));
+          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_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,
@@ -28139,10 +28654,56 @@ 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);
       }
 
+      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];
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          min_value = _mp_arg(4),
+          max_value = _mp_arg(5),
+          constant_case_ratio = _mp_arg(6);
+        CImg<doubleT>(ptrd,siz,1,1,1,true) = CImg<doubleT>(ptrs,siz,1,1,1,true).
+          get_normalize(min_value,max_value,constant_case_ratio);
+        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;
@@ -28212,6 +28773,25 @@ namespace cimg_library {
         return p?cimg::abs(val):(val!=0);
       }
 
+      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) {
         const bool print_string = (bool)mp.opcode[4];
         cimg_pragma_openmp(critical(mp_vector_print))
@@ -28249,20 +28829,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();
       }
@@ -28302,8 +28901,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();
       }
 
@@ -28333,6 +28932,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],
+          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),
+          boundary_conditions = (int)_mp_arg(14);
+        const double
+          *const ptrs = &_mp_arg(2) + 1,
+          *const ptrw = &_mp_arg(7) + 1;
+        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();
+      }
+
 #define _cimg_mp_vfunc(func) \
       const longT sizd = (longT)mp.opcode[2];\
       const unsigned int nbargs = (unsigned int)(mp.opcode[3] - 4)/2; \
@@ -31297,7 +31919,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;
@@ -31317,7 +31939,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) {
@@ -31325,7 +31947,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);
@@ -31338,7 +31960,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);
@@ -33265,6 +33887,61 @@ namespace cimg_library {
       return (+*this).rand(val_min,val_max);
     }
 
+    //! Fill image with random values following specified distribution and range.
+    /**
+       \param val_min Minimal authorized random value.
+       \param val_max Maximal authorized random value.
+       \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 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(),
+        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);
+
+      // Compute inverse cdf.
+      CImg<tfloat> cdf = pdf.get_max((t)0).cumulate(), icdf(prec);
+      unsigned int k = 0;
+      tfloat p = 0;
+      cdf*=(prec - 1)/cdf.back();
+      cimg_forX(icdf,x) {
+        while (k<siz && (!(p = cdf[k]) || p<x)) ++k;
+        if (k>=siz) { while (x<icdf.width()) icdf[x++] = val_max; break; }
+        icdf[x] = val_min + k*delta_over_siz1;
+      }
+
+      // Generate random numbers.
+      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)prec,&rng),
+            ind = _ind==prec?0:_ind;
+          _data[off] = (T)icdf[ind];
+        }
+        cimg::srand(rng);
+      }
+      return *this;
+    }
+
+    //! Fill image with random values following specified distribution and range \newinstance.
+    template<typename t>
+    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.
     /**
        \param y Rounding precision.
@@ -33301,13 +33978,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)) {
@@ -33318,7 +33995,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;
@@ -33335,7 +34012,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;
@@ -33344,7 +34021,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(); }
@@ -33356,7 +34033,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;
@@ -33384,8 +34061,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;
@@ -33405,8 +34082,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.
@@ -33810,7 +34487,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
@@ -33822,14 +34499,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) "
@@ -33842,7 +34519,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;
@@ -33868,7 +34545,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);
           }
@@ -33903,7 +34580,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);
@@ -33948,7 +34625,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);
             }
@@ -33982,7 +34659,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; }
               }
@@ -34005,7 +34682,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;
@@ -34023,7 +34700,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);
             }
           }
@@ -34046,7 +34723,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);
@@ -34066,7 +34743,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; }
               }
@@ -40454,6 +41131,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;
     }
@@ -44474,7 +45152,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)<patch_penalization)
+      if (!allow_identity && cimg::hypot((float)x1 - x2,(float)y1 - y2)<patch_penalization)
         return cimg::type<float>::inf();
       const T *p1 = img1.data(x1*psizec,y1), *p2 = img2.data(x2*psizec,y2);
       const unsigned int psizewc = psizew*psizec;
@@ -55977,15 +56655,21 @@ namespace cimg_library {
       bool endian = false;
       unsigned int header_size;
       cimg::fread(&header_size,1,nfile_header);
-      if (!header_size)
-        throw CImgIOException(_cimg_instance
-                              "load_analyze(): Invalid zero-size header in file '%s'.",
-                              cimg_instance,
-                              filename?filename:"(FILE*)");
       if (header_size>=4096) { endian = true; cimg::invert_endianness(header_size); }
+      if (header_size<128)
+        throw CImgIOException(_cimg_instance
+                              "load_analyze(): Invalid header size (%u) specified in file '%s'.",
+                              cimg_instance,
+                              header_size,filename?filename:"(FILE*)");
 
       unsigned char *const header = new unsigned char[header_size];
-      cimg::fread(header + 4,header_size - 4,nfile_header);
+      const size_t header_size_read = cimg::fread(header + 4,header_size - 4,nfile_header);
+      if (header_size_read!=header_size - 4)
+        throw CImgIOException(_cimg_instance
+                              "load_analyze(): Cannot read header (of size %u) in file '%s'.",
+                              cimg_instance,
+                              header_size,filename?filename:"(FILE*)");
+
       if (!file && nfile_header!=nfile) cimg::fclose(nfile_header);
       if (endian) {
         cimg::invert_endianness((short*)(header + 40),5);
@@ -65090,7 +65774,7 @@ namespace cimg_library {
     CImg<charT> __display() const {
       CImg<charT> res, str;
       cimglist_for(*this,l) {
-        CImg<charT>::string((char*)_data[l]).move_to(str);
+        CImg<charT>::string((char*)_data[l]._data).move_to(str);
         if (l!=width() - 1) {
           str.resize(str._width + 1,1,1,1,0);
           str[str._width - 2] = ',';