OpenCVSharp 4.5 PCA 主成分分析

用 OpenCVSharp 4.5 跑一遍 OpenCV 官方教程

原 OpenCV 官方教程链接：Introduction to Principal Component Analysis (PCA)

核心函数：PCA

using System;
using OpenCvSharp;

namespace ConsoleApp1
{
    class tutorial42 : ITutorial
    {
        public void Run()
        {
            Mat src = Cv2.ImRead(@"I:\csharp\images\pca_test1.jpg");
            // Check if image is loaded successfully
            if (src.Empty())
            {
                Console.WriteLine("Problem loading image!!!");
                return;
            }
            Cv2.ImShow("src", src);
            // Convert image to grayscale
            Mat gray = new Mat();
            Cv2.CvtColor(src, gray, ColorConversionCodes.BGR2GRAY);
            // Convert image to binary
            Mat bw = new Mat();
            Cv2.Threshold(gray, bw, 50, 255, ThresholdTypes.Binary | ThresholdTypes.Otsu);
            // Find all the contours in the thresholded image
            Point[][] contours;
            HierarchyIndex[] hierarchy;
            Cv2.FindContours(bw, out contours, out hierarchy, RetrievalModes.List, ContourApproximationModes.ApproxNone);
            for (int i = 0; i < contours.Length; i++)
            {
                // Calculate the area of each contour
                double area = Cv2.ContourArea(contours[i]);
                // Ignore contours that are too small or too large
                if (area < 1e2 || 1e5 < area) continue;
                // Draw each contour only for visualisation purposes
                Cv2.DrawContours(src, contours, i, new Scalar(0, 0, 255), 2);
                // Find the orientation of each shape
                getOrientation(contours[i], src);
            }
            Cv2.ImShow("output", src);
            Cv2.WaitKey();
            return;

        }

        private void drawAxis(Mat img, Point p, Point q, Scalar colour, double scale = 0.2)
        {
            double angle = Math.Atan2((double)p.Y - q.Y, (double)p.X - q.X); // angle in radians

            double hypotenuse = Math.Sqrt((double)(p.Y - q.Y) * (p.Y - q.Y) + (p.X - q.X) * (p.X - q.X));
            // Here we lengthen the arrow by a factor of scale
            q.X = (int)(p.X - scale * hypotenuse * Math.Cos(angle));
            q.Y = (int)(p.Y - scale * hypotenuse * Math.Sin(angle));
            Cv2.Line(img, p, q, colour, 1, LineTypes.AntiAlias);
            // create the arrow hooks
            p.X = (int)(q.X + 9 * Math.Cos(angle + Cv2.PI / 4));
            p.Y = (int)(q.Y + 9 * Math.Sin(angle + Cv2.PI / 4));

            Cv2.Line(img, p, q, colour, 1, LineTypes.AntiAlias);
            p.X = (int)(q.X + 9 * Math.Cos(angle - Cv2.PI / 4));
            p.Y = (int)(q.Y + 9 * Math.Sin(angle - Cv2.PI / 4));
            Cv2.Line(img, p, q, colour, 1, LineTypes.AntiAlias);
        }
        private double getOrientation(Point[] pts, Mat img)
        {
            //Construct a buffer used by the pca analysis
            int sz = (int)pts.Length;
            Mat data_pts = new Mat(sz, 2, MatType.CV_64F);

            for (int i = 0; i < data_pts.Rows; i++)
            {
                data_pts.At<double>(i, 0) = pts[i].X;
                data_pts.At<double>(i, 1) = pts[i].Y;
            }
            //Perform PCA analysis
            PCA pca_analysis = new PCA(data_pts, new Mat(), PCA.Flags.DataAsRow);
            //Store the center of the object
            Point cntr = new Point((int)(pca_analysis.Mean.At<double>(0, 0)),
                              (int)(pca_analysis.Mean.At<double>(0, 1)));

            //Store the eigenvalues and eigenvectors
            Point2d[] eigen_vecs = new Point2d[2];
            double[] eigen_val = new double[2];
            for (int i = 0; i < 2; i++)
            {
                eigen_vecs[i] = new Point2d(pca_analysis.Eigenvectors.At<double>(i, 0),
                                        pca_analysis.Eigenvectors.At<double>(i, 1));
                eigen_val[i] = pca_analysis.Eigenvalues.At<double>(i);
            }
            // Draw the principal components
            Cv2.Circle(img, cntr, 3, new Scalar(255, 0, 255), 2);
            Point px = new Point((int)(0.02 * eigen_vecs[0].X * eigen_val[0]), (int)(0.02 * eigen_vecs[0].Y * eigen_val[0]));
            Point py = new Point((int)(0.02 * eigen_vecs[1].X * eigen_val[1]), (int)(0.02 * eigen_vecs[1].Y * eigen_val[1]));
            Point p1 = cntr + px;
            Point p2 = cntr - py;

            drawAxis(img, cntr, p1, new Scalar(0, 255, 0), 1);
            drawAxis(img, cntr, p2, new Scalar(255, 255, 0), 5);
            double angle = Math.Atan2(eigen_vecs[0].Y, eigen_vecs[0].X); // orientation in radians
            return angle;
        }

    }
}