detect_ignition/ellipse_fit.m

   1 function [e ] = ellipse_fit( data,ci )\r
   2 %UNTITLED Summary of this function goes here\r
   3 %   Detailed explanation goes here\r
   4 % Calculate the eigenvectors and eigenvalues\r
   5 covariance = cov(data);\r
   6 [eigenvec, eigenval ] = eig(covariance);\r
   7 \r
   8 % Get the index of the largest eigenvector\r
   9 [largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));\r
  10 largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);\r
  11 \r
  12 % Get the largest eigenvalue\r
  13 largest_eigenval = max(max(eigenval));\r
  14 \r
  15 % Get the smallest eigenvector and eigenvalue\r
  16 if(largest_eigenvec_ind_c == 1)\r
  17     smallest_eigenval = max(eigenval(:,2));\r
  18     smallest_eigenvec = eigenvec(:,2);\r
  19 else\r
  20     smallest_eigenval = max(eigenval(:,1));\r
  21     smallest_eigenvec = eigenvec(1,:);\r
  22 end\r
  23 \r
  24 % Calculate the angle between the x-axis and the largest eigenvector\r
  25 angle = atan2(largest_eigenvec(2), largest_eigenvec(1));\r
  26 \r
  27 % This angle is between -pi and pi.\r
  28 % Let's shift it such that the angle is between 0 and 2pi\r
  29 if(angle < 0)\r
  30     angle = angle + 2*pi;\r
  31 end\r
  32 \r
  33 % Get the coordinates of the data mean\r
  34 avg = mean(data);\r
  35 \r
  36 % Get the 95% confidence interval error ellipse\r
  37 %chisquare_val = 2.4477;\r
  38 %chisquare_val = 2.2;\r
  39 chisquare_val = ci;\r
  40 theta_grid = linspace(0,2*pi);\r
  41 phi = angle;\r
  42 X0=avg(1);\r
  43 Y0=avg(2);\r
  44 a=chisquare_val*sqrt(largest_eigenval);\r
  45 b=chisquare_val*sqrt(smallest_eigenval);\r
  46 \r
  47 % the ellipse in x and y coordinates \r
  48 ellipse_x_r  = a*cos( theta_grid );\r
  49 ellipse_y_r  = b*sin( theta_grid );\r
  50 \r
  51 %Define a rotation matrix\r
  52 R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];\r
  53 \r
  54 %let's rotate the ellipse to some angle phi\r
  55 r_ellipse = [ellipse_x_r;ellipse_y_r]' * R;\r
  56 \r
  57 % Draw the error ellipse\r
  58 % axis square\r
  59 figure(1)\r
  60 plot(r_ellipse(:,1) + X0,r_ellipse(:,2) + Y0,'-')\r
  61 hold on;\r
  62 \r
  63 % Plot the original data\r
  64 plot(data(:,1), data(:,2), '.');\r
  65 mindata = min(min(data));\r
  66 maxdata = max(max(data));\r
  67 %xlim([mindata-3, maxdata+3]);\r
  68 %ylim([mindata-3, maxdata+3]);\r
  69 hold on;\r
  70 \r
  71 % Plot the eigenvectors\r
  72 quiver(X0, Y0, largest_eigenvec(1)*sqrt(largest_eigenval), largest_eigenvec(2)*sqrt(largest_eigenval), '-m', 'LineWidth',2);\r
  73 quiver(X0, Y0, smallest_eigenvec(1)*sqrt(smallest_eigenval), smallest_eigenvec(2)*sqrt(smallest_eigenval), '-g', 'LineWidth',2);\r
  74 hold on;\r
  75 \r
  76 % Set the axis labels\r
  77 hXLabel = xlabel('x');\r
  78 hYLabel = ylabel('y');\r
  79 \r
  80 e.v = eigenvec;\r
  81 e.d = eigenval;\r
  82 e.center = [X0,Y0];\r
  83 e.ellipse = r_ellipse;\r
  84 \r
  85 end\r
  86 \r