surf.net 之关键点检测

//还未创建描述子哦

//尽量将计算量减少到最少，但是速度还是比不上CPP ,640X480竟然平均600MS，如果c++应该快1倍

public class SurfFeature : ICloneable
    {
        public SurfFeature(PointD pos, int lap, int size, double dir, double hessian, double[] discriptors)
        {
            this.pt = pos;
            this.laplacian = lap;
            this.size = size;
            this.dir = dir;
            this.hessian = hessian;
            this.discriptors = discriptors; 
        }


        public PointD pt;
        public int laplacian;
        public int size;
        public double dir;
        public double hessian;
        public double[] discriptors;


        public SurfFeature clone
        {
            get { return new SurfFeature(pt == null ? null : pt.clone, laplacian, size, dir, hessian, discriptors == null ? null : (double[])discriptors.Clone()); }
        }


        public object Clone()
        {
            return clone;
        }
    }

[Serializable]

    public struct SurfBox : ICloneable
    {
        public int top_left, top_right, bottom_left, bottom_right;
        public float weight;

        public SurfBox(int top_left, int top_right, int bottom_left, int bottom_right, float weight)
        {
            this.top_left = top_left;
            this.top_right = top_right;
            this.bottom_left = bottom_left;
            this.bottom_right = bottom_right;
            this.weight = weight;
        }

        public SurfBox clone
        {
            get { return new SurfBox(top_left, top_right, bottom_left, bottom_right, weight); }
        }
        public object Clone()
        {
            return clone;
        }
        public float Hear(ImageMapF iimg)
        {
            return Hear(iimg, 0);
        }
        public float Hear(ImageMapF iimg,int pos)
        {
            return (iimg[pos + top_left] + iimg[pos + bottom_right] - iimg[pos + top_right] - iimg[pos + bottom_left]) * weight; 
        }
        public static float Hears(ImageMapF iimg, int pos, SurfBox[] areas)
        {
            float sum = 0f;
            int n = areas.Length;
            for (int i = 0; i < n; i++)
                sum += (iimg[pos + areas[i].top_left] + iimg[pos + areas[i].bottom_right] - iimg[pos + areas[i].top_right] - iimg[pos + areas[i].bottom_left]) * areas[i].weight;            
            return sum;
        }

        public static SurfBox[][] Boxes(int[,] src, int oldSize, int imgWidth, int imgHeight, params int[] newSizes)
        {
            int n = src.GetLength(0);
            int m = newSizes.Length;
            int dx1, dy1, dx2, dy2;


            SurfBox[][] dst = new SurfBox[m][];
            for (int sz = 0; sz < m; sz++)
            {
                dst[sz] = new SurfBox[n];
                int newSize = newSizes[sz];
                float ratio = (float)newSize / (float)oldSize;


                for (int k = 0; k < n; k++)
                {
                    dx1 = Math.Max(0, (int)Math.Round(ratio * src[k, 0]));
                    dy1 = Math.Max(0, (int)Math.Round(ratio * src[k, 1]));
                    dx2 = Math.Min(imgWidth - 1, (int)Math.Round(ratio * src[k, 2]));
                    dy2 = Math.Min(imgHeight - 1, (int)Math.Round(ratio * src[k, 3]));


                    dst[sz][k] = new SurfBox();
                    dst[sz][k].top_left = dy1 * imgWidth + dx1;
                    dst[sz][k].top_right = dy1 * imgWidth + dx2;
                    dst[sz][k].bottom_left = dy2 * imgWidth + dx1;
                    dst[sz][k].bottom_right = dy2 * imgWidth + dx2;
                    dst[sz][k].weight = (float)src[k, 4] / ((float)(dx2 - dx1) * (dy2 - dy1));
                }
            }
            return dst;
        }

    }

//integral

public static ImageMapF Integral(Bitmap image, float alpha)
        {
            ImageMapF pic = new ImageMapF(image.Width, image.Height);
            float rowsum = 0;
            BitmapData dataIn = image.LockBits(new Rectangle(0, 0, image.Width, image.Height), ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);
            unsafe
            {
                byte* pIn = (byte*)(dataIn.Scan0.ToPointer());
                for (int y = 0; y < dataIn.Height; y++)
                {
                    rowsum = 0;
                    for (int x = 0; x < dataIn.Width; x++)
                    {
                        int cb = (byte)(pIn[0]);
                        int cg = (byte)(pIn[1]);
                        int cr = (byte)(pIn[2]);


                        rowsum += (float)(cr + cb + cg) / alpha;
                        if (y == 0)
                            pic[0, x] = rowsum;
                        else
                            pic[y, x] = rowsum + pic[y - 1, x];


                        pIn += 3;
                    }
                    pIn += dataIn.Stride - dataIn.Width * 3;
                }
            }
            image.UnlockBits(dataIn);
            return pic;
        }

    public class ImageMapF
    {
        public ImageMapF() { }
        public ImageMapF(int width, int height)
        {
            this.width = width;
            this.height = height;
            datas = new float[width * height];
        }
        public ImageMapF(int width, int height, float value)
        {
            this.width = width;
            this.height = height;
            int area = width * height;
            datas = new float[area];
            for (int i = 0; i < area; i++)
                datas[i] = value;
        }
        public ImageMapF(int width, int height, float[] datas)
        {
            this.datas = datas;
            this.width = width;
            this.height = height;
        }


        int width, height;
        float[] datas;
        
        public int Width
        {
            get { return width; }
            set { width = value; }
        }
        public int Height
        {
            get { return height; }
            set { height = value; }
        }
        public float[] Datas
        {
            get { return datas; }
            set { datas = value; }
        }
        
        public float this[int index]
        {
            get { return datas[index]; }
            set { datas[index] = value; }
        }
        public float this[int y, int x]
        {
            get { return datas[x + y * width]; }
            set { datas[x + y * width] = value; }
        }


        public float BoxIntegral(int x, int y, int rows, int cols)
        {
            int r1 = Math.Min(y, height) - 1;
            int c1 = Math.Min(x, width) - 1;
            int r2 = Math.Min(y + rows, height) - 1;
            int c2 = Math.Min(x + cols, width) - 1;


            float A = 0, B = 0, C = 0, D = 0;
            if (r1 >= 0 && c1 >= 0) A = datas[r1 * width + c1];
            if (r1 >= 0 && c2 >= 0) B = datas[r1 * width + c2];
            if (r2 >= 0 && c1 >= 0) C = datas[r2 * width + c1];
            if (r2 >= 0 && c2 >= 0) D = datas[r2 * width + c2];


            return Math.Max(0, A - B - C + D);
        }
    }

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ImageProcessing.Base.Imaging;
using ImageProcessing.Graph;
using ImageProcessing.Generic;


namespace Surf.Base.Net.FastSurf
{
    public class SurfOctave
    {
        public static bool InterpolateKeypoint(float[,] mat, int dx, int dy, int dscale, SurfFeature feature)
        {
            bool solve_ok = true;
            Matrix A = new Matrix(3, 3);
            Matrix X = new Matrix(3, 1);


            X[0, 0] = -(mat[1, 5] - mat[1, 3]) / 2;  /* Negative 1st deriv with respect to x */
            X[1, 0] = -(mat[1, 7] - mat[1, 1]) / 2;  /* Negative 1st deriv with respect to y */
            X[2, 0] = -(mat[2, 4] - mat[0, 4]) / 2;  /* Negative 1st deriv with respect to s */


            A[0, 0] = mat[1, 3] - 2 * mat[1, 4] + mat[1, 5];            /* 2nd deriv x, x */
            A[0, 1] = (mat[1, 8] - mat[1, 6] - mat[1, 2] + mat[1, 0]) / 4; /* 2nd deriv x, y */
            A[0, 2] = (mat[2, 5] - mat[2, 3] - mat[0, 5] + mat[0, 3]) / 4; /* 2nd deriv x, s */
            A[1, 0] = A[0, 1];                                    /* 2nd deriv y, x */
            A[1, 1] = mat[1, 1] - 2 * mat[1, 4] + mat[1, 7];            /* 2nd deriv y, y */
            A[1, 2] = (mat[2, 7] - mat[2, 1] - mat[0, 7] + mat[0, 1]) / 4; /* 2nd deriv y, s */
            A[2, 0] = A[0, 2];                                    /* 2nd deriv s, x */
            A[2, 1] = A[1, 2];                                    /* 2nd deriv s, y */
            A[2, 2] = mat[0, 4] - 2 * mat[1, 4] + mat[2, 4];            /* 2nd deriv s, s */


            A.SolveLinear(X);


            if (X[0, 0] > 1 || X[0, 0] < -1 || X[1, 0] > 1 || X[1, 0] < -1 || X[2, 0] > 1 || X[2, 0] < -1)
                solve_ok = false;
            else
            {
                feature.pt.X += X[0, 0] * dx;
                feature.pt.Y += X[1, 0] * dy;
                feature.size = (int)Math.Round(feature.size + X[2, 0] * dscale);
            }
            return solve_ok;
        }


