1. 程式人生 > >《數字影象處理原理與實踐(MATLAB版)》一書之程式碼Part3

《數字影象處理原理與實踐(MATLAB版)》一書之程式碼Part3

本文系《數字影象處理原理與實踐(MATLAB版)》一書之程式碼系列的Part3,輯錄該書第135至第184頁之程式碼,供有需要讀者下載研究使用。程式碼執行結果請參見原書配圖。

-------------------------------------------

P139

original = imread('snowflakes.png');
figure, imshow(original);
se = strel('disk',5);
afterOpening = imopen(original,se);
figure, imshow(afterOpening,[]);

P140

originalBW = imread('circles.png');
imshow(originalBW);
se = strel('disk',10);
closeBW = imclose(originalBW,se);
figure, imshow(closeBW)

P144

bw = imread('bw.bmp');
shape1 = [0 0 0 0 1
            0 0 0 1 1
            0 0 1 1 1
            0 1 1 1 1
            1 1 1 1 1];
shape2 = [1 1 0 0 0
            1 0 0 0 0
            0 0 0 0 0
            0 0 0 0 0
            0 0 0 0 0];
bw2 = bwhitmiss(bw, shape1, shape2);
imshow(bw2)

P146-1

I = imread('letter2.jpg');
I = im2bw(I);
I1 = bwmorph(I, 'thin',inf);
figure(1), imshow(I1);

P146-2

I = imread('letter2.jpg');
I = im2bw(I);
I2 = bwmorph(I, 'skel',inf);
figure(2), imshow(I2);

P153

I = imread('lena.jpg');
I = rgb2gray(I);
BW1 = edge(I, 'roberts');
BW2 = edge(I, 'sobel');
BW3 = edge(I, 'prewitt');
figure
subplot(2,2,1),imshow(I),title('original')
subplot(2,2,2),imshow(BW1),title('roberts')
subplot(2,2,3),imshow(BW2),title('sobel')
subplot(2,2,4),imshow(BW3),title('prewitt')

P157

I = imread('einstein.bmp');
I = rgb2gray(I);
N = [1, 2, 1
     0, 0, 0
     -1,-2,-1];
edge_n = imfilter(I,N,'symmetric','conv');
imwrite(edge_n, 'edge_n.jpg');

P160

I = rgb2gray(imread('lena.jpg'));
M = [1,1,1
        1,-8,1
        1,1,1];
img=imfilter(I,M);
[x,y]=size(I);
img2 = img;
for i = 2:x-1
        for j = 2:y-1
            a = [img(i,j+1),img(i,j-1),img(i+1,j+1),img(i+1,j-1), ...
                 img(i-1,j+1),img(i-1,j-1),img(i+1,j),img(i-1,j)];
            if ( (max(a)-min(a))>64 && max(a)>img(i,j) && min(a)<img(i,j))
                img2(i,j)=255;
            else
                img2(i,j)=0;
            end
        end
end

P165

I = imread('lena.jpg');
IMG = rgb2gray(I);
Edge_LoG = edge(IMG, 'log');
imshow(Edge_LoG);
figure
subplot(1,2,1), imshow(IMG);
subplot(1,2,2), imshow(Edge_LoG);

P167
I = double(rgb2gray(imread('lena.jpg')));
figure, imshow(uint8(I))
DoG=fspecial('gaussian',5,0.8)-fspecial('gaussian',5,0.6);
ImageDoG=imfilter(I,DoG,'symmetric','conv');
figure, imshow(ImageDoG)
% threshold = 2
proc_Img1 = ImageDoG;
proc_Img1(find(proc_Img1 < 2))=0;
figure, imshow(proc_Img1)
% threshold = 3
proc_Img2 = ImageDoG;
proc_Img2(find(proc_Img2 < 3))=0;
figure, imshow(proc_Img2)

P172

img = edge(I, 'canny',[0.032,0.08], 3);

P183

RGB= imread('building.jpg');
I = rgb2gray(RGB);
BW = edge(I, 'canny');

[H, T, R]=hough(BW, 'RhoResolution',0.5,'ThetaResolution',0.5);
figure, imshow(imadjust(mat2gray(H)), 'XData', T, ...
'YData', R, 'InitialMagnification', 'fit');
xlabel('\theta'), ylabel('\rho');
axis on; axis normal; hold on;
colormap(hot);
peaks = houghpeaks(H, 15);
figure, imshow(BW);
hold on;
lines = houghlines(BW, T, R, peaks, 'FillGap',25, 'MinLength',15);

