# HG changeset patch
# User Hartmut Gimpel <hg_code@gmx.de>
# Date 1424865582 -3600
#      Wed Feb 25 12:59:42 2015 +0100
# Node ID a5c0e6e3b4031c0036228fcd114457f718bd3667
# Parent  e4ced9d8db353e8bf2b2e5b6c71365455762c6b2
add 4 ellipse properties to regionprops.m (bug #44100)

diff -r e4ced9d8db35 -r a5c0e6e3b403 inst/regionprops.m
--- a/inst/regionprops.m	Mon Feb 23 20:29:17 2015 +0200
+++ b/inst/regionprops.m	Wed Feb 25 12:59:42 2015 +0100
@@ -44,6 +44,10 @@
 ## @item "Centroid"
 ## The center coordinate of the object.
 ##
+## @item "Eccentricity"
+## The eccentricity of the ellipse that has the same normalized 
+## second central moments as the object (value between 0 and 1).
+##
 ## @item "EulerNumber"
 ## @itemx "euler_number"
 ## The Euler number of the object (see @code{bweuler} for details).
@@ -63,6 +67,11 @@
 ## @item "Image"
 ## An image with the same size as the bounding box that contains the original pixels.
 ##
+## @item "MajorAxisLength"
+## @itemx "major_axis_length"
+## The length of the major axis of the ellipse that has the same 
+## normalized second central moments as the object.
+##
 ## @item "MaxIntensity"
 ## @itemx "max_intensity"
 ## The maximum intensity inside the object.
@@ -75,6 +84,16 @@
 ## @itemx "min_intensity"
 ## The minimum intensity inside the object.
 ##
+## @item "MinorAxisLength"
+## @itemx "minor_axis_length"
+## The length of the minor axis of the ellipse that has the same 
+## normalized second central moments as the object.
+##
+## @item "Orientation"
+## The angle between the x-axis and the major axis of the ellipse that  
+## has the same normalized second central moments as the object 
+## (value in degrees between -90 and 90).
+##
 ## @item "Perimeter"
 ## The length of the boundary of the object.
 ##
@@ -145,7 +164,7 @@
                "Centroid", "PixelIdxList", "FilledArea", "PixelList",...
                "FilledImage", "Image", "MaxIntensity", "MinIntensity",...
                "WeightedCentroid", "MeanIntensity", "PixelValues",...
-               "Orientation"};
+               "Orientation", "MajorAxisLength", "MinorAxisLength","Eccentricity"};
 
   if (ismember ("basic", properties))
     properties = union (properties, {"Area", "Centroid", "BoundingBox"});
@@ -290,7 +309,67 @@
           retval (k).PixelValues = I(L == k)(:);
         endfor
 
-      case "orientation"
+     case "majoraxislength" 
+        if (N > 2)
+          warning ("regionprops: skipping majoraxislength for Nd image");
+          break
+        endif
+        
+        for k = 1:num_labels
+          [Y, X] = find (L == k);        
+          
+          if (numel(Y) > 1)
+            # calculate (centralised) second moment of region
+            C = cov ([X(:), Y(:)], 1);    # option 1 for normalisation with n instead of n-1
+            C = C + 1/12.*eye(size(C,1)); # centralised second moment of 1 pixel is 1/12
+            lambda = eig (C);
+            retval (k).MajorAxisLength =4 * sqrt(max(lambda));            
+          else
+            retval (k).MajorAxisLength = 1;
+          endif
+        endfor     
+        
+     case "minoraxislength" 
+         if (N > 2)
+          warning ("regionprops: skipping minoraxislength for Nd image");
+          break
+        endif
+        
+        for k = 1:num_labels
+          [Y, X] = find (L == k);
+          if (numel(Y) > 1)
+            # see case majoraxislength for explanations
+            C = cov ([X(:), Y(:)], 1);         
+            C = C + 1/12.*eye(size(C,1)); 
+            lambda = eig (C);            
+            retval (k).MinorAxisLength =4 * sqrt(min(lambda));
+          else
+            retval (k).MinorAxisLength = 1;
+          endif
+        endfor    
+        
+     case "eccentricity"
+         if (N > 2)
+          warning ("regionprops: skipping eccentricity for Nd image");
+          break
+        endif
+        
+        for k = 1:num_labels
+          [Y, X] = find (L == k);
+          if (numel(Y) > 1)
+            # see case majoraxislength for explanations
+            C = cov ([X(:), Y(:)], 1);        
+            C = C + 1/12.*eye(size(C,1)); 
+            lambda = eig (C);
+            major =  2 * sqrt(max(lambda));
+            minor =  2 * sqrt(min(lambda));
+            retval (k).Eccentricity = sqrt(1-(minor/major)^2);
+          else
+            retval (k).Eccentricity = 0; # a circle has 0 eccentricity
+          endif
+        endfor     
+        
+     case "orientation"
         if (N > 2)
           warning ("regionprops: skipping orientation for Nd image");
           break
