1 // Copyright (C) 2004 Stefan van der Walt <stefan@sun.ac.za> 2 // All rights reserved. 3 // 4 // Redistribution and use in source and binary forms, with or without 5 // modification, are permitted provided that the following conditions are met: 6 // 7 // 1 Redistributions of source code must retain the above copyright notice, 8 // this list of conditions and the following disclaimer. 9 // 2 Redistributions in binary form must reproduce the above copyright 10 // notice, this list of conditions and the following disclaimer in the 11 // documentation and/or other materials provided with the distribution. 12 // 13 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS'' 14 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 15 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 16 // ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR 17 // ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 18 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 19 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 20 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 21 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 22 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 23 24 #include <octave/oct.h> 25 26 DEFUN_DLD(hough_line, args, , "\ 27 -*- texinfo -*-\n\ 28 @deftypefn {Loadable Function} {[@var{H}, @var{R}] =} hough_line(@var{I}, @var{angles})\n\ 29 Calculate the straight line Hough transform of a binary image @var{I}.\n\ 30 \n\ 31 The angles are given in radians and default to -pi/2:pi/2.\n\ 32 \n\ 33 @var{H} is the resulting Hough transform, and @var{R} is the radial distances.\n\ 34 \n\ 35 The algorithm is described in\n\ 36 Digital Image Processing by Gonzales & Woods (2nd ed., p. 587)\n\ 37 \n\ 38 For a Matlab compatible Hough transform see hough.m\n\ 39 @end deftypefn\n\ 40 ") 41 { 42 octave_value_list retval; 43 const int nargin = args.length (); 44 const bool DEF_THETA = (nargin == 1); 45 46 if (nargin < 1 || nargin > 2) 47 print_usage (); 48 49 const Matrix I = args (0).matrix_value (); 50 const ColumnVector thetas = (DEF_THETA) ? ColumnVector (Range (-M_PI/2.0, M_PI/2.0, M_PI/180.0).matrix_value ()) 51 : ColumnVector (args (1).vector_value ()); 52 53 const int r = I.rows (); 54 const int c = I.columns (); 55 const int thetas_length = thetas.numel (); 56 57 const double diag_length = sqrt ((r-1)*(r-1) + (c-1)*(c-1)); 58 const int nr_bins = 2 * (int)ceil (diag_length) + 1; 59 RowVector bins = RowVector (Range(1, nr_bins).matrix_value ()) - ceil (nr_bins/2.0); 60 const int bins_length = bins.numel (); 61 62 Matrix J (bins_length, thetas_length, 0.0); 63 64 for (int i = 0; i < thetas_length; i++) 65 { 66 const double theta = thetas (i); 67 68 const double cT = cos (theta); 69 const double sT = sin (theta); 70 for (int x = 0; x < r; x++) 71 { 72 for (int y = 0; y < c; y++) 73 { 74 if (I(x, y) == 1) 75 { 76 const int rho = (int)floor (cT*x + sT*y + 0.5); 77 const int bin = (int)(rho - bins (0)); 78 if ((bin > 0) && (bin < bins_length)) 79 J (bin, i)++; 80 } 81 } 82 } 83 } 84 85 retval.append (J); 86 retval.append (bins); 87 return retval; 88 } 89 90 /* 91 %!test 92 %! I = zeros(100, 100); 93 %! I(1,1) = 1; I(100,100) = 1; I(1,100) = 1; I(100, 1) = 1; I(50,50) = 1; 94 %! [J, R] = houghtf(I); J = J / max(J(:)); 95 %! assert(size(J) == [length(R) 181]); 96 %! 97 98 %!demo 99 %! I = zeros(100, 150); 100 %! I(30,:) = 1; I(:, 65) = 1; I(35:45, 35:50) = 1; 101 %! for i = 1:90, I(i,i) = 1;endfor 102 %! I = imnoise(I, 'salt & pepper'); 103 %! imshow(I); 104 %! J = houghtf(I); J = J / max(J(:)); 105 %! imshow(J, bone(128), 'truesize'); 106 107 */ 108