max_len = 0;
for k=1:length(lines)
xy = [lines(k).point1; lines(k).point2];
   plot(xy(:,1),xy(:,2),'LineWidth',3,'Color','b');
   plot(xy(1,1),xy(1,2),'x','LineWidth',3,'Color','yellow');
   plot(xy(2,1),xy(2,2),'x','LineWidth',3,'Color','red');

   len = norm(lines(k).point1 - lines(k).point2);
   if ( len > max_len)
      max_len = len;
      xy_long = xy;
   end
end

-------------------------------------------

原書184頁,我曾提到,MATLAB中影象處理工具箱提供的Hough直線檢測的效果不是特別好,並推薦了一個我個人認為實現得更好的版本。但由於本人並非該程式碼之原作者,遂未將其收錄於書中。近來有讀者留言,希望可以獲得這部分程式碼,所以我特將其釋出到此部落格上,供有需要的讀者參考學習。

程式碼原作者為Tao Peng,請尊重原作者權利。

%  Author:  Tao Peng  
%           Department of Mechanical Engineering  
%           University of Maryland, College Park, Maryland 20742, USA  
%           [email protected]  
%  Version: alpha       Revision: Dec. 02, 2005

使用時可參考下面這兩個例子:
% EXAMPLE #1:  
rawimg = imread('TestHT_Chkbd1.bmp');  
fltr4img = [1 2 3 2 1; 2 3 4 3 2; 3 4 6 4 3; 2 3 4 3 2; 1 2 3 2 1];  
fltr4img = fltr4img / sum(fltr4img(:));  
imgfltrd = filter2( fltr4img , rawimg );  
  
  
[accum, axis_rho, axis_theta, lineprm, lineseg] = Hough_Grd(imgfltrd, 6, 0.02);  
  
figure(1); imagesc(axis_theta*(180/pi), axis_rho, accum); axis xy;  
xlabel('Theta (degree)'); ylabel('Pho (pixels)');  
title('Accumulation Array from Hough Transform');  
figure(2); imagesc(rawimg); colormap('gray'); axis image;  
DrawLines_2Ends(lineseg);  
title('Raw Image with Line Segments Detected'); 

