diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc
index 368a8edde..03a292c11 100644
--- a/rtengine/demosaic_algos.cc
+++ b/rtengine/demosaic_algos.cc
@@ -92,7 +92,7 @@ void RawImageSource::eahd_demosaic ()
 
     // end of cielab preparation
 
-    const JaggedArray<float>
+    JaggedArray<float>
     rh (W, 3), gh (W, 4), bh (W, 3),
     rv (W, 3), gv (W, 4), bv (W, 3),
     lLh (W, 3), lah (W, 3), lbh (W, 3),
@@ -497,7 +497,7 @@ void RawImageSource::hphd_demosaic ()
         plistener->setProgress (0.0);
     }
 
-    const JaggedArray<float> hpmap (W, H, true);
+    JaggedArray<float> hpmap (W, H, true);
 
 #ifdef _OPENMP
     #pragma omp parallel
@@ -3943,7 +3943,7 @@ void RawImageSource::rcd_demosaic()
         plistener->setProgressStr(Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), "rcd"));
         plistener->setProgress(0);
     }
-    
+
     int width = W, height = H;
 
     std::vector<float> cfa(width * height);
@@ -3963,15 +3963,15 @@ void RawImageSource::rcd_demosaic()
         plistener->setProgress(0.05);
     }
     // ------------------------------------------------------------------------
-/* RT    
+/* RT
     int row, col, indx, c;
-*/    
+*/
     int w1 = width, w2 = 2 * width, w3 = 3 * width, w4 = 4 * width;
 
     //Tolerance to avoid dividing by zero
     static const float eps = 1e-5, epssq = 1e-10;
 
-/* RT    
+/* RT
     //Gradients
     float N_Grad, E_Grad, W_Grad, S_Grad, NW_Grad, NE_Grad, SW_Grad, SE_Grad;
 
@@ -3981,7 +3981,7 @@ void RawImageSource::rcd_demosaic()
     //Directional discrimination
     //float V_Stat, H_Stat, P_Stat, Q_Stat;
     float VH_Central_Value, VH_Neighbour_Value, PQ_Central_Value, PQ_Neighbour_Value;
-*/    
+*/
     float *VH_Dir, *PQ_Dir;
 
     //Low pass filter
