smart-bow-infrared/Assets/InfraredProject/o0/o0IdentifyEdge.cs

using System;
using System.Collections.Generic;
using System.Threading.Tasks;

namespace o0.Project
{
    public static partial class Extension
    {
        static private Matrix IdentifyEdgeVariance(this Matrix mat)
        {
            var conSize = 3;
            var conArea = conSize * conSize;
            var convolution = new Matrix(new Geometry2D.Vector<int>(mat.Size.x - (conSize - 1), mat.Size.y - (conSize - 1)), Tiling: true);
            var pixelC = convolution.Element;
            var convolutionWidth = convolution.Size.x;
            Parallel.For(0, convolution.Size.x * convolution.Size.y, i => {
                var iPx = i % convolutionWidth;
                var iPy = i / convolutionWidth;
                var mean = 0.0f;
                //var max = 0.0f;
                //var min = 1.0f;
                List<float> rList = new List<float>();
                for (var jX = 0; jX < conSize; ++jX)
                {
                    for (var jY = 0; jY < conSize; ++jY)
                    {
                        var jIndex = (iPx + jX) + (iPy + jY) * mat.Size.x;
                        var r = mat.Element[jIndex];
                        mean += r;
                        rList.Add(r);
                        /*
                        if (max < r)
                            max = r;
                        if (min > r)
                            min = r;/**/
                    }
                }
                mean /= conArea;
                var variance = 0.0f;
                foreach (var j in rList)
                {
                    var d = (j - mean);// / mean;
                    variance += d * d;
                }
                variance /= conArea;
                variance /= mean;
                /*
                foreach (var j in (Geometry2D.Vector<int>.One * 3).Range())
                {
                    var jP = iP + j;
                    var jIndex = jP.x + jP.y * Size.x;
                    rList.Add(pixelLighted[jIndex].r);
                }/**/
                //pixelC[i] = new UnityEngine.Color(Num.Extension.Variance(rList), 0, 0);
                pixelC[i] = variance;
            });
            return convolution.Map(ignoreExtremum: 10);
        }
        static private Matrix IdentifyEdgeGradientX(this Matrix mat)
        {
            var conSize = 3;
            var conArea = conSize * conSize;
            var convolution = new Matrix(new Geometry2D.Vector<int>(mat.Size.x - (conSize - 1), mat.Size.y - (conSize - 1)), Tiling: true);
            var pixelC = convolution.Element;
            var convolutionWidth = convolution.Size.x;
            Parallel.For(0, convolution.Size.x * convolution.Size.y, i => {
                var iPx = i % convolutionWidth;
                var iPy = i / convolutionWidth;

                pixelC[i] = MathF.Abs(mat.Element[(iPx + 0) + (iPy + 1) * mat.Size.x] - mat.Element[(iPx + 2) + (iPy + 1) * mat.Size.x]);
            });
            return convolution.Map(ignoreExtremum: 10);
        }
        static private Matrix IdentifyEdgeGradientY(this Matrix mat)
        {
            var conSize = 3;
            var conArea = conSize * conSize;
            var convolution = new Matrix(new Geometry2D.Vector<int>(mat.Size.x - (conSize - 1), mat.Size.y - (conSize - 1)), Tiling: true);
            var pixelC = convolution.Element;
            var convolutionWidth = convolution.Size.x;
            Parallel.For(0, convolution.Size.x * convolution.Size.y, i => {
                var iPx = i % convolutionWidth;
                var iPy = i / convolutionWidth;

                pixelC[i] = MathF.Abs(mat.Element[(iPx + 1) + (iPy + 0) * mat.Size.x] - mat.Element[(iPx + 1) + (iPy + 2) * mat.Size.x]);
            });
            return convolution.Map(ignoreExtremum: 10);
        }
        static private Matrix IdentifyEdgeGradient(this Matrix mat)
        {
            var conSize = 3;
            var conArea = conSize * conSize;
            var convolution = new Matrix(new Geometry2D.Vector<int>(mat.Size.x - (conSize - 1), mat.Size.y - (conSize - 1)), Tiling: true);
            var pixelC = convolution.Element;
            var convolutionWidth = convolution.Size.x;
            Parallel.For(0, convolution.Size.x * convolution.Size.y, i => {
                var iPx = i % convolutionWidth;
                var iPy = i / convolutionWidth;

                var gX = mat.Element[(iPx + 2) + (iPy + 1) * mat.Size.x] * 2
                - mat.Element[(iPx + 0) + (iPy + 1) * mat.Size.x] * 2;

                var gY = mat.Element[(iPx + 1) + (iPy + 2) * mat.Size.x] * 2
                - mat.Element[(iPx + 1) + (iPy + 0) * mat.Size.x] * 2;

                float root2 = 1.4142135623730950488016887242097f;

                var gXY = //x0y0向x2y2方向的斜率
                mat.Element[(iPx + 2) + (iPy + 2) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 0) * mat.Size.x];
                gXY *= root2;

                var gYX = //x2y0向x0y2方向的斜率
                + mat.Element[(iPx + 0) + (iPy + 2) * mat.Size.x]
                - mat.Element[(iPx + 2) + (iPy + 0) * mat.Size.x];
                gYX *= root2;

                pixelC[i] = o0.Max(MathF.Abs(gX), MathF.Abs(gY), MathF.Abs(gXY), MathF.Abs(gYX));
            });
            return convolution.Map(ignoreExtremum: 10);
        }