@@ -299,49 +378,24 @@
         for k = 1:num_labels
           [Y, X] = find (L == k);
           if (numel (Y) > 1)
-            C = cov ([X(:), Y(:)]);
+            # see case majoraxislength for explanations
+            C = cov ([X(:), Y(:)], 1);
+            C = C + 1/12.*eye(size(C,1)); 
             [V, lambda] = eig (C);
             [max_val, max_idx] = max (diag (lambda));
-            v = V (:, max_idx);
-            retval (k).Orientation = 180 - 180 * atan2 (v (2), v (1)) / pi;
+            max_vec = V(:, max_idx);
+            retval (k).Orientation = -(180/pi)  * atan (max_vec (2) / max_vec (1));      
           else
-            retval (k).Orientation = 0; # XXX: What does the other brand do?
+            retval (k).Orientation = 0; 
           endif
         endfor
 
-      %{
-      case "majoraxislength"
-        for k = 1:num_labels
-          [Y, X] = find (L == k);
-          if (numel (Y) > 1)
-            C = cov ([X(:), Y(:)]);
-            lambda = eig (C);
-            retval (k).MajorAxisLength = (max (lambda));
-          else
-            retval (k).MajorAxisLength = 1;
-          endif
-        endfor
-
-      case "minoraxislength"
-        for k = 1:num_labels
-          [Y, X] = find (L == k);
-          if (numel (Y) > 1)
-            C = cov ([X(:), Y(:)]);
-            lambda = eig (C);
-            retval (k).MinorAxisLength = (min (lambda));
-          else
-            retval (k).MinorAxisLength = 1;
-          endif
-        endfor
-      %}
-
       #case "extrema"
       #case "convexarea"
       #case "convexhull"
       #case "solidity"
       #case "conveximage"
       #case "subarrayidx"
-      #case "eccentricity"
       #case "equivdiameter"
 
       otherwise
@@ -426,4 +480,48 @@
 %! assert (props(1).WeightedCentroid(1), x, 10*eps)
 %! assert (props(1).WeightedCentroid(2), 2, 10*eps)
 
-%!assert (size (regionprops ([1 0 0; 0 0 2], "Area")), [2, 1])
\ No newline at end of file
+%!assert (size (regionprops ([1 0 0; 0 0 2], "Area")), [2, 1])
+
+%!test
+%! a = eye(4);
+%! t = regionprops(a, 'majoraxislength');
+%! assert(t.MajorAxisLength, 6.4291,  1e-3);
+%! b = ones(5);
+%! t = regionprops(b, 'majoraxislength');
+%! assert(t.MajorAxisLength, 5.7735 , 1e-3);
+%! c = [0 0 1; 0 1 1; 1 1 0];
+%! t = regionprops(c, 'majoraxislength');
+%! assert(t.MajorAxisLength, 4.1633 , 1e-3);
+
+%!test
+%! a = eye(4);
+%! t = regionprops(a, 'minoraxislength');
+%! assert(t.MinorAxisLength, 1.1547 , 1e-3);
+%! b = ones(5);
+%! t = regionprops(b, 'minoraxislength');
+%! assert(t.MinorAxisLength, 5.7735 , 1e-3);
+%! c = [0 0 1; 0 1 1; 1 1 0];
+%! t = regionprops(c, 'minoraxislength');
+%! assert(t.MinorAxisLength, 1.8037 , 1e-3);
+
+%!test
+%! a = eye(4);
+%! t = regionprops(a, 'eccentricity');
+%! assert(t.Eccentricity, 0.98374 , 1e-3);
+%! b = ones(5);
+%! t = regionprops(b, 'eccentricity');
+%! assert(t.Eccentricity, 0 , 1e-3);
+%! c = [0 0 1; 0 1 1; 1 1 0];
+%! t = regionprops(c, 'eccentricity');
+%! assert(t.Eccentricity, 0.90128 , 1e-3);
+
+%!test
+%! a = eye(4);
+%! t = regionprops(a, 'orientation');
+%! assert(t.Orientation, -45 , 1e-3);
+%! b = ones(5);
+%! t = regionprops(b, 'orientation');
+%! assert(t.Orientation, 0 , 1e-3);
+%! c = [0 0 1; 0 1 1; 1 1 0];
+%! t = regionprops(c, 'orientation');
+%! assert(t.Orientation, 45 , 1e-3);
