工业视觉_71:连通域轮廓的Fourier变换与小波变换

 

/*
    工业视觉_71:连通域轮廓的Fourier变换与小波变换
                
    未来20年里,机器换人,智能使用机器人,改造传统工厂成智慧工厂,将成为最火的行业之一.
    工业视觉,目标很明确:"快,准,稳"三个字. 快:开发快,运行速度快;准:高精度;稳:稳健可靠
    使用高级语言做工程主要优势在:已经有丰富的数据结构和成熟的类型库,如List,Dictionary,Lambda,Accord,...
    所以,目前"机器换人"项目大多采用工控电脑,搭建 Windows7+VS2019+EMGU(或AForge+Accord),这一方案见效最快.
    Halcon,Emgu的视觉库都很强大,而AForge+Accord库更全面,更丰富(如:数学,人工智能,机器学习).
                
        许多产品,为了识别与匹配,划分出连通域后,或者运用AI机器人学习,或者深入开展数据挖掘,以提取连通域或轮廓的一系列特征.
    将边缘上的点变换成相对于Blob,contour中心的半径集合,分析半径集,基本的统计有:min,max,avg,std,...,area, perimeter,...
    原先处理边缘上全部的点,后来只选择0,1,2,3,...,359这些角度上的点也够了.
    更深入的分析,可运用富立叶变换与小波变换,对其频谱特征分析与利用.                
    富立叶变换与小波变换都有它们的逆变换.本文也对各自的逆变换进行了绘制.
    富立叶变换中,忽略R<80的数据,逆变换图像仍然失真很小.
    小波变换中,删除后一半数据,逆变换图像仍然失真很小.
    这些实验,说明了经变换后特征向量的数量可以大大缩小.

    本文系作者在"安吉八塔机器人公司"产品线研发中的拓展与随笔.
    [已经发布,链接:...

        ---------     编撰:    项道德(微信:daode1212),2021-07-21
             
    */



//一,基本图像与绘制工具准备:##########################################################

pictureBox1.Image = Image.FromFile("2D156_.png");//黑色背景的图片,单一工件对象呈白色
//MSG("马上呈现: ...");//自定义函数"MSG()"即"MessageBox.Show()"
Bitmap bmp = (Bitmap)(pictureBox1.Image);
int ww = bmp.Width, hh = bmp.Height;
string txtAll = "";
Graphics g = Graphics.FromImage(bmp);
SolidBrush bh0 = new SolidBrush(Color.FromArgb(0,111,  111));
SolidBrush bh1 = new SolidBrush(Color.FromArgb(161, 0, 161));
SolidBrush bh2 = new SolidBrush(Color.FromArgb(111, 111, 0));
SolidBrush bh3 = new SolidBrush(Color.FromArgb(211, 211, 211));
Pen pen0 = new Pen(Color.FromArgb(182,122,182), 1);
Pen pen1 = new Pen(Color.FromArgb(255, 122, 0), 1);
Pen pen2 = new Pen(Color.FromArgb(0, 110, 255), 1);
Pen pen3 = new Pen(Color.FromArgb(0, 255, 110), 1);


//一,轮廓处理:###################################################################

var Gry = new Accord.Imaging.Filters.Grayscale(0.2, 0.5, 0.1).Apply(bmp);//灰度化
//pictureBox1.Image = Gry;
//MSG("位图处理呈现:灰度化");

var Dila = new Accord.Imaging.Filters.Dilation3x3().Apply(Gry);//[白色点]膨胀
//pictureBox1.Image = Dila;
//MSG("位图处理呈现:膨胀");

var fhs = new Accord.Imaging.Filters.FillHoles().Apply(Dila); //填补黑色的小洞
//pictureBox1.Image = fhs;
//MSG("位图处理呈现:填补小洞");

var blr = new Accord.Imaging.Filters.FastBoxBlur(13,13).Apply(fhs);//平滑化
//pictureBox1.Image = blr;
//MSG("位图处理呈现:平滑化");

Bitmap blbs = new Accord.Imaging.Filters.BlobsFiltering(250, 250, 600, 600).Apply(blr);//过滤出合适大小的连通域
//pictureBox1.Image = blbs;
//MSG("位图处理呈现:BlobsFiltering");

Bitmap tsd = new Accord.Imaging.Filters.Threshold(32).Apply(blbs);
//pictureBox1.Image = tsd;
//MSG("位图处理呈现:二值化");