        ////////////////梯度大小图///梯度方向图 0-360度 数学4象限
        static public (Matrix, Matrix) IdentifyEdge(this Matrix mat)//这个效果最好
        {
            var conSize = 3;
            var conArea = conSize * conSize;
            var convolution = new Matrix(new Geometry2D.Vector<int>(mat.Size.x - (conSize - 1), mat.Size.y - (conSize - 1)), Tiling: true);
            var pixelC = convolution.Element;
            var dirMat = new Matrix(convolution.Size, Tiling: true);
            var pixelDir = dirMat.Element;
            var convolutionWidth = convolution.Size.x;
            Parallel.For(0, convolution.Size.x * convolution.Size.y, i => {
                var iPx = i % convolutionWidth;
                var iPy = i / convolutionWidth;

                var sum = 0.0f;
                //var max = 0.0f;
                //var min = 1.0f;
                for (var jX = 0; jX < conSize; ++jX)
                    for (var jY = 0; jY < conSize; ++jY)
                    {
                        var v = mat.Element[(iPx + jX) + (iPy + jY) * mat.Size.x];
                        sum += v;
                        /*
                        if (max < v)
                            max = v;
                        if (min > v)
                            min = v;/**/
                    }
                if (sum == 0.0f)
                    return;
                /*
                if (min == 0.0f)
                    return;/**/

                var gX = mat.Element[(iPx + 2) + (iPy + 0) * mat.Size.x]
                + mat.Element[(iPx + 2) + (iPy + 1) * mat.Size.x] * 2
                + mat.Element[(iPx + 2) + (iPy + 2) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 1) * mat.Size.x] * 2
                - mat.Element[(iPx + 0) + (iPy + 2) * mat.Size.x];

                var gY = mat.Element[(iPx + 0) + (iPy + 2) * mat.Size.x]
                + mat.Element[(iPx + 1) + (iPy + 2) * mat.Size.x] * 2
                + mat.Element[(iPx + 2) + (iPy + 2) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 1) + (iPy + 0) * mat.Size.x] * 2
                - mat.Element[(iPx + 2) + (iPy + 0) * mat.Size.x];

                float root2 = 1.4142135623730950488016887242097f;

                var gXY = mat.Element[(iPx + 1) + (iPy + 2) * mat.Size.x] //x0y0向x2y2方向的斜率
                + mat.Element[(iPx + 2) + (iPy + 2) * mat.Size.x]
                + mat.Element[(iPx + 2) + (iPy + 1) * mat.Size.x]
                - mat.Element[(iPx + 1) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 1) * mat.Size.x];
                gXY *= root2;

