function [M,Mh] = RadialLines(N,L)
%RADIAL LINES  Generates a radial-line mask
%
%   [M,MH] = RADIAL LINES(N,L) creates a mask of size N by N
%   containing L radial lines. The center of the mask is shifted to
%   coincide with the ordering of the coefficients generated by the
%   Fourier transform.
%
%   This function is based on LineMask.m by Justin Romberg, which
%   is part of the l1-magic code.
%
%   References
%
%   [l1magic] E.J. Candes and J. Romberg, l1-Magic, 2005
%     http://www.l1-magic.org/ 

%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
%   http://www.cs.ubc.ca/labs/scl/sparco
%   $Id: RadialLines.m 1040 2008-06-26 20:29:02Z ewout78 $


% Generate spectral lines
M = sparse(N,N); k = floor((N-1)/2);
for i=1:L
   a  = pi*(i-1) / L;
   dx = cos(a);
   dy = sin(a);

   if (abs(dx) > abs(dy))
      % Full width
      x = -k:k;
      y = round(x*dy/dx);
   else
      % Full height
      y = -k:k;
      x = round(y*dx/dy);
   end

   x = x + floor(N/2);
   y = y + floor(N/2) + 1;
   M(y+(x * N)) = 1;
end

% Half of the radial lines, excluding origin
k = floor(N/2);
Mh            = M;
Mh(k+2:N,:)   = 0;
Mh(k+1,k+1:N) = 0;

% Shift mask to correspond to Fourier coefficients
M  = ifftshift(M);
Mh = ifftshift(Mh);
