function modgraph(shift,grow_hrt,grow_lng,res)
%
%  modgraph(shift,grow_hrt,grow_lng,res)
%
%  Modgraph will display a synthetic 2-D slice model
%  of a human chest using some of the input arguments
%  needed by PROJECT.  SHIFT allows a horizontal translation
%  while GROW_HRT allows a dilation of the "heart", and
%  GROW_LNG allows a dialation of the "lungs" and "chest wall".
%  The model is contained within a square box with vertices at (0,0)
%  and (1,1).  SHIFT is limited to +/-0.04, GROW_HRT
%  is limited from 0 to 0.04, and GROW_LNG is limited from
%  0 to 0.08.  RES is the rectangular sampling
%  grid.  A suggestion for RES is 0.004.
%
%  SEE project.m, projgraf.m
%
%  Developed by Timothy D. Dorney
%               Rice University
%               February, 1999
%               tdorney@ieee.org
%
%  Coded using MATLAB 5.X.X with NO additional toolboxes.
%
%	REVISION HISTORY
%
%	VERSION 1.0.0		FEB. 15, 1999	TIM DORNEY
%       VERSION 1.0.1		APR. 12, 1999	TIM DORNEY
%		Modified the sampling grid to match the
%		one defined in the back-projection algorithm.
%
%
%  NOTE:  MODGRAPH WILL GIVE ONLY APPROXIMATE ANSWERS!!!!
%         MODGRAPH IS FAST, BUT FOR BEST RESULTS USE PROJGRAF.M
%

if (nargin ~= 4)
        disp('MODGRAPH requires 4 input arguments!')
	return;
end
 
if ((shift < -0.04) | (shift > 0.04))
	shift = 0;
	disp('SHIFT limited to +/- 0.04.  The limit was exceeded.  SHIFT set to 0!');
end
if ((grow_hrt < 0) | (grow_hrt > 0.04))
	grow_hrt = 0;
	disp('GROW_HRT limited from 0 to 0.04.  The limit was exceeded.  GROW_HRT set to 0!');
end
if ((grow_lng < 0) | (grow_lng > 0.08))
	grow_lng = 0;
	disp('GROW_LNG limited from 0 to 0.08.  The limit was exceeded.  GROW_LNG set to 0!');
end

% DENSITY		REGION		PART			GROWTH
DA = 1.00;		% A		center spine circle     N
DB = 0.99;		% B             left spine ellipse	N
DC = 0.98;		% C		right spine ellipse	N
DD = 0.5;		% D		center heart circle	Symmetric
DE = 0.6;		% E		left lung ellipse	Asymmetric left
DF = 0.61;		% F		right lung ellipse	Asymmetric right
DG = 0.4;		% G		outer skin ellipse	Asymmetric vertical

colormap(gray(255));	xres = res;	yres = res;
xxres = [xres/2:xres:1]; yyres = [yres/2:yres:1];
yy = (fliplr(yyres))'*ones(size(yyres)); xx = (xxres'*ones(size(xxres)))';
map = zeros(length(xx),length(yy));

% REGION G
A = 2.25;	B = 1.95;	C = 0.03;	D = (sqrt(C+grow_lng/10) - sqrt(C))*B;
[e,f] = find( ( ((xx-0.5-shift)./A).^2 + ((yy-0.35-D)./B).^2 ) <= C+grow_lng/10);
map((f-1)*length(yy)+e) = DG;

% REGION A
[e,f] = find (((xx-0.5-shift).^2 + (yy-0.1).^2) <= 0.0009);
map((f-1)*length(yy)+e) = DA;

% REGION B
A = 2;	B = 1;
[e,f] = find (( ((xx-0.415-shift)./A).^2 + ((yy-0.1)./B).^2) <= 0.000625);
map((f-1)*length(yy)+e) = DB;

% REGION C
A = 1;	B = 2;
[e,f] = find (( ((xx-0.56-shift)./A).^2 + ((yy-0.1)./B).^2) <= 0.000625);
map((f-1)*length(yy)+e) = DC;

% REGION D
[e,f] = find (((xx-0.5-shift).^2 + (yy-0.25).^2) <= 0.003+grow_hrt/8);
map((f-1)*length(yy)+e) = DD;

xxres = [-0.5+xres/2:xres:0.5];	yyres = [-0.5+yres/2:yres:0.5];
yy = (fliplr(yyres))'*ones(size(yyres));	xx = (xxres'*ones(size(xxres)))';
% map = zeros(length(xx),length(yy));

% REGION E
A = 2.2;	B = 1.3;	C = 0.008;	D = (sqrt(C+grow_lng/20) - sqrt(C))*B;	ANG = pi/3;
xxp = xx*cos(ANG)+yy*sin(ANG);	yyp = -xx*sin(ANG)+yy*cos(ANG);
[e,f] = find( ( (xxp./A).^2 + (yyp./B).^2 )  <= C+grow_lng/20);
e = e + round(0.08/yres - D/yres);
f = f - round(0.192/xres - shift/xres);
map((f-1)*length(yy)+e) = DE;

% REGION F
A = 2.2;	B = 1.3;	C = 0.008;	D = (sqrt(C+grow_lng/20) - sqrt(C))*B;	ANG = 2*pi/3;
xxp = xx*cos(ANG)+yy*sin(ANG);	yyp = -xx*sin(ANG)+yy*cos(ANG);
[e,f] = find( ( (xxp./A).^2 + (yyp./B).^2 )  <= C+grow_lng/20);
e = e + round(0.08/yres - D/yres);
f = f + round(0.192/xres + shift/xres);
map((f-1)*length(yy)+e) = DF;

clf
disp('THE IMAGE SIZE IS:')
size(map)
imagesc([xres/2:xres:1],[yres/2:yres:1],map);
axis('image');