函式Houg_Grd()程式碼實現:

    function [accum, axis_rho, axis_theta, varargout] = ...  
        Hough_Grd(img, varargin)  
    %Detect lines in a grayscale image  
    %  
    %  [accum, axis_rho, axis_theta, lineprm, lineseg, dbg_label] =  
    %      Hough_Grd(img, grdthres, detsens)  
    %  Hough transform for line detection based on image's gradient field.  
    %  NOTE:    Operates on grayscale images, NOT B/W bitmaps.  
    %           NO loops involved in the implementation of Hough transform,  
    %               which makes the operation fast.  
    %           Able to detect the two ends of line segments.  
    %  
    %  INPUT: (img, grdthres, detsens)  
    %  img:         A 2-D grayscale image (NO B/W bitmap)  
    %  grdthres:    (Optional, default is 8, must be non-negative)  
    %               The algorithm is based on the gradient field of the  
    %               input image. A thresholding on the gradient magnitude  
    %               is performed before the voting process of the Hough  
    %               transform to remove the 'uniform intensity' (sort-of)  
    %               image background from the voting process. In other words,  
    %               pixels with gradient magnitudes smaller than 'grdthres'  
    %               are considered not belong to any line.  
    %  detsens:     (Optional, default is 0.08)  
    %               A value that controls the sensitivity of line detection.  
    %               The range of the value is from 0 to 1, open ends.  
    %               Although the physical meaning of this parameter is not  
    %               obvious, it controls the algorithm in a fuzzy manner.  
    %               The SMALLER the value is, the MORE features in the image  
    %               will be considered as lines.  
    %  
    %  OUTPUT: [accum, axis_rho, axis_theta, lineprm, lineseg, dbg_label]  
    %  accum:       The result accumulation array from the Hough transform.  
    %               The horizontal dimension is 'theta' and the vertical  
    %               dimension is 'rho'.  
    %  axis_rho:    Vector that contains the 'rho' values for the rows in  
    %               the accumulation array.  
    %  axis_theta:  Vector that contains the 'theta' values for the columns  
    %               in the accumulation array.  
    %  lineprm:     (Optional)  
    %               Parameters (rho, theta) of the lines detected. Is a N-by-2  
    %               matrix with each row contains the parameters (rho, theta)   
    %               for a line. The definitions of 'rho' and 'theta' are as  
    %               the following:  
    %               'rho' is the perpendicular distance from the line to the  
    %               origin of the image. 'rho' can be negative since 'theta'  
    %               is constrained to [0, pi]. The unit of 'rho' is in pixels.  
    %               'theta' is the sweep angle from axis X (i.e. horizontal  
    %               axis of the image) to the direction that is perpendicular  
    %               to the line. The range of 'theta' is [0, pi].  
    %  lineseg:     (Optional)  
    %               Parameters (x1, x2, y1, y2) of line segments detected.  
    %               Lines given by (rho, theta) are infinite lines. This  
    %               function tries to detect line segments in the raw image  
    %               and output them as 'lineseg', which is a Ns-by-4 matrix  
    %               with each row contains the parameters (x1, x2, y1, y2)  
    %               that defines the two ends of a line segment.  
    %  dbg_label:   (Optional, for debug purpose)  
    %               Labelled segmentation map from the search of peaks in  
    %               the accumulation array  
    %  
    %  
    %  BUG REPORT:  
    %  This is an alpha version. Please send your bug reports, comments and  
    %  suggestions to [email protected] . Thanks.  
      
    %  Reserved for extension  
    %  accumres:    (Optional, default is [4, 1], minimum is [2, 0.5])  
    %               The desired resolutions in 'rho' and 'theta'. Is a  
    %               2-element vector in following format:  
    %               [resolution in 'rho' (pixels),  
    %               resolution in 'theta' (degrees)]  
      
      
      
      
    %%%%%%%% Arguments and parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
      
    % Validation of arguments  
    if ndims(img) ~= 2 || ~isnumeric(img),  
        error('Hough_Grd: ''img'' has to be 2 dimensional');  
    end  
    if ~all(size(img) >= 16),  
        error('Hough_Grd: ''img'' has to be larger than 16-by-16');  
    end  
      
      
    % Parameters (default values)  
    prm_grdthres = 8;  
    prm_accumres = [4, 1];  
    prm_detsens = 0.08;  
      
      
    func_compu_lineprm = true;  
    prm_fltraccum = true;  
      
      
    % Validation of arguments  
    vap_grdthres = 1;  
    if nargin > vap_grdthres,  
        if isnumeric(varargin{vap_grdthres}) && ...  
                varargin{vap_grdthres}(1) >= 0,  
            prm_grdthres = varargin{vap_grdthres}(1);  
        else  
            error(['Hough_Grd: ''grdthres'' has to be ', ...  
                'a non-negative number']);  
        end  
    end  
    %{  
    vap_accumres = 3;  
    if nargin > vap_accumres,  
        if numel(varargin{vap_accumres}) == 2 && ...  
                isnumeric(varargin{vap_accumres}) && ...  
                ( varargin{vap_accumres}(1) >= 2 && ...  
                varargin{vap_accumres}(2) >= 0.5 ),  
            prm_accumres = varargin{vap_accumres};  
        else  
            error(['Hough_Grd: ''accumres'' has to be a two-element ', ...  
                'vector and no smaller than [2, 0.5]']);  
        end  
    end  
    %}  
    vap_detsens = 2;  
    if nargin > vap_detsens,  
        if isnumeric(varargin{vap_detsens}) && ...  
                varargin{vap_detsens}(1) > 0 && varargin{vap_detsens}(1) < 1,  
            prm_detsens = varargin{vap_detsens};  
        else  
            error('Hough_Grd: ''detsens'' has to be in the range (0, 1)');  
        end  
    end  
      
      
    func_compu_lineprm = ( nargout > 3 );  
      
      
    % Size of the accumulation array  
    imgsize = size(img);  
    coef_rhorng = [ -imgsize(2), sqrt(sum(imgsize.^2)) ];  
    coef_thetarng = [-pi/18, pi+pi/18];  
      
      
    prm_accumsize = [ ...  
        round( (coef_rhorng(2)-coef_rhorng(1)) * (2/prm_accumres(1)) ) , ...  
        round( (coef_thetarng(2)-coef_thetarng(1)) * (180/pi) * ...  
        (2/prm_accumres(2)) ) ];  
      
      
    % Default filter for the accumulation array  
    prm_acmfltr_R = 4;  
    prm_acmfltr_w = [1 2 4 8];  
      
      
    fltr4accum = ones(2 * prm_acmfltr_R - 1) * prm_acmfltr_w(1);  
    for k = 2 : prm_acmfltr_R,  
        fltr4accum(k:(2*prm_acmfltr_R-k), k:(2*prm_acmfltr_R-k)) = ...  
            prm_acmfltr_w(k);  
    end  
    fltr4accum = fltr4accum / sum(fltr4accum(:));  
      
      
    % Parameters for the algorithm using repeated  
    % thresholding and segmentation  
    prm_lp_dthres = min([ 0.1, (0.1+prm_detsens)/2 ]);  
    prm_lp_maxthres = 0.8;  
      
      
    % Parameters for the algorithm using local maximum filter  
    prm_useaoi = false;  
    prm_aoiminsize = [8, 8];  
      
      
    prm_fltrLM_R = 4;  
    prm_fltrLM_s = 1.3;  
    prm_fltrLM_r = ceil( prm_fltrLM_R * 0.6 );  
    prm_fltrLM_npix = 6;  
      
      
    % Reserved parameters  
    dbg_on = false;      % debug information  
    dbg_bfigno = 5;  
    if nargout > 5,  dbg_on = true;  end  
      
      
      
      
    %%%%%%%% Building accumulation array %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
      
      
    % Compute the gradient and the magnitude of gradient  
    img = double(img);  
    [grdx, grdy] = gradient(img);  
    grdmag = sqrt(grdx.^2 + grdy.^2);  
      
      
    % Clear the margins of the gradient field  
    prm_grdfldmgn = 4;  
    grdmag([1:prm_grdfldmgn, (end-prm_grdfldmgn+1):end], :) = 0;  
    grdmag(:, [1:prm_grdfldmgn, (end-prm_grdfldmgn+1):end]) = 0;  
      
      
    % Get the linear indices, as well as the subscripts, of the pixels  
    % whose gradient magnitudes are larger than the given threshold  
    grdmasklin = find(grdmag > prm_grdthres);  
    [grdmask_IdxI, grdmask_IdxJ] = ind2sub(imgsize, grdmasklin);  
      
      
    % Compute (the line parameter) 'theta' for all voting pixels  
    grdphs_vot = atan2( grdy(grdmasklin), grdx(grdmasklin) );  
      
      
    % -- Regulate 'grdphs_vot' from [-pi, pi] to [0, pi]  
    grdphs_vot = grdphs_vot + pi * (grdphs_vot < 0);  
      
      
    % Compute the 'theta'-subscript (to the accumulation array)  
    % for all voting pixels  
    coef_subtheta = prm_accumsize(2) / ...  
        (coef_thetarng(2) - coef_thetarng(1)) * (1 - 1e-6);  
    sub_theta = ceil( (grdphs_vot - coef_thetarng(1)) * coef_subtheta );  
      
      
    % 'theta' vector for the accumulation array  
    axis_theta = (coef_thetarng(1) + 0.5 / coef_subtheta) : ...  
        (1 / coef_subtheta) : coef_thetarng(2);  
      
      
    % Compute the 'rho' values for all voting pixels  
    % rho = (J - 0.5) * cos(theta) + (I - 0.5) * sin(theta)  
    rho_vot = (grdmask_IdxJ - 0.5) .* cos(grdphs_vot) + ...  
        (grdmask_IdxI - 0.5) .* sin(grdphs_vot);  
      
      
    % Compute the 'rho'-subscript (to the accumulation array)  
    % for all voting pixels  
    coef_subrho = prm_accumsize(1) / ...  
        (coef_rhorng(2) - coef_rhorng(1)) * (1 - 1e-6);  
    sub_rho = ceil( (rho_vot - coef_rhorng(1)) * coef_subrho );  
      
      
    % 'rho' vector for the accumulation array  
    axis_rho = (coef_rhorng(1) + 0.5 / coef_subrho) : ...  
        (1 / coef_subrho) : coef_rhorng(2);  
      
      
    % Build the accumulation array, using gradient magnitude as weight  
    accum = accumarray( sub2ind(prm_accumsize, sub_rho, sub_theta), ...  
        grdmag(grdmasklin) );  
    accum = [ accum ; ...  
        zeros(prm_accumsize(1) * prm_accumsize(2) - numel(accum), 1) ];  
    accum = reshape( accum, prm_accumsize );  
      
      
      
      
    %%%%%%%% Locating peaks in the accumulation array %%%%%%%%%%%%%%%%%%%  
      
      
    % Stop if no need to estimate the parameters of the lines  
    if ~func_compu_lineprm,  
        return;  
    end  
      
      
    % Smooth the accumulation array  
    if prm_fltraccum,  
        accum = filter2( fltr4accum, accum );  
    end  
      
      
    % Find the maximum value in the accumulation array  
    accum_max = max(accum(:));  
      
      
      
      
    %------- Algorithm 1: Local maximum filter (begin) -------------  
    % Build the local maximum filter  
    fltr4LM = zeros(2 * prm_fltrLM_R + 1);  
      
      
    [mesh4fLM_x, mesh4fLM_y] = meshgrid(-prm_fltrLM_R : prm_fltrLM_R);  
    mesh4fLM_r = sqrt( mesh4fLM_x.^2 + mesh4fLM_y.^2 );  
    fltr4LM_mask = ...  
        ( mesh4fLM_r > prm_fltrLM_r & mesh4fLM_r <= prm_fltrLM_R );  
    fltr4LM = fltr4LM - ...  
        fltr4LM_mask * (prm_fltrLM_s / sum(fltr4LM_mask(:)));  
      
      
    if prm_fltrLM_R >= 4,  
        fltr4LM_mask = ( mesh4fLM_r < (prm_fltrLM_r - 1) );  
    else  
        fltr4LM_mask = ( mesh4fLM_r < prm_fltrLM_r );  
    end  
    fltr4LM = fltr4LM + fltr4LM_mask / sum(fltr4LM_mask(:));  
      
      
    % Select a number of Areas-Of-Interest from the accumulation array  
    if prm_useaoi,  
        % Thresholding and segmentation  
        accummask = ( accum > (accum_max * prm_detsens) );  
        [accumlabel, accum_nRgn] = bwlabel( accummask, 8 );  
      
      
        % Select AOIs from segmented regions  
        accumAOI = ones(0,4);  
        for k = 1 : accum_nRgn,  
            accumrgn_lin = find( accumlabel == k );  
            [accumrgn_IdxI, accumrgn_IdxJ] = ...  
                ind2sub( size(accumlabel), accumrgn_lin );  
            rgn_top = min( accumrgn_IdxI );  
            rgn_bottom = max( accumrgn_IdxI );  
            rgn_left = min( accumrgn_IdxJ );  
            rgn_right = max( accumrgn_IdxJ );          
            % The AOIs selected must satisfy a minimum size  
            if ( (rgn_bottom - rgn_top + 1) >= prm_aoiminsize(1) || ...  
                    (rgn_right - rgn_left + 1) >= prm_aoiminsize(2) ),  
                accumAOI = [ accumAOI; ...  
                    max([ 1, rgn_top - prm_fltrLM_R ]), ...  
                    min([ size(accum,1), rgn_bottom + prm_fltrLM_R ]), ...  
                    max([ 1, rgn_left - prm_fltrLM_R ]), ...  
                    min([ size(accum,2), rgn_right + prm_fltrLM_R ]) ];  
            end  
        end  
    else  
        % Whole accumulation array as the one AOI  
        accumAOI = [1, size(accum,1), 1, size(accum,2)];  
    end  
      
      
    % **** Debug code (begin)  
    if dbg_on && prm_useaoi,  
        dbg_accumLM = zeros(size(accum));  
    end  
    % **** Debug code (end)  
      
      
    % For each of the AOIs selected, locate the local maxima  
    lineprm = zeros(0,2);  
    for k = 1 : size(accumAOI, 1),  
        aoi = accumAOI(k,:);    % just for referencing convenience  
      
      
        % Apply the local maxima filter  
        candLM = conv2( accum(aoi(1):aoi(2), aoi(3):aoi(4)) , ...  
            fltr4LM , 'same' );  
      
      
        % Thresholding of 'candLM' & 'accum'  
        if prm_useaoi,  
            candLM_mask = ( candLM > (max(candLM(:))*prm_detsens) & ...  
                accummask(aoi(1):aoi(2),aoi(3):aoi(4)) );  
        else  
            candLM_mask = ( candLM > (max(candLM(:))*prm_detsens) & ...  
                accum(aoi(1):aoi(2),aoi(3):aoi(4)) > (accum_max*prm_detsens) );  
        end  
      
      
        % Clear the margins of 'candLM_mask'  
        candLM_mask([1:prm_fltrLM_R, (end-prm_fltrLM_R+1):end], :) = 0;  
        candLM_mask(:, [1:prm_fltrLM_R, (end-prm_fltrLM_R+1):end]) = 0;  
      
      
        % **** Debug code (begin)  
        if dbg_on && prm_useaoi,  
            dbg_accumLM(aoi(1):aoi(2), aoi(3):aoi(4)) = ...  
                dbg_accumLM(aoi(1):aoi(2), aoi(3):aoi(4)) + candLM;  
        end  
        % **** Debug code (end)  
      
      
        % Group the local maxima candidates by adjacency, compute the  
        % centroid position for each group and take that as the center  
        % of one circle detected  
        [candLM_label, candLM_nRgn] = bwlabel( candLM_mask, 8 );  
      
      
        for ilabel = 1 : candLM_nRgn,  
            % Indices (to current AOI) of the pixels in the group  
            candgrp_masklin = find( candLM_label == ilabel );  
            [candgrp_IdxI, candgrp_IdxJ] = ...  
                ind2sub( size(candLM_label) , candgrp_masklin );  
      
      
            % Indices (to 'accum') of the pixels in the group  
            candgrp_IdxI = candgrp_IdxI + ( aoi(1) - 1 );  
            candgrp_IdxJ = candgrp_IdxJ + ( aoi(3) - 1 );  
            candgrp_idx2acm = ...  
                sub2ind( size(accum) , candgrp_IdxI , candgrp_IdxJ );  
      
      
            % Minimum number of qulified pixels in the group  
            if sum(numel(candgrp_masklin)) < prm_fltrLM_npix,  
                continue;  
            end  
      
      
            % Compute the centroid position  
            candgrp_acmsum = sum( accum(candgrp_idx2acm) );  
            cc_rho = sum( candgrp_IdxI .* accum(candgrp_idx2acm) ) / ...  
                candgrp_acmsum;  
            cc_theta = sum( candgrp_IdxJ .* accum(candgrp_idx2acm) ) / ...  
                candgrp_acmsum;  
            lineprm = [lineprm; cc_rho, cc_theta];  
        end  
    end  
      
      
    % **** Debug code (begin)  
    if dbg_on,  
        figure(dbg_bfigno);  
        if prm_useaoi,  
            imagesc(dbg_accumLM); axis xy;  
        else  
            imagesc(candLM); axis xy;  
        end  
    end  
    % **** Debug code (end)  
    %------- Algorithm 1: Local maximum filter (end) ---------------  
      
      
    %------- Algo 2: Repeated thresholding & segmentation (begin) --  
    %{  
    % Locate the peaks (in pixel coordinates) in the accumulation array  
    if dbg_on,  
        [lineprm, dbg_label] = RecursSegment( accum , ...  
            accum_max * [prm_detsens, prm_lp_dthres, prm_lp_maxthres] );  
    else  
        lineprm = RecursSegment( accum , ...  
            accum_max * [prm_detsens, prm_lp_dthres, prm_lp_maxthres] );  
    end  
      
      
    % **** Debug code (begin)  
    if dbg_on,  
        figure(dbg_bfigno);  
        imagesc(dbg_label); axis xy; axis equal;  
        hold on;  
        plot(lineprm(:,2), lineprm(:,1), ...  
            'w+', 'LineWidth', 2, 'MarkerSize', 6);  
        hold off;  
    end  
    % **** Debug code (end)  
    %}  
    %------- Algo 2: Repeated thresholding & segmentation (end) ----  
      
      
      
      
    % Convert 'lineprm' from pixel coordinates to (rho, theta)  
    lineprm = [ (lineprm(:,1) - 0.5)/coef_subrho + coef_rhorng(1), ...  
        (lineprm(:,2) - 0.5)/coef_subtheta + coef_thetarng(1) ];  
      
      
    % Output 'lineprm'  
    varargout{1} = lineprm;  
    if nargout <= 4,  
        return;  
    end  
      
      
      
      
    %%%%%%%% Locating the line segments in the raw image %%%%%%%%%%%%%%%%  
      
      
    % Parameters for locating the line segments  
    prm_ls_tTol = [-pi/7.5, pi/7.5];  
    prm_ls_pTol = [-3.5, 3.5];  
    prm_ls_minlen = 10;  
      
      
    % Locate the two ends for all lines detected  
    rgnmask = logical(zeros(imgsize));  
    lineseg = zeros(0, 4);  
      
      
    for k = 1 : size(lineprm, 1),  
        % Compute the 'rho' values for all pixels voted,  
        % assuming they contribute to the current line  
        rho2_vot = (grdmask_IdxJ - 0.5) * cos( lineprm(k,2) ) + ...  
            (grdmask_IdxI - 0.5) .* sin( lineprm(k,2) );  
      
      
        % Find the pixels that belong to the line  
        PixOnLn_votmask = ( grdphs_vot > (lineprm(k,2) + prm_ls_tTol(1)) & ...  
            grdphs_vot < (lineprm(k,2) + prm_ls_tTol(2)) & ...  
            rho2_vot > (lineprm(k,1) + prm_ls_pTol(1)) & ...  
            rho2_vot < (lineprm(k,1) + prm_ls_pTol(2)) );  
      
      
        PixOnLn_lin = grdmasklin( find(PixOnLn_votmask) );  
        if isempty(PixOnLn_lin),  
            continue;  
        end  
        [PixOnLn_IdxI, PixOnLn_IdxJ] = ind2sub(imgsize, PixOnLn_lin);  
      
      
        % Find the axis-aligned bounding box for these pixels  
        bndbox = [ min(PixOnLn_IdxI), max(PixOnLn_IdxI), ...  
            min(PixOnLn_IdxJ), max(PixOnLn_IdxJ) ];  
      
      
        % Seperate the line segments by adjacencies  
        rgnmask( bndbox(1):bndbox(2), bndbox(3):bndbox(4) ) = 0;  
        rgnmask( PixOnLn_lin ) = 1;  
        [rgnlabel, nrgn] = bwlabel( ...  
            rgnmask(bndbox(1):bndbox(2), bndbox(3):bndbox(4)), 8 );  
      
      
        for k_rgn = 1 : nrgn,  
            % Get the linear indices of pixels that belong to one segment  
            seg_lin = find( rgnlabel == k_rgn );  
            [seg_IdxI, seg_IdxJ] = ind2sub( size(rgnlabel), seg_lin );  
      
      
            % Ignore if the line segment is too short  
            segspan = [ min(seg_IdxI) + bndbox(1) - 1, ...  
                max(seg_IdxI) + bndbox(1) - 1, ...  
                min(seg_IdxJ) + bndbox(3) - 1, ...  
                max(seg_IdxJ) + bndbox(3) - 1 ];  
            if (segspan(2) - segspan(1) + 1) < prm_ls_minlen && ...  
                    (segspan(4) - segspan(3) + 1) < prm_ls_minlen,  
                continue;  
            end  
      
      
            % Get the [x1 x2 y1 y2] structure for the line SEGMENT  
            if lineprm(k,2) > pi/4 && lineprm(k,2) < 3*pi/4,  
                % ls_base_x = 0;  
                ls_base_y = lineprm(k,1) / sin(lineprm(k,2));  
                lineseg = [ lineseg; segspan(3) - 0.5, segspan(4) - 0.5, ...  
                    ls_base_y - (segspan(3) - 0.5) / tan(lineprm(k,2)), ...  
                    ls_base_y - (segspan(4) - 0.5) / tan(lineprm(k,2)) ];  
            else  
                % ls_base_y = 0;  
                ls_base_x = lineprm(k,1) / cos(lineprm(k,2));  
                lineseg = [ lineseg; ...  
                    ls_base_x - (segspan(1) - 0.5) * tan(lineprm(k,2)), ...  
                    ls_base_x - (segspan(2) - 0.5) * tan(lineprm(k,2)), ...  
                    segspan(1) - 0.5, segspan(2) - 0.5 ];  
            end  
        end  
    end  
      
      
    % Output 'lineseg'  
    varargout{2} = lineseg;  
    if nargout > 5,  
        varargout{3} = dbg_label;  
    end  
      
      
      
      
      
      
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    %%%%%%%% Recursive thresholding and segmentation ********************  
      
      
    function [lineprm, varargout] = RecursSegment(accum, struct_thres)  
    % 'struct_thres' contains [thres, deltathres, maxthres]  
    % 'lineprm' is in pixel coordinate (w.r.t. 'accum')  
      
      
    % Parameters  
    prm_as_minpixn = 3;  
    prm_as_maxsize = [12, 12];  
      
      
    % Thresholding  
    accummask = ( accum > struct_thres(1) );  
      
      
    % Segmentation and locating the centroids of individual regions  
    [accumlabel, accum_nRgn] = bwlabel( accummask, 8 );  
      
      
    % Segmentation label (for debug purpose)  
    func_seglbl = ( nargout > 1 );  
    if func_seglbl,  
        seglbl_lblshft = 4;  
        seglabel = accumlabel;  
    end  
      
      
    lineprm = zeros(0, 2);  
    for k = 1 : accum_nRgn,  
        % Linear indices of the pixels in one connected component  
        acmrgn_lin = find( accumlabel == k );  
        if numel(acmrgn_lin) < prm_as_minpixn,  
            continue;  
        end  
        % Subscripts of the pixels in one connected component  
        [acmrgn_IdxPho, acmrgn_IdxTheta] = ...  
            ind2sub( size(accumlabel), acmrgn_lin );  
      
      
        % Further segmentation if the connected region is too big, or  
        % computing the centroid of the region  
        % -- Axis-aligned bounding box for the region  
        bndbox = [ min(acmrgn_IdxPho), max(acmrgn_IdxPho), ...  
            min(acmrgn_IdxTheta), max(acmrgn_IdxTheta) ];  
      
      
        % -- Further segmentation  
        bAddCentrdOfRgn = true;  
        if ( (bndbox(2) - bndbox(1) + 1) > prm_as_maxsize(1) || ...  
                (bndbox(4) - bndbox(3) + 1) > prm_as_maxsize(2) ) && ...  
                (struct_thres(1) + struct_thres(2)) <= struct_thres(3),  
            if func_seglbl,  
                [lineprm_sub, seglabel_sub] = RecursSegment( ...  
                    accum(bndbox(1):bndbox(2), bndbox(3):bndbox(4)), ...  
                    [struct_thres(1) + struct_thres(2), struct_thres(2:3)] );  
                seglabel(bndbox(1):bndbox(2), bndbox(3):bndbox(4)) = ...  
                    seglabel(bndbox(1):bndbox(2), bndbox(3):bndbox(4)) + ...  
                    seglabel_sub + (seglabel_sub > 0) * seglbl_lblshft;  
            else  
                lineprm_sub = RecursSegment( ...  
                    accum(bndbox(1):bndbox(2), bndbox(3):bndbox(4)), ...  
                    [struct_thres(1) + struct_thres(2), struct_thres(2:3)] );  
            end  
            if ~isempty(lineprm_sub),  
                lineprm = [lineprm; lineprm_sub(:,1) + bndbox(1) - 1, ...  
                    lineprm_sub(:,2) + bndbox(3) - 1 ];  
                bAddCentrdOfRgn = false;  
            end  
        end  
      
      
        % -- Computing the centroid of the whole region  
        if bAddCentrdOfRgn,  
            acmrgn_acmsum = sum( accum(acmrgn_lin) );  
            lp_IdxPho = sum( acmrgn_IdxPho .* accum(acmrgn_lin) ) / ...  
                acmrgn_acmsum;  
            lp_IdxTheta = sum( acmrgn_IdxTheta .* accum(acmrgn_lin) ) / ...  
                acmrgn_acmsum;  
            lineprm = [ lineprm; lp_IdxPho, lp_IdxTheta ];  
        end  
    end  
      
      
    % Output the segmentation label  
    if func_seglbl,  
        varargout{1} = seglabel;  
    end  
      
      
    function DrawLines_2Ends(lineseg, varargin)  
      
    hold on;  
    for k = 1 : size(lineseg, 1),  
        % The image origin defined in function '[...] = Hough_Grd(...)' is  
        % different from what is defined in Matlab, off by (0.5, 0.5).  
        if nargin > 1,  
            plot(lineseg(k,1:2)+0.5, lineseg(k,3:4)+0.5, varargin{1});  
        else  
            plot(lineseg(k,1:2)+0.5, lineseg(k,3:4)+0.5, 'LineWidth', 2);  
        end  
    end  
    hold off;  
      
      
    function DrawLines_Polar(imgsize, lineprm, varargin)  
      
    hold on;  
    line = zeros(2,2);  
    for k = 1 : size(lineprm, 1),  
        if lineprm(k,2) > pi/4 && lineprm(k,2) < 3*pi/4,  
            line(1,1) = 0;  
            line(1,2) = lineprm(k,1) / sin(lineprm(k,2));  
            line(2,1) = imgsize(2);  
            line(2,2) = line(1,2) - line(2,1) / tan(lineprm(k,2));  
        else  
            line(1,2) = 0;  
            line(1,1) = lineprm(k,1) / cos(lineprm(k,2));  
            line(2,2) = imgsize(1);  
            line(2,1) = line(1,1) - line(2,2) * tan(lineprm(k,2));  
        end  
        % The image origin defined in function '[...] = Hough_Grd(...)' is  
        % different from what is defined in Matlab, off by (0.5, 0.5).  
        line = line + 0.5;  
        % Draw lines using 'plot'  
        if nargin > 2,  
            plot(line(:,1), line(:,2), varargin{1});  
        else  
            plot(line(:,1), line(:,2));  
        end  
    end  
    hold off;  

(程式碼釋出未完,請待後續...)

20170601164758578