Advertisement
Guest User

grain counter

a guest
Nov 10th, 2014
1,462
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C# 27.59 KB | None | 0 0
  1. using System;
  2. using System.Collections.Concurrent;
  3. using System.Collections.Generic;
  4. using System.Diagnostics;
  5. using System.Linq;
  6. using System.Text;
  7. using System.Threading.Tasks;
  8. using OpenCvSharp;
  9.  
  10. namespace GrainTest2
  11. {
  12.     class Program
  13.     {
  14.         static void Main(string[] args)
  15.         {
  16.             System.Threading.Thread.CurrentThread.CurrentCulture = System.Globalization.CultureInfo.InvariantCulture;
  17.  
  18.             // hardcoded filenames and grain counts
  19.             string[] files = new[]{"sourceA.jpg", "sourceB.jpg", "sourceC.jpg","sourceD.jpg", "sourceE.jpg", "sourceF.jpg","sourceG.jpg", "sourceH.jpg", "sourceI.jpg", "sourceJ.jpg",};
  20.             int[] expectedGrains = new[]{3, 5, 12, 25, 50, 83, 120, 150, 151, 200,};
  21.  
  22.             // if any command line arguments are given, use those instead of the hardcoded filenames
  23.             // and ignore the grain counts
  24.             if (args.Any())
  25.             {
  26.                 files = args;
  27.                 expectedGrains = Enumerable.Repeat(0, 100).ToArray();
  28.             }
  29.  
  30.             for(int fileNo = 0;fileNo<files.Length;fileNo+=1)
  31.             {
  32.                 Stopwatch watch = Stopwatch.StartNew();
  33.                 string file = files[fileNo];
  34.  
  35.                 using (CvMat source = new CvMat(file))
  36.                 using (CvMat hsv = new CvMat(source.Rows, source.Cols, MatrixType.U8C3))
  37.                 using (CvMat hue = new CvMat(source.Rows, source.Cols, MatrixType.U8C1))
  38.                 using (CvMat saturation = new CvMat(source.Rows, source.Cols, MatrixType.U8C1))
  39.                 using (CvMat value = new CvMat(source.Rows, source.Cols, MatrixType.U8C1))
  40.                 using (CvMat foreground = new CvMat(source.Rows, source.Cols, MatrixType.U8C1))
  41.                 using (CvMat shadows = new CvMat(source.Rows, source.Cols, MatrixType.U8C1))
  42.                 using (CvMat erodedForeground = new CvMat(source.Rows, source.Cols, MatrixType.U8C1))
  43.                 using (CvMat display = new CvMat(source.Rows, source.Cols, MatrixType.U8C3))
  44.                 {
  45.                    
  46.                     // convert image to hsv
  47.                     source.CvtColor(hsv, ColorConversion.RgbToHsv);
  48.                     hsv.Split(hue, saturation, value, null);
  49.                     saturation.SaveImage("01-saturation.png");
  50.                     value.SaveImage("02-value.png");
  51.  
  52.  
  53.                     // separate foreground from background using the saturation channel and Otsu
  54.                     saturation.Threshold(foreground, 0, 255, ThresholdType.Otsu);
  55.                     foreground.SaveImage("03-otsu.png");
  56.  
  57.  
  58.                     // remove very small holes in the foreground
  59.                     using (CvMat negative = new CvMat(foreground.Rows, foreground.Cols, MatrixType.U8C1))
  60.                     using (CvMat filledHoles = new CvMat(foreground.Rows, foreground.Cols, MatrixType.U8C1))
  61.                     {
  62.                         filledHoles.SetZero();
  63.                         foreground.Not(negative);
  64.  
  65.                         using (CvMemStorage storage = new CvMemStorage())
  66.                         using (CvContourScanner scanner = new CvContourScanner(negative, storage, CvContour.SizeOf, ContourRetrieval.List,ContourChain.ApproxNone))
  67.                         {
  68.                             foreach (CvSeq<CvPoint> c in scanner.Where(c=>Math.Abs(c.ContourArea())<20))
  69.                             {
  70.                                 foreground.DrawContours(c, CvColor.White, CvColor.White, 0, -1);
  71.                                 filledHoles.DrawContours(c, CvColor.White, CvColor.White, 0, -1);
  72.                             }
  73.                         }
  74.                         filledHoles.SaveImage("04-filled-holes.png");
  75.                         foreground.SaveImage("04-background-without-holes.png");
  76.                     }
  77.  
  78.  
  79.                     // reomve the grain shadows from the foreground, using a fixed threshold for the value channel
  80.                     value.Threshold(shadows, 140, 255, ThresholdType.Binary);
  81.                     shadows.And(foreground, foreground);
  82.                     shadows.SaveImage("05-grain-shadows.png");
  83.  
  84.                     foreground.SaveImage("06-background-without-shadows.png");
  85.  
  86.                     display.SetZero();
  87.                     display.Merge(saturation,saturation,null,null);
  88.                     display.Scale(display,0.5);
  89.  
  90.                     // list of recognized grains
  91.                     List<Grain> grains = new List<Grain>();
  92.  
  93.                     Random rand = new Random(4); // determined by fair dice throw, guaranteed to be random
  94.  
  95.                     // repeat until we have found all grains (to a maximum of 10000)
  96.                     for (int numIterations = 0; numIterations < 10000; numIterations++ )
  97.                     {
  98.                         // erode the image of the remaining foreground pixels, only big blobs can be grains
  99.                         foreground.Erode(erodedForeground,null,7);
  100.  
  101.                         // pick a number of starting points to fit grains
  102.                         List<CvPoint> startPoints = new List<CvPoint>();
  103.                         using (CvMemStorage storage = new CvMemStorage())
  104.                         using (CvContourScanner scanner = new CvContourScanner(erodedForeground, storage, CvContour.SizeOf, ContourRetrieval.List, ContourChain.ApproxNone))
  105.                         {
  106.                             if (!scanner.Any()) break; // no grains left, finished!
  107.  
  108.                             // search for grains within the biggest blob first (this is arbitrary)
  109.                             var biggestBlob = scanner.OrderByDescending(c => c.Count()).First();
  110.  
  111.                             // pick 10 random edge pixels
  112.                             for (int i = 0; i < 10; i++)
  113.                             {
  114.                                 startPoints.Add(biggestBlob.ElementAt(rand.Next(biggestBlob.Count())).Value);
  115.                             }
  116.                         }
  117.  
  118.                         // for each starting point, try to fit a grain there
  119.                         ConcurrentBag<Grain> candidates = new ConcurrentBag<Grain>();
  120.                         Parallel.ForEach(startPoints, point =>
  121.                         {
  122.                             Grain candidate = new Grain(point);
  123.                             candidate.Fit(foreground);
  124.                             candidates.Add(candidate);
  125.                         });
  126.                        
  127.                         Grain grain = candidates
  128.                             .OrderByDescending(g=>g.Converged) // we don't want grains where the iterative fit did not finish
  129.                             .ThenBy(g=>g.IsTooSmall) // we don't want tiny grains
  130.                             .ThenByDescending(g => g.CircumferenceRatio) // we want grains that have many edge pixels close to the fitted elipse
  131.                             .ThenBy(g => g.MeanSquaredError)
  132.                             .First(); // we only want the best fit among the 10 candidates
  133.  
  134.                         // count the number of foreground pixels this grain has
  135.                         grain.CountPixel(foreground);
  136.  
  137.                         // remove the grain from the foreground
  138.                         grain.Draw(foreground,CvColor.Black);
  139.  
  140.                         // add the grain to the colection fo found grains
  141.                         grains.Add(grain);
  142.                         grain.Index = grains.Count;
  143.                        
  144.                         // draw the grain for visualisation
  145.                         grain.Draw(display, CvColor.Random());
  146.                         grain.DrawContour(display, CvColor.Random());
  147.                         grain.DrawEllipse(display, CvColor.Random());
  148.  
  149.                         //display.SaveImage("10-foundGrains.png");
  150.                     }
  151.  
  152.                     // throw away really bad grains
  153.                     grains = grains.Where(g => g.PixelRatio >= 0.73).ToList();
  154.  
  155.                     // estimate the average grain size, ignoring outliers
  156.                     double avgGrainSize =
  157.                         grains.OrderBy(g => g.NumPixel).Skip(grains.Count/10).Take(grains.Count*9/10).Average(g => g.NumPixel);
  158.  
  159.                     //ignore the estimated grain size, use a fixed size
  160.                     avgGrainSize = 2530;
  161.  
  162.                     // count the number of grains, using the average grain size
  163.                     double numGrains = grains.Sum(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize));
  164.  
  165.                     // get some statistics
  166.                     double avgWidth = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) == 1).Average(g => g.Width);
  167.                     double avgHeight = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) == 1).Average(g => g.Height);
  168.                     double avgPixelRatio = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) == 1).Average(g => g.PixelRatio);
  169.  
  170.                     int numUndersized = grains.Count(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1);
  171.                     int numOversized = grains.Count(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1);
  172.  
  173.                     double avgWidthUndersized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).Select(g=>g.Width).DefaultIfEmpty(0).Average();
  174.                     double avgHeightUndersized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).Select(g => g.Height).DefaultIfEmpty(0).Average();
  175.                     double avgGrainSizeUndersized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).Select(g => g.NumPixel).DefaultIfEmpty(0).Average();
  176.                     double avgPixelRatioUndersized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).Select(g => g.PixelRatio).DefaultIfEmpty(0).Average();
  177.  
  178.                     double avgWidthOversized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).Select(g => g.Width).DefaultIfEmpty(0).Average();
  179.                     double avgHeightOversized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).Select(g => g.Height).DefaultIfEmpty(0).Average();
  180.                     double avgGrainSizeOversized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).Select(g => g.NumPixel).DefaultIfEmpty(0).Average();
  181.                     double avgPixelRatioOversized = grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).Select(g => g.PixelRatio).DefaultIfEmpty(0).Average();
  182.  
  183.  
  184.                     Console.WriteLine("===============================");
  185.                     Console.WriteLine("Grains: {0}|{1:0.} of {2} (e{3}), size {4:0.}px, {5:0.}x{6:0.}  {7:0.000}  undersized:{8}  oversized:{9}   {10:0.0} minutes  {11:0.0} s per grain",grains.Count,numGrains,expectedGrains[fileNo],expectedGrains[fileNo]-numGrains,avgGrainSize,avgWidth,avgHeight, avgPixelRatio,numUndersized,numOversized,watch.Elapsed.TotalMinutes, watch.Elapsed.TotalSeconds/grains.Count);
  186.  
  187.                     Console.WriteLine("undersized:  {0:0.}px  {1:0.}x{2:0.}", avgGrainSizeUndersized, avgWidthUndersized, avgHeightUndersized);
  188.                     foreach (Grain grain in grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) < 1).OrderBy(g=>g.NumPixel))
  189.                     {
  190.                         Console.WriteLine("          :  {0:0.}px  {1:0.}x{2:0.}, {3:0.000}", grain.NumPixel, grain.Width, grain.Height,grain.PixelRatio);
  191.                     }
  192.                     Console.WriteLine("oversized :  {0:0.}px  {1:0.}x{2:0.}", avgGrainSizeOversized, avgWidthOversized, avgHeightOversized);
  193.                     foreach (Grain grain in grains.Where(g => Math.Round(g.NumPixel * 1.0 / avgGrainSize) > 1).OrderBy(g=>g.NumPixel))
  194.                     {
  195.                         Console.WriteLine("          :  {0:0.}px  {1:0.}x{2:0.} {3:0.000}", grain.NumPixel, grain.Width, grain.Height, grain.PixelRatio);
  196.                     }
  197.  
  198.  
  199.                     // draw the description for each grain
  200.                     foreach (Grain grain in grains)
  201.                     {
  202.                         grain.DrawText(avgGrainSize, display, CvColor.Black);
  203.                     }
  204.  
  205.                     display.SaveImage("10-foundGrains.png");
  206.                     display.SaveImage("X-" + file + "-foundgrains.png");
  207.                 }
  208.             }
  209.         }
  210.     }
  211.  
  212.  
  213.  
  214.     public class Grain
  215.     {
  216.         private const int MIN_WIDTH = 70;
  217.         private const int MAX_WIDTH = 130;
  218.         private const int MIN_HEIGHT = 20;
  219.         private const int MAX_HEIGHT = 35;
  220.  
  221.         private static CvFont font01 = new CvFont(FontFace.HersheyPlain, 0.5, 1);
  222.         private Random random = new Random(4); // determined by fair dice throw; guaranteed to be random
  223.  
  224.  
  225.         /// <summary> center of grain </summary>
  226.         public CvPoint2D32f Position { get; private set; }
  227.         /// <summary> Width of grain (always bigger than height)</summary>
  228.         public float Width { get; private set; }
  229.         /// <summary> Height of grain (always smaller than width)</summary>
  230.         public float Height { get; private set; }
  231.  
  232.         public float MinorRadius { get { return this.Height / 2; } }
  233.         public float MajorRadius { get { return this.Width / 2; } }
  234.         public double Angle { get; private set; }
  235.         public double AngleRad { get { return this.Angle * Math.PI / 180; } }
  236.  
  237.         public int Index { get; set; }
  238.         public bool Converged { get; private set; }
  239.         public int NumIterations { get; private set; }
  240.         public double CircumferenceRatio { get; private set; }
  241.         public int NumPixel { get; private set; }
  242.         public List<EllipsePoint> EdgePoints { get; private set; }
  243.         public double MeanSquaredError { get; private set; }
  244.         public double PixelRatio { get { return this.NumPixel / (Math.PI * this.MajorRadius * this.MinorRadius); } }
  245.         public bool IsTooSmall { get { return this.Width < MIN_WIDTH || this.Height < MIN_HEIGHT; } }
  246.  
  247.         public Grain(CvPoint2D32f position)
  248.         {
  249.             this.Position = position;
  250.             this.Angle = 0;
  251.             this.Width = 10;
  252.             this.Height = 10;
  253.             this.MeanSquaredError = double.MaxValue;
  254.         }
  255.  
  256.         /// <summary>  fit a single rice grain of elipsoid shape </summary>
  257.         public void Fit(CvMat img)
  258.         {
  259.             // distance between the sampled points on the elipse circumference in degree
  260.             int angularResolution = 1;
  261.  
  262.             // how many times did the fitted ellipse not change significantly?
  263.             int numConverged = 0;
  264.  
  265.             // number of iterations for this fit
  266.             int numIterations;
  267.  
  268.             // repeat until the fitted ellipse does not change anymore, or the maximum number of iterations is reached
  269.             for (numIterations = 0; numIterations < 100 && !this.Converged; numIterations++)
  270.             {
  271.                 // points on an ideal ellipse
  272.                 CvPoint[] points;
  273.                 Cv.Ellipse2Poly(this.Position, new CvSize2D32f(MajorRadius, MinorRadius), Convert.ToInt32(this.Angle), 0, 359, out points,
  274.                                 angularResolution);
  275.  
  276.                 // points on the edge of foregroudn to background, that are close to the elipse
  277.                 CvPoint?[] edgePoints = new CvPoint?[points.Length];
  278.  
  279.                 // remeber if the previous pixel in a given direction was foreground or background
  280.                 bool[] prevPixelWasForeground = new bool[points.Length];
  281.  
  282.                 // when the first edge pixel is found, this value is updated
  283.                 double firstEdgePixelOffset = 200;
  284.  
  285.                 // from the center of the elipse towards the outside:
  286.                 for (float offset = -this.MajorRadius + 1; offset < firstEdgePixelOffset + 20; offset++)
  287.                 {
  288.                     // draw an ellipse with the given offset
  289.                     Cv.Ellipse2Poly(this.Position, new CvSize2D32f(MajorRadius + offset, MinorRadius + (offset > 0 ? offset : MinorRadius / MajorRadius * offset)), Convert.ToInt32(this.Angle), 0,
  290.                                     359, out points, angularResolution);
  291.  
  292.                     // for each angle
  293.                     Parallel.For(0, points.Length, i =>
  294.                     {
  295.                         if (edgePoints[i].HasValue) return; // edge for this angle already found
  296.  
  297.                         // check if the current pixel is foreground
  298.                         bool foreground = points[i].X < 0 || points[i].Y < 0 || points[i].X >= img.Cols || points[i].Y >= img.Rows
  299.                                               ? false // pixel outside of image borders is always background
  300.                                               : img.Get2D(points[i].Y, points[i].X).Val0 > 0;
  301.  
  302.  
  303.                         if (prevPixelWasForeground[i] && !foreground)
  304.                         {
  305.                             // found edge pixel!
  306.                             edgePoints[i] = points[i];
  307.  
  308.                             // if this is the first edge pixel we found, remember its offset. the other pixels cannot be too far away, so we can stop searching soon
  309.                             if (offset < firstEdgePixelOffset && offset > 0) firstEdgePixelOffset = offset;
  310.                         }
  311.  
  312.                         prevPixelWasForeground[i] = foreground;
  313.                     });
  314.                 }
  315.  
  316.                 // estimate the distance of each found edge pixel from the ideal elipse
  317.                 // this is a hack, since the actual equations for estimating point-ellipse distnaces are complicated
  318.                 Cv.Ellipse2Poly(this.Position, new CvSize2D32f(MajorRadius, MinorRadius), Convert.ToInt32(this.Angle), 0, 360,
  319.                                 out points, angularResolution);
  320.                 var pointswithDistance =
  321.                     edgePoints.Select((p, i) => p.HasValue ? new EllipsePoint(p.Value, points[i], this.Position) : null)
  322.                               .Where(p => p != null).ToList();
  323.  
  324.                 if (pointswithDistance.Count == 0)
  325.                 {
  326.                     Console.WriteLine("no points found! should never happen! ");
  327.                     break;
  328.                 }
  329.  
  330.                 // throw away all outliers that are too far outside the current ellipse
  331.                 double medianSignedDistance = pointswithDistance.OrderBy(p => p.SignedDistance).ElementAt(pointswithDistance.Count / 2).SignedDistance;
  332.                 var goodPoints = pointswithDistance.Where(p => p.SignedDistance < medianSignedDistance + 15).ToList();
  333.  
  334.                 // do a sort of ransack fit with the inlier points to find a new better ellipse
  335.                 CvBox2D bestfit = ellipseRansack(goodPoints);
  336.  
  337.                 // check if the fit has converged
  338.                 if (Math.Abs(this.Angle - bestfit.Angle) < 3 && // angle has not changed much (<3°)
  339.                     Math.Abs(this.Position.X - bestfit.Center.X) < 3 && // position has not changed much (<3 pixel)
  340.                     Math.Abs(this.Position.Y - bestfit.Center.Y) < 3)
  341.                 {
  342.                     numConverged++;
  343.                 }
  344.                 else
  345.                 {
  346.                     numConverged = 0;
  347.                 }
  348.  
  349.                 if (numConverged > 2)
  350.                 {
  351.                     this.Converged = true;
  352.                 }
  353.  
  354.                 //Console.WriteLine("Iteration {0}, delta {1:0.000} {2:0.000} {3:0.000}    {4:0.000}-{5:0.000} {6:0.000}-{7:0.000} {8:0.000}-{9:0.000}",
  355.                 //  numIterations, Math.Abs(this.Angle - bestfit.Angle), Math.Abs(this.Position.X - bestfit.Center.X), Math.Abs(this.Position.Y - bestfit.Center.Y), this.Angle, bestfit.Angle, this.Position.X, bestfit.Center.X, this.Position.Y, bestfit.Center.Y);
  356.  
  357.                 double msr = goodPoints.Sum(p => p.Distance * p.Distance) / goodPoints.Count;
  358.  
  359.                 // for drawing the polygon, filter the edge points more strongly
  360.                 if (goodPoints.Count(p => p.SignedDistance < 5) > goodPoints.Count / 2)
  361.                     goodPoints = goodPoints.Where(p => p.SignedDistance < 5).ToList();
  362.                 double cutoff = goodPoints.Select(p => p.Distance).OrderBy(d => d).ElementAt(goodPoints.Count * 9 / 10);
  363.                 goodPoints = goodPoints.Where(p => p.SignedDistance <= cutoff + 1).ToList();
  364.  
  365.                 int numCertainEdgePoints = goodPoints.Count(p => p.SignedDistance > -2);
  366.                 this.CircumferenceRatio = numCertainEdgePoints * 1.0 / points.Count();
  367.  
  368.                 this.Angle = bestfit.Angle;
  369.                 this.Position = bestfit.Center;
  370.                 this.Width = bestfit.Size.Width;
  371.                 this.Height = bestfit.Size.Height;
  372.                 this.EdgePoints = goodPoints;
  373.                 this.MeanSquaredError = msr;
  374.  
  375.             }
  376.             this.NumIterations = numIterations;
  377.             //Console.WriteLine("Grain found after {0,3} iterations, size={1,3:0.}x{2,3:0.}   pixel={3,5}    edgePoints={4,3}   msr={5,2:0.00000}", numIterations, this.Width,
  378.             //                        this.Height, this.NumPixel, this.EdgePoints.Count, this.MeanSquaredError);
  379.         }
  380.  
  381.         /// <summary> a sort of ransakc fit to find the best ellipse for the given points </summary>
  382.         private CvBox2D ellipseRansack(List<EllipsePoint> points)
  383.         {
  384.             using (CvMemStorage storage = new CvMemStorage(0))
  385.             {
  386.                 // calculate minimum bounding rectangle
  387.                 CvSeq<CvPoint> fullPointSeq = CvSeq<CvPoint>.FromArray(points.Select(p => p.Point), SeqType.EltypePoint, storage);
  388.                 var boundingRect = fullPointSeq.MinAreaRect2();
  389.  
  390.                 // the initial candidate is the previously found ellipse
  391.                 CvBox2D bestEllipse = new CvBox2D(this.Position, new CvSize2D32f(this.Width, this.Height), (float)this.Angle);
  392.                 double bestError = calculateEllipseError(points, bestEllipse);
  393.  
  394.                 Queue<EllipsePoint> permutation = new Queue<EllipsePoint>();
  395.                 if (points.Count >= 5) for (int i = -2; i < 20; i++)
  396.                     {
  397.                         CvBox2D ellipse;
  398.                         if (i == -2)
  399.                         {
  400.                             // first, try the ellipse described by the boundingg rect
  401.                             ellipse = boundingRect;
  402.                         }
  403.                         else if (i == -1)
  404.                         {
  405.                             // then, try the best-fit ellipsethrough all points
  406.                             ellipse = fullPointSeq.FitEllipse2();
  407.                         }
  408.                         else
  409.                         {
  410.                             // then, repeatedly fit an ellipse through a random sample of points
  411.  
  412.                             // pick some random points
  413.                             if (permutation.Count < 5) permutation = new Queue<EllipsePoint>(permutation.Concat(points.OrderBy(p => random.Next())));
  414.                             CvSeq<CvPoint> pointSeq = CvSeq<CvPoint>.FromArray(permutation.Take(10).Select(p => p.Point), SeqType.EltypePoint, storage);
  415.                             for (int j = 0; j < pointSeq.Count(); j++) permutation.Dequeue();
  416.  
  417.                             // fit an ellipse through these points
  418.                             ellipse = pointSeq.FitEllipse2();
  419.                         }
  420.  
  421.                         // assure that the width is greater than the height
  422.                         ellipse = NormalizeEllipse(ellipse);
  423.  
  424.                         // if the ellipse is too big for agrain, shrink it
  425.                         ellipse = rightSize(ellipse, points.Where(p => isOnEllipse(p.Point, ellipse, 10, 10)).ToList());
  426.  
  427.                         // sometimes the ellipse given by FitEllipse2 is totally off
  428.                         if (boundingRect.Center.DistanceTo(ellipse.Center) > Math.Max(boundingRect.Size.Width, boundingRect.Size.Height) * 2)
  429.                         {
  430.                             // ignore this bad fit
  431.                             continue;
  432.                         }
  433.  
  434.                         // estimate the error
  435.                         double error = calculateEllipseError(points, ellipse);
  436.  
  437.                         if (error < bestError)
  438.                         {
  439.                             // found a better ellipse!
  440.                             bestError = error;
  441.                             bestEllipse = ellipse;
  442.                         }
  443.                     }
  444.  
  445.                 return bestEllipse;
  446.             }
  447.         }
  448.  
  449.         /// <summary> The proper thing to do would be to use the actual distance of each point to the elipse.
  450.         /// However that formula is complicated, so ...  </summary>
  451.         private double calculateEllipseError(List<EllipsePoint> points, CvBox2D ellipse)
  452.         {
  453.             const double toleranceInner = 5;
  454.             const double toleranceOuter = 10;
  455.             int numWrongPoints = points.Count(p => !isOnEllipse(p.Point, ellipse, toleranceInner, toleranceOuter));
  456.             double ratioWrongPoints = numWrongPoints * 1.0 / points.Count;
  457.  
  458.             int numTotallyWrongPoints = points.Count(p => !isOnEllipse(p.Point, ellipse, 10, 20));
  459.             double ratioTotallyWrongPoints = numTotallyWrongPoints * 1.0 / points.Count;
  460.  
  461.             // this pseudo-distance is biased towards deviations on the major axis
  462.             double pseudoDistance = Math.Sqrt(points.Sum(p => Math.Abs(1 - ellipseMetric(p.Point, ellipse))) / points.Count);
  463.  
  464.             // primarily take the number of points far from the elipse border as an error metric.
  465.             // use pseudo-distance to break ties between elipses with the same number of wrong points
  466.             return ratioWrongPoints * 1000  + ratioTotallyWrongPoints+ pseudoDistance / 1000;
  467.         }
  468.  
  469.  
  470.         /// <summary> shrink an ellipse if it is larger than the maximum grain dimensions </summary>
  471.         private static CvBox2D rightSize(CvBox2D ellipse, List<EllipsePoint> points)
  472.         {
  473.             if (ellipse.Size.Width < MAX_WIDTH && ellipse.Size.Height < MAX_HEIGHT) return ellipse;
  474.  
  475.             // elipse is bigger than the maximum grain size
  476.             // resize it so it fits, while keeping one edge of the bounding rectangle constant
  477.  
  478.             double desiredWidth = Math.Max(10, Math.Min(MAX_WIDTH, ellipse.Size.Width));
  479.             double desiredHeight = Math.Max(10, Math.Min(MAX_HEIGHT, ellipse.Size.Height));
  480.  
  481.             CvPoint2D32f average = points.Average();
  482.  
  483.             // get the corners of the surrounding bounding box
  484.             var corners = ellipse.BoxPoints().ToList();
  485.  
  486.             // find the corner that is closest to the center of mass of the points
  487.             int i0 = ellipse.BoxPoints().Select((point, index) => new { point, index }).OrderBy(p => p.point.DistanceTo(average)).First().index;
  488.             CvPoint p0 = corners[i0];
  489.  
  490.             // find the two corners that are neighbouring this one
  491.             CvPoint p1 = corners[(i0 + 1) % 4];
  492.             CvPoint p2 = corners[(i0 + 3) % 4];
  493.  
  494.             // p1 is the next corner along the major axis (widht), p2 is the next corner along the minor axis (height)
  495.             if (p0.DistanceTo(p1) < p0.DistanceTo(p2))
  496.             {
  497.                 CvPoint swap = p1;
  498.                 p1 = p2;
  499.                 p2 = swap;
  500.             }
  501.  
  502.             // calculate the three other corners with the desired widht and height
  503.  
  504.             CvPoint2D32f edge1 = (p1 - p0);
  505.             CvPoint2D32f edge2 = p2 - p0;
  506.             double edge1Length = Math.Max(0.0001, p0.DistanceTo(p1));
  507.             double edge2Length = Math.Max(0.0001, p0.DistanceTo(p2));
  508.  
  509.             CvPoint2D32f newCenter = (CvPoint2D32f)p0 + edge1 * (desiredWidth / edge1Length) + edge2 * (desiredHeight / edge2Length);
  510.  
  511.             CvBox2D smallEllipse = new CvBox2D(newCenter, new CvSize2D32f((float)desiredWidth, (float)desiredHeight), ellipse.Angle);
  512.  
  513.             return smallEllipse;
  514.         }
  515.  
  516.         /// <summary> assure that the width of the elipse is the major axis, and the height is the minor axis.
  517.         /// Swap widht/height and rotate by 90° otherwise  </summary>
  518.         private static CvBox2D NormalizeEllipse(CvBox2D ellipse)
  519.         {
  520.             if (ellipse.Size.Width < ellipse.Size.Height)
  521.             {
  522.                 ellipse = new CvBox2D(ellipse.Center, new CvSize2D32f(ellipse.Size.Height, ellipse.Size.Width), (ellipse.Angle + 90 + 360) % 360);
  523.             }
  524.             return ellipse;
  525.         }
  526.  
  527.         /// <summary> greater than 1 for points outside ellipse, smaller than 1 for points inside ellipse </summary>
  528.         private static double ellipseMetric(CvPoint p, CvBox2D ellipse)
  529.         {
  530.             double theta = ellipse.Angle * Math.PI / 180;
  531.             double u = Math.Cos(theta) * (p.X - ellipse.Center.X) + Math.Sin(theta) * (p.Y - ellipse.Center.Y);
  532.             double v = -Math.Sin(theta) * (p.X - ellipse.Center.X) + Math.Cos(theta) * (p.Y - ellipse.Center.Y);
  533.  
  534.             return u * u / (ellipse.Size.Width * ellipse.Size.Width / 4) + v * v / (ellipse.Size.Height * ellipse.Size.Height / 4);
  535.         }
  536.  
  537.         /// <summary> Is the point on the ellipseBorder, within a certain tolerance </summary>
  538.         private static bool isOnEllipse(CvPoint p, CvBox2D ellipse, double toleranceInner, double toleranceOuter)
  539.         {
  540.             double theta = ellipse.Angle * Math.PI / 180;
  541.             double u = Math.Cos(theta) * (p.X - ellipse.Center.X) + Math.Sin(theta) * (p.Y - ellipse.Center.Y);
  542.             double v = -Math.Sin(theta) * (p.X - ellipse.Center.X) + Math.Cos(theta) * (p.Y - ellipse.Center.Y);
  543.  
  544.             double innerEllipseMajor = (ellipse.Size.Width - toleranceInner) / 2;
  545.             double innerEllipseMinor = (ellipse.Size.Height - toleranceInner) / 2;
  546.             double outerEllipseMajor = (ellipse.Size.Width + toleranceOuter) / 2;
  547.             double outerEllipseMinor = (ellipse.Size.Height + toleranceOuter) / 2;
  548.  
  549.             double inside = u * u / (innerEllipseMajor * innerEllipseMajor) + v * v / (innerEllipseMinor * innerEllipseMinor);
  550.             double outside = u * u / (outerEllipseMajor * outerEllipseMajor) + v * v / (outerEllipseMinor * outerEllipseMinor);
  551.             return inside >= 1 && outside <= 1;
  552.         }
  553.  
  554.  
  555.         /// <summary> count the number of foreground pixels for this grain </summary>
  556.         public int CountPixel(CvMat img)
  557.         {
  558.             // todo: this is an incredibly inefficient way to count, allocating a new image with the size of the input each time
  559.             using (CvMat mask = new CvMat(img.Rows, img.Cols, MatrixType.U8C1))
  560.             {
  561.                 mask.SetZero();
  562.                 mask.FillPoly(new CvPoint[][] { this.EdgePoints.Select(p => p.Point).ToArray() }, CvColor.White);
  563.                 mask.And(img, mask);
  564.                 this.NumPixel = mask.CountNonZero();
  565.             }
  566.             return this.NumPixel;
  567.         }
  568.  
  569.         /// <summary> draw the recognized shape of the grain </summary>
  570.         public void Draw(CvMat img, CvColor color)
  571.         {
  572.             img.FillPoly(new CvPoint[][] { this.EdgePoints.Select(p => p.Point).ToArray() }, color);
  573.         }
  574.  
  575.         /// <summary> draw the contours of the grain </summary>
  576.         public void DrawContour(CvMat img, CvColor color)
  577.         {
  578.             img.DrawPolyLine(new CvPoint[][] { this.EdgePoints.Select(p => p.Point).ToArray() }, true, color);
  579.         }
  580.  
  581.         /// <summary> draw the best-fit ellipse of the grain </summary>
  582.         public void DrawEllipse(CvMat img, CvColor color)
  583.         {
  584.             img.DrawEllipse(this.Position, new CvSize2D32f(this.MajorRadius, this.MinorRadius), this.Angle, 0, 360, color, 1);
  585.         }
  586.  
  587.         /// <summary> print the grain index and the number of pixels divided by the average grain size</summary>
  588.         public void DrawText(double averageGrainSize, CvMat img, CvColor color)
  589.         {
  590.             img.PutText(String.Format("{0}|{1:0.0}", this.Index, this.NumPixel / averageGrainSize), this.Position + new CvPoint2D32f(-5, 10), font01, color);
  591.         }
  592.  
  593.     }
  594.  
  595.     public static class Extensions
  596.     {
  597.         public static double DistanceTo(this CvPoint2D32f p1, CvPoint2D32f p2)
  598.         {
  599.             return Math.Sqrt(Math.Pow((p2.X - p1.X), 2.0) + Math.Pow((p2.Y - p1.Y), 2.0));
  600.         }
  601.         public static double DistanceTo(this CvPoint p1, CvPoint2D32f p2)
  602.         {
  603.             return Math.Sqrt(Math.Pow((p2.X - p1.X), 2.0) + Math.Pow((p2.Y - p1.Y), 2.0));
  604.         }
  605.         public static CvPoint2D32f Average(this IEnumerable<EllipsePoint> points)
  606.         {
  607.             float x = 0, y = 0;
  608.             int i = 0;
  609.             foreach (var point in points)
  610.             {
  611.                 x += point.Point.X;
  612.                 y += point.Point.Y;
  613.                 i++;
  614.             }
  615.             return new CvPoint2D32f(x / i, y / i);
  616.         }
  617.  
  618.         public static double Median<T>(this IEnumerable<T> list, Func<T, double> selector)
  619.         {
  620.             var sorted = list.Select(selector).OrderBy(d => d).ToList();
  621.             return sorted.ElementAt(sorted.Count / 2);
  622.  
  623.         }
  624.  
  625.     }
  626.  
  627.     /// <summary> Point with signed distance to a given ellipse  
  628.     /// </summary>
  629.     public class EllipsePoint
  630.     {
  631.         public CvPoint Point { get; private set; }
  632.         public double SignedDistance { get; private set; }
  633.         public double Distance { get { return Math.Abs(this.SignedDistance); } }
  634.  
  635.         public EllipsePoint(CvPoint point, CvPoint2D32f closestPointOnElipse, CvPoint2D32f ellipseCenter)
  636.         {
  637.             this.Point = point;
  638.  
  639.             double distance = point.DistanceTo(closestPointOnElipse);
  640.             bool isInner = point.DistanceTo(ellipseCenter) < closestPointOnElipse.DistanceTo(ellipseCenter);
  641.             this.SignedDistance = isInner ? -distance : distance;
  642.         }
  643.     }
  644. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement