I don't have the rep to comment on Dima's answer, but the code in the question does convert the $\theta$ to radians already.

This is my implementation of the routine for detecting lines:

[ theta,rho,height,width,Accumulator] = SHTanalyse( m,delta_theta,delta_rho )
%% Get dimensions of input image
[height, width] = size(m);

% Create array for theta 
theta = -90:delta_theta:90;
costheta = cosd(theta); % Create lookup tables to optimize computation of cos theta and sin theta
sintheta = sind(theta);

% Find bounds on rho
rhomax = round(sqrt(height^2 + width^2));

% Create array for rho
rho = -rhomax:delta_rho:rhomax;

% Get edgel coordinates
[yEdges,xEdges] =  find(m);

% Initialize accumulator
Accumulator = zeros(numel(rho),numel(theta));

% Voting process
for k = 1:numel(yEdges) % Loop through edge pixels
    for t = 1:numel(theta)
        % Compute the corresponding rho
        rhvalue = round(xEdges(k).*costheta(t) + yEdges(k).*sintheta(t));
        % Increment accumulator
        Accumulator(rhvalue + round(0.5*numel(rho)) + 1,theta(t) + 90 + 1) = Accumulator(rhvalue + round(0.5*numel(rho)) + 1,theta(t) + 90 + 1) + 1;
    end

You can compare its performance with the MATLAB's own Hough function to figure out the discrepancy.

I don't have the rep to comment on Dima's answer, but the code in the question does convert the $\theta$ to radians already.

This is my implementation of the routine for detecting lines:

[ theta,rho,height,width,Accumulator] = SHTanalyse( m,delta_theta,delta_rho )
%% Get dimensions of input image
[height, width] = size(m);

% Create array for theta 
theta = -90:delta_theta:90;
costheta = cosd(theta); % Create lookup tables to optimize computation of cos theta and sin theta
sintheta = sind(theta);

% Find bounds on rho
rhomax = round(sqrt(height^2 + width^2));

% Create array for rho
rho = -rhomax:delta_rho:rhomax;

% Get edgel coordinates
[yEdges,xEdges] =  find(m);

% Initialize the accumulator
Accumulator = zeros(numel(rho),numel(theta));

% Voting process
for k = 1:numel(yEdges) % Loop through edge pixels
    for t = 1:numel(theta)
        % Compute the corresponding rho
        rhvalue = round(xEdges(k).*costheta(t) + yEdges(k).*sintheta(t));
        % Increment accumulator
        Accumulator(rhvalue + round(0.5*numel(rho)) + 1,theta(t) + 90 + 1) = Accumulator(rhvalue + round(0.5*numel(rho)) + 1,theta(t) + 90 + 1) + 1;
    end
end    

You can compare its performance with the MATLAB's own Hough function to figure out the discrepancy.