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

using o0.Geometry2D.Float;
using Org.BouncyCastle.Crypto.Macs;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Threading.Tasks;
using UnityEngine;

namespace o0.Project
{
    public static partial class Extension
    {
        static public (Matrix, Matrix) zimIdentifyEdge55GradientFull(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 factor0 = new float[] { 1, 1.5f, 2 };
                var factor1 = new float[] { 0.94868329805051379959966806332982f, 1.3416407864998738178455042012388f, 1.7888543819998317571273389349851f };
                var factor2 = new float[] { 0.70710678118654752440084436210485f, 1.0606601717798212866012665431573f, 1.4142135623730950488016887242097f };
                //float root2 = 1.4142135623730950488016887242097f;       // 根号2
                //float root5d4 = 1.1180339887498948482045868343656f;     // 根号5/4
                //float root10d3 = 1.0540925533894597773329645148109f;    // 根号10/9

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

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

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

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

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

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

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

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

                float[] gra = new float[] { MathF.Abs(g0), MathF.Abs(g26), MathF.Abs(g45), MathF.Abs(g63), MathF.Abs(g90), MathF.Abs(g116), MathF.Abs(g135), MathF.Abs(g153) };
                var maxIndex = o0.MaxIndex(gra);
                pixelC[i] = gra[maxIndex];
                switch (maxIndex)
                {
                    case 0:
                        pixelDir[i] = g0 > 0 ? 0 : 180;
                        break;
                    case 1:
                        pixelDir[i] = g26 > 0 ? 26.565f : 206.565f;
                        break;
                    case 2:
                        pixelDir[i] = g45 > 0 ? 45 : 225;
                        break;
                    case 3:
                        pixelDir[i] = g63 > 0 ? 63.435f : 243.435f;
                        break;
                    case 4:
                        pixelDir[i] = g90 > 0 ? 90 : 270;
                        break;
                    case 5:
                        pixelDir[i] = g116 > 0 ? 116.565f : 296.565f;
                        break;
                    case 6:
                        pixelDir[i] = g135 > 0 ? 135 : 315;
                        break;
                    case 7:
                        pixelDir[i] = g153 > 0 ? 153.435f : 333.435f;
                        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 / 25f + 1f;//+1是为了防止除数太小，结果不合理地过大
            });
            return (convolution.Map(ignoreExtremum: 10), dirMat);
        }

