Quantification of movement directionality of EB3-GFP "comet tails"
The following script can be copied into a MATLAB script editor:
%% presenting the images for analysis
image1=imread('3.tif');
image2=imread('4.tif');
I1 = im2double(image1);
I2 = im2double(image2);
figure, subplot(2,1,1),imshow(I1), xlabel('image 1'),subplot(2,1,2), imshow(I2),xlabel('image2');
%% limiting analysis to ROI:
figure,imshow(I1),title('choose an area of interest');
rect=getrect
shaft=I1(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3));
shaft2=I2(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3));
figure,subplot(2,1,1),imshow(shaft),xlabel(' axonal shaft of image 1'),subplot(2,1,2),imshow(shaft2),xlabel(' axonal shaft of image 2')
%% adjusting 2 images to saturate 2 percent of pixels+gamma enhancment:
adjusted = imadjust(shaft,stretchlim(shaft),[],1.2);
adjusted2 = imadjust(shaft2,stretchlim(shaft),[],1.2);
figure, subplot(2,1,1),imshow(adjusted),xlabel('0 sec'),subplot(2,1,2),imshow(adjusted2),xlabel('5 sec')
%% creating a padded matrice in order not to exceed "adjusted" dimensions
dimension=size (adjusted)
pad1 =zeros(dimension(1)+100,dimension(2)+100);
pad2 =zeros(dimension(1)+200,dimension(2)+200);
pad1(51:end-50,51:end-50)=pad1(51:end-50,51:end-50)+adjusted;
pad2(101:end-100,101:end-100)=pad2(101:end-100,101:end-100)+adjusted2;
figure,subplot (2,1,1),imshow(pad1),xlabel ('padded image 1'),subplot (2,1,2),imshow(pad2),xlabel ('padded image 2'),
%% converting the padded matrice into binary images
level = graythresh(pad1);
level2= graythresh(pad2);
binary = im2bw(pad1,0.5);
binary2 = im2bw(pad2,0.5);
figure,subplot(2,1,1),imshow(binary),xlabel('this is the first frame converted into a binary image'),subplot(2,1,2),imshow(binary2),xlabel('this is the second frame converted into a binary image'),
% %% eroding the pictures to reduce noise:
% SE=strel('diamond',1);
% clean1=imerode(binary,SE);
% clean2=imerode(binary2,SE);
% figure,subplot(2,1,1),imshow(clean1),subplot(2,1,2),imshow(clean2);
%% creating a labeled matrix and presenting each candidate EB3 as a discrete component:
[labeled,num]=bwlabel(binary,4);
num
show = label2rgb(labeled, 'jet', 'k', 'shuffle');
figure,imshow(show),xlabel('every labeled object is marked differently'),title('candidate signals')
region=regionprops(labeled,'All')
%% Now i shall eliminate any signals smaller than 4 pixels (noise) yielding relevat signals
Real_EB3=binary;
for j=(1:num)
if region(j).Area<4
Real_EB3(region(j).BoundingBox(2):region(j).BoundingBox(2)+region(j).BoundingBox(4),region(j).BoundingBox(1):region(j).BoundingBox(1)+region(j).BoundingBox(3))=0;
elseif region(j).Area>100
Real_EB3(region(j).BoundingBox(2):region(j).BoundingBox(2)+region(j).BoundingBox(4),region(j).BoundingBox(1):region(j).BoundingBox(1)+region(j).BoundingBox(3))=0;
end
end
%% recreating the regionprops array and displaying the cleaned image
Real_EB3=bwlabel(Real_EB3,4);
region2=regionprops(Real_EB3,'All')
show2=label2rgb(Real_EB3, 'jet', 'k', 'shuffle');
figure,imshow(show2),xlabel('these are the discrete EB3'),colormap(jet)
%% checking the new number of EB3s:
[number,num]=bwlabel(Real_EB3,4);
num
%% now i will create vectors to contain the movement info
angle=(1:num);
x=(1:num);
y=(1:num);
%% here i shall obtain the vector for each movement
show2=label2rgb(Real_EB3, 'jet', 'k', 'shuffle');
figure,imshow(show2),xlabel('lines represent movement vectors')
size1=size(Real_EB3)
for i=1:num
if region2(i).BoundingBox(3)>0
smallwin=pad1(ceil(region2(i).BoundingBox(2)):ceil(region2(i).BoundingBox(2))+ceil(region2(i).BoundingBox(4)),ceil(region2(i).BoundingBox(1)):ceil(region2(i).BoundingBox(1))+ceil(region2(i).BoundingBox(3)));
bigwin=pad2(ceil(region2(i).BoundingBox(2))-4.2+50:ceil(region2(i).BoundingBox(2))+ceil(region2(i).BoundingBox(4))+4.2+50,ceil(region2(i).BoundingBox(1))-4.2+50:ceil(region2(i).BoundingBox(1))+ceil(region2(i).BoundingBox(3))+4.2+50);
%correlating the two windows
cor=imfilter(bigwin,smallwin);
% finding in the correlation window the regional maxima
high=max(cor(:));
[row,col]=find(cor==high);
%finding center of mass for eb3 (the small window)
center2(1)=(region2(i).Centroid(1));
center2(2)=(region2(i).Centroid(2));
%determining the xy coordinates for movement vector
beginX=center2(1);
beginY=center2(2);
endX= col-5+ region2(i).BoundingBox(1);
endY= row-5+ region2(i).BoundingBox(2);
endX=min(endX);
endY=min(endY);
line([beginX endX],[beginY endY],'color','y','LineWidth',5)
angle(i)=atan2((endY-beginY),(endX-beginX))*57.296;
x(i)=(endX-beginX);
y(i)=(endY-beginY);
end
end
%% analizing
% plotting a histogram of angles
snum=num2str(num);
STD=round(std(angle));
sSTD=num2str(STD);
figure, hist(angle,20),title('Histogram of angles'), xlabel('angles'),ylabel('N');
text(-90,11,[' n=' snum ],'EdgeColor','red')
set (gca,'xtick',[-180 -90 0 90 180]);
% plotting the sum of vectors
radians=angle*pi/180;
figure, rose(radians);
%%