Merge pull request #825 from volgunin/development

Speed optimization: Replace multi-dimentional with jagged arrays in IntegralImage.cs

Author: cesarsouza

Sources/Accord.Imaging/AForge.Imaging/IntegralImage.cs

@@ -11,6 +11,7 @@ namespace Accord.Imaging
     using System;
     using System.Drawing;
     using System.Drawing.Imaging;
+    using Accord.Math;
 
     /// <summary>
     /// Integral image.
@@ -52,12 +53,20 @@ namespace Accord.Imaging
     public class IntegralImage : ICloneable
     {
         /// <summary>
-        /// Intergral image's array.
+        /// Integral image's array.
+        /// </summary>
+        /// 
+        /// <remarks>See remarks to <see cref="Matrix"/> property.</remarks>
+        /// 
+        protected readonly uint[][] matrix = null;
+
+        /// <summary>
+        /// Integral image's array.
         /// </summary>
         /// 
         /// <remarks>See remarks to <see cref="InternalData"/> property.</remarks>
         /// 
-        protected readonly uint[,] integralImage = null;
+        protected uint[,] integralImage = null;
 
         // image's width and height
         private readonly int width;
@@ -79,6 +88,23 @@ public int Height
             get { return height; }
         }
 
+        /// <summary>
+        /// Provides access to internal array keeping integral image data.
+        /// </summary>
+        /// 
+        /// <remarks>
+        /// <para><note>The array should be accessed by [y][x] indexing.</note></para>
+        /// 
+        /// <para><note>The array's size is [<see cref="Height"/>+1, <see cref="Width"/>+1]. The first
+        /// row and column are filled with zeros, what is done for more efficient calculation of
+        /// rectangles' sums.</note></para>
+        /// </remarks>
+        /// 
+        public uint[][] Matrix
+        {
+            get { return matrix; }
+        }
+
         /// <summary>
         /// Provides access to internal array keeping integral image data.
         /// </summary>
@@ -91,9 +117,25 @@ public int Height
         /// rectangles' sums.</note></para>
         /// </remarks>
         /// 
+        [Obsolete("Please use Matrix property instead.")]
         public uint[,] InternalData
         {
-            get { return integralImage; }
+            get
+            {
+                if (integralImage == null)
+                {
+                    integralImage = new uint[height + 1, width + 1];
+                    for (int y = 0; y <= height; y++)
+                    {
+                        for (int x = 0; x <= width; x++)
+                        {
+                            integralImage[y, x] = matrix[y][x];
+                        }
+                    }
+                }
+
+                return integralImage;
+            }
         }
 
         /// <summary>
@@ -103,15 +145,15 @@ public int Height
         /// <param name="width">Image width.</param>
         /// <param name="height">Image height.</param>
         /// 
-        /// <remarks>The constractor is protected, what makes it imposible to instantiate this
+        /// <remarks>The constructor is protected, what makes it impossible to instantiate this
         /// class directly. To create an instance of this class <see cref="FromBitmap(Bitmap)"/> or
         /// <see cref="FromBitmap(BitmapData)"/> method should be used.</remarks>
         ///
         protected IntegralImage(int width, int height)
         {
             this.width = width;
             this.height = height;
-            integralImage = new uint[height + 1, width + 1];
+            this.matrix = Jagged.Zeros<uint>(height + 1, width + 1);
         }
 
         /// <summary>
@@ -129,7 +171,7 @@ public static IntegralImage FromBitmap(Bitmap image)
             // check image format
             if (image.PixelFormat != PixelFormat.Format8bppIndexed)
             {
-                throw new UnsupportedImageFormatException("Source image can be graysclae (8 bpp indexed) image only.");
+                throw new UnsupportedImageFormatException("Source image can be grayscale (8 bpp indexed) image only.");
             }
 
             // lock source image
@@ -174,7 +216,7 @@ public static IntegralImage FromBitmap(UnmanagedImage image)
             // check image format
             if (image.PixelFormat != PixelFormat.Format8bppIndexed)
             {
-                throw new ArgumentException("Source image can be graysclae (8 bpp indexed) image only.");
+                throw new ArgumentException("Source image can be grayscale (8 bpp indexed) image only.");
             }
 
             // get source image size
@@ -184,7 +226,7 @@ public static IntegralImage FromBitmap(UnmanagedImage image)
 
             // create integral image
             var im = new IntegralImage(width, height);
-            uint[,] integralImage = im.integralImage;
+            uint[][] matrix = im.matrix;
 
             // do the job
             unsafe
@@ -203,7 +245,7 @@ public static IntegralImage FromBitmap(UnmanagedImage image)
 
                         rowSum += *src;
 
-                        integralImage[y, x] = rowSum + integralImage[y - 1, x];
+                        matrix[y][x] = rowSum + matrix[y - 1][x];
                     }
                     src += offset;
                 }
@@ -240,7 +282,7 @@ public uint GetRectangleSum(int x1, int y1, int x2, int y2)
             if (x2 > width) x2 = width;
             if (y2 > height) y2 = height;
 
-            return integralImage[y2, x2] + integralImage[y1, x1] - integralImage[y2, x1] - integralImage[y1, x2];
+            return matrix[y2][x2] + matrix[y1][x1] - matrix[y2][x1] - matrix[y1][x2];
         }
 
         /// <summary>
@@ -256,7 +298,7 @@ public uint GetRectangleSum(int x1, int y1, int x2, int y2)
         /// <remarks><para>The method calculates horizontal wavelet, which is a difference
         /// of two horizontally adjacent boxes' sums, i.e. <b>A-B</b>. A is the sum of rectangle with coordinates
         /// (x, y-radius, x+radius-1, y+radius-1). B is the sum of rectangle with coordinates
-        /// (x-radius, y-radius, x-1, y+radiys-1).</para></remarks>
+        /// (x-radius, y-radius, x-1, y+radius-1).</para></remarks>
         ///
         public int GetHaarXWavelet(int x, int y, int radius)
         {
@@ -313,7 +355,7 @@ public uint GetRectangleSumUnsafe(int x1, int y1, int x2, int y2)
             x2++;
             y2++;
 
-            return integralImage[y2, x2] + integralImage[y1, x1] - integralImage[y2, x1] - integralImage[y1, x2];
+            return matrix[y2][x2] + matrix[y1][x1] - matrix[y2][x1] - matrix[y1][x2];
         }
 
         /// <summary>
@@ -385,7 +427,7 @@ public float GetRectangleMean(int x1, int y1, int x2, int y2)
             if (y2 > height) y2 = height;
 
             // return sum divided by actual rectangles size
-            return (float)((double)(integralImage[y2, x2] + integralImage[y1, x1] - integralImage[y2, x1] - integralImage[y1, x2]) /
+            return (float)((double)(matrix[y2][x2] + matrix[y1][x1] - matrix[y2][x1] - matrix[y1][x2]) /
                 (double)((x2 - x1) * (y2 - y1)));
         }
 
@@ -408,7 +450,7 @@ public float GetRectangleMeanUnsafe(int x1, int y1, int x2, int y2)
             y2++;
 
             // return sum divided by actual rectangles size
-            return (float)((double)(integralImage[y2, x2] + integralImage[y1, x1] - integralImage[y2, x1] - integralImage[y1, x2]) /
+            return (float)((double)(matrix[y2][x2] + matrix[y1][x1] - matrix[y2][x1] - matrix[y1][x2]) /
                 (double)((x2 - x1) * (y2 - y1)));
         }
 

Sources/Accord.Imaging/Interest Points/FREAK/FastRetinaKeypointDescriptor.cs

@@ -345,10 +345,10 @@ private unsafe int mean(double kx, double ky, int scale, int orientation, int po
             int x_right = (int)(xf + radius + 1.5);  // integral image is 1px wider
             int y_bottom = (int)(yf + radius + 1.5); // integral image is 1px higher
 
-            ret_val = (int)Integral.InternalData[y_bottom, x_right]; // bottom right corner
-            ret_val -= (int)Integral.InternalData[y_bottom, x_left];
-            ret_val += (int)Integral.InternalData[y_top, x_left];
-            ret_val -= (int)Integral.InternalData[y_top, x_right];
+            ret_val = (int)Integral.Matrix[y_bottom][x_right]; // bottom right corner
+            ret_val -= (int)Integral.Matrix[y_bottom][x_left];
+            ret_val += (int)Integral.Matrix[y_top][x_left];
+            ret_val -= (int)Integral.Matrix[y_top][x_right];
             ret_val = ret_val / ((x_right - x_left) * (y_bottom - y_top));
             return ret_val;
         }