@@ -4001,7 +4001,7 @@ void RawImageSource::rcd_demosaic()
         for (int col = 4, indx = row * width + col; col < width - 4; col++, indx++ ) {
             //Calculate h/v local discrimination
             float V_Stat = max(epssq, - 18.0f  *  cfa[indx] * cfa[indx - w1] - 18.0f * cfa[indx] * cfa[indx + w1] - 36.0f * cfa[indx] * cfa[indx - w2] - 36.0f * cfa[indx] * cfa[indx + w2] + 18.0f * cfa[indx] * cfa[indx - w3] + 18.0f * cfa[indx] * cfa[indx + w3] - 2.0f * cfa[indx] * cfa[indx - w4] - 2.0f * cfa[indx] * cfa[indx + w4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx - w1] * cfa[indx + w1] - 12.0f * cfa[indx - w1] * cfa[indx - w2] + 24.0f * cfa[indx - w1] * cfa[indx + w2] - 38.0f * cfa[indx - w1] * cfa[indx - w3] + 16.0f * cfa[indx - w1] * cfa[indx + w3] + 12.0f * cfa[indx - w1] * cfa[indx - w4] - 6.0f * cfa[indx - w1] * cfa[indx + w4] + 46.0f * cfa[indx - w1] * cfa[indx - w1] + 24.0f * cfa[indx + w1] * cfa[indx - w2] - 12.0f * cfa[indx + w1] * cfa[indx + w2] + 16.0f * cfa[indx + w1] * cfa[indx - w3] - 38.0f * cfa[indx + w1] * cfa[indx + w3] - 6.0f * cfa[indx + w1] * cfa[indx - w4] + 12.0f * cfa[indx + w1] * cfa[indx + w4] + 46.0f * cfa[indx + w1] * cfa[indx + w1] + 14.0f * cfa[indx - w2] * cfa[indx + w2] - 12.0f * cfa[indx - w2] * cfa[indx + w3] - 2.0f * cfa[indx - w2] * cfa[indx - w4] + 2.0f * cfa[indx - w2] * cfa[indx + w4] + 11.0f * cfa[indx - w2] * cfa[indx - w2] - 12.0f * cfa[indx + w2] * cfa[indx - w3] + 2.0f * cfa[indx + w2] * cfa[indx - w4] - 2.0f * cfa[indx + w2] * cfa[indx + w4] + 11.0f * cfa[indx + w2] * cfa[indx + w2] + 2.0f * cfa[indx - w3] * cfa[indx + w3] - 6.0f * cfa[indx - w3] * cfa[indx - w4] + 10.0f * cfa[indx - w3] * cfa[indx - w3] - 6.0f * cfa[indx + w3] * cfa[indx + w4] + 10.0f * cfa[indx + w3] * cfa[indx + w3] + 1.0f * cfa[indx - w4] * cfa[indx - w4] + 1.0f * cfa[indx + w4] * cfa[indx + w4]);
-            
+
             float H_Stat = max(epssq, - 18.0f  *  cfa[indx] * cfa[indx -  1] - 18.0f * cfa[indx] * cfa[indx +  1] - 36.0f * cfa[indx] * cfa[indx -  2] - 36.0f * cfa[indx] * cfa[indx +  2] + 18.0f * cfa[indx] * cfa[indx -  3] + 18.0f * cfa[indx] * cfa[indx +  3] - 2.0f * cfa[indx] * cfa[indx -  4] - 2.0f * cfa[indx] * cfa[indx +  4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx -  1] * cfa[indx +  1] - 12.0f * cfa[indx -  1] * cfa[indx -  2] + 24.0f * cfa[indx -  1] * cfa[indx +  2] - 38.0f * cfa[indx -  1] * cfa[indx -  3] + 16.0f * cfa[indx -  1] * cfa[indx +  3] + 12.0f * cfa[indx -  1] * cfa[indx -  4] - 6.0f * cfa[indx -  1] * cfa[indx +  4] + 46.0f * cfa[indx -  1] * cfa[indx -  1] + 24.0f * cfa[indx +  1] * cfa[indx -  2] - 12.0f * cfa[indx +  1] * cfa[indx +  2] + 16.0f * cfa[indx +  1] * cfa[indx -  3] - 38.0f * cfa[indx +  1] * cfa[indx +  3] - 6.0f * cfa[indx +  1] * cfa[indx -  4] + 12.0f * cfa[indx +  1] * cfa[indx +  4] + 46.0f * cfa[indx +  1] * cfa[indx +  1] + 14.0f * cfa[indx -  2] * cfa[indx +  2] - 12.0f * cfa[indx -  2] * cfa[indx +  3] - 2.0f * cfa[indx -  2] * cfa[indx -  4] + 2.0f * cfa[indx -  2] * cfa[indx +  4] + 11.0f * cfa[indx -  2] * cfa[indx -  2] - 12.0f * cfa[indx +  2] * cfa[indx -  3] + 2.0f * cfa[indx +  2] * cfa[indx -  4] - 2.0f * cfa[indx +  2] * cfa[indx +  4] + 11.0f * cfa[indx +  2] * cfa[indx +  2] + 2.0f * cfa[indx -  3] * cfa[indx +  3] - 6.0f * cfa[indx -  3] * cfa[indx -  4] + 10.0f * cfa[indx -  3] * cfa[indx -  3] - 6.0f * cfa[indx +  3] * cfa[indx +  4] + 10.0f * cfa[indx +  3] * cfa[indx +  3] + 1.0f * cfa[indx -  4] * cfa[indx -  4] + 1.0f * cfa[indx +  4] * cfa[indx +  4]);
 
             VH_Dir[indx] = V_Stat / (V_Stat + H_Stat);
@@ -4023,7 +4023,7 @@ void RawImageSource::rcd_demosaic()
 
 #ifdef _OPENMP
     #pragma omp parallel for
-#endif    
+#endif
     for ( int row = 2; row < height - 2; row++ ) {
         for ( int col = 2 + (FC( row, 0 ) & 1), indx = row * width + col; col < width - 2; col += 2, indx += 2 ) {
 
@@ -4084,7 +4084,7 @@ void RawImageSource::rcd_demosaic()
         plistener->setProgress(0.5);
     }
     // ------------------------------------------------------------------------
-    
+
     /**
     * STEP 4: Populate the red and blue channels
     */
@@ -4154,7 +4154,7 @@ void RawImageSource::rcd_demosaic()
         plistener->setProgress(0.825);
     }
     // -------------------------------------------------------------------------
-        
+
     // Step 4.3: Populate the red and blue channels at green CFA positions
 #ifdef _OPENMP
     #pragma omp parallel for
@@ -4165,7 +4165,7 @@ void RawImageSource::rcd_demosaic()
             // Refined vertical and horizontal local discrimination
             float VH_Central_Value   = VH_Dir[indx];
             float VH_Neighbourhood_Value = 0.25f * ( VH_Dir[indx - w1 - 1] + VH_Dir[indx - w1 + 1] + VH_Dir[indx + w1 - 1] + VH_Dir[indx + w1 + 1] );
-            
+
             float VH_Disc = ( fabs( 0.5f - VH_Central_Value ) < fabs( 0.5f - VH_Neighbourhood_Value ) ) ? VH_Neighbourhood_Value : VH_Central_Value;
 
             for ( int c = 0; c <= 2; c += 2 ) {
diff --git a/rtengine/jaggedarray.h b/rtengine/jaggedarray.h
index 153ed2902..01da776a6 100644
--- a/rtengine/jaggedarray.h
+++ b/rtengine/jaggedarray.h
@@ -19,6 +19,8 @@
  */
 #pragma once
 
+#include <cstring>
+
 #include "noncopyable.h"
 
 namespace rtengine
@@ -26,60 +28,46 @@ namespace rtengine
 
 // These emulate a jagged array, but use only 2 allocations instead of 1 + H.
 
-template<class T>
-inline T** const allocJaggedArray (const int W, const int H, const bool initZero = false)
-{
-    T** const a = new T*[H];
-    a[0] = new T[H * W];
-
-    for (int i = 1; i < H; ++i) {
-        a[i] = a[i - 1] + W;
-    }
-
-    if (initZero) {
-        std::memset(a[0], 0, sizeof(T) * W * H);
-    }
-
-    return a;
-}
-
-template<class T>
-inline void freeJaggedArray (T** const a)
-{
-    delete [] a[0];
-    delete [] a;
-}
-
-template<class T>
+template<typename T>
 class JaggedArray :
     public NonCopyable
 {
 public:
-    JaggedArray (const int W, const int H, const bool initZero = false)
+    JaggedArray(std::size_t width, std::size_t height, bool init_zero = false) :
+        array(
+            [width, height, init_zero]() -> T**
+            {
+                T** const res = new T*[height];
+                res[0] = new T[height * width];
+
+                for (std::size_t i = 1; i < height; ++i) {
+                    res[i] = res[i - 1] + width;
+                }
+
+                if (init_zero) {
+                    std::memset(res[0], 0, sizeof(T) * width * height);
+                }
+
+                return res;
+            }()
+        )
     {
-        a = allocJaggedArray<T> (W, H, initZero);
-    }
-    ~JaggedArray ()
-    {
-        if (a) {
-            freeJaggedArray<T> (a);
-            a = nullptr;
-        }
     }
 
-    operator T** const () const
+    ~JaggedArray ()
     {
-        return a;
+        delete[] array[0];
+        delete[] array;
+    }
+
+    operator T** ()
+    {
+        return array;
     }
 
 private:
-    T** a;
+    T** const array;
 
 };
 
-// Declared but not defined to prevent
-// explicitly freeing a JaggedArray<T> implicitly cast to T**.
-template<class T>
-void freeJaggedArray (JaggedArray<T>&);
-
 } // rtengine
diff --git a/rtengine/shmap.cc b/rtengine/shmap.cc
index b34aac7af..f819bf1bd 100644
--- a/rtengine/shmap.cc
+++ b/rtengine/shmap.cc
@@ -124,7 +124,7 @@ void SHMap::update (Imagefloat* img, double radius, double lumi[3], bool hq, int
         rangefn[lutSize - 1] = 1e-15f;
 
         // We need one temporary buffer
-        const JaggedArray<float> buffer (W, H);
+        JaggedArray<float> buffer (W, H);
 
         // the final result has to be in map
         // for an even number of levels that means: map => buffer, buffer => map
@@ -255,7 +255,7 @@ void SHMap::updateL (float** L, double radius, bool hq, int skip)
         //printf("lut=%d rf5=%f rfm=%f\n thre=%f",lutSize, rangefn[5],rangefn[lutSize-10],thresh );
 
         // We need one temporary buffer
-        const JaggedArray<float> buffer (W, H);
+        JaggedArray<float> buffer (W, H);
 
         // the final result has to be in map
         // for an even number of levels that means: map => buffer, buffer => map