1 /* 2 Copyright (c) 2019- Ferhat Kurtulmuş 3 Boost Software License - Version 1.0 - August 17th, 2003 4 Permission is hereby granted, free of charge, to any person or organization 5 obtaining a copy of the software and accompanying documentation covered by 6 this license (the "Software") to use, reproduce, display, distribute, 7 execute, and transmit the Software, and to prepare derivative works of the 8 Software, and to permit third-parties to whom the Software is furnished to 9 do so, all subject to the following: 10 The copyright notices in the Software and this entire statement, including 11 the above license grant, this restriction and the following disclaimer, 12 must be included in all copies of the Software, in whole or in part, and 13 all derivative works of the Software, unless such copies or derivative 14 works are solely in the form of machine-executable object code generated by 15 a source language processor. 16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 17 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 18 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT 19 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE 20 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, 21 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 22 DEALINGS IN THE SOFTWARE. 23 */ 24 25 module measure.ellipsefit; 26 27 import std.stdio; 28 import std.math; 29 import std.range; 30 import std.algorithm; 31 import std.algorithm.searching; 32 33 import measure.types; 34 import measure.regionprops; 35 import measure.eigen; 36 37 // gives exact results of matlab 38 Ellipse ellipseFit(XYList xylist) 39 { 40 const double meanX = xylist.xs.mean; 41 const double meanY = xylist.ys.mean; 42 43 auto cov = covNx2(xylist.xs, xylist.ys); 44 45 auto eigSolver = new Eig(cov); 46 auto V = eigSolver.getV(); 47 auto lambda = eigSolver.getD(); 48 49 auto lambda_d = lambda.diag(); 50 foreach(ref val; lambda_d) val = 4*sqrt(val); 51 52 const double minorAxisLength = minElement(lambda_d); 53 //const auto minor_idx = minIndex(lambda_d); 54 const double majorAxisLength = maxElement(lambda_d); 55 auto major_idx = maxIndex(lambda_d); 56 57 58 double[2] major_vec = [V[1, major_idx], V[0, major_idx]] ; 59 60 double orientation; 61 if (majorAxisLength == minorAxisLength) 62 orientation = 0; 63 else if (major_vec[1] == 0) 64 orientation = 90; 65 else 66 orientation = -(180/PI) * atan (major_vec[0] / major_vec[1]); 67 68 return Ellipse(orientation, meanX + 1, meanY + 1, majorAxisLength, minorAxisLength); 69 } 70 71 Mat2D!double covNx2(U)(U[] xs, U[] ys) 72 { 73 const auto npoints = xs.length; 74 const auto _npoints = ys.length; 75 76 assert(npoints == _npoints, "array sizes are mismatch!"); 77 78 auto _cov = Mat2D!double(npoints, 2); 79 80 const double meanX = xs.mean; 81 const double meanY = ys.mean; 82 83 foreach(i; 0..npoints) 84 { 85 _cov[i, 0] = xs[i] - meanX; 86 _cov[i, 1] = ys[i] - meanY; 87 } 88 89 auto rescov = _cov.transpose() * _cov; 90 double[2] d; d[0..$] = 1.0/12; 91 92 auto covv = Mat2D!double.diagFromArray(d); 93 94 for(size_t i = 0; i < covv.rows; i++) 95 for(size_t j = 0; j < covv.cols; j++) 96 { 97 rescov[i, j] = rescov[i, j] * (1.0/cast(double)npoints); 98 rescov[i, j] = rescov[i, j] + covv[i, j]; 99 } 100 return rescov; 101 } 102 103 // gives exact results of matlab and 1 ms faster :) 104 Ellipse ellipseFit2(Region rg) @nogc 105 { 106 const double center_x = rg.m10/rg.m00; 107 const double center_y = rg.m01/rg.m00; 108 109 // central moments 110 const double a = rg.m20/rg.m00 - center_x^^2 + 1.0/12; 111 const double b = 2*(rg.m11/rg.m00 - center_x*center_y); 112 const double c = rg.m02/rg.m00 - center_y^^2 + 1.0/12; 113 114 const double theta = 0.5*atan(b/(a-c)) + (a<c)*PI/2; 115 116 const double axis1 = sqrt(8*(a+c-sqrt(b^^2+(a-c)^^2)))/2; 117 const double axis2 = sqrt(8*(a+c+sqrt(b^^2+(a-c)^^2)))/2; 118 119 double orientation; 120 if (axis1 == axis2) 121 orientation = 0; 122 else 123 { 124 orientation = -(180/PI) * theta; 125 if (abs(orientation) > 90) 126 orientation = 180 - abs(orientation); 127 } 128 129 return Ellipse(orientation, rg.bBox.x + center_x + 1, rg.bBox.y + center_y + 1, 2*axis2, 2*axis1); 130 }