Sources/Accord.Imaging/Interest Points/SURF/ResponseLayer.cs

@@ -29,6 +29,7 @@ namespace Accord.Imaging
     using System.Collections;
     using System.Collections.Generic;
     using Accord.Imaging;
+    using Accord.Math;
 
     /// <summary>
     ///   Response filter.
@@ -53,7 +54,7 @@ internal class ResponseLayerCollection : IEnumerable<ResponseLayer[]>
         private int step;
         private int octaves;
 
-        private static readonly int[,] map = 
+        private static readonly int[,] map =
         {
             { 0,  1,  2,  3 },
             { 1,  3,  4,  5 },
@@ -268,13 +269,13 @@ internal class ResponseLayer
         ///   Gets the responses computed from the filter.
         /// </summary>
         /// 
-        public float[,] Responses { get; private set; }
+        public float[][] Responses { get; private set; }
 
         /// <summary>
         ///   Gets the Laplacian computed from the filter.
         /// </summary>
         /// 
-        public int[,] Laplacian { get; private set; }
+        public int[][] Laplacian { get; private set; }
 
 
         /// <summary>
@@ -288,8 +289,8 @@ public ResponseLayer(int width, int height, int step, int filter)
             this.Step = step;
             this.Size = filter;
 
-            this.Responses = new float[height, width];
-            this.Laplacian = new int[height, width];
+            this.Responses = Jagged.Zeros<float>(height, width);
+            this.Laplacian = Jagged.Zeros<int>(height, width);
         }
 
         /// <summary>
@@ -301,8 +302,8 @@ public void Update(int width, int height, int step, int filter)
         {
             if (height > Height || width > Width)
             {
-                this.Responses = new float[height, width];
-                this.Laplacian = new int[height, width];
+                this.Responses = Jagged.Zeros<float>(height, width);
+                this.Laplacian = Jagged.Zeros<int>(height, width);
             }
 
             this.Width = width;
@@ -351,8 +352,8 @@ public void Compute(IntegralImage image)
                     Dxy *= inv / 255f;
 
                     // Get the determinant of Hessian response & Laplacian sign
-                    Responses[y, x] = (Dxx * Dyy) - (0.9f * 0.9f * Dxy * Dxy);
-                    Laplacian[y, x] = (Dxx + Dyy) >= 0 ? 1 : 0;
+                    Responses[y][x] = (Dxx * Dyy) - (0.9f * 0.9f * Dxy * Dxy);
+                    Laplacian[y][x] = (Dxx + Dyy) >= 0 ? 1 : 0;
                 }
             }
         }

Sources/Accord.Imaging/Interest Points/SURF/SpeededUpRobustFeaturesDetector.cs

@@ -342,7 +342,7 @@ private List<SpeededUpRobustFeaturePoint> processImage(UnmanagedImage image)
                         int mscale = mid.Width / top.Width;
                         int bscale = bot.Width / top.Width;
 
-                        double currentValue = mid.Responses[y * mscale, x * mscale];
+                        double currentValue = mid.Responses[y * mscale][x * mscale];
 
                         // for each windows' row
                         for (int i = -r; (currentValue >= threshold) && (i <= r); i++)
@@ -354,9 +354,9 @@ private List<SpeededUpRobustFeaturePoint> processImage(UnmanagedImage image)
                                 int xj = x + j;
 
                                 // for each response layer