        public const int haar_sizeo = 9;
        public const int haar_size_inc = 6;
        public static int[,] dx_s = { { 0, 2, 3, 7, 1 }, { 3, 2, 6, 7, -2 }, { 6, 2, 9, 7, 1 } };
        public static int[,] dy_s = { { 2, 0, 7, 3, 1 }, { 2, 3, 7, 6, -2 }, { 2, 6, 7, 9, 1 } };
        public static int[,] dxy_s = { { 1, 1, 4, 4, 1 }, { 5, 1, 8, 4, -1 }, { 1, 5, 4, 8, -1 }, { 5, 5, 8, 8, 1 } };
        public static List<SurfFeature> GetKeypoints(ImageMapF iimg, int octaves, int layers, float thresh)
        {
            List<SurfFeature> fts = new List<SurfFeature>();
            int margin, size, size0, size1;
            SurfBox[][] xBoxes, yBoxes, xyBoxes;
            SurfFeature ft;
            int step = 1;
            ImageMapF[][] spaces = new ImageMapF[octaves][];
            for (int i = 0; i < octaves; i++)
            {
                spaces[i] = new ImageMapF[layers + 2];
                for (int j = 0; j < layers + 2; j++)
                    spaces[i][j] = new ImageMapF(iimg.Width / step, iimg.Height / step, -1f);
                step *= 2;
            }
            step = 1;
            for (int octave = 0; octave < octaves; octave++)
            {
                for (int layer = 0; layer < layers + 2; layer++)
                {
                    if (layer != 0 && layer != layers + 1)
                    {
                        size = (haar_sizeo + haar_size_inc * layer) << octave;
                        size0 = (haar_sizeo + haar_size_inc * (layer - 1)) << octave;
                        size1 = (haar_sizeo + haar_size_inc * (layer + 1)) << octave;


                        margin = (((step + size1 / 2) + 1) / 2) * 2; //取偶数


                        xBoxes = SurfBox.Boxes(dx_s, haar_sizeo, iimg.Width, iimg.Height, size0, size, size1);
                        yBoxes = SurfBox.Boxes(dy_s, haar_sizeo, iimg.Width, iimg.Height, size0, size, size1);
                        xyBoxes = SurfBox.Boxes(dxy_s, haar_sizeo, iimg.Width, iimg.Height, size0, size, size1);


                        for (int y = margin; y < iimg.Height - 1 - margin; y += step)
                            for (int x = margin; x < iimg.Width - 1 - margin; x += step)
                                if (TryGetFeature(iimg, spaces, octave, layer, step, x, y, thresh, xBoxes, yBoxes, xyBoxes, out ft, size0, size, size1))
                                    fts.Add(ft);
                    }
                }
                step *= 2;
            }
            return fts;
        }


        /// <summary>
        /// 判断(x,y)是否是脚点
        /// </summary>
        static bool TryGetFeature(ImageMapF iimg,ImageMapF[][] spaces, int octave, int layer, int step, int x, int y,
            float thresh, SurfBox[][] xBoxes, SurfBox[][] yBoxes, SurfBox[][] xyBoxes, out SurfFeature result, params int[] sizes)
        {
            result = null;
            float dx, dy, dxy;
            float det, db;
            byte trc;


            int x0 = x / step;
            int y0 = y / step;
            int idx = (x - sizes[1] / 2 + (y - sizes[1] / 2) * iimg.Width);


            dx = SurfBox.Hears(iimg, idx, xBoxes[1]); //获取中心窗口的权值
            dy = SurfBox.Hears(iimg, idx, yBoxes[1]);
            dxy = SurfBox.Hears(iimg, idx, xyBoxes[1]);
            det = dx * dy - 0.81f * dxy * dxy;


            if (det >= thresh) //如果大于阈值,判断和邻域两个layer之间的大小关系
            {
                trc = (byte)((dx + dy) >= 0 ? 1: 0);
                float[,] mat = new float[3, 9];


                bool ex = true;
                for (int i = -1; i <= 1 && ex; i++) //layer
                {
                    int size = sizes[i + 1]; //获取窗口尺寸
                    int rd = size / 2;       //窗口半径


                    int nx, ny;
                    for (int j = -1; j <= 1 && ex; j++) //Y
                    {
                        ny = y + j * step - rd;
                        for (int k = -1; k <= 1 && ex; k++) //X
                        {
                            if (spaces[octave][layer + i][y0 + j, x0 + k] == -1)
                            {
                                nx = x + k * step - rd;
                                idx = nx + ny * iimg.Width;


                                dx = SurfBox.Hears(iimg, idx, xBoxes[i + 1]);
                                dy = SurfBox.Hears(iimg, idx, yBoxes[i + 1]);
                                dxy = SurfBox.Hears(iimg, idx, xyBoxes[i + 1]);


                                db = dx * dy - 0.81f * dxy * dxy;
                                spaces[octave][layer + i][y0 + j, x0 + k] = db;
                            }
                            else                            
                                db = spaces[octave][layer + i][y0 + j, x0 + k];
                            if (i != 0 && j != 0 && k != 0)
                                if (db >= det) ex = false;


                            mat[i + 1, 4 + j * 3 + k] = db;
                        }
                    }
                }


                if (ex)
                {
                    SurfFeature ft = new SurfFeature(new PointD(x, y), trc, sizes[1], 0, det, null);
                    if (InterpolateKeypoint(mat, step, step, sizes[1] - sizes[0], ft))
                    {
                        result = ft;
                        return true;
                    }
                }
            }
            return false;
        }
    }
}