# HG changeset patch
# User Hartmut Gimpel <hg_code@gmx.de>
# Date 1486418265 -3600
#      Mon Feb 06 22:57:45 2017 +0100
# Node ID 9f28da48350e3f7b0de24dc15cad2546f487944f
# Parent  ec1a16084bfc2849d050b1dd8c959543faf60c42
hough.m: new m-file wrapper around hough_line.cc, Matlab compatible, with tests

diff -r ec1a16084bfc -r 9f28da48350e NEWS
--- a/NEWS	Sun Jan 29 16:06:35 2017 +0000
+++ b/NEWS	Mon Feb 06 22:57:45 2017 +0100
@@ -8,6 +8,10 @@
     single if the input was class single (they were previously
     returning a double for such cases).
 
+ ** The following functions are new:
+
+      hough
+
 
  Summary of important user-visible changes for image 2.6.1 (2016/10/21):
 -------------------------------------------------------------------------
diff -r ec1a16084bfc -r 9f28da48350e inst/hough.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/hough.m	Mon Feb 06 22:57:45 2017 +0100
@@ -0,0 +1,280 @@
+## Copyright (C) 2017 Hartmut Gimpel <hg_code@gmx.de>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {[@var{H}, @var{theta}, @var{rho}] =} hough (@var{BW})
+## @deftypefnx  {Function File} {[@var{H}, @var{theta}, @var{rho}] =} hough (@var{BW}, @var{property}, @var{value}, @dots{})
+## Compute the Hough transform to find lines in a binary image @var{BW}.
+##
+## The input image @var{BW} needs to be logical or only have values of 0s and 1s.
+##
+## The resulting Hough transform matrix @var{H} (accumulator array) is 2D. Its rows
+## correspond to the distance values @var{rho} and its columns to the angle
+## values @var{theta}. Points of high value in @var{H} correspond to present lines 
+## in the given image.
+##
+## The distance @var{rho} is measured with respect to the image origin. The angle @var{theta}
+## is measured clockwise to the vertical axis (this is Matlab compatible).
+##
+## The following @var{property} values are possible:
+## 
+## @table
+## @item "Theta"
+## A vector of angle values. The Hough transform will be calculated at 
+## those Theta angle values. The angles are given in degrees. 
+## This defaults to [-90:89].
+##
+## @item "ThetaResolution" (not recommended in Matlab)
+## A scalar value to specify the Theta angles for the Hough transform in a different way. 
+## This will result in Theta = [-90:ThetaResolution:90], with the +90° angle excluded.
+##
+## @item "RhoResolution"
+## A scalar value to specify the Resolution of the @var{rho} distance values.
+## This is currently not implemented in Octave. It defaults to the value 1.
+##
+## @end table
+##
+## @seealso{hough_line, hough_circle, immaximas}
+## @end deftypefn
+
+%% TODO: correct help text in hough_line.cc (Radian statt Grad!)
+
+function [H, theta, rho] = hough (bw, varargin)
+    
+  ## check input image, try to do obvious conversion to logical  
+  if (! (ismatrix (bw) & islogical(bw)))
+    if (all (all (logical (bw) == bw)))
+      bw = logical (bw);
+    else
+      error ("hough: BW must be a binary image");
+    endif
+  endif
+  
+  ## set default parameters:
+  theta = [-90:1:89];
+  theta_res = 1;
+  
+  ## process property/value pairs
+  if (rem (numel (varargin), 2) != 0)
+    error ("hough: PROPERTY/VALUE arguments must occur in pairs");
+  endif
+  
+  for idx = 1:2:(numel (varargin))
+    switch (tolower (varargin{idx}))
+
+      case "rhoresolution"
+          warning ("hough: parameter RHORESOLUTION is not implemented. A default value of 1 will be used instead.");
+      
+      case "thetaresolution"
+        theta_res = varargin{idx+1};
+        if (! (isreal (theta_res) & isscalar (theta_res) & (theta_res > 0) & (theta_res < 180)))
+          error ("hough: value THETARESOLUTION must be between 0 and 180");
+        endif
+        
+      case "theta"
+        theta = varargin{idx+1};
+        if (! (isreal (theta) & isvector (theta)))
+          error ("hough: values THETA must be a vector of real numbers");
+        endif
+        
+     endswitch
+   endfor
+
+   ## set theta values
+   if (theta_res != 1)
+     theta = [-90:theta_res:90];
+     theta = theta(theta != 90);    # exclude +90 degrees
+   endif
+   
+   ## Matlab's hough.m function measures the angle theta clockwise to the
+   ## vertical axis, but Octave's hough_line.cc measure this angle (as usual)
+   ## counter-clockwise the the horizontal axis. Octave's hough_line also
+   ## need radians, Matlab uses degrees. So we translate this.
+   theta_oct = (-theta+90) *  (pi/180); 
+   
+   # eventually call hough_line.cc to do the real work
+   [H, rho] = hough_line(bw, theta_oct);
+  
+endfunction
+
+%!shared BW0, BW1, BW2, BWx, BWy
+%!
+%! BW0 = false (5);
+%! BW0(2,2) = true;
+%!
+%! BW1 = zeros (100, 100);
+%! BW1(1,1) = 1;
+%! BW1(100,100) = 1;
+%! BW1(1,100) = 1;
+%! BW1(100, 1) = 1;
+%! BW1(50,50) = 1;
+%!
+%! n = 100;
+%! BW2 = false (n);
+%! a = 50;    # line starts at left side at row a
+%! b = 3;      # slope of line is 1:b
+%! for column = 1:n
+%!   if rem (column, b) == 0
+%!     row = a - column/b;
+%!     BW2(row, column) = true;
+%!   endif;
+%! endfor;
+%!
+%! BWx = false (10);
+%! BWx(:,5) = true;
+%!
+%! BWy = false (10);
+%! BWy(5,:) = true;
+
+%!test
+%! [H, T, R] = hough (BW1);
+%! assert (size (H) == [283 180]);
+
+%!test
+%! [H, T, R] = hough (BW1, 'Theta', [-90 0 45 79]);
+%! assert (size (H) == [283 4]);
+
+%!test
+%! [H, T, R] = hough (BW1, 'ThetaResolution', 0.5);
+%! assert (size (H) == [283 360]);
+
+%!error <must be a binary image>
+%! I = 5.* ones  (10);
+%! I(3,3) = 2;
+%! [H, T, R] = hough (I);
+
+%!error <must occur in pairs>
+%! [H, T, R] = hough (BW0, 'Theta');
+
+%!error <must be a vector of real numbers>
+%! [H, T, R] = hough (BW0, 'Theta', ones (10));
+
+%!error <must be a vector of real numbers>
+%! [H, T, R] = hough (BW0, 'Theta', [5 -i 7]);
+
+%!warning <default value of 1>
+%! [H, T, R] = hough (BW0, 'RhoResolution', 0.5);
+
+%!test
+%! [H, theta, rho] = hough (BW2);
+%! H_max = max (H(:));
+%! H_size = size (H);
+%! [~, max_idx_lin] = max (H(:));
+%! [max_row, max_column] = ind2sub (size (H), max_idx_lin);
+%! theta_max = theta(max_column);
+%! rho_max = rho(max_row);
+%! assert (H_max , 33);
+%! assert (H_size, [283 180]);
+%! assert (max_row, 188);
+%! assert (max_column, 163);
+%! assert (theta_max, 72);
+%! assert (rho_max, 46);
+
+%!test
+%! [H, theta, rho] = hough (BW2, 'Theta', [65:1:75]);
+%! H_max = max (H(:));
+%! H_size = size (H);
+%! [~, max_idx_lin] = max (H(:));
+%! [max_row, max_column] = ind2sub (size (H), max_idx_lin);
+%! theta_max = theta(max_column);
+%! rho_max = rho(max_row);
+%! assert (H_max , 33);
+%! assert (H_size, [283 11]);
+%! assert (max_row, 188);
+%! assert (max_column, 8);
+%! assert (theta_max, 72);
+%! assert (rho_max, 46);
+
+%!test
+%! [H, theta, rho] = hough (BW2, 'Theta', [-90:0.5:89.5]);
+%! H_max = max (H(:));
+%! H_size = size (H);
+%! [~, max_idx_lin] = max (H(:));
+%! [max_row, max_column] = ind2sub (size (H), max_idx_lin);
+%! theta_max = theta(max_column);
+%! rho_max = rho(max_row);
+%! assert (H_max , 33);
+%! assert (H_size, [283 360]);
+%! assert (max_row, 188);
+%! assert (max_column, 324);
+%! assert (theta_max, 71.5);
+%! assert (rho_max, 46);
+
+%!test
+%! [H, theta, rho] = hough (BW2, 'ThetaResolution', 0.5);
+%! H_max = max (H(:));
+%! H_size = size (H);
+%! [~, max_idx_lin] = max (H(:));
+%! [max_row, max_column] = ind2sub (size (H), max_idx_lin);
+%! theta_max = theta(max_column);
+%! rho_max = rho(max_row);
+%! assert (H_max , 33);
+%! assert (H_size, [283 360]);
+%! assert (max_row, 188);
+%! assert (max_column, 324);
+%! assert (theta_max, 71.5);
+%! assert (rho_max, 46);
+
+%!test
+%! [H, theta, rho] = hough (BWx);
+%! H_max = max (H(:));
+%! [~, max_idx_lin] = max (H(:));
+%! [max_row, max_column] = ind2sub (size (H), max_idx_lin);
+%! theta_max = theta(max_column);
+%! rho_max = rho(max_row);
+%! assert (H_max , 10);
+%! assert (max_column, 88);
+%! assert (theta_max, -3);
+%! assert (rho_max, 4);
+
+%!xtest    # those ones are currently not Matlab compatible
+%! [H, theta, rho] = hough (BWx);
+%! H_size = size (H);
+%! [~, max_idx_lin] = max (H(:));
+%! [max_row, max_column] = ind2sub (size (H), max_idx_lin);
+%! assert (H_size, [27 180]);
+%! assert (max_row, 18);
+
+%!test
+%! [H, theta, rho] = hough (BWy);
+%! H_max = max (H(:));
+%! [~, max_idx_lin] = max (H(:));
+%! [max_row, max_column] = ind2sub (size (H), max_idx_lin);
+%! theta_max = theta(max_column);
+%! rho_max = rho(max_row);
+%! assert (H_max , 10);
+%! assert (max_column, 1);
+%! assert (theta_max, -90);
+%! assert (rho_max, -4);
+
+%!xtest    # those ones are currently not Matlab compatible
+%! [H, theta, rho] = hough (BWy);
+%! H_size = size (H);
+%! [~, max_idx_lin] = max (H(:));
+%! [max_row, max_column] = ind2sub (size (H), max_idx_lin);
+%! assert (H_size, [27 180]);
+%! assert (max_row, 11);
+
+%!demo
+%! BW = zeros (100, 150);
+%! BW(30,:) = 1; BW(:, 65) = 1; BW(35:45, 35:50) = 1;
+%! for i = 1:90, BW(i,i) = 1; endfor
+%! BW = imnoise (BW, 'salt & pepper');
+%! figure, imshow (BW); title ('BW');
+%! [H, theta, rho] = hough (BW); H = H / max (H(:));
+%! figure, imshow (H, 'XData', theta,'YData', rho);
+%! title ('hough transform of BW');
+%! axis on; xlabel('angle \theta [degrees]'); 
+%! ylabel('distance \rho to origin [pixels]');
diff -r ec1a16084bfc -r 9f28da48350e src/hough_line.cc
--- a/src/hough_line.cc	Sun Jan 29 16:06:35 2017 +0000
+++ b/src/hough_line.cc	Mon Feb 06 22:57:45 2017 +0100
@@ -28,12 +28,14 @@
 @deftypefn {Loadable Function} {[@var{H}, @var{R}] =} hough_line(@var{I}, @var{angles})\n\
 Calculate the straight line Hough transform of a binary image @var{I}.\n\
 \n\
-The angles are given in degrees and defaults to -90:90.\n\
+The angles are given in radians and default to -pi/2:pi/2.\n\
 \n\
 @var{H} is the resulting Hough transform, and @var{R} is the radial distances.\n\
 \n\
 The algorithm is described in\n\
 Digital Image Processing by Gonzales & Woods (2nd ed., p. 587)\n\
+\n\
+For a Matlab compatible Hough transform see hough.m\n\
 @end deftypefn\n\
 ")
 {