-                                if (top.Responses[yi, xj] >= currentValue ||
-                                    bot.Responses[yi * bscale, xj * bscale] >= currentValue || ((i != 0 || j != 0) &&
-                                    mid.Responses[yi * mscale, xj * mscale] >= currentValue))
+                                if (top.Responses[yi][xj] >= currentValue ||
+                                    bot.Responses[yi * bscale][xj * bscale] >= currentValue || ((i != 0 || j != 0) &&
+                                    mid.Responses[yi * mscale][xj * mscale] >= currentValue))
                                 {
                                     currentValue = 0;
                                     break;
@@ -378,7 +378,7 @@ private List<SpeededUpRobustFeaturePoint> processImage(UnmanagedImage image)
                                     (x + offset[0]) * tstep,
                                     (y + offset[1]) * tstep,
                                     0.133333333 * (mid.Size + offset[2] * mstep),
-                                    mid.Laplacian[y * mscale, x * mscale]));
+                                    mid.Laplacian[y * mscale][x * mscale]));
                             }
                         }
 
@@ -490,9 +490,9 @@ private static double[] interpolate(int y, int x, ResponseLayer top, ResponseLay
             int xm1 = x - 1, ym1 = y - 1;
 
             // Compute first order scale-space derivatives
-            double dx = (mid.Responses[y * ms, xp1 * ms] - mid.Responses[y * ms, xm1 * ms]) / 2f;
-            double dy = (mid.Responses[yp1 * ms, x * ms] - mid.Responses[ym1 * ms, x * ms]) / 2f;
-            double ds = (top.Responses[y, x] - bot.Responses[y * bs, x * bs]) / 2f;
+            double dx = (mid.Responses[y * ms][xp1 * ms] - mid.Responses[y * ms][xm1 * ms]) / 2f;
+            double dy = (mid.Responses[yp1 * ms][x * ms] - mid.Responses[ym1 * ms][x * ms]) / 2f;
+            double ds = (top.Responses[y][x] - bot.Responses[y * bs][x * bs]) / 2f;
 
             double[] d =
             {
@@ -502,14 +502,14 @@ private static double[] interpolate(int y, int x, ResponseLayer top, ResponseLay
             };
 
             // Compute Hessian
-            double v = mid.Responses[y * ms, x * ms] * 2.0;
-            double dxx = (mid.Responses[y * ms, xp1 * ms] + mid.Responses[y * ms, xm1 * ms] - v);
-            double dyy = (mid.Responses[yp1 * ms, x * ms] + mid.Responses[ym1 * ms, x * ms] - v);
-            double dxs = (top.Responses[y, xp1] - top.Responses[y, x - 1] - bot.Responses[y * bs, xp1 * bs] + bot.Responses[y * bs, xm1 * bs]) / 4f;
-            double dys = (top.Responses[yp1, x] - top.Responses[y - 1, x] - bot.Responses[yp1 * bs, x * bs] + bot.Responses[ym1 * bs, x * bs]) / 4f;
-            double dss = (top.Responses[y, x] + bot.Responses[y * ms, x * ms] - v);
-            double dxy = (mid.Responses[yp1 * ms, xp1 * ms] - mid.Responses[yp1 * ms, xm1 * ms]
-                - mid.Responses[ym1 * ms, xp1 * ms] + mid.Responses[ym1 * ms, xm1 * ms]) / 4f;
+            double v = mid.Responses[y * ms][x * ms] * 2.0;
+            double dxx = (mid.Responses[y * ms][xp1 * ms] + mid.Responses[y * ms][xm1 * ms] - v);
+            double dyy = (mid.Responses[yp1 * ms][x * ms] + mid.Responses[ym1 * ms][x * ms] - v);
+            double dxs = (top.Responses[y][xp1] - top.Responses[y][x - 1] - bot.Responses[y * bs][xp1 * bs] + bot.Responses[y * bs][xm1 * bs]) / 4f;
+            double dys = (top.Responses[yp1][x] - top.Responses[y - 1][x] - bot.Responses[yp1 * bs][x * bs] + bot.Responses[ym1 * bs][x * bs]) / 4f;
+            double dss = (top.Responses[y][x] + bot.Responses[y * ms][x * ms] - v);
+            double dxy = (mid.Responses[yp1 * ms][xp1 * ms] - mid.Responses[yp1 * ms][xm1 * ms]
+                - mid.Responses[ym1 * ms][xp1 * ms] + mid.Responses[ym1 * ms][xm1 * ms]) / 4f;
 
             double[,] H =
             {