%%loading
cd 'f:/LAB/LAB2011/17january112mMCaQd/Series089/' %assign dir with files
files=dir('*.tif'); %enumerate all files
%nr = 60;
nr=numel(files); % count all files
img=zeros(512,512,nr,'uint8');%create nr images
for i=1:nr;
img(:,:,i)=imread(files(i).name);
end
filt_img=zeros(512,512,nr,'uint8');
for i=1:nr;
filt_img(:,:,i)=medfilt2(img(:,:,i),[6,6]);
end
for i=1:nr;
imshow(filt_img(:,:,i));
pause(0.01);
end
%for i=1:nr;                         show images
%imshow(filt_img(:,:,i));
%pause(0.05);
%end
%for i=1:nr;          was used to determine center point with cursor
%    imshow(filt_img(:,:,902))
% white mask with black hole in center (hole diam is determined from patch)
% ->> start mask code
 %determine radius l for mask from center line medfilt2
%% radius
 mrad=zeros(512,'double');
for m =1:nr;  %determine patch edge
   for l=0:250;
           intensity=filt_img(256-l,256,m); % determine intensity point 1
           intensity1=filt_img(256-l,257,m); % determine intensity point 2
           if (intensity < 120)  && (intensity1  < 120) % compare with threshold
              mrad(m) = l; %write intensity array
              break % leave if edge is reached
           end 
   end
end

limit_img=zeros(512,512,'uint8'); %create a blank mask image filled black 
%used to create a circle of radius l pixels

% <<- end mask code

for i=1:nr; % substract mask from median filtered images
    
     for ii=1:512;
      for jj=1:512;
        limit_img(ii,jj)=255; %make white
        r=sqrt((ii-256)*(ii-256)+(jj-256)*(jj-256)); %determine distance of pixel from center
           if (r < (mrad(i)+1)) limit_img(ii,jj)=0; %check if within lipid patch
           end
      end
     end
filt_img(:,:,i)=filt_img(:,:,i)-limit_img; % subtract masks
imshow(filt_img(:,:,i)); %show image for checking
pause(0.01);
end
for i=1:nr; % bw convert and show masked images
imshow(im2bw(filt_img(:,:,i),0.6)); %show black/white image for checking 
pause(0.01);
end


%% display
bw_img=zeros(512,512,nr,'uint8'); %create bw images
for i=1:nr;
bw_img(:,:,i)=im2bw(filt_img(:,:,i),0.6); %fill with filtered and masked images
end
mass=zeros(nr,1,'double'); %define sum variable array

for k=1:nr;
mass(k)=sum(sum(bw_img(:,:,k)));% get sum of intensities of all image pixels
end
x=97.8698/512; % length per pixel
area=mass*x*x; % area of the patch (bw image: mass is nr of pixels)
rr=sqrt(mass/3.141519)*x; % derive radius from no of pixels and unit length
rr=tsmovavg(rr,'s',15,1) %
y=1474.77/902; % total time devided by no of frames = time betw 2 frames
t=zeros(1,nr); % create zero filled array for time assignment
for i=1:nr; % array desribing timeline, t(1) = 0s ->>
t(i)=(i-1)*y;
end  %  <<-

%% plotting area
plot(t,area); % plot area vs time

%%plotting velocities
% velocities for three points
vv=zeros(nr,1,'double'); % create zero filled array for velocities
for i=4:nr-3
vv(i)=(rr(i+1)-rr(i-1))/(2*y); %velocity over three points
end
plot(t,vv) % plots velocity vs time

%% end here

% velocities for five points 
vv2=zeros(nr,1,'double'); % create zero filled array for velocities

%cut 10 frames from end
for i=4:nr-10 
vv2(i)=((rr(i+1)-rr(i-1))+rr(i+2)-rr(i-2))/(4*y);
end
plot(t,vv2) % plots velocity vs time

%% write images to files

for i=1:nr
imwrite(filt_img(:,:,i),strcat('fseries',int2str(i),'.tif')); %,'tiff');    
end