% returns a sXs grayscale scene that is:
%       (1) rotated about the center by phi radians, then
%       (2) shifted by (x0,y0)
%       (3) sheared along x by sx

function [Im] = shepplogan(s,x0,y0,phi,sx)

  % create hi res image at hisXhis resolution
  fac = 4;
  his = fac*s;

  Im = zeros(s,s);
  hiIm = zeros(his,his);

  index = 1:his;
  dsin = (index - his/2)*sin(phi);
  dcos = (index - his/2)*cos(phi);

  cx = [0 0 0.22 -0.22 0 0 0 -0.08 0 0.06]*his/2;
  cy = [0 -0.0184 0 0 0.35 0.1 -0.1 -0.605 -0.605 -0.605]*his/2;
  maj = [0.92 0.874 0.31 0.41 0.25 0.046 0.046 0.046 0.023 0.046]*his/2;
  min = [0.69 0.6624 0.11 0.16 0.21 0.046 0.046 0.023 0.023 0.023]*his/2;
  rot = [90 90 72 108 90 0 0 0 0 90]*pi/180;
  ref = [1 -0.98 -0.02 -0.02 0.01 0.01 0.01 0.01 0.01 0.01];
  rcos = cos(rot); rsin = sin(rot);

  for i=1:his,
      for j=1:his,
	for k=1:10,
	    x = rcos(k)*(dcos(i)-dsin(j)+sx*dcos(j)+cx(k)) ...
		+ rsin(k)*(dsin(i)+dcos(j)+sx*dsin(j)+cy(k));
	    y = - rsin(k)*(dcos(i)-dsin(j)+sx*dcos(j)+cx(k)) ...
		+ rcos(k)*(dsin(i)+dcos(j)+sx*dsin(j)+cy(k));
	    d = sqrt( x^2/maj(k)^2 + y^2/min(k)^2 );
	    if (d <= 1), hiIm(i,j) = hiIm(i,j)+ref(k); end;
	end;
      end;
  end;

  pts = find(hiIm>0.5);
  hiIm(pts) = hiIm(pts) - 0.95;
  hiIm = hiIm'/0.05;

  % low-pass filter hi res image to nyquist of lo res
  hlow = sinc((-3*fac:3*fac)/fac)/fac;
  hlow = hlow.*hamming(length(hlow))';
  filtHiIm = conv2(hlow,hlow,hiIm);

  % sample hi res image to get low res
  shiftx = round(x0*fac);
  xi = shiftx+1:fac:shiftx+his;
  xi = mod(xi-1,his)+1;

  shifty = round(y0*fac);
  yi = shifty+1:fac:shifty+his;
  yi = mod(yi-1,his)+1;
	
  Im = filtHiIm(xi,yi);

return;