diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc
index aa1b90278..5f5e07139 100644
--- a/rtengine/demosaic_algos.cc
+++ b/rtengine/demosaic_algos.cc
@@ -21,6 +21,7 @@
 
 #include "rawimagesource.h"
 #include "rawimagesource_i.h"
+#include "jaggedarray.h"
 #include "median.h"
 #include "rawimage.h"
 #include "mytime.h"
@@ -94,40 +95,12 @@ void RawImageSource::eahd_demosaic ()
 
     // end of cielab preparation
 
-    float* rh[3];
-    float* gh[4];
-    float* bh[3];
-    float* rv[3];
-    float* gv[4];
-    float* bv[3];
-    float* lLh[3];
-    float* lah[3];
-    float* lbh[3];
-    float* lLv[3];
-    float* lav[3];
-    float* lbv[3];
-    float* homh[3];
-    float* homv[3];
-
-    for (int i = 0; i < 4; i++) {
-        gh[i] = new float[W];
-        gv[i] = new float[W];
-    }
-
-    for (int i = 0; i < 3; i++) {
-        rh[i] = new float[W];
-        bh[i] = new float[W];
-        rv[i] = new float[W];
-        bv[i] = new float[W];
-        lLh[i] = new float[W];
-        lah[i] = new float[W];
-        lbh[i] = new float[W];
-        lLv[i] = new float[W];
-        lav[i] = new float[W];
-        lbv[i] = new float[W];
-        homh[i] = new float[W];
-        homv[i] = new float[W];
-    }
+    const 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),
+            lLv (W, 3), lav (W, 3), lbv (W, 3),
+            homh (W, 3), homv (W, 3);
 
     // interpolate first two lines
     interpolate_row_g (gh[0], gv[0], 0);
@@ -311,27 +284,9 @@ void RawImageSource::eahd_demosaic ()
             }
         }
 
-    freeArray2<float>(rh, 3);
-    freeArray2<float>(gh, 4);
-    freeArray2<float>(bh, 3);
-    freeArray2<float>(rv, 3);
-    freeArray2<float>(gv, 4);
-    freeArray2<float>(bv, 3);
-    freeArray2<float>(lLh, 3);
-    freeArray2<float>(lah, 3);
-    freeArray2<float>(lbh, 3);
-    freeArray2<float>(homh, 3);
-    freeArray2<float>(lLv, 3);
-    freeArray2<float>(lav, 3);
-    freeArray2<float>(lbv, 3);
-    freeArray2<float>(homv, 3);
-
     // Interpolate R and B
     for (int i = 0; i < H; i++) {
-        if (i == 0)
-            // rm, gm, bm must be recovered
-            //interpolate_row_rb_mul_pp (red, blue, NULL, green[i], green[i+1], i, rm, gm, bm, 0, W, 1);
-        {
+        if (i == 0) {
             interpolate_row_rb_mul_pp (red[i], blue[i], NULL, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1);
         } else if (i == H - 1) {
             interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], NULL, i, 1.0, 1.0, 1.0, 0, W, 1);
@@ -545,7 +500,7 @@ void RawImageSource::hphd_demosaic ()
         plistener->setProgress (0.0);
     }
 
-    float** hpmap = allocArray< float >(W, H, true);
+    const JaggedArray<float> hpmap (W, H, true);
 
 #ifdef _OPENMP
     #pragma omp parallel
@@ -590,17 +545,13 @@ void RawImageSource::hphd_demosaic ()
 #endif
 
     hphd_green (hpmap);
-    freeArray<float>(hpmap, H);//TODO: seems to cause sigabrt ???  why???
 
     if (plistener) {
         plistener->setProgress (0.66);
     }
 
     for (int i = 0; i < H; i++) {
-        if (i == 0)
-            // rm, gm, bm must be recovered
-            //interpolate_row_rb_mul_pp (red, blue, NULL, green[i], green[i+1], i, rm, gm, bm, 0, W, 1);
-        {
+        if (i == 0) {
             interpolate_row_rb_mul_pp (red[i], blue[i], NULL, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1);
         } else if (i == H - 1) {
             interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], NULL, i, 1.0, 1.0, 1.0, 0, W, 1);
diff --git a/rtengine/hlmultipliers.cc b/rtengine/hlmultipliers.cc
index 8f7c5ef85..35e1d3e87 100644
--- a/rtengine/hlmultipliers.cc
+++ b/rtengine/hlmultipliers.cc
@@ -20,6 +20,7 @@
 #include <cstdio>
 #include "rawimagesource.h"
 #include "rawimagesource_i.h"
+#include "jaggedarray.h"
 #include "../rtgui/options.h"
 
 namespace rtengine
@@ -321,7 +322,7 @@ void RawImageSource::updateHLRecoveryMap_ColorPropagation ()
     int** rec[3];
 
     for (int i = 0; i < 3; i++) {
-        rec[i] = allocArray<int> (dw, dh);
+        rec[i] = allocJaggedArray<int> (dw, dh);
     }
 
     float* reds[HR_SCALE];
@@ -333,10 +334,10 @@ void RawImageSource::updateHLRecoveryMap_ColorPropagation ()
     }
 
     if (needhr) {
-        freeArray<char>(needhr, H);
+        freeJaggedArray<char>(needhr);
     }
 