var Eros = new Accord.Imaging.Filters.Erosion3x3().Apply(tsd);//腐蚀:黑色点扩大
//pictureBox1.Image = Eros;
//MSG("位图处理呈现:腐蚀");

Bitmap ced = new Accord.Imaging.Filters.CannyEdgeDetector().Apply(Eros);
pictureBox1.Image = ced;
//MSG("位图处理呈现:Canny边缘化");

//二, Canny边缘化后的白色点转换:##############################################################
//1,白色点集 --->> List<Point> LP
List<Point> LP = new List<Point>();
BitmapData data = ced.LockBits(
                new System.Drawing.Rectangle(0, 0, ced.Width, ced.Height),
                ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
System.IntPtr Scan0 = data.Scan0;
int stride = data.Stride;
unsafe
{
    byte* p = (byte*)Scan0;
    int offset = stride - 3 * ced.Width;

    //双层循环读取数据:
    for (int y = 0; y < ced.Height; y++)
    {
        for (int x = 0; x < ced.Width; x++)
        {
            if (0.30 * p[2] + 0.59 * p[1] + 0.11 * p[0] > 128)
            {
                LP.Add(new Point(x,y));
            }                        
            p += 3;
        }
        p += offset;
    }
}
ced.UnlockBits(data);

//2,计算中心:
var Org = new Point((int)LP.Average(o=>o.X),(int)LP.Average(o=>o.Y));

//3,生成360度扫描整数角度上的点:
Dictionary<int, float> DP = new Dictionary<int, float>();//小集合:360度扫描,半径
foreach (var p in LP)
{
    float dx = p.X - Org.X;
    float dy = p.Y - Org.Y;
    int A = (int)(180 + Math.Atan2(dy, dx) * toD)%360;
    float R = (float)Math.Sqrt(dx*dx+dy*dy);              
    if(!DP.ContainsKey((int)A)) DP.Add((int)A, R);
}

//4,按角度排序:
var AR = DP.Select(o => o).OrderBy(o => o.Key);//小集合
//绘制半径分布图:
g = Graphics.FromImage(bmp);
            
foreach(var m in AR){
    int a = m.Key;
    g.DrawLine(pen1, 75 + 3 * a, hh/3, 75 + 3 * a, hh/3 - m.Value);
}

//三, 复数域的Fourier变换: ##############################################################
double[] real0 = new double[AR.Count()], imag0 = new double[AR.Count()],wh = new double[AR.Count()];
for (int i = 0; i < AR.Count(); i++)
{
    real0[i] = AR.ElementAt(i).Key;
    imag0[i] = AR.ElementAt(i).Value;
    wh[i] = AR.ElementAt(i).Value;
}
//复数域.快速富立叶变换.正向:
Accord.Math.Transforms.FourierTransform2.FFT(real0, imag0, Accord.Math.FourierTransform.Direction.Forward);

//数据监视:
//txtAll += "\r\n"+ String.Join(",", real0);
//txtAll += "\r\n"+ String.Join(",", imag0);
//ScvWrite2Array(real0, imag0, "button71.scv");//SCV中可用图表工具观察.

//逆变换:
int cnt = AR.Count();
double[] real2 = new double[cnt], imag2 = new double[cnt];
for (int i = 0; i < cnt; i++)
{
    real2[i] = real0[i];
    imag2[i] = imag0[i];
    if (imag0[i]*imag0[i]<6400) { 
        real2[i] = 0; imag2[i] = 0;//忽略R<80的数据
    }
}

//复数域.快速富立叶变换.逆向:
Accord.Math.Transforms.FourierTransform2.FFT(real2, imag2, Accord.Math.FourierTransform.Direction.Backward);

for (int i = 0; i < cnt; i++)
{
    float A = (float)real2[i];
    float R = (float)imag2[i] ;
    double x = 300 - R * Math.Cos(A/toD);
    double y = 400 - R * Math.Sin(A/toD);
    g.DrawLine(pen2, 300, 400, (float)x, (float)y);//仿真图
    g.DrawLine(pen3, 75 + 3 * A, hh, 75 + 3 * A, hh - R);//频谱图
}

//四, 运用小波变换[实数域]处理半径集:  ##############################################################
//Accord.Math.Wavelets.Haar WH = new Accord.Math.Wavelets.Haar(1);
Accord.Math.Wavelets.CDF97 WH = new Accord.Math.Wavelets.CDF97(1);
WH.Forward(wh);//小波.正变换
//Accord.Math.Wavelets.Haar.FWT(wh);

double[] wh2=new double[cnt];
//wh.CopyTo(wh2,0);

for (int i = 0; i < cnt; i++)
{
    float A = (float)i;
    float R = (float)wh[i] / 2;
    g.DrawLine(pen0, 450 + 3 * A, hh - 310, 450 + 3 * A, hh - 310 - R);//频谱图
    //截取前部[低频部分]:
    if (i < cnt / 2)
    {
        wh2[i] = wh[i];
    }
    else
    {
        wh2[i] = 0;
    }
}

WH.Backward(wh2);//小波.逆变换
//Accord.Math.Wavelets.Haar.IWT(wh2);

for (int i = 0; i < cnt; i++)
{
    float A = (float)i;
    float R = (float)wh2[i];
    double x = 1100 + R * Math.Cos(PI+A / toD);
    double y = 500 + R * Math.Sin(PI+A / toD);
    g.DrawLine(pen3, 1100, 500, (float)x, (float)y);//仿真图
}

//五,连通域处理:  ##############################################################
           
//创建非托管位图:
UnmanagedImage mdImage = UnmanagedImage.FromManagedImage(Eros);

//处理二进制位图:          
BlobCounter bCounter = new BlobCounter(); //创建连通计数器
bCounter.ProcessImage(mdImage);
Blob[] blobs = bCounter.GetObjectsInformation();
MSG("blobs.Count()=" + blobs.Count());

//凸多边形:
var grahamScan = new Accord.Math.Geometry.GrahamConvexHull();//凸多边形识别器
Dictionary<int, List<IntPoint>> hulls = new Dictionary<int, List<IntPoint>>();//凸多边形容器

//各连通域处理:
float MaxArea = blobs.Max(o => o.Area);// 0.09f * bmp.Width * bmp.Height;
for (int i = 1; i < blobs.Length; i++)
{

    //int a = 0;
    if (blobs[i].Area < MaxArea)
    {
        continue;
    }
    if (blobs[i].Area >= MaxArea)
    {
        List<IntPoint> edgePts = bCounter.GetBlobsEdgePoints(blobs[i]);

        //BlobsEdgePoints是无序的,绘制后一片线段:
        //Point[] pt2 = PointsListToArray(edgePts);
        //g.DrawPolygon(pen1, pt2);

        //以点画出BlobsEdgePoints:
        //List<IntPoint> hullPts = grahamScan.FindHull(edgePts);
        //hulls.Add(blobs[i].ID, hullPts);//添加到凸多边形集合中                   
        foreach (var m in edgePts)
        {
            g.FillEllipse(bh1, m.X - 4, m.Y - 4, 8, 8);
        }

        //List<IntPoint> --> Point[]:
        //Point[] EP3 = hulls[blobs[i].ID].Select(o => new Point(o.X, o.Y)).ToArray();

        //g.DrawPolygon(pen1, EP3);//绘制空心凸多边形
        //g.FillPolygon(bh2, EP3);//绘制实心凸多边形

        //绘制内接空心四边形或三角形:
        //List<IntPoint> QP = Accord.Math.Geometry.PointsCloud.FindQuadrilateralCorners(edgePts);
        //Point[] QP2 = QP.Select(o => new Point(o.X, o.Y)).ToArray();
        //g.DrawPolygon(pen2, QP2);//绘制空心四边形

        //获取点集的最小外接矩形(PointF[]:四个角点上的坐标):
        //PointF[] minRect = getMinRect(edgePts);
        //PointF[] minRect2 = minRect.Select(o => new PointF(o.X + blobs[i].CenterOfGravity.X, o.Y + blobs[i].CenterOfGravity.Y)).ToArray();
        //g.DrawPolygon(pen3, minRect2);//绘制空心四边形

        string txt = "@" + i + "\r\nedgePts.Count=" + edgePts.Count();// +" \r\nhullPts.Count=" + hullPts.Count();
        txt += "\r\nCenterOfGravity(x,y):" + (int)blobs[i].CenterOfGravity.X + "," + (int)blobs[i].CenterOfGravity.Y;
        txt += "\r\nArea:" + blobs[i].Area;
        txt += "\r\nFullness:" + blobs[i].Fullness;
        g.DrawString(txt, new Font("", 12), bh3, blobs[i].Rectangle.Left, blobs[i].Rectangle.Top);
        txtAll += "\r\n\r\n" + txt;
    }
}
g.Dispose();

//显示结果:
textBox1.Text = txtAll;
pictureBox1.Image = bmp;