        static public (Matrix, Matrix) zimIdentifyEdgeGradient45(this Matrix mat, int conSizeLevel = 2)
        {
            var conSize = 4 * conSizeLevel + 1;        // conSize一定 = 4 * n + 1

            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;

                float root2 = 1.4142135623730950488016887242097f;       // 根号2
                float root5d4 = 1.1180339887498948482045868343656f;     // 根号5/4
                //float root10d3 = 1.0540925533894597773329645148109f;    // 根号10/9

                var g0 =
                  mat.Element[(iPx + 4 * conSizeLevel) + (iPy + 1 * conSizeLevel) * mat.Size.x] * 1f
                + mat.Element[(iPx + 4 * conSizeLevel) + (iPy + 2 * conSizeLevel) * mat.Size.x] * 1.5f
                + mat.Element[(iPx + 4 * conSizeLevel) + (iPy + 3 * conSizeLevel) * mat.Size.x] * 1f
                - mat.Element[(iPx + 0 * conSizeLevel) + (iPy + 3 * conSizeLevel) * mat.Size.x] * 1f
                - mat.Element[(iPx + 0 * conSizeLevel) + (iPy + 2 * conSizeLevel) * mat.Size.x] * 1.5f
                - mat.Element[(iPx + 0 * conSizeLevel) + (iPy + 1 * conSizeLevel) * mat.Size.x] * 1f;

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

                g26 /= root5d4;

                var g45 =
                  mat.Element[(iPx + 4 * conSizeLevel) + (iPy + 3 * conSizeLevel) * mat.Size.x] * 1f
                + mat.Element[(iPx + 4 * conSizeLevel) + (iPy + 4 * conSizeLevel) * mat.Size.x] * 1.5f
                + mat.Element[(iPx + 3 * conSizeLevel) + (iPy + 4 * conSizeLevel) * mat.Size.x] * 1f
                - mat.Element[(iPx + 0 * conSizeLevel) + (iPy + 1 * conSizeLevel) * mat.Size.x] * 1f
                - mat.Element[(iPx + 0 * conSizeLevel) + (iPy + 0 * conSizeLevel) * mat.Size.x] * 1.5f
                - mat.Element[(iPx + 1 * conSizeLevel) + (iPy + 0 * conSizeLevel) * mat.Size.x] * 1f;

                g45 /= root2;

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

                g63 /= root5d4;

                var g90 =
                  mat.Element[(iPx + 3 * conSizeLevel) + (iPy + 4 * conSizeLevel) * mat.Size.x] * 1f
                + mat.Element[(iPx + 2 * conSizeLevel) + (iPy + 4 * conSizeLevel) * mat.Size.x] * 1.5f
                + mat.Element[(iPx + 1 * conSizeLevel) + (iPy + 4 * conSizeLevel) * mat.Size.x] * 1f
                - mat.Element[(iPx + 1 * conSizeLevel) + (iPy + 0 * conSizeLevel) * mat.Size.x] * 1f
                - mat.Element[(iPx + 2 * conSizeLevel) + (iPy + 0 * conSizeLevel) * mat.Size.x] * 1.5f
                - mat.Element[(iPx + 3 * conSizeLevel) + (iPy + 0 * conSizeLevel) * mat.Size.x] * 1f;

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

                g116 /= root5d4;

                var g135 =
                  mat.Element[(iPx + 1 * conSizeLevel) + (iPy + 4 * conSizeLevel) * mat.Size.x] * 1f
                + mat.Element[(iPx + 0 * conSizeLevel) + (iPy + 4 * conSizeLevel) * mat.Size.x] * 1.5f
                + mat.Element[(iPx + 0 * conSizeLevel) + (iPy + 3 * conSizeLevel) * mat.Size.x] * 1f
                - mat.Element[(iPx + 3 * conSizeLevel) + (iPy + 0 * conSizeLevel) * mat.Size.x] * 1f
                - mat.Element[(iPx + 4 * conSizeLevel) + (iPy + 0 * conSizeLevel) * mat.Size.x] * 1.5f
                - mat.Element[(iPx + 4 * conSizeLevel) + (iPy + 1 * conSizeLevel) * mat.Size.x] * 1f;

                g135 /= root2;

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

                g153 /= root5d4;

                float[] gra = new float[] { MathF.Abs(g0), MathF.Abs(g26), MathF.Abs(g45), MathF.Abs(g63), MathF.Abs(g90), MathF.Abs(g116), MathF.Abs(g135), MathF.Abs(g153) };
                var maxIndex = o0.MaxIndex(gra);
                pixelC[i] = gra[maxIndex];
                switch (maxIndex)
                {
                    case 0:
                        pixelDir[i] = g0 > 0 ? 0 : 180;
                        break;
                    case 1:
                        pixelDir[i] = g26 > 0 ? 26.565f : 206.565f;
                        break;
                    case 2:
                        pixelDir[i] = g45 > 0 ? 45 : 225;
                        break;
                    case 3:
                        pixelDir[i] = g63 > 0 ? 63.435f : 243.435f;
                        break;
                    case 4:
                        pixelDir[i] = g90 > 0 ? 90 : 270;
                        break;
                    case 5:
                        pixelDir[i] = g116 > 0 ? 116.565f : 296.565f;
                        break;
                    case 6:
                        pixelDir[i] = g135 > 0 ? 135 : 315;
                        break;
                    case 7:
                        pixelDir[i] = g153 > 0 ? 153.435f : 333.435f;
                        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 / 25f + 1f;//+1是为了防止除数太小，结果不合理地过大
            });
            return (convolution.Map(ignoreExtremum: 10), dirMat);
        }

        // conSize 应 ≥ 5 
        static public (Matrix, Matrix) zimIdentifyEdgeGradientAny(this Matrix mat, int conSize = 7)
        {
            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;

            var sampleCount = conSize + conSize - 6;
            (Vector a, Vector b)[] samples = new (Vector, Vector)[sampleCount];
            int[] angles = new int[sampleCount];
            float[] scales = new float[sampleCount];

            int ax = conSize - 2, ay = (conSize - 1) / 2, bx, by;
            for (int i = 0; i < samples.Length; i++)
            {
                bx = conSize - ax - 1;
                by = conSize - ay - 1;
                samples[i] = (new Vector(ax, ay), new Vector(bx, by));
                var dir = new Vector(ax - bx, ay - by);
                scales[i] = dir.Length;
                if (ax > bx)
                    angles[i] = (int)(Math.Atan(dir.y / dir.x) * 180 / Math.PI);
                else if (ax < bx)
                    angles[i] = 180 + (int)(Math.Atan(dir.y / dir.x) * 180 / Math.PI);
                else
                    angles[i] = 90;

                if (ax == 1)
                    ay--;
                else if (ay == conSize - 2)
                    ax--;
                else
                    ay++;
            }
            //Debug.Log(angles.ToJson() + "\r\n" + samples.ToJson() +"\r\n" + scales.ToJson());

            Parallel.For(0, convolution.Size.x * convolution.Size.y, i => {
                var iPx = i % convolutionWidth;
                var iPy = i / convolutionWidth;

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

                var gra = new float[sampleCount];
                //var gAvg = new Vector();
                for (int j = 0; j < samples.Length; j++)
                {
                    var a = samples[j].a;
                    var b = samples[j].b;
                    float aAvg = 0, bAvg = 0;
                    for (int m = -1; m <= 1; m++)       //  取像素点周围9格平均值
                    {
                        for (int n = -1; n <= 1; n++)
                        {
                            aAvg += mat.Element[(iPx + (int)a.x + m) + (iPy + (int)a.y + n) * mat.Size.x];
                            bAvg += mat.Element[(iPx + (int)b.x + m) + (iPy + (int)b.y + n) * mat.Size.x];
                        }
                    }
                    aAvg /= 9;
                    bAvg /= 9;

                    gra[j] = (aAvg - bAvg) / scales[j];
                    //gAvg += (a - b) * gra[j] / scales[j];
                }
                //gAvg /= samples.Length;

                int maxIndex = 0;
                float gMax = gra[0];
                for (int k = 1; k < gra.Length; k++)
                {
                    var gAbs = Math.Abs(gra[k]);
                    if (gMax < gAbs)
                    {
                        maxIndex = k;
                        gMax = gAbs;
                    }
                }
                pixelC[i] = gMax / (avg + 1f);  // +1是防止结果过大

                pixelDir[i] = gra[maxIndex] > 0 ? angles[maxIndex] : angles[maxIndex] + 180;
            });
            return (convolution.Map(ignoreExtremum: 10), dirMat);
        }

        static public (Matrix, Matrix) zimIdentifyEdgeGradientGroup(this Matrix mat, int 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 convList = new List<(float, Vector)>();
                for (var jX = 0; jX < conSize; ++jX)
                    for (var jY = 0; jY < conSize; ++jY)
                        convList.Add((mat.Element[(iPx + jX) + (iPy + jY) * mat.Size.x], new Vector(jX, jY))); ;

                convList.Sort((a, b) => a.Item1.CompareTo(b.Item1));

                var vMean0 = Vector.Zero;  // 暗区
                var vMean1 = Vector.Zero;  // 明区
                for (int k = 0; k < convList.Count / 2; k++)
                {
                    vMean0 += convList[k].Item2;
                    pixelC[i] -= convList[k].Item1;
                }
                vMean0 /= convList.Count / 2;
                for (int k = convList.Count / 2; k < convList.Count; k++)
                {
                    vMean1 += convList[k].Item2;
                    pixelC[i] += convList[k].Item1;
                }
                vMean1 /= convList.Count - convList.Count / 2;
                pixelC[i] /= sum / convList.Count + 1f;
                
                var dir = vMean1 - vMean0;
                //pixelDir[i] = new Vector(dir.y, -dir.x);
                pixelDir[i] = dir.y > 0 ? (float)Math.Atan((double)(dir.y / dir.x)) : -(float)Math.Atan((double)(dir.y / dir.x));
                pixelDir[i] *= 180 / (float)Math.PI;
            });
            return (convolution.Map(ignoreExtremum: 10), dirMat);
        }

        
    }
}