%	to read a physiology file and look at sort
%	and correct image series
%	rev 1 	11/22/99	does the phys sort internally


	more off;
	npix = 128*128;

	fnin = input('gimme infile path=', 's');
	fnout = input('gimme outfile path=', 's');

	foo = input('gimme baseline frames [img1 img2] =');
	img1 = foo(1);  img2 = foo(2);
	foo = input('gimme [nframes/slc islc nslc TR (ms)]=');
	nframes = foo(1);  islc = foo(2);  nslc = foo(3); TR = .001*foo(4);
	flag = input('ECG (1), resp(2) or both(3)=');
        imgo = (islc-1)*nframes;

%  get the physiology data
 
phaser = zeros(1,nframes);
phasec = phaser;
FSAMP = 40;             % sampling rate for phys data/s
fname = input('gimme physiology file =','s');
fid = fopen(fname, 'r');
[dat nf] = fscanf(fid, '%g %g\n');
fclose(fid);
n=nf/2;
ecg = -dat(1:2:nf);     % invert polarity
resp = dat(2:2:nf);

dt = 1./FSAMP;          % in seconds

timei = input('gimme number of seconds to ignore [1] =');
if(isempty(timei))
  timei = 1;
end

% find the ecg trigs

nt = 0;
for (j=2:n)
  if (ecg(j)>=0 & ecg(j-1)<0)  % first only
    nt = nt + 1;
    trig(nt) = j*dt;
  end
end
trig(nt+1) = 2*trig(nt) -trig(nt-1);
averr = (trig(nt) - trig(1))/nt;
avehr = 60/averr;               % beats/min
fprintf('average HR = %d beats/min\n', fix(avehr));

% calc ecg phase

time0 = timei + (islc-1)*TR/nslc;
tfoo = input('gimme fudge(s)=');
time0 = time0 + tfoo;

m = 1;
for ifr=0:nframes-1;
  time = time0 + ifr*TR;
  while(time>trig(m))
    m = m+1;
  end
  t2=trig(m);
  t1=trig(m-1);
  phasec(ifr+1) = (time - t1)*2*pi/(t2 - t1);
end
 
%  now do the resp.  Calculate histogram and generate transfer function
%  from the integral
 
respx = max(resp);
respn = min(resp);
resp = 100*(resp - respn)/(respx - respn);
 
rbin = zeros(1,101);  % 1% increments
for (j=1:n)
  k = 1 + round(resp(j));
  rbin(k) = rbin(k) + 1;
end
 
pdf = cumsum(rbin);
pdf = 2*pi*pdf/max(pdf);
 
for ifr=0:nframes-1;
  time = time0 + ifr*TR;
  isamp = round(time/dt);
  phaser(ifr+1) = pdf(1 + round(resp(isamp)));
end

%   fit timeseries to low order sine function
	
	nfrb = img2-img1+1;
	scl = 2/nfrb;
	c1c = cos(phasec);
	s1c = sin(phasec);
	c2c = cos(2*phasec);
	s2c = sin(2*phasec);
	a0=zeros(npix,1);
	a1c=a0;
	b1c=a0;
	a2c=a0;
	b2c=a0;
	c1r = cos(phaser);
	s1r = sin(phaser);
	c2r = cos(2*phaser);
	s2r = sin(2*phaser);
	a1r=a0;
	b1r=a0;
	a2r=a0;
	b2r=a0;
	for img = img1:img2
	  fnm = sprintf('%s.%03d',fnin, img+imgo);
	  fid = fopen(fnm, 'r');
	  if(fid==-1) 
	    fprintf('U turkey!  poo poo input file...\n'); 
	    return; 
	  end;
	  fprintf('read %s\n ', fnm);
          buf = fread(fid, 'short');
	  fclose(fid);
          npix=length(buf);
	  a0 = a0+ buf;
	  a1c = a1c + buf*c1c(img);
	  b1c = b1c + buf*s1c(img);
	  a2c = a2c + buf*c2c(img);
	  b2c = b2c + buf*s2c(img);
	  a1r = a1r + buf*c1r(img);
	  b1r = b1r + buf*s1r(img);
	  a2r = a2r + buf*c2r(img);
	  b2r = b2r + buf*s2r(img);
        end
	a0 = a0/nfrb;
	a1c = (a1c - a0*sum(c1c(img1:img2)))*scl;
	b1c = (b1c - a0*sum(s1c(img1:img2)))*scl;
	a2c = (a2c - a0*sum(c2c(img1:img2)))*scl;
	b2c = (b2c - a0*sum(s2c(img1:img2)))*scl;
	a1r = (a1r - a0*sum(c1r(img1:img2)))*scl;
	b1r = (b1r - a0*sum(s1r(img1:img2)))*scl;
	a2r = (a2r - a0*sum(c2r(img1:img2)))*scl;
	b2r = (b2r - a0*sum(s2r(img1:img2)))*scl;

%  apply correction

	foo = input('gimme frames to correct [img1 img2] =');
	img1 = foo(1);  img2 = foo(2);
	for img = img1:img2
	  fnm = sprintf('%s.%03d',fnin, img+imgo);
	  fid = fopen(fnm, 'r');
	  if(fid==-1) 
	    fprintf('U turkey!  poo poo input file...\n'); 
	    return; 
	  end;
	  fprintf('read %s  ', fnm);
          buf = fread(fid, 'short');
	  fclose(fid);
	  yc = a1c*c1c(img) + b1c*s1c(img) + a2c*c2c(img) + b2c*s2c(img);
	  yr = a1r*c1r(img) + b1r*s1r(img) + a2r*c2r(img) + b2r*s2r(img);
	  if(flag==1)
	    buf = buf - yc; 
          elseif(flag==2)
	    buf = buf - yr; 
	  else
	    buf = buf - yc - yr;
	  end
	  fnm = sprintf('%s.%03d',fnout, img+imgo);
	  fprintf('write %s\n', fnm);
	  fid = fopen(fnm, 'w');
	  fwrite(fid, buf, 'short');
	  fclose(fid);
	end

	more on;
