From: Subject: Date: Tue, 11 Dec 2007 15:46:24 -0800 MIME-Version: 1.0 Content-Type: text/html; charset="Windows-1252" Content-Transfer-Encoding: quoted-printable Content-Location: http://www.stanford.edu/class/ee263/matlab/line_pixel_length.m X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2900.3198

function L=3Dline_pixel_length(d,theta,n)

% image reconstruction from line measurements (tomography via =
least-squares)
%
% given a grid of n by n square pixels, and a line over that grid,
% this function computes the length of line that goes over each pixel
%
% INPUTS
% d:     displacement of line,
%        distance of line from center of image,
%        measured in pixel lengths (and orthogonally to line)
% theta: angle of line,
%        measured in radians from x-axis
% n:     image size is n by n
%
% OUTPUTS
% L:     matrix of size n by n (same as image),
%        length of the line over each pixel (most entries are zero)
%
% expects an angle theta in [0,pi]
% (but will work at least for angles in [-pi/4,pi])


  % for angle in [pi/4,3*pi/4],
  % flip along diagonal (transpose) and call recursively
if theta>pi/4 & theta<3/4*pi
  L=3Dline_pixel_length(d,pi/2-theta,n)';
  return
end

  % for angle in [3*pi/4,pi],
  % redefine line to go in opposite direction
if theta>pi/2
  d=3D-d;
  theta=3Dtheta-pi;
end

  % for angle in [-pi/4,0],
  % flip along x-axis (up/down) and call recursively
if theta<0
  L=3Dflipud(line_pixel_length(-d,-theta,n));
  return
end

if theta>pi/2 | theta<0
  disp('invalid angle')
  return
end

L=3Dzeros(n,n);

ct=3Dcos(theta);
st=3Dsin(theta);

x0=3Dn/2-d*st;
y0=3Dn/2+d*ct;

y=3Dy0-x0*st/ct;
jy=3Dceil(y);
dy=3Drem(y+n,1);

for jx=3D1:n
  dynext=3Ddy+st/ct;
  if dynext<1
    if jy>=3D1 & jy<=3Dn, L(n+1-jy,jx)=3D1/ct; end
    dy=3Ddynext;
  else
    if jy>=3D1 & jy<=3Dn, L(n+1-jy,jx)=3D(1-dy)/st; end
    if jy+1>=3D1 & jy+1<=3Dn, =
L(n+1-(jy+1),jx)=3D(dynext-1)/st; end
    dy=3Ddynext-1;
    jy=3Djy+1;
  end
end