                var gYX = mat.Element[(iPx + 1) + (iPy + 2) * mat.Size.x] //x2y0向x0y2方向的斜率
                + mat.Element[(iPx + 0) + (iPy + 2) * mat.Size.x]
                + mat.Element[(iPx + 0) + (iPy + 1) * mat.Size.x]
                - mat.Element[(iPx + 1) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 2) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 2) + (iPy + 1) * mat.Size.x];
                gYX *= root2;

                float[] gra = new float[] { MathF.Abs(gX), MathF.Abs(gY), MathF.Abs(gXY), MathF.Abs(gYX) };
                var maxIndex = o0.MaxIndex(gra);
                pixelC[i] = gra[maxIndex];
                switch (maxIndex)
                {
                    case 0:
                        pixelDir[i] = gX > 0 ? 0 : 180;
                        break;
                    case 1:
                        pixelDir[i] = gY > 0 ? 90 : 270;
                        break;
                    case 2:
                        pixelDir[i] = gXY > 0 ? 45 : 225;
                        break;
                    case 3:
                        pixelDir[i] = gYX > 0 ? 135 : 315;
                        break;
                }
                //pixelC[i] = MathF.Abs(gX) + MathF.Abs(gY) + MathF.Abs(gXY) + MathF.Abs(gYX);
                /*
                if (pixelC[i] / 9f < 1.0f / 256f)
                {
                    pixelC[i] = 0;
                    return;
                }/**/
                /*
                if (sum < 100f/256f)
                {
                    pixelC[i] = 0;
                    return;
                }/**/
                //pixelC[i] = o0.Max(MathF.Abs(gX), MathF.Abs(gY), MathF.Abs(gXY), MathF.Abs(gYX));
                //pixelC[i] /= MathF.Pow(sum / 9, 0.1f);
                //pixelC[i] /= min + 0.5f;
                //pixelC[i] /= min + 0.1f;
                pixelC[i] /= sum / 9f + 1f;//+1是为了防止除数太小，结果不合理地过大
            });
            return (convolution.Map(ignoreExtremum: 10), dirMat);
        }
        static public (Matrix, Matrix) IdentifyEdge55GradientFull(this Matrix mat)
        {
            var conSize = 5;
            var conArea = conSize * conSize;
            var convolution = new Matrix(new Geometry2D.Vector<int>(mat.Size.x - (conSize - 1), mat.Size.y - (conSize - 1)), Tiling: true);
            var pixelC = convolution.Element;
            var dirMat = new Matrix(convolution.Size, Tiling: true);
            var pixelDir = dirMat.Element;
            var convolutionWidth = convolution.Size.x;
            Parallel.For(0, convolution.Size.x * convolution.Size.y, i => {
                var iPx = i % convolutionWidth;
                var iPy = i / convolutionWidth;

                var sum = 0.0f;
                for (var jX = 0; jX < conSize; ++jX)
                    for (var jY = 0; jY < conSize; ++jY)
                        sum += mat.Element[(iPx + jX) + (iPy + jY) * mat.Size.x];
                if (sum == 0.0f)
                    return;

                var gX = mat.Element[(iPx + 4) + (iPy + 1) * mat.Size.x]     //梯度0度
                + mat.Element[(iPx + 4) + (iPy + 2) * mat.Size.x] * 2
                + mat.Element[(iPx + 4) + (iPy + 3) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 1) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 2) * mat.Size.x] * 2
                - mat.Element[(iPx + 0) + (iPy + 3) * mat.Size.x];

                var gY = mat.Element[(iPx + 1) + (iPy + 4) * mat.Size.x]     //梯度90度
                + mat.Element[(iPx + 2) + (iPy + 4) * mat.Size.x] * 2
                + mat.Element[(iPx + 3) + (iPy + 4) * mat.Size.x]
                - mat.Element[(iPx + 1) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 2) + (iPy + 0) * mat.Size.x] * 2
                - mat.Element[(iPx + 3) + (iPy + 0) * mat.Size.x];

                float root2 = 1.4142135623730950488016887242097f;

                var gXY = mat.Element[(iPx + 3) + (iPy + 4) * mat.Size.x] //x0y0向x2y2方向的斜率, 45度
                + mat.Element[(iPx + 4) + (iPy + 4) * mat.Size.x]
                + mat.Element[(iPx + 4) + (iPy + 3) * mat.Size.x]
                - mat.Element[(iPx + 1) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 0) + (iPy + 1) * mat.Size.x];
                gXY *= root2;

                var gYX = mat.Element[(iPx + 1) + (iPy + 4) * mat.Size.x] //x2y0向x0y2方向的斜率, 135度
                + mat.Element[(iPx + 0) + (iPy + 4) * mat.Size.x]
                + mat.Element[(iPx + 0) + (iPy + 3) * mat.Size.x]
                - mat.Element[(iPx + 3) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 4) + (iPy + 0) * mat.Size.x]
                - mat.Element[(iPx + 4) + (iPy + 1) * mat.Size.x];
                gYX *= root2;

                float[] gra = new float[] { MathF.Abs(gX), MathF.Abs(gY), MathF.Abs(gXY), MathF.Abs(gYX) };
                var maxIndex = o0.MaxIndex(gra);
                pixelC[i] = gra[maxIndex];
                switch (maxIndex)
                {
                    case 0:
                        pixelDir[i] = gX > 0 ? 0 : 180;
                        break;
                    case 1:
                        pixelDir[i] = gY > 0 ? 90 : 270;
                        break;
                    case 2:
                        pixelDir[i] = gXY > 0 ? 45 : 225;
                        break;
                    case 3:
                        pixelDir[i] = gYX > 0 ? 135 : 315;
                        break;
                }
                //pixelC[i] = MathF.Abs(gX) + MathF.Abs(gY) + MathF.Abs(gXY) + MathF.Abs(gYX);
                /*
                if (pixelC[i] / 9f < 1.0f / 256f)
                {
                    pixelC[i] = 0;
                    return;
                }/**/
                /*
                if (sum < 100f/256f)
                {
                    pixelC[i] = 0;
                    return;
                }/**/
                //pixelC[i] = o0.Max(MathF.Abs(gX), MathF.Abs(gY), MathF.Abs(gXY), MathF.Abs(gYX));
                //pixelC[i] /= MathF.Pow(sum / 9, 0.1f);
                //pixelC[i] /= min + 0.5f;
                //pixelC[i] /= min + 0.1f;
                pixelC[i] /= sum / 9f + 1f;//+1是为了防止除数太小，结果不合理地过大
            });
            return (convolution.Map(ignoreExtremum: 10), dirMat);
        }

    }
}