%radon demo
%
%this computer artifical data for use is testing
%The analysis program is below
%
function []=radon_demo_pd(void)
%%code for generating simulated data
spatial_freq=1;
npoints=64*spatial_freq;%number of points in a line
nlines=64*8*2*spatial_freq;%number of lines 
zz=zeros(nlines,npoints);
f0=30;
phi0=(pi/25);%-0.3571;
phi1=1;
f_phi=pi;
for i=1:npoints
    for j=1:nlines
        phi=phi0*i;%
        f=f0+phi1*sin(((f_phi*2*pi*j)/spatial_freq)/1000);
        zz(j,i)=.5+.5*sin((2*pi*f*j/spatial_freq)/1000+phi);
    end
end
figure(1)
subplot(211)
imagesc(zz')%plots the artificial 'linescan'
colormap gray
[thetasz32,the_tz32,spread_radon32]=GetVelocityRadonFig_demo(zz,npoints);
subplot(212)
plot(the_tz32,thetasz32,'b')% plots the angle at any given time point
hold on
xlabel('time,A.U.')
ylabel('angle, degrees')

%% data = zz

%____________________________________________________

% enter the data file below
		%data =zz if you use the fictional data above.
		%An example of a reading a TIFF id given below
data = imread('XT_20X_3X_HLFS_30mA_R12_175um.tif');
imagesc(data(1:1000,5:34));
dx=1.0e-6; %pixel length
dt=2.05e-6; %pixel clock
%INPUTS
%data - the matrix of time X space data 
%windowsize - number of lines to use in estimating velocity; must be a multiple of 4
data=double(data(1:29000,5:34));
%function [thetas,the_t,spread_matrix]=GetVelocityRadonFig_demo(data,windowsize);
windowsize=32% 45 ms of fNyquistof 11Hz
%%OUTPUTS
%thetas - the time varying angle of the space-time image
%the_t - time pointsof the angle estimates (in lines)
%spreadmatrix - matrix of variances as a function of angles at each time point
stepsize=.25*windowsize;
nlines=size(data,1);
npoints=size(data,2);
nsteps=floor(nlines/stepsize)-3;
%find the edges
angles=(0:179);
angles_fine=-2:.25:2;

spread_matrix=zeros(nsteps,length(angles));
spread_matrix_fine=zeros(nsteps,length(angles_fine));
thetas=zeros(nsteps,1);

hold_matrix=ones(windowsize,npoints);
blank_matrix=ones(nsteps,length(angles));
the_t=NaN*ones(nsteps,1);

for k=1:nsteps
    the_t(k)=1+(k-1)*stepsize+windowsize/2;
    data_hold=data(1+(k-1)*stepsize:(k-1)*stepsize+windowsize,:);
    data_hold=data_hold-mean(data_hold(:))*hold_matrix;%subtract the mean
    radon_hold=radon(data_hold,angles);%radon transform
    spread_matrix(k,:)=var(radon_hold);%take variance
    [m the_theta]=max(spread_matrix(k,:));%find max variace
    thetas(k)=angles(the_theta);     
    radon_hold_fine=radon(data_hold,thetas(k)+angles_fine);%re-do radon with finer increments around first estiamte of the maximum
    spread_matrix_fine(k,:)=var(radon_hold_fine);
    [m the_theta]=max(spread_matrix_fine(k,:));
    thetas(k)=thetas(k)+angles_fine(the_theta);
end
thetas=thetas-90; %rotate
v=(dx/dt)*cot(pi*thetas/180);
plot(the_t, v)