-    needhr = allocArray<char> (W, H);
+    needhr = allocJaggedArray<char> (W, H);
 
     for (int i = 0; i < dh; i++) {
         for (int j = 0; j < HR_SCALE; j++) {
@@ -413,14 +414,14 @@ void RawImageSource::updateHLRecoveryMap_ColorPropagation ()
     hlmultipliers (rec, max_3, dh, dw);
 
     if (hrmap[0] != NULL) {
-        freeArray<float> (hrmap[0], dh);
-        freeArray<float> (hrmap[1], dh);
-        freeArray<float> (hrmap[2], dh);
+        freeJaggedArray<float> (hrmap[0]);
+        freeJaggedArray<float> (hrmap[1]);
+        freeJaggedArray<float> (hrmap[2]);
     }
 
-    hrmap[0] = allocArray<float> (dw, dh);
-    hrmap[1] = allocArray<float> (dw, dh);
-    hrmap[2] = allocArray<float> (dw, dh);
+    hrmap[0] = allocJaggedArray<float> (dw, dh);
+    hrmap[1] = allocJaggedArray<float> (dw, dh);
+    hrmap[2] = allocJaggedArray<float> (dw, dh);
 
     for (int i = 0; i < dh; i++)
         for (int j = 0; j < dw; j++) {
@@ -431,9 +432,9 @@ void RawImageSource::updateHLRecoveryMap_ColorPropagation ()
 
     delete ds;
 
-    freeArray<int> (rec[0], dh);
-    freeArray<int> (rec[1], dh);
-    freeArray<int> (rec[2], dh);
+    freeJaggedArray<int> (rec[0]);
+    freeJaggedArray<int> (rec[1]);
+    freeJaggedArray<int> (rec[2]);
 }
 
 }
diff --git a/rtengine/jaggedarray.h b/rtengine/jaggedarray.h
new file mode 100644
index 000000000..86d556c7f
--- /dev/null
+++ b/rtengine/jaggedarray.h
@@ -0,0 +1,70 @@
+#ifndef JAGGEDARRAY_H
+#define JAGGEDARRAY_H
+
+namespace rtengine
+{
+
+// These emulate a jagged array, but use only 2 allocations instead of 1 + W.
+
+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>
+class JaggedArray
+{
+public:
+    JaggedArray (const int W, const int H, const bool initZero = false)
+    {
+        a = allocJaggedArray<T> (W, H, initZero);
+    }
+    ~JaggedArray ()
+    {
+        if (a) {
+            freeJaggedArray<T> (a);
+            a = nullptr;
+        }
+    }
+
+    JaggedArray (const JaggedArray&) = delete;
+    JaggedArray& operator= (const JaggedArray&) = delete;
+
+public:
+    operator T** const () const
+    {
+        return a;
+    }
+
+private:
+    T** a;
+
+};
+
+// Declared but not defined to prevent
+// explicitly freeing a JaggedArray<T> implicitly cast to T**.
+template<class T>
+void freeJaggedArray (JaggedArray<T>&);
+
+} // rtengine
+
+#endif // JAGGEDARRAY_H
diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc
index 0d8c9abd3..d1d4da54f 100644
--- a/rtengine/rawimagesource.cc
+++ b/rtengine/rawimagesource.cc
@@ -22,6 +22,7 @@
 #include "rtengine.h"
 #include "rawimagesource.h"
 #include "rawimagesource_i.h"
+#include "jaggedarray.h"
 #include "median.h"
 #include "rawimage.h"
 #include "mytime.h"
@@ -480,15 +481,11 @@ RawImageSource::~RawImageSource ()
 
     if (hrmap[0] != NULL) {
         int dh = H / HR_SCALE;
-        freeArray<float>(hrmap[0], dh);
-        freeArray<float>(hrmap[1], dh);
-        freeArray<float>(hrmap[2], dh);
+        freeJaggedArray<float>(hrmap[0]);
+        freeJaggedArray<float>(hrmap[1]);
+        freeJaggedArray<float>(hrmap[2]);
     }
 
-    //if (needhr)
-    //    freeArray<char>(needhr, H);
-    //if (hpmap)
-    //    freeArray<char>(hpmap, H);
     if (camProfile) {
         cmsCloseProfile (camProfile);
     }
diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h
index a1bac4363..bef4ad4e2 100644
--- a/rtengine/rawimagesource.h
+++ b/rtengine/rawimagesource.h
@@ -31,37 +31,6 @@
 namespace rtengine
 {
 
-// these two functions "simulate" and jagged array, but just use two allocs
-template<class T> T** allocArray (int W, int H, bool initZero = false)
-{
-
-    T** t = new T*[H];
-    t[0] = new T[H * W];
-
-    if (initZero) {
-        memset(t[0], 0, sizeof(T)*W * H);
-    }
-
-    for (int i = 1; i < H; i++) {
-        t[i] = t[i - 1] + W;
-    }
-
-    return t;
-}
-
-template<class T> void freeArray (T** a, int H)
-{
-
-    delete [] a[0];
-    delete [] a;
-}
-
-
-template<class T> void freeArray2 (T** a, int H)
-{
-    delete [] a[0];
-}
-
 class RawImageSource : public ImageSource
 {
 
diff --git a/rtengine/shmap.cc b/rtengine/shmap.cc
index 4f93c80d3..1783a6bc8 100644
--- a/rtengine/shmap.cc
+++ b/rtengine/shmap.cc
@@ -21,6 +21,7 @@
 #include "rtengine.h"
 #include "rt_math.h"
 #include "rawimagesource.h"
+#include "jaggedarray.h"
 #undef THREAD_PRIORITY_NORMAL
 #include "opthelper.h"
 
@@ -114,7 +115,7 @@ void SHMap::update (Imagefloat* img, double radius, double lumi[3], bool hq, int
         rangefn[lutSize - 1] = 1e-15f;
 
         // We need one temporary buffer
-        float ** buffer = allocArray<float> (W, H);
+        const 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
@@ -157,23 +158,6 @@ void SHMap::update (Imagefloat* img, double radius, double lumi[3], bool hq, int
         }
 
         dirpyr_shmap(dirpyrlo[indx], dirpyrlo[1 - indx], W, H, rangefn, level, scale );
-
-        freeArray<float>(buffer, H);
-
-        //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-        /*
-                // anti-alias filtering the result
-        #ifdef _OPENMP
-        #pragma omp for
-        #endif
-                for (int i=0; i<H; i++)
-                    for (int j=0; j<W; j++)
-                        if (i>0 && j>0 && i<H-1 && j<W-1)
-                            map[i][j] = (buffer[i-1][j-1]+buffer[i-1][j]+buffer[i-1][j+1]+buffer[i][j-1]+buffer[i][j]+buffer[i][j+1]+buffer[i+1][j-1]+buffer[i+1][j]+buffer[i+1][j+1])/9;
-                        else
-                            map[i][j] = buffer[i][j];
-        */
-
     }
 
     // update average, minimum, maximum
@@ -262,7 +246,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
-        float ** buffer = allocArray<float> (W, H);
+        const 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
@@ -306,23 +290,6 @@ void SHMap::updateL (float** L, double radius, bool hq, int skip)
         }
 
         dirpyr_shmap(dirpyrlo[indx], dirpyrlo[1 - indx], W, H, rangefn, level, scale );
-
-        freeArray<float>(buffer, H);
-
-        //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-        /*
-                // anti-alias filtering the result
-        #ifdef _OPENMP
-        #pragma omp for
-        #endif
-                for (int i=0; i<H; i++)
-                    for (int j=0; j<W; j++)
-                        if (i>0 && j>0 && i<H-1 && j<W-1)
-                            map[i][j] = (buffer[i-1][j-1]+buffer[i-1][j]+buffer[i-1][j+1]+buffer[i][j-1]+buffer[i][j]+buffer[i][j+1]+buffer[i+1][j-1]+buffer[i+1][j]+buffer[i+1][j+1])/9;
-                        else
-                            map[i][j] = buffer[i][j];
-        */
-
     }
 
     // update average, minimum, maximum