SHARE
TWEET

Untitled

a guest Jul 17th, 2017 74 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. %----------------------------------------------------------------------------------------------
  2. %A script to recontruct  a microtubule image, and find the bending modes
  3. %Authors: David Apigo, Arooj Aslam, John Palmieri
  4. %Date: September 5,
  5. %Version Number: 1.2
  6. %----------------------------------------------------------------------------------------------
  7.  
  8. %create file called snakes
  9. clearvars -except snakes
  10.  
  11. %% Image Import
  12. %==============
  13. ii = 0;
  14. %input('How many microns per pixel: ');  %100x Mag has 15.3 pix/micron
  15. u = .06;
  16.  
  17. %User Input
  18. %==========
  19. display('Select the file for the original images')
  20.  
  21. %path and filename for original image
  22. [filename_orig,path_orig] = uigetfile('*.tif','Select Original Image Stack');
  23.  
  24. frame_num = input('How many images are you processing: ');
  25. begin = input('Number of the first image ');
  26. increment = input('How many images do you want to increment by? ');
  27. end_image = input('Number of the last image ');
  28. pt_spacing = input('What is the point spacing along the MT? ');
  29.  
  30. %initializing index
  31. index = begin;
  32.  
  33. %Inputs all images into a 3D matrix
  34. for ii = 1:((end_image-begin)/increment)+1
  35.     I_orig(:,:,ii) = imread(strcat(path_orig, filename_orig),'Index',index);
  36.     index =  index + 1;
  37. end
  38.  
  39. %% Setting up or initializing parameters for the images that will be used later in the code
  40. %==========================================================================================
  41.  
  42. %Initializing Certain Parameters
  43. %===============================
  44. %spline smoothing parameter parameter = .9;
  45. sz = 10;
  46.  
  47. %initializing count for gauss fit errors
  48. catch_count = 0;
  49.  
  50. %Pads the array with zeroes and gets image information for skeletonized image
  51. I_orig_padded = padarray(I_orig,[sz sz]);
  52.  
  53. % I_skel = padarray(I_skel_unpadded,[sz sz]);
  54.  
  55. %Width and height of the image in pixels
  56. [width,height] = size(I_orig_padded);
  57.  
  58. %Initializing the size of the reconstructed microtubule matrix to be big
  59. %enough to include the recostructed pixels from ther largest skeletonized image
  60. previous_frame = 1;
  61.  
  62. %number of coordinates ("beads")
  63. N = 0;
  64.  
  65. count = 1;
  66.  
  67. for recon_size = 1:size(snakes,1)
  68.     current_frame = snakes(recon_size,1);
  69.     if previous_frame ~= current_frame
  70.         N(count)= snakes(recon_size-1,2)+1;  %plus 1 b/c the count starts at zero in data file
  71.         count = count+1;
  72.     end
  73.     previous_frame = current_frame;
  74. end
  75.  
  76. %puts in count of points for last frame and finds max N
  77. %======================================================
  78. %number of beads in each frame
  79. N = [N snakes(size(snakes,1),2)];
  80.  
  81. %max number of beads across all frames
  82. N_max = max(N);
  83.  
  84. %Generate coord matrix using max N row dimension
  85. coord_reconstructed = zeros(N_max,2,length(N)); %[number of pixels, x-y cols, number of images]
  86.  
  87. %% Filtering Original Image
  88. %===========================
  89. sigma = 2;
  90. w = 20;
  91.  
  92. filter_image  = input('Does this image need to be filtered? ','s');
  93.  
  94. if filter_image == 'y' || filter_image == 'Y'
  95.     filter_flag = 0;
  96. else
  97.     filter_flag =1;
  98.     I_filtered = I_orig;
  99. end
  100.  
  101. while filter_flag == 0
  102.    
  103.     %First smoothing function does a scanning average of the neighboring pixels
  104.     %using a mask the size of w by w
  105.     I_boxavg = imboxfilt3(I_orig,(2*w+1));
  106.    
  107.     %Second smoothing function uses a gaussian
  108.     %convolution kernel
  109.     I_gaussfilt = imgaussfilt3(I_orig,sigma);%gauss_smooth(I_orig,w);
  110.     %I_gaussfilt = imgaussfilt3(I_orig);%gauss_smooth(I_orig,w);
  111.    
  112.     %Made uint8 a double to deal with neg values
  113.     I_filtered = abs((I_gaussfilt) - (I_boxavg));
  114.     imshow(I_filtered(:,:,1),[])
  115.    
  116.     % subplot(221)
  117.     % imshow(I_boxavg(:,:,1),[])
  118.     % subplot(222)
  119.     % imshow(I_gaussfilt(:,:,1),[])
  120.     % subplot(223)
  121.    
  122.    
  123.     f = input('Is this filter satisfactory? Y/N ','s');
  124.    
  125.     if f == 'N' || f == 'n'
  126.        
  127.         %used as the region around the MT for image rotation
  128.         w = input('What is the window size for box filter?');
  129.        
  130.         %used as the region around the MT for image rotation
  131.         sigma = input('Standard deviation for Gaussian Filter?');
  132.        
  133.     end
  134.    
  135.     if f == 'Y' || f == 'y'
  136.         filter_flag = 1;
  137.     end
  138.    
  139. end
  140.  
  141. %% Filament Reconstruction
  142. %Looks at one image at a time and:
  143. %Finds the fit using the skeletonized image to find the tangential angles
  144. %Rotates small region of the bandpassed image around the points taken from the polyfit function
  145. %Fits a gaussian to the intensity profile in each column of the rotated image
  146. %Generates a reconstructed MT using the max of the gaussian as the center of the MT
  147. k =0;
  148. N_previous = 0;
  149.  
  150. for filament_recon=1:size(I_orig,3) %size (I_skel,3) returns the number of images, so this runs for each image in the tif file
  151.     clear x;
  152.     clear y;
  153.     %Pads each image matrix with zeroes to deal with points near the edge
  154.     I_new_filtered = padarray(I_filtered(:,:,filament_recon),[sz sz]);
  155.    
  156.     %Get points from snake
  157.     x = snakes((N_previous+1):(N(filament_recon)+ N_previous),3)+10;
  158.     y = snakes((N_previous+1):(N(filament_recon)+ N_previous),4)+10;
  159.    
  160.     N_previous = N(filament_recon)+ N_previous;
  161.    
  162.     %xx = [x(1):pt_spacing:N(filament_recon)]';
  163.    
  164.     %Fit a spline to the skeleton and evaluate its derivative
  165.     spline = fit(x,y,'smoothingspline');
  166.     %     plot(spline,x,y,'-')
  167.     y_spline = spline(x);
  168.     % %     hold on
  169.     %      plot(xx,y_spline,'-')
  170.     %     hold off
  171.     differential_spline = differentiate(spline,x);
  172.    
  173.     %Go through each point in the differential to
  174.     %find the angle of rotation
  175.     theta = double(rot_angle(differential_spline,x));
  176.    
  177.    
  178.     %Grab a section (+-w) of the image around x0 to rotate by theta(x0)
  179.     %Take the rotated section and find the new x,y coord
  180.     for k = 1:length(x)%1:size(x)
  181.        
  182.         try
  183.             %zooming into a wxw section around the point of interest
  184.             %row_range = y(k)-(w):y(k)+(w);   %following the points of the skel image
  185.             %col_range = x(k)-(sz):x(k)+(sz);
  186.             row_range = round(y_spline(k))-(sz):round(y_spline(k))+(sz);   %following the points of the spline
  187.             col_range = round(x(k))-(sz):round(x(k))+(sz);
  188.            
  189.            
  190.            
  191.             I = I_new_filtered(row_range,col_range);
  192.             I_rot = imagerot(I,(theta(k)));
  193.            
  194.             %Use for checking rotation of MT
  195.             %          figure(2);
  196.             %          imshow(I_rot,[]);
  197.             %
  198.             %
  199.             %Find max intensity of the middle column of the rotated image
  200.             %section
  201.             I_midpoint = (size(I,1)+1)/2;
  202.             rows_rot = 1:size(I_rot,1); %number of rows
  203.             intensity = double(I_rot(:,I_midpoint)); %gaussfit needs Y data to be a double
  204.             gaussfit = fit(rows_rot.',intensity,'gauss1');
  205.             gauss_centroid = gaussfit.b1;  %x-coord of the peak, gives the y value of the rotated image section
  206.             %         gauss_centroid_matrix(k) = gaussfit.b1;  %x-coord of the peak, gives the y value of the rotated image section
  207.             gauss_centroid_matrix(k,filament_recon) = gaussfit.b1;  %x-coord of the peak, gives the y value of the rotated image section
  208.            
  209.             % Use for plotting the gaussian fit
  210.             %          figure(5)
  211.             %          plot(rows_rot,intensity,'-d')
  212.             %          hold on
  213.             %          plot(gaussfit,rows_rot,intensity)
  214.             %          hold off
  215.             %          pause
  216.             % %
  217.            
  218.             % rotate the wxw section back
  219.             coord_unrotated = [cos(-theta(k)) sin(-theta(k));-sin(-theta(k)) cos(-theta(k))]*[0;gauss_centroid-I_midpoint];
  220.             coord_unrotated_gauss_peak(k,filament_recon) = coord_unrotated(2,1);
  221.             %translate coordinate back to its position in the larger image using skel image coorinates
  222.            
  223.             %         row_reconstructed(k) = coord_unrotated(2,1)+(y(k));  %row number of final position
  224.             %         col_reconstructed(k) = coord_unrotated(1,1)+x(k);  %column number of final position
  225.             %
  226.             %coord_reconstructed(k,1,i) = x;  %uses the x values from the
  227.             %spline fit
  228.             coord_reconstructed(k,1,filament_recon) = coord_unrotated(1,1)+round(x(k));  %column number (x-coord) of final position
  229.             coord_reconstructed(k,2,filament_recon) = coord_unrotated(2,1)+(round(y_spline(k)));  %row number (y-coord) of final position
  230.            
  231.            
  232.         catch exception
  233.             catch_count = catch_count+1;
  234.             continue   %if the gauss fit does not converge, then return to the beginning of the loop
  235.         end
  236.     end
  237.    
  238.     %Generate a concatenated [max number of pixels x 2 x number of images] matrix of coords for all images
  239.     %coord_reconstructed(:,:,filament_recon)  = [col_reconstructed' row_reconstructed'];
  240.    
  241.     %     %Plot skel image, its fit, and the final reconconstructed MT
  242.     %      f = figure(filament_recon);
  243.     %      hold on
  244.     %      plot(spline,x,y, 'b');     %Plot spline fit over skeletonized image
  245.     % %      print(f,'MT_fit.tif')
  246.     %
  247.     %   %plot reconstructed MT on original image
  248.     figure(6)
  249.     imshow(I_new_filtered(:,:),[])
  250.     hold on
  251.     coord_reconstructed(coord_reconstructed == 0) = NaN;
  252.     plot(coord_reconstructed(:,1,filament_recon),coord_reconstructed(:,2,filament_recon),'r-');
  253.     hold off
  254. end
  255.  
  256. %% Use for saving coords of each frame so that it can be used in origin
  257.  
  258. Colx = 1;
  259. Coly = 2;
  260. for frame = 1:size(I_orig,3)
  261.     if frame == 1
  262.         excel_data(:,[Colx Coly]) = coord_reconstructed(:,[1 2],frame);
  263.     else
  264.         Colx = Colx +2;
  265.         Coly = Coly +2;
  266.         excel_data(:,[(Colx) (Coly)]) = coord_reconstructed(:,[1 2],frame);
  267.     end
  268.    
  269.    
  270. end
  271.  
  272. %% Filament Tracking - Fourier Analysis
  273.  
  274. %User Input
  275. mode_num = input('What mode number do you want to plot up to?');
  276. %initializing variables
  277. count = 1;
  278. frame = 1;
  279. a_k = 0;
  280. k = 0;
  281. n = 0;
  282.  
  283.  
  284. %Calculate fourier coefficients
  285. for frame_fourier = 1:size(I_skel,3) %number of frames
  286.     %Grab coordinates of frame j
  287.     %row = find(~isnan(coord_reconstructed(:,1,j)),1,'last');
  288.     coord_xy_pixels = coord_reconstructed(:,:,frame_fourier);
  289.     coord_xy = (1/15.3).*coord_reconstructed(:,:,frame_fourier); %converting from pixels to microns 15.3 pix/um
  290.    
  291.     %Get the angle and segment length of each segment
  292.     %Calculate the inverse fourier transform for each segment and save
  293.     %in a vector
  294.    
  295.     %If there are zeroes in the coord_xy variable, they are ignored by
  296.     %forcing the last non-zero entry to be the last xy_coord for that
  297.     %frame
  298.     if ~any(coord_xy)
  299.         end_point = find(coord_xy == 0,1,'first')-1;
  300.     else
  301.         end_point =  find(coord_xy(:,1),1,'last');
  302.     end
  303.    
  304.     for k = 1:end_point-1
  305.         %Calculate angle of segment k
  306.         delta_x = coord_xy(k+1,1) - coord_xy(k,1);
  307.         delta_y = coord_xy(k+1,2) - coord_xy(k,2);
  308.         slope = (delta_y/delta_x);
  309.         theta_s(k) = atan(slope);
  310.         theta_s(isnan(theta_s)) = 0;
  311.        
  312.         %Calculate segment k length
  313.         delta_s(k) = sqrt(delta_x^2 + delta_y^2); %discrete ds
  314.         count = k;
  315.     end
  316.     %Sum all the segments for all modes to get a vector of a1..an
  317.     L = sum(delta_s);
  318.    
  319.     %for n = 1:mode_num
  320.     for k = 1:end_point-1
  321.         %Calculate s_mid, the position you are at along the MT
  322.         if k == 1
  323.             s_mid(k) = delta_s(k);
  324.         else
  325.             %s_mid(k) = sum(delta_s(k-1))+ .5.*delta_s(k);
  326.             s_mid(k) = s_mid(k-1)+0.5.*delta_s(k-1)+0.5.*delta_s(k);
  327.            
  328.         end
  329.         %                     a_k(k)= theta_s(k)*delta_s(k)*cos((n*pi()*s_mid(k))/L) %a_n vector for each k
  330.     end
  331.     %summation to get the final coeff based on all the segments, a_1 a_2 a_3, etc
  332.     %coeff = [coeff_mode_1 coeff_mode_2 coeff_mode_3...] for
  333.     %one image
  334.    
  335.     for n = 1:mode_num
  336.         for k = 1:end_point-1
  337.             a_k = a_k + theta_s(k)*delta_s(k)*exp(i*((n*pi()*s_mid(k))/L));
  338.         end
  339.         a_n(n) = sqrt(2/L)*a_k;
  340.         a_k = 0;
  341.     end
  342.     %all the modes for all the frames
  343.     %coeff_frame = each row corresponds to each frame, each column
  344.     %corresponds to the coeff for each mode (col 1 is mode 1, col 2
  345.     %is mode 2..)
  346.     coeff_frame(frame,:) = a_n;
  347.     coeff_frame_stats_on_rows(:,frame) = a_n;
  348.     frame = frame + 1;
  349.     magnitude_an(frame,:) = abs(a_n);
  350.    
  351. end
  352.  
  353. figure
  354.  
  355. %Plotting Results
  356.  
  357.  
  358. for plot_mode = 1:mode_num
  359.     subplot(1,mode_num,plot_mode); %divides the current figure into an 1-by-n grid and creates axes in the position specified by the fram
  360.     plot(1:size(magnitude_an,1),magnitude_an(:,plot_mode),'.')
  361.     % linkaxes(mode_plot,'y')
  362.     ylim([-3 3]);
  363. end
  364.  
  365. %plotting only reals
  366. for i = 1:mode_num
  367.     subplot(1,mode_num,i); %divides the current figure into an 1-by-n grid and creates axes in the position specified by the fram
  368.     plot(1:size(coeff_frame,1),coeff_frame(:,i),'.')
  369.     % linkaxes(mode_plot,'y')
  370.     ylim([-3 3]);
  371. end
  372.  
  373. %Plot variance of the amplitudes
  374. %Get the squared  difference between the same coeff between each frame and the average of the coeff over all frames, then sum over all frames, and divide by Frames - 1
  375. var_coeff = var(magnitude_an,0,1);
  376. figure (4)
  377. plot(1:n,var_coeff)
RAW Paste Data
Top