Optimization pass

source/dcv/features/corner/harris.d

@@ -9,7 +9,10 @@ License: $(LINK3 http://www.boost.org/LICENSE_1_0.txt, Boost Software License -
 */
 module dcv.features.corner.harris;
 
+import std.parallelism : parallel, taskPool, TaskPool;
+
 import mir.ndslice;
+import mir.ndslice.algorithm : ndEach, Yes;
 
 import dcv.core.utils : emptySlice;
 import dcv.imgproc.filter : calcPartialDerivatives;
@@ -31,12 +34,23 @@ Returns:
  */
 Slice!(2, OutputType*) harrisCorners(InputType, OutputType = InputType)(Slice!(2,
         InputType*) image, in uint winSize = 3, in float k = 0.64f, in float gauss = 0.84f, Slice!(2,
-        OutputType*) prealloc = emptySlice!(2, OutputType))
+        OutputType*) prealloc = emptySlice!(2, OutputType), TaskPool pool = taskPool)
+in
 {
-
-    HarrisDetector det;
-    det.k = k;
-    return calcCorners(image, winSize, gauss, prealloc, det);
+    assert(!image.empty, "Empty image given.");
+    assert(winSize % 2 != 0, "Kernel window size has to be odd.");
+    assert(gauss > 0.0, "Gaussian sigma value has to be greater than 0.");
+    assert(k > 0.0, "K value has to be greater than 0.");
+}
+body
+{
+    if (prealloc.shape != image.shape)
+    {
+        prealloc = uninitializedSlice!OutputType(image.shape);
+    }
+    HarrisDetector detector;
+    detector.k = k;
+    return calcCorners(image, winSize, gauss, prealloc, detector, pool);
 }
 
 /**
@@ -55,10 +69,22 @@ Returns:
  */
 Slice!(2, OutputType*) shiTomasiCorners(InputType, OutputType = InputType)(Slice!(2,
         InputType*) image, in uint winSize = 3, in float gauss = 0.84f, Slice!(2,
-        OutputType*) prealloc = emptySlice!(2, OutputType))
+        OutputType*) prealloc = emptySlice!(2, OutputType), TaskPool pool = taskPool)
+in
 {
-    ShiTomasiDetector det;
-    return calcCorners(image, winSize, gauss, prealloc, det);
+    assert(!image.empty, "Empty image given.");
+    assert(winSize % 2 != 0, "Kernel window size has to be odd.");
+    assert(gauss > 0.0, "Gaussian sigma value has to be greater than 0.");
+}
+body
+{
+    if (prealloc.shape != image.shape)
+    {
+        prealloc = uninitializedSlice!OutputType(image.shape);
+    }
+
+    ShiTomasiDetector detector;
+    return calcCorners(image, winSize, gauss, prealloc, detector, pool);
 }
 
 unittest
@@ -115,15 +141,15 @@ struct HarrisDetector
 {
     float k;
 
-    float opCall(float r1, float r2, float r3)
+    @nogc nothrow float opCall(float r1, float r2, float r3)
     {
         return (((r1 * r1) - (r2 * r3)) - k * ((r1 + r3) * r1 + r3));
     }
 }
 
 struct ShiTomasiDetector
 {
-    float opCall(float r1, float r2, float r3)
+    @nogc nothrow float opCall(float r1, float r2, float r3)
     {
         import std.math : sqrt;
 
@@ -132,64 +158,37 @@ struct ShiTomasiDetector
 }
 
 Slice!(2, OutputType*) calcCorners(Detector, InputType, OutputType)(Slice!(2, InputType*) image,
-        uint winSize, float gaussSigma, Slice!(2, OutputType*) prealloc, Detector detector)
+        uint winSize, float gaussSigma, Slice!(2, OutputType*) prealloc, Detector detector, TaskPool pool)
 {
+    // TODO: implement gaussian weighting!
 
-    import std.math : exp, PI;
-    import std.array : uninitializedArray;
-    import std.range : iota;
-    import std.algorithm.iteration : reduce, each;
-    import std.algorithm.comparison : equal;
+    Slice!(2, InputType*) fx, fy;
+    calcPartialDerivatives(image, fx, fy);
 
-    assert(!image.empty);
+    auto windowPack = assumeSameStructure!("corners", "fx", "fy")(prealloc, fx, fy).windows(winSize, winSize);
 
-    if (!prealloc.shape[].equal(image.shape[]))
+    foreach (windowRow; pool.parallel(windowPack))
     {
-        prealloc = uninitializedArray!(OutputType[])(image.shape[].reduce!"a*b").sliced(image.shape);
+        windowRow.ndEach!(win => calcCornersImpl(win, detector));
     }
-    prealloc[] = cast(OutputType)0;
 
-    auto rows = image.length!0;
-    auto cols = image.length!1;
+    return prealloc;
+}
 
-    Slice!(2, InputType*) fx, fy;
-    calcPartialDerivatives(image, fx, fy);
+@nogc nothrow
+void calcCornersImpl(Window, Detector)(ref Window window, ref Detector detector)
+{
+    float r1 = 0.0f, r2 = 0.0f, r3 = 0.0f;
+    float winSqr = float(window.length!0);
+    winSqr *= winSqr;
 
-    auto winSqr = winSize ^^ 2;
-    auto winHalf = winSize / 2;
+    window.ndEach!((pack) { auto gx = pack.fx; auto gy = pack.fy; r1 += gx * gx; r2 += gx * gy; r3 += gy * gy;  });
 
-    float R; // Score value
-    float w, gx, gy, r1, r2, r3;
-    float gaussMul = 1.0f / (2.0f * PI * gaussSigma);
-    float gaussDel = 2.0f * (gaussSigma ^^ 2);
+    r1 = (r1 / winSqr) * 0.5;
+    r2 /= winSqr;
+    r3 = (r3 / winSqr) * 0.5;
 
-    foreach (i; winHalf.iota(rows - winHalf))
-    {
-        foreach (j; winHalf.iota(cols - winHalf))
-        {
-            r1 = 0.;
-            r2 = 0.;
-            r3 = 0.;
-            for (int cr = cast(int)(i - winHalf); cr < i + winHalf; cr++)
-            {
-                for (int cc = cast(int)(j - winHalf); cc < j + winHalf; cc++)
-                {
-                    w = 1.0f; //gaussMul * exp(-((cast(float)cr - i)^^2 + (cast(float)cc - j)^^2) / gaussDel);
-
-                    gx = fx[cr, cc];
-                    gy = fy[cr, cc];
-                    r1 += w * (gx * gx);
-                    r2 += w * (gx * gy);
-                    r3 += w * (gy * gy);
-                }
-            }
-            r1 = (r1 / winSqr) * 0.5;
-            r2 /= winSqr;
-            r3 = (r3 / winSqr) * 0.5;
-            R = detector(r1, r2, r3);
-            if (R > 0)
-                prealloc[i, j] = cast(OutputType)R;
-        }
-    }
-    return prealloc;
+    auto r = detector(r1, r2, r3);
+    if (r > 0)
+        window[$ / 2, $ / 2].corners = r;
 }

source/dcv/features/utils.d

@@ -49,35 +49,27 @@ return:
 Dynamic array of size_t[2], as in array of 2D points, of corner reponses 
 which fit the given criteria.
 */
-size_t[2][] extractCorners(T)(Slice!(2, T*) cornerResponse, int count = -1, T threshold = cast(T)0)@trusted pure nothrow 
-        if (isNumeric!T)
+pure nothrow auto extractCorners(T)(Slice!(2, T*) cornerResponse, int count = -1, T threshold = 0)
+    if (isNumeric!T)
 {
-    import std.algorithm : sort, map, min;
+    import std.algorithm.sorting : sort;
+    import std.algorithm.iteration : map, filter;
+    import std.range : take;
     import std.array : array;
-    import std.range : take, iota;
-
-    size_t[2][float] features;
 
     if (cornerResponse.empty)
     {
         return null;
     }
 
-    auto rows = cornerResponse.length!0;
-    auto cols = cornerResponse.length!1;
-
-    foreach (r; rows.iota)
-    {
-        foreach (c; cols.iota)
-        {
-            auto res = cornerResponse[r, c];
-            if (res > threshold)
-                features[res] = [r, c];
-        }
-    }
-
-    return sort!"a > b"(features.keys).take(count > 0 ? min(count, features.length) : features.length)
-        .map!(a => features[a]).array;
+    return assumeSameStructure!("indices", "value")(indexSlice(cornerResponse.shape), cornerResponse)
+        .byElement
+        .filter!(p => p.value > threshold)
+        .array
+        .sort!((a, b) => a.value > b.value)
+        .take(count)
+        .map!(p => p.indices)
+        .array;
 }
 
 unittest