Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Added various unorganized matlab files to the repo.

  • Loading branch information...
commit 03308c4e0d160a051fce43e4ef7662c878c8c548 0 parents
ajdecon authored
Showing with 7,295 additions and 0 deletions.
  1. +646 −0 CircularHough_Grd.m
  2. +30 −0 DrawCircle.m
  3. +37 −0 HoughList.m
  4. +20 −0 HoughMovie.m
  5. +28 −0 HoughTracking.m
  6. +28 −0 HoughTracking2.m
  7. +177 −0 Watershed3D_Final.m
  8. +33 −0 assign_fields.m
  9. +32 −0 average_intensity.m
  10. +15 −0 avg_array.m
  11. +10 −0 avg_array_inc.m
  12. +21 −0 avg_frames.m
  13. +25 −0 bersen.m
  14. +33 −0 bondlengths.m
  15. +63 −0 bounded_array.m
  16. +56 −0 boxtracking.m
  17. +155 −0 bpass.m
  18. +82 −0 canny_edge.m
  19. +146 −0 cefit_4C.m
  20. +20 −0 cell2num.m
  21. +22 −0 centerpos.m
  22. +36 −0 check_dist.m
  23. +143 −0 cntrd.m
  24. +33 −0 compare_labels.m
  25. +34 −0 compare_stacks.m
  26. +34 −0 compare_stacks_labels.m
  27. +20 −0 connect_circles.m
  28. +55 −0 convert_tifs.m
  29. +67 −0 create_help_menu.m
  30. +16 −0 d2dgauss.m
  31. +35 −0 demo_thresh_box.m
  32. +12 −0 devices_volfrac.m
  33. +2 −0  dgauss.m
  34. +141 −0 dna_boxtracking.m
  35. +82 −0 dna_pairing.m
  36. +77 −0 dna_particle_tracking.m
  37. +82 −0 dpt2.m
  38. +25 −0 drawline.m
  39. +78 −0 dt_percent_threshold.m
  40. +163 −0 dumbbell_pos.m
  41. +22 −0 extract_circles_1img.m
  42. +24 −0 fast_tracking.m
  43. +27 −0 fixed_ints.m
  44. +3 −0  gauss.m
  45. +13 −0 gen_gauss_signal.m
  46. +129 −0 group_backbone.m
  47. +17 −0 histstack.m
  48. +48 −0 houghtrackmovie.m
  49. +22 −0 info2label.m
  50. +35 −0 intensity_at_points.m
  51. +32 −0 label_and_time_stamp.m
  52. +17 −0 lentrk.m
  53. +6 −0 loadrgbtiffstack.m
  54. +6 −0 loadtiffstack.m
  55. +25 −0 local_otsu.m
  56. +83 −0 local_percent_threshold.m
  57. +31 −0 local_threshold.m
  58. +91 −0 lpt3d.m
  59. +12 −0 makesplinederiv.m
  60. +96 −0 maketables.m
  61. +24 −0 mean_thresh.m
  62. +47 −0 meanplot.m
  63. +24 −0 median_thresh.m
  64. +28 −0 msd.m
  65. +163 −0 msd_a.m
  66. +180 −0 msd_d.m
  67. +12 −0 myfunc.m
  68. +40 −0 myview.m
  69. +22 −0 niblack.m
  70. +19 −0 norm_ratio.m
  71. +40 −0 orient_batch.m
  72. +47 −0 orient_corr.m
  73. +27 −0 parse_noran.m
  74. +9 −0 particle_averages.m
  75. +6 −0 particle_split.m
  76. +16 −0 perc_thresh.m
  77. +102 −0 pkfnd.m
  78. +44 −0 pretrack.m
  79. +45 −0 pretrack2.m
  80. +39 −0 random_particles.m
  81. +57 −0 read_gdf.m
  82. +41 −0 read_noran.m
  83. +85 −0 ring_section.m
  84. +26 −0 sauvola.m
  85. +31 −0 showmovie.m
  86. +195 −0 slice_view_2d.m
  87. +14 −0 slopes.m
  88. +19 −0 sortvtime.m
  89. +34 −0 stack_intensity.m
  90. +16 −0 stackmovie.m
  91. +15 −0 stackmovie_labels.m
  92. +9 −0 temp.m
  93. +22 −0 test.m
  94. +23 −0 test_function.m
  95. +23 −0 test_function2.m
  96. +23 −0 test_function3.m
  97. +1,112 −0 track.m
  98. +116 −0 track_intensities.m
  99. +29 −0 units2tex.m
  100. +38 −0 vector2str.m
  101. +200 −0 vol3d.m
  102. +136 −0 vol3dtool.m
  103. +225 −0 volume_browser.m
  104. +50 −0 watershed_pos.m
  105. +113 −0 watershed_tracking.m
  106. +77 −0 watertracking2.m
  107. +27 −0 welch_t_test.m
  108. +26 −0 x_section.m
  109. +26 −0 y_section.m
646 CircularHough_Grd.m
@@ -0,0 +1,646 @@
+function [accum, varargout] = CircularHough_Grd(img, radrange, varargin)
+%Detect circular shapes in a grayscale image. Resolve their center
+%positions and radii.
+%
+% [accum, circen, cirrad, dbg_LMmask] = CircularHough_Grd(
+% img, radrange, grdthres, fltr4LM_R, multirad, fltr4accum)
+% Circular Hough transform based on the gradient field of an image.
+% NOTE: Operates on grayscale images, NOT B/W bitmaps.
+% NO loops in the implementation of Circular Hough transform,
+% which means faster operation but at the same time larger
+% memory consumption.
+%
+%%%%%%%% INPUT: (img, radrange, grdthres, fltr4LM_R, multirad, fltr4accum)
+%
+% img: A 2-D grayscale image (NO B/W bitmap)
+%
+% radrange: The possible minimum and maximum radii of the circles
+% to be searched, in the format of
+% [minimum_radius , maximum_radius] (unit: pixels)
+% **NOTE**: A smaller range saves computational time and
+% memory.
+%
+% grdthres: (Optional, default is 10, must be non-negative)
+% The algorithm is based on the gradient field of the
+% input image. A thresholding on the gradient magnitude
+% is performed before the voting process of the Circular
+% Hough transform to remove the 'uniform intensity'
+% (sort-of) image background from the voting process.
+% In other words, pixels with gradient magnitudes smaller
+% than 'grdthres' are NOT considered in the computation.
+% **NOTE**: The default parameter value is chosen for
+% images with a maximum intensity close to 255. For cases
+% with dramatically different maximum intensities, e.g.
+% 10-bit bitmaps in stead of the assumed 8-bit, the default
+% value can NOT be used. A value of 4% to 10% of the maximum
+% intensity may work for general cases.
+%
+% fltr4LM_R: (Optional, default is 8, minimum is 3)
+% The radius of the filter used in the search of local
+% maxima in the accumulation array. To detect circles whose
+% shapes are less perfect, the radius of the filter needs
+% to be set larger.
+%
+% multirad: (Optional, default is 0.5)
+% In case of concentric circles, multiple radii may be
+% detected corresponding to a single center position. This
+% argument sets the tolerance of picking up the likely
+% radii values. It ranges from 0.1 to 1, where 0.1
+% corresponds to the largest tolerance, meaning more radii
+% values will be detected, and 1 corresponds to the smallest
+% tolerance, in which case only the "principal" radius will
+% be picked up.
+%
+% fltr4accum: (Optional. A default filter will be used if not given)
+% Filter used to smooth the accumulation array. Depending
+% on the image and the parameter settings, the accumulation
+% array built has different noise level and noise pattern
+% (e.g. noise frequencies). The filter should be set to an
+% appropriately size such that it's able to suppress the
+% dominant noise frequency.
+%
+%%%%%%%% OUTPUT: [accum, circen, cirrad, dbg_LMmask]
+%
+% accum: The result accumulation array from the Circular Hough
+% transform. The accumulation array has the same dimension
+% as the input image.
+%
+% circen: (Optional)
+% Center positions of the circles detected. Is a N-by-2
+% matrix with each row contains the (x, y) positions
+% of a circle. For concentric circles (with the same center
+% position), say k of them, the same center position will
+% appear k times in the matrix.
+%
+% cirrad: (Optional)
+% Estimated radii of the circles detected. Is a N-by-1
+% column vector with a one-to-one correspondance to the
+% output 'circen'. A value 0 for the radius indicates a
+% failed detection of the circle's radius.
+%
+% dbg_LMmask: (Optional, for debugging purpose)
+% Mask from the search of local maxima in the accumulation
+% array.
+%
+%%%%%%%%% EXAMPLE #0:
+% rawimg = imread('TestImg_CHT_a2.bmp');
+% tic;
+% [accum, circen, cirrad] = CircularHough_Grd(rawimg, [15 60]);
+% toc;
+% figure(1); imagesc(accum); axis image;
+% title('Accumulation Array from Circular Hough Transform');
+% figure(2); imagesc(rawimg); colormap('gray'); axis image;
+% hold on;
+% plot(circen(:,1), circen(:,2), 'r+');
+% for k = 1 : size(circen, 1),
+% DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, 'b-');
+% end
+% hold off;
+% title(['Raw Image with Circles Detected ', ...
+% '(center positions and radii marked)']);
+% figure(3); surf(accum, 'EdgeColor', 'none'); axis ij;
+% title('3-D View of the Accumulation Array');
+%
+% COMMENTS ON EXAMPLE #0:
+% Kind of an easy case to handle. To detect circles in the image whose
+% radii range from 15 to 60. Default values for arguments 'grdthres',
+% 'fltr4LM_R', 'multirad' and 'fltr4accum' are used.
+%
+%%%%%%%%% EXAMPLE #1:
+% rawimg = imread('TestImg_CHT_a3.bmp');
+% tic;
+% [accum, circen, cirrad] = CircularHough_Grd(rawimg, [15 60], 10, 20);
+% toc;
+% figure(1); imagesc(accum); axis image;
+% title('Accumulation Array from Circular Hough Transform');
+% figure(2); imagesc(rawimg); colormap('gray'); axis image;
+% hold on;
+% plot(circen(:,1), circen(:,2), 'r+');
+% for k = 1 : size(circen, 1),
+% DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, 'b-');
+% end
+% hold off;
+% title(['Raw Image with Circles Detected ', ...
+% '(center positions and radii marked)']);
+% figure(3); surf(accum, 'EdgeColor', 'none'); axis ij;
+% title('3-D View of the Accumulation Array');
+%
+% COMMENTS ON EXAMPLE #1:
+% The shapes in the raw image are not very good circles. As a result,
+% the profile of the peaks in the accumulation array are kind of
+% 'stumpy', which can be seen clearly from the 3-D view of the
+% accumulation array. (As a comparison, please see the sharp peaks in
+% the accumulation array in example #0) To extract the peak positions
+% nicely, a value of 20 (default is 8) is used for argument 'fltr4LM_R',
+% which is the radius of the filter used in the search of peaks.
+%
+%%%%%%%%% EXAMPLE #2:
+% rawimg = imread('TestImg_CHT_b3.bmp');
+% fltr4img = [1 1 1 1 1; 1 2 2 2 1; 1 2 4 2 1; 1 2 2 2 1; 1 1 1 1 1];
+% fltr4img = fltr4img / sum(fltr4img(:));
+% imgfltrd = filter2( fltr4img , rawimg );
+% tic;
+% [accum, circen, cirrad] = CircularHough_Grd(imgfltrd, [15 80], 8, 10);
+% toc;
+% figure(1); imagesc(accum); axis image;
+% title('Accumulation Array from Circular Hough Transform');
+% figure(2); imagesc(rawimg); colormap('gray'); axis image;
+% hold on;
+% plot(circen(:,1), circen(:,2), 'r+');
+% for k = 1 : size(circen, 1),
+% DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, 'b-');
+% end
+% hold off;
+% title(['Raw Image with Circles Detected ', ...
+% '(center positions and radii marked)']);
+%
+% COMMENTS ON EXAMPLE #2:
+% The circles in the raw image have small scale irregularities along
+% the edges, which could lead to an accumulation array that is bad for
+% local maxima detection. A 5-by-5 filter is used to smooth out the
+% small scale irregularities. A blurred image is actually good for the
+% algorithm implemented here which is based on the image's gradient
+% field.
+%
+%%%%%%%%% EXAMPLE #3:
+% rawimg = imread('TestImg_CHT_c3.bmp');
+% fltr4img = [1 1 1 1 1; 1 2 2 2 1; 1 2 4 2 1; 1 2 2 2 1; 1 1 1 1 1];
+% fltr4img = fltr4img / sum(fltr4img(:));
+% imgfltrd = filter2( fltr4img , rawimg );
+% tic;
+% [accum, circen, cirrad] = ...
+% CircularHough_Grd(imgfltrd, [15 105], 8, 10, 0.7);
+% toc;
+% figure(1); imagesc(accum); axis image;
+% figure(2); imagesc(rawimg); colormap('gray'); axis image;
+% hold on;
+% plot(circen(:,1), circen(:,2), 'r+');
+% for k = 1 : size(circen, 1),
+% DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, 'b-');
+% end
+% hold off;
+% title(['Raw Image with Circles Detected ', ...
+% '(center positions and radii marked)']);
+%
+% COMMENTS ON EXAMPLE #3:
+% Similar to example #2, a filtering before circle detection works for
+% noisy image too. 'multirad' is set to 0.7 to eliminate the false
+% detections of the circles' radii.
+%
+%%%%%%%%% BUG REPORT:
+% This is a beta version. Please send your bug reports, comments and
+% suggestions to pengtao@glue.umd.edu . Thanks.
+%
+%
+%%%%%%%%% INTERNAL PARAMETERS:
+% The INPUT arguments are just part of the parameters that are used by
+% the circle detection algorithm implemented here. Variables in the code
+% with a prefix 'prm_' in the name are the parameters that control the
+% judging criteria and the behavior of the algorithm. Default values for
+% these parameters can hardly work for all circumstances. Therefore, at
+% occasions, the values of these INTERNAL PARAMETERS (parameters that
+% are NOT exposed as input arguments) need to be fine-tuned to make
+% the circle detection work as expected.
+% The following example shows how changing an internal parameter could
+% influence the detection result.
+% 1. Change the value of the internal parameter 'prm_LM_LoBndRa' to 0.4
+% (default is 0.2)
+% 2. Run the following matlab code:
+% fltr4accum = [1 2 1; 2 6 2; 1 2 1];
+% fltr4accum = fltr4accum / sum(fltr4accum(:));
+% rawimg = imread('Frame_0_0022_portion.jpg');
+% tic;
+% [accum, circen] = CircularHough_Grd(rawimg, ...
+% [4 14], 10, 4, 0.5, fltr4accum);
+% toc;
+% figure(1); imagesc(accum); axis image;
+% title('Accumulation Array from Circular Hough Transform');
+% figure(2); imagesc(rawimg); colormap('gray'); axis image;
+% hold on; plot(circen(:,1), circen(:,2), 'r+'); hold off;
+% title('Raw Image with Circles Detected (center positions marked)');
+% 3. See how different values of the parameter 'prm_LM_LoBndRa' could
+% influence the result.
+
+% Author: Tao Peng
+% Department of Mechanical Engineering
+% University of Maryland, College Park, Maryland 20742, USA
+% pengtao@glue.umd.edu
+% Version: Beta Revision: Mar. 07, 2007
+
+
+%%%%%%%% Arguments and parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Validation of arguments
+if ndims(img) ~= 2 || ~isnumeric(img),
+ error('CircularHough_Grd: ''img'' has to be 2 dimensional');
+end
+if ~all(size(img) >= 32),
+ error('CircularHough_Grd: ''img'' has to be larger than 32-by-32');
+end
+
+if numel(radrange) ~= 2 || ~isnumeric(radrange),
+ error(['CircularHough_Grd: ''radrange'' has to be ', ...
+ 'a two-element vector']);
+end
+prm_r_range = sort(max( [0,0;radrange(1),radrange(2)] ));
+
+% Parameters (default values)
+prm_grdthres = 10;
+prm_fltrLM_R = 8;
+prm_multirad = 0.5;
+func_compu_cen = true;
+func_compu_radii = true;
+
+% Validation of arguments
+vap_grdthres = 1;
+if nargin > (1 + vap_grdthres),
+ if isnumeric(varargin{vap_grdthres}) && ...
+ varargin{vap_grdthres}(1) >= 0,
+ prm_grdthres = varargin{vap_grdthres}(1);
+ else
+ error(['CircularHough_Grd: ''grdthres'' has to be ', ...
+ 'a non-negative number']);
+ end
+end
+
+vap_fltr4LM = 2; % filter for the search of local maxima
+if nargin > (1 + vap_fltr4LM),
+ if isnumeric(varargin{vap_fltr4LM}) && varargin{vap_fltr4LM}(1) >= 3,
+ prm_fltrLM_R = varargin{vap_fltr4LM}(1);
+ else
+ error(['CircularHough_Grd: ''fltr4LM_R'' has to be ', ...
+ 'larger than or equal to 3']);
+ end
+end
+
+vap_multirad = 3;
+if nargin > (1 + vap_multirad),
+ if isnumeric(varargin{vap_multirad}) && ...
+ varargin{vap_multirad}(1) >= 0.1 && ...
+ varargin{vap_multirad}(1) <= 1,
+ prm_multirad = varargin{vap_multirad}(1);
+ else
+ error(['CircularHough_Grd: ''multirad'' has to be ', ...
+ 'within the range [0.1, 1]']);
+ end
+end
+
+vap_fltr4accum = 4; % filter for smoothing the accumulation array
+if nargin > (1 + vap_fltr4accum),
+ if isnumeric(varargin{vap_fltr4accum}) && ...
+ ndims(varargin{vap_fltr4accum}) == 2 && ...
+ all(size(varargin{vap_fltr4accum}) >= 3),
+ fltr4accum = varargin{vap_fltr4accum};
+ else
+ error(['CircularHough_Grd: ''fltr4accum'' has to be ', ...
+ 'a 2-D matrix with a minimum size of 3-by-3']);
+ end
+else
+ % Default filter (5-by-5)
+ fltr4accum = ones(5,5);
+ fltr4accum(2:4,2:4) = 2;
+ fltr4accum(3,3) = 6;
+end
+
+func_compu_cen = ( nargout > 1 );
+func_compu_radii = ( nargout > 2 );
+
+% Reserved parameters
+dbg_on = false; % debug information
+dbg_bfigno = 4;
+if nargout > 3, dbg_on = true; end
+
+
+%%%%%%%% Building accumulation array %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Convert the image to single if it is not of
+% class float (single or double)
+img_is_double = isa(img, 'double');
+if ~(img_is_double || isa(img, 'single')),
+ imgf = single(img);
+end
+
+% Compute the gradient and the magnitude of gradient
+if img_is_double,
+ [grdx, grdy] = gradient(img);
+else
+ [grdx, grdy] = gradient(imgf);
+end
+grdmag = sqrt(grdx.^2 + grdy.^2);
+
+% Get the linear indices, as well as the subscripts, of the pixels
+% whose gradient magnitudes are larger than the given threshold
+grdmasklin = find(grdmag > prm_grdthres);
+[grdmask_IdxI, grdmask_IdxJ] = ind2sub(size(grdmag), grdmasklin);
+
+% Compute the linear indices (as well as the subscripts) of
+% all the votings to the accumulation array.
+% The Matlab function 'accumarray' accepts only double variable,
+% so all indices are forced into double at this point.
+% A row in matrix 'lin2accum_aJ' contains the J indices (into the
+% accumulation array) of all the votings that are introduced by a
+% same pixel in the image. Similarly with matrix 'lin2accum_aI'.
+rr_4linaccum = double( prm_r_range );
+linaccum_dr = [ (-rr_4linaccum(2) + 0.5) : -rr_4linaccum(1) , ...
+ (rr_4linaccum(1) + 0.5) : rr_4linaccum(2) ];
+
+lin2accum_aJ = floor( ...
+ double(grdx(grdmasklin)./grdmag(grdmasklin)) * linaccum_dr + ...
+ repmat( double(grdmask_IdxJ)+0.5 , [1,length(linaccum_dr)] ) ...
+);
+lin2accum_aI = floor( ...
+ double(grdy(grdmasklin)./grdmag(grdmasklin)) * linaccum_dr + ...
+ repmat( double(grdmask_IdxI)+0.5 , [1,length(linaccum_dr)] ) ...
+);
+
+% Clip the votings that are out of the accumulation array
+mask_valid_aJaI = ...
+ lin2accum_aJ > 0 & lin2accum_aJ < (size(grdmag,2) + 1) & ...
+ lin2accum_aI > 0 & lin2accum_aI < (size(grdmag,1) + 1);
+
+mask_valid_aJaI_reverse = ~ mask_valid_aJaI;
+lin2accum_aJ = lin2accum_aJ .* mask_valid_aJaI + mask_valid_aJaI_reverse;
+lin2accum_aI = lin2accum_aI .* mask_valid_aJaI + mask_valid_aJaI_reverse;
+clear mask_valid_aJaI_reverse;
+
+% Linear indices (of the votings) into the accumulation array
+lin2accum = sub2ind( size(grdmag), lin2accum_aI, lin2accum_aJ );
+
+lin2accum_size = size( lin2accum );
+lin2accum = reshape( lin2accum, [numel(lin2accum),1] );
+clear lin2accum_aI lin2accum_aJ;
+
+% Weights of the votings, currently using the gradient maginitudes
+% but in fact any scheme can be used (application dependent)
+weight4accum = ...
+ repmat( double(grdmag(grdmasklin)) , [lin2accum_size(2),1] ) .* ...
+ mask_valid_aJaI(:);
+clear mask_valid_aJaI;
+
+% Build the accumulation array using Matlab function 'accumarray'
+accum = accumarray( lin2accum , weight4accum );
+accum = [ accum ; zeros( numel(grdmag) - numel(accum) , 1 ) ];
+accum = reshape( accum, size(grdmag) );
+
+
+%%%%%%%% Locating local maxima in the accumulation array %%%%%%%%%%%%
+
+% Stop if no need to locate the center positions of circles
+if ~func_compu_cen,
+ return;
+end
+clear lin2accum weight4accum;
+
+% Parameters to locate the local maxima in the accumulation array
+% -- Segmentation of 'accum' before locating LM
+prm_useaoi = true;
+prm_aoithres_s = 2;
+prm_aoiminsize = floor(min([ min(size(accum)) * 0.25, ...
+ prm_r_range(2) * 1.5 ]));
+
+% -- Filter for searching for local maxima
+prm_fltrLM_s = 1.35;
+prm_fltrLM_r = ceil( prm_fltrLM_R * 0.6 );
+prm_fltrLM_npix = max([ 6, ceil((prm_fltrLM_R/2)^1.8) ]);
+
+% -- Lower bound of the intensity of local maxima
+prm_LM_LoBndRa = 0.2; % minimum ratio of LM to the max of 'accum'
+
+% Smooth the accumulation array
+fltr4accum = fltr4accum / sum(fltr4accum(:));
+accum = filter2( fltr4accum, accum );
+
+% Select a number of Areas-Of-Interest from the accumulation array
+if prm_useaoi,
+ % Threshold value for 'accum'
+ prm_llm_thres1 = prm_grdthres * prm_aoithres_s;
+
+ % Thresholding over the accumulation array
+ accummask = ( accum > prm_llm_thres1 );
+
+ % Segmentation over the mask
+ [accumlabel, accum_nRgn] = bwlabel( accummask, 8 );
+
+ % Select AOIs from segmented regions
+ accumAOI = ones(0,4);
+ for k = 1 : accum_nRgn,
+ accumrgn_lin = find( accumlabel == k );
+ [accumrgn_IdxI, accumrgn_IdxJ] = ...
+ ind2sub( size(accumlabel), accumrgn_lin );
+ rgn_top = min( accumrgn_IdxI );
+ rgn_bottom = max( accumrgn_IdxI );
+ rgn_left = min( accumrgn_IdxJ );
+ rgn_right = max( accumrgn_IdxJ );
+ % The AOIs selected must satisfy a minimum size
+ if ( (rgn_right - rgn_left + 1) >= prm_aoiminsize && ...
+ (rgn_bottom - rgn_top + 1) >= prm_aoiminsize ),
+ accumAOI = [ accumAOI; ...
+ rgn_top, rgn_bottom, rgn_left, rgn_right ];
+ end
+ end
+else
+ % Whole accumulation array as the one AOI
+ accumAOI = [1, size(accum,1), 1, size(accum,2)];
+end
+
+% Thresholding of 'accum' by a lower bound
+prm_LM_LoBnd = max(accum(:)) * prm_LM_LoBndRa;
+
+% Build the filter for searching for local maxima
+fltr4LM = zeros(2 * prm_fltrLM_R + 1);
+
+[mesh4fLM_x, mesh4fLM_y] = meshgrid(-prm_fltrLM_R : prm_fltrLM_R);
+mesh4fLM_r = sqrt( mesh4fLM_x.^2 + mesh4fLM_y.^2 );
+fltr4LM_mask = ...
+ ( mesh4fLM_r > prm_fltrLM_r & mesh4fLM_r <= prm_fltrLM_R );
+fltr4LM = fltr4LM - ...
+ fltr4LM_mask * (prm_fltrLM_s / sum(fltr4LM_mask(:)));
+
+if prm_fltrLM_R >= 4,
+ fltr4LM_mask = ( mesh4fLM_r < (prm_fltrLM_r - 1) );
+else
+ fltr4LM_mask = ( mesh4fLM_r < prm_fltrLM_r );
+end
+fltr4LM = fltr4LM + fltr4LM_mask / sum(fltr4LM_mask(:));
+
+% **** Debug code (begin)
+if dbg_on,
+ dbg_LMmask = zeros(size(accum));
+end
+% **** Debug code (end)
+
+% For each of the AOIs selected, locate the local maxima
+circen = zeros(0,2);
+for k = 1 : size(accumAOI, 1),
+ aoi = accumAOI(k,:); % just for referencing convenience
+
+ % Thresholding of 'accum' by a lower bound
+ accumaoi_LBMask = ...
+ ( accum(aoi(1):aoi(2), aoi(3):aoi(4)) > prm_LM_LoBnd );
+
+ % Apply the local maxima filter
+ candLM = conv2( accum(aoi(1):aoi(2), aoi(3):aoi(4)) , ...
+ fltr4LM , 'same' );
+ candLM_mask = ( candLM > 0 );
+
+ % Clear the margins of 'candLM_mask'
+ candLM_mask([1:prm_fltrLM_R, (end-prm_fltrLM_R+1):end], :) = 0;
+ candLM_mask(:, [1:prm_fltrLM_R, (end-prm_fltrLM_R+1):end]) = 0;
+
+ % **** Debug code (begin)
+ if dbg_on,
+ dbg_LMmask(aoi(1):aoi(2), aoi(3):aoi(4)) = ...
+ dbg_LMmask(aoi(1):aoi(2), aoi(3):aoi(4)) + ...
+ accumaoi_LBMask + 2 * candLM_mask;
+ end
+ % **** Debug code (end)
+
+ % Group the local maxima candidates by adjacency, compute the
+ % centroid position for each group and take that as the center
+ % of one circle detected
+ [candLM_label, candLM_nRgn] = bwlabel( candLM_mask, 8 );
+
+ for ilabel = 1 : candLM_nRgn,
+ % Indices (to current AOI) of the pixels in the group
+ candgrp_masklin = find( candLM_label == ilabel );
+ [candgrp_IdxI, candgrp_IdxJ] = ...
+ ind2sub( size(candLM_label) , candgrp_masklin );
+
+ % Indices (to 'accum') of the pixels in the group
+ candgrp_IdxI = candgrp_IdxI + ( aoi(1) - 1 );
+ candgrp_IdxJ = candgrp_IdxJ + ( aoi(3) - 1 );
+ candgrp_idx2acm = ...
+ sub2ind( size(accum) , candgrp_IdxI , candgrp_IdxJ );
+
+ % Minimum number of qulified pixels in the group
+ if sum(accumaoi_LBMask(candgrp_masklin)) < prm_fltrLM_npix,
+ continue;
+ end
+
+ % Compute the centroid position
+ candgrp_acmsum = sum( accum(candgrp_idx2acm) );
+ cc_x = sum( candgrp_IdxJ .* accum(candgrp_idx2acm) ) / ...
+ candgrp_acmsum;
+ cc_y = sum( candgrp_IdxI .* accum(candgrp_idx2acm) ) / ...
+ candgrp_acmsum;
+ circen = [circen; cc_x, cc_y];
+ end
+end
+
+% **** Debug code (begin)
+if dbg_on,
+ figure(dbg_bfigno); imagesc(dbg_LMmask); axis image;
+ title('Generated map of local maxima');
+ if size(accumAOI, 1) == 1,
+ figure(dbg_bfigno+1);
+ surf(candLM, 'EdgeColor', 'none'); axis ij;
+ title('Accumulation array after local maximum filtering');
+ end
+end
+% **** Debug code (end)
+
+
+%%%%%%%% Estimation of the Radii of Circles %%%%%%%%%%%%
+
+% Stop if no need to estimate the radii of circles
+if ~func_compu_radii,
+ varargout{1} = circen;
+ return;
+end
+
+% Parameters for the estimation of the radii of circles
+fltr4SgnCv = [2 1 1];
+fltr4SgnCv = fltr4SgnCv / sum(fltr4SgnCv);
+
+% Find circle's radius using its signature curve
+cirrad = zeros( size(circen,1), 1 );
+
+for k = 1 : size(circen,1),
+ % Neighborhood region of the circle for building the sgn. curve
+ circen_round = round( circen(k,:) );
+ SCvR_I0 = circen_round(2) - prm_r_range(2) - 1;
+ if SCvR_I0 < 1,
+ SCvR_I0 = 1;
+ end
+ SCvR_I1 = circen_round(2) + prm_r_range(2) + 1;
+ if SCvR_I1 > size(grdx,1),
+ SCvR_I1 = size(grdx,1);
+ end
+ SCvR_J0 = circen_round(1) - prm_r_range(2) - 1;
+ if SCvR_J0 < 1,
+ SCvR_J0 = 1;
+ end
+ SCvR_J1 = circen_round(1) + prm_r_range(2) + 1;
+ if SCvR_J1 > size(grdx,2),
+ SCvR_J1 = size(grdx,2);
+ end
+
+ % Build the sgn. curve
+ SgnCvMat_dx = repmat( (SCvR_J0:SCvR_J1) - circen(k,1) , ...
+ [SCvR_I1 - SCvR_I0 + 1 , 1] );
+ SgnCvMat_dy = repmat( (SCvR_I0:SCvR_I1)' - circen(k,2) , ...
+ [1 , SCvR_J1 - SCvR_J0 + 1] );
+ SgnCvMat_r = sqrt( SgnCvMat_dx .^2 + SgnCvMat_dy .^2 );
+ SgnCvMat_rp1 = round(SgnCvMat_r) + 1;
+
+ f4SgnCv = abs( ...
+ double(grdx(SCvR_I0:SCvR_I1, SCvR_J0:SCvR_J1)) .* SgnCvMat_dx + ...
+ double(grdy(SCvR_I0:SCvR_I1, SCvR_J0:SCvR_J1)) .* SgnCvMat_dy ...
+ ) ./ SgnCvMat_r;
+ SgnCv = accumarray( SgnCvMat_rp1(:) , f4SgnCv(:) );
+
+ SgnCv_Cnt = accumarray( SgnCvMat_rp1(:) , ones(numel(f4SgnCv),1) );
+ SgnCv_Cnt = SgnCv_Cnt + (SgnCv_Cnt == 0);
+ SgnCv = SgnCv ./ SgnCv_Cnt;
+
+ % Suppress the undesired entries in the sgn. curve
+ % -- Radii that correspond to short arcs
+ SgnCv = SgnCv .* ( SgnCv_Cnt >= (pi/4 * [0:(numel(SgnCv_Cnt)-1)]') );
+ % -- Radii that are out of the given range
+ SgnCv( 1 : (round(prm_r_range(1))+1) ) = 0;
+ SgnCv( (round(prm_r_range(2))+1) : end ) = 0;
+
+ % Get rid of the zero radius entry in the array
+ SgnCv = SgnCv(2:end);
+ % Smooth the sgn. curve
+ SgnCv = filtfilt( fltr4SgnCv , [1] , SgnCv );
+
+ % Get the maximum value in the sgn. curve
+ SgnCv_max = max(SgnCv);
+ if SgnCv_max <= 0,
+ cirrad(k) = 0;
+ continue;
+ end
+
+ % Find the local maxima in sgn. curve by 1st order derivatives
+ % -- Mark the ascending edges in the sgn. curve as 1s and
+ % -- descending edges as 0s
+ SgnCv_AscEdg = ( SgnCv(2:end) - SgnCv(1:(end-1)) ) > 0;
+ % -- Mark the transition (ascending to descending) regions
+ SgnCv_LMmask = [ 0; 0; SgnCv_AscEdg(1:(end-2)) ] & (~SgnCv_AscEdg);
+ SgnCv_LMmask = SgnCv_LMmask & [ SgnCv_LMmask(2:end) ; 0 ];
+
+ % Incorporate the minimum value requirement
+ SgnCv_LMmask = SgnCv_LMmask & ...
+ ( SgnCv(1:(end-1)) >= (prm_multirad * SgnCv_max) );
+ % Get the positions of the peaks
+ SgnCv_LMPos = sort( find(SgnCv_LMmask) );
+
+ % Save the detected radii
+ if isempty(SgnCv_LMPos),
+ cirrad(k) = 0;
+ else
+ cirrad(k) = SgnCv_LMPos(end);
+ for i_radii = (length(SgnCv_LMPos) - 1) : -1 : 1,
+ circen = [ circen; circen(k,:) ];
+ cirrad = [ cirrad; SgnCv_LMPos(i_radii) ];
+ end
+ end
+end
+
+% Output
+varargout{1} = circen;
+varargout{2} = cirrad;
+if nargout > 3,
+ varargout{3} = dbg_LMmask;
+end
30 DrawCircle.m
@@ -0,0 +1,30 @@
+function DrawCircle(x, y, r, nseg, S)
+% Draw a circle on the current figure using ploylines
+%
+% DrawCircle(x, y, r, nseg, S)
+% A simple function for drawing a circle on graph.
+%
+% INPUT: (x, y, r, nseg, S)
+% x, y: Center of the circle
+% r: Radius of the circle
+% nseg: Number of segments for the circle
+% S: Colors, plot symbols and line types
+%
+% OUTPUT: None
+%
+% BUG REPORT:
+% Please send your bug reports, comments and suggestions to
+% pengtao@glue.umd.edu . Thanks.
+
+% Author: Tao Peng
+% Department of Mechanical Engineering
+% University of Maryland, College Park, Maryland 20742, USA
+% pengtao@glue.umd.edu
+% Version: alpha Revision: Jan. 10, 2006
+
+
+theta = 0 : (2 * pi / nseg) : (2 * pi);
+pline_x = r * cos(theta) + x;
+pline_y = r * sin(theta) + y;
+
+plot(pline_x, pline_y, S);
37 HoughList.m
@@ -0,0 +1,37 @@
+function [poslist,cirrad]=HoughList(img,radrange,plot)
+
+
+ [accum,circen,cirrad]=CircularHough_Grd(img(:,:),radrange);
+
+ % First, remove duplicates.
+ k=0;
+ rem_ind = [];
+ for i=1:size(circen,1)
+ x0=circen(i,1);
+ y0=circen(i,2);
+ for j=i:size(circen,1)
+ x1=circen(j,1);
+ y1=circen(j,2);
+ dist=sqrt((x1-x0)^2 + (y1-y0)^2);
+ if ((dist<min(radrange)) && (i~=j))
+ k=k+1;
+ rem_ind(k)=j;
+ end
+ end
+ end
+ circen(rem_ind,:)=[];
+ cirrad(rem_ind)=[];
+
+ if (plot>0)
+ imagesc(img);
+ colormap('gray');axis image;hold on;
+ for k=1:size(circen,1),DrawCircle(circen(k,1),circen(k,2),cirrad(k),32,'b-');end;
+ for i=1:size(circen,1)
+ h=text(circen(i,1),circen(i,2),num2str(i));
+ set(h,'Color','g');
+ set(h,'FontWeight','bold');
+ end
+ hold off;
+ end
+
+poslist=circen;
20 HoughMovie.m
@@ -0,0 +1,20 @@
+function poslist=HoughList(img,radrange)
+
+numim=size(img,3);
+
+for m=1:numim
+ [accum,circen,cirrad]=CircularHough_Grd(img(:,:,m),radrange);
+ imagesc(img(:,:,m));
+ colormap('gray');axis image;hold on;plot(circen(:,1),circen(:,2),'r+');
+ for k=1:size(circen,1),DrawCircle(circen(k,1),circen(k,2),cirrad(k),32,'b-');end;hold off;
+
+ for x=1:size(circen,1)
+ for y=1:size(circen,2)
+ cirarray(x,y,m)=circen(x,y);
+ end
+ end
+
+ pause(0.05);
+end
+
+poslist=cirarray;
28 HoughTracking.m
@@ -0,0 +1,28 @@
+function [positions,tracklist]=HoughTracking(img,radrange,timerange,plot)
+
+interval=5;
+for t=timerange
+ img0=img(:,:,t);
+ poslist=HoughList(img0,radrange,plot);
+ pause(0.05);
+ for k=1:size(poslist,1)
+ for m=1:size(poslist,2)
+ positions(k,m,t)=poslist(k,m);
+ end
+ end
+ if (mod(t,interval)==0)
+ t
+ end
+end
+
+k=0;
+for i=1:size(positions,3)
+ for j=1:size(positions,1)
+ if (positions(j,1,i)~=0)
+ k=k+1;
+ tracklist(k,1)=positions(j,1,i);
+ tracklist(k,2)=positions(j,2,i);
+ tracklist(k,3)=i;
+ end
+ end
+end
28 HoughTracking2.m
@@ -0,0 +1,28 @@
+function [positions,tracklist]=HoughTracking2(img,radrange,timerange,plot)
+
+interval=5;
+for t=timerange
+ img0=img(:,:,t);
+ [poslist,cirrad]=HoughList(img0,radrange,plot);
+ pause(0.05);
+ for k=1:size(poslist,1)
+ for m=1:size(poslist,2)
+ positions(k,m,t)=poslist(k,m);
+ end
+ end
+ if (mod(t,interval)==0)
+ t
+ end
+end
+
+k=0;
+for i=1:size(positions,3)
+ for j=1:size(positions,1)
+ if (positions(j,1,i)~=0)
+ k=k+1;
+ tracklist(k,1)=positions(j,1,i);
+ tracklist(k,2)=positions(j,2,i);
+ tracklist(k,3)=i;
+ end
+ end
+end
177 Watershed3D_Final.m
@@ -0,0 +1,177 @@
+SKIZdepth=input('number of images in sequence: '); %61
+SKIZheight=input('height of images: '); %312
+SKIZwidth=input('width of images: '); %346
+THRESH=.2; %Threshold amount
+FudgeFactor=2;
+scrsz = get(0,'ScreenSize');
+stack=zeros([SKIZheight SKIZwidth SKIZdepth]);
+BWstack=zeros([SKIZheight SKIZwidth SKIZdepth]);
+
+
+
+for count=1:SKIZdepth
+ stack(:,:,count) = imread('BEADZ_INV.tif',count);%Probably not necessary
+ BWstack(:,:,count) = im2bw(stack(:,:,count),THRESH);
+end
+
+InvDst=bwdist(~BWstack);
+InvDst2=bwdist(BWstack);
+
+Dst=-InvDst;
+Dst(~BWstack) = -Inf;
+Wtr = watershed(Dst,26);
+WtrTrnsfrm = Wtr;
+
+%Find a way to automatically detect background pixellation!!!!!
+for w=1:SKIZwidth
+ for h=1:SKIZheight
+ for d=1:SKIZdepth
+ if Wtr(h,w,d)==FudgeFactor
+ WtrTrnsfrm(h,w,d)=0;
+ end
+ end
+ d=1;
+ end
+h=1;
+end
+
+
+%%
+%Makes Binary 3D matrix of Wtr, for Distance Transform
+%Then takes a distance transform of total 3D backgroun
+BWwtr=WtrTrnsfrm;
+
+figure('Position', [scrsz(4)/2 scrsz(3) scrsz(3)/2 scrsz(4)/2]),
+for d=1:SKIZdepth
+ imshow(label2rgb(uint8(WtrTrnsfrm(:,:,d)),'jet'));
+ pause(.1);
+end
+title ('Segmented Granules');
+
+figure,
+for d=1:SKIZdepth
+ BWwtr(:,:,d) = im2bw(WtrTrnsfrm(:,:,d));
+end
+
+BWwtrDst=bwdist(BWwtr); %used for background subtraction at end
+figure('Position', [scrsz(4)/2 scrsz(3) scrsz(3)/2 scrsz(4)/2]),
+imshow(label2rgb(uint8(BWwtr(:,:,d))))
+title ('Global Distance Transform')
+
+
+
+%%
+WtrSegment=WtrTrnsfrm;
+TotalGranules=max(max((max(Wtr))));
+
+GRANS=zeros([SKIZheight SKIZwidth SKIZdepth TotalGranules]);
+for Granule=1:TotalGranules
+%%
+%This will take only a single granule from the Watertransform, then make it
+%binary and take its distance transform, stored in GRANS
+for w=1:SKIZwidth
+ for h=1:SKIZheight
+ for d=1:SKIZdepth
+ if WtrTrnsfrm(h,w,d)==Granule
+ WtrSegment(h,w,d)=1;
+ end
+ if WtrTrnsfrm(h,w,d)~=Granule
+ WtrSegment(h,w,d)=0;
+ end
+ end
+ d=1;
+ end
+h=1;
+end
+%%
+GRANS(:,:,:,Granule)=bwdist(WtrSegment);
+end
+%% Subtracts individual distance transform from total distance transform
+Voronoi=GRANS;
+Voronoi2=GRANS;
+for i=1:TotalGranules
+ Voronoi(:,:,:,i) = GRANS(:,:,:,i) - BWwtrDst;
+ Voronoi2(:,:,:,i) = GRANS(:,:,:,i) - InvDst2;
+ Voronoi(:,:,:,i)=255-Voronoi(:,:,:,i) - InvDst*333;%makes lowest hill into highest hill, making granule center brightest
+ Voronoi2(:,:,:,i)=255-Voronoi2(:,:,:,i);
+ %Voronoi(:,:,:,i)=imextendedmax(Voronoi(:,:,:,i),26)*i;
+end
+
+figure('Position', [scrsz(4)/2 scrsz(1) scrsz(3)/2 scrsz(4)/2]),
+for k=1:SKIZdepth
+imshow(label2rgb((uint8(Voronoi(:,:,k,10))),'jet'))
+pause(.1);
+end
+title ('Local Void Region');
+PureVoronoi=zeros([SKIZheight SKIZwidth SKIZdepth]);
+
+maxVal = 0;
+for j=1:SKIZwidth %interrogates and eliminates non 250 pixels
+ for i=1:SKIZheight
+
+ for k=1:SKIZdepth
+ maxVal=0;
+ for count=1:TotalGranules
+ if Voronoi(i,j,k,count)> maxVal;
+ maxVal=Voronoi(i,j,k,count);
+ PureVoronoi(i,j,k)= count;
+ end
+ end
+
+ end
+ end
+end
+
+
+figure('Position', [scrsz(2)/2 scrsz(1) scrsz(3)/2 scrsz(4)/2]),
+for k=1:SKIZdepth
+imshow(label2rgb((uint8(PureVoronoi(:,:,k))),'jet'))
+pause(.1);
+end
+title ('Voronoi Void Region');
+
+%So far I have made the background transform, called BWwtrDst and labelled
+%each individual granules background and stored this transform in GRANS
+%which is a 4D array, where the last column is the granule ID from the
+%watershed transform at the beginning. Voronoi subtracts the Individual
+%granules background from the total background, determining the voronoi
+%map. To get the standard Vornoi volume, I need to find the pixels that are
+%closest (255) to the granule and eliminate everything else, then squeeze
+%into a 3D array from the 4D array and make an isosurface.
+%% Imaging in 3D the voronoi volumes
+image_num = 3;
+image(PureVoronoi(:,:,image_num));
+x=xlim;
+y=ylim;
+contourslice(PureVoronoi,[],[],image_num)
+axis ij
+xlim(x)
+ylim(y)
+daspect([1,1,1])
+colormap(jet)
+pause(1);
+%layered contour images
+phandles = contourslice(PureVoronoi,[],[],[1,2,3,4,5],8);
+view(3); axis tight
+set(phandles,'LineWidth',2)
+pause(1);
+%Isosurface
+Ls = smooth3(PureVoronoi);
+hiso = patch(isosurface(Ls,500),...
+ 'FaceColor',[1,.75,.65],...
+ 'EdgeColor','red');
+hcap = patch(isocaps(Ls,100),...
+ 'FaceColor','interp',...
+ 'EdgeColor','none');
+colormap(jet)
+view(45,30)
+axis tight
+
+daspect([1,1,.4])
+%better rendering
+lightangle(305,30);
+set(gcf,'Renderer','zbuffer'); lighting phong
+isonormals(Ls,hiso)
+set(hcap,'AmbientStrength',.6)
+set(hiso,'SpecularColorReflectance',0,'SpecularExponent',50)
+view(215,30);
33 assign_fields.m
@@ -0,0 +1,33 @@
+function param=assign_fields(param,newparam)
+% Copy fields from structure "newparam" to like-named fields of
+% structure "param"
+%
+% Written by: E. Rietsch: October 15, 2006
+% Last updated:
+%
+% param=assign_fields(param,newparam)
+% INPUT
+% param structure with fields representing default parameters
+% newparam structure with fields representing new parameters
+% OUTPUT
+% param input structure with updated fields
+
+% Check if the fields of newparam are a subset of those of "param"
+if isempty(newparam)
+ return
+end
+
+fields=fieldnames(param);
+newfields=fieldnames(newparam);
+bool=ismember(newfields,fields);
+if ~all(bool)
+ disp(' The following fieldes are probably misspelled:')
+ disp([' ',cell2str(newfields(~bool),',')])
+ disp(' Possible fields are:')
+ disp([' ',cell2str(fields,',')])
+ error('Abnormal termination')
+end
+
+for ii=1:length(newfields)
+ param.(newfields{ii})=newparam.(newfields{ii});
+end
32 average_intensity.m
@@ -0,0 +1,32 @@
+function result=average_intensity(img,mask,r)
+
+[x0,y0]=centerpos(mask);
+
+sizex=size(img,2);
+sizey=size(img,1);
+
+pxcount=0.0;
+intensity=0.0;
+
+for i=1:sizex
+ for j=1:sizey
+
+ distance=sqrt( (i-x0)^2 + (j-y0)^2 );
+ if (distance <= r)
+ pxcount = pxcount + 1;
+ end
+ end
+end
+
+for i=1:sizex
+ for j=1:sizey
+
+ distance=sqrt( (i-x0)^2 + (j-y0)^2 );
+ if (distance <= r)
+ intensity = intensity + img(j,i)/pxcount;
+ end
+ end
+end
+
+
+result = intensity;
15 avg_array.m
@@ -0,0 +1,15 @@
+function result=avg_array(data,n)
+
+length=size(data,2);
+last=1;
+num=1;
+
+while (last<length)
+ last=last+n;
+ if (last<length)
+ newdata(num)=mean(data(1,(last-n):last));
+ num=num+1;
+ end
+end
+
+result=newdata;
10 avg_array_inc.m
@@ -0,0 +1,10 @@
+function result=avg_array_inc(data,num)
+
+length=size(data);
+n=num/2;
+
+for i=(1+n):(length-n)
+ newdata(i)=mean(data( (i-n):(i+n) ));
+end
+
+result=newdata;
21 avg_frames.m
@@ -0,0 +1,21 @@
+function result = avg_frames(img,n)
+
+sizeX=size(img,2);
+sizeY=size(img,1);
+length=size(img,3);
+
+last_frame=1;
+framenum=1;
+
+while (last_frame < length)
+ last_frame=last_frame+n;
+ for x=1:sizeX
+ for y=1:sizeY
+
+ newimg(y,x,framenum)=mean(img(y,x,(last_frame-n):last_frame));
+
+ end
+ end
+end
+
+result = newimg;
25 bersen.m
@@ -0,0 +1,25 @@
+function final_result = bersen(img,boxsize,ct)
+
+% Implements Bersen contrast threshold.
+[Ny,Nx]=size(img);
+
+timg=zeros(Ny,Nx);
+w=round((boxsize-1)/2);
+
+for y=(2+w):(Ny-w-1)
+ for x=(2+w):(Nx-w-1)
+
+ %pixel=img(y,x);
+ box=img( (y-w):(y+w), (x-w):(x+w) );
+ contrast=max(max(double(box))) - min(min(double(box)));
+ if contrast>ct
+ timg(y,x)=1;
+ end;
+ end;
+
+ if mod(y,20)==0
+ y
+ end;
+end;
+
+final_result=timg;
33 bondlengths.m
@@ -0,0 +1,33 @@
+function final_result = bondlengths(db)
+
+num_circles = size(db,1);
+
+used=zeros(num_circles,2);
+
+k=0;
+for i=1:num_circles
+ flag=0;
+
+ x0=db(i,1); y0=db(i,2); x1=db(i,4); y1=db(i,5);
+
+ % Prevent double counting.
+ for j=1:numcircles
+ xc=used(j,1); yc=used(j,2);
+ if x0==xc && y0==yc
+ flag=1;
+ end
+ end
+
+ if flag==1
+ continue;
+ end
+
+ if x1==0 || y1==0
+ continue;
+ end
+
+ k=k+1;
+ bonds(k)=sqrt( (x1-x0)^2 + (y1-y0)^2);
+ used(i,1)=x0; used(i,2)=y0;
+
+
63 bounded_array.m
@@ -0,0 +1,63 @@
+function [x,ier]=bounded_array(x,lowerx,upperx)
+% Modify array "x" in such a way that its entries do not exceed "upperx"
+% and are not below "lowerx".
+%
+% Written by: E. Rietsch: October 21, 2006
+% Last updated:
+%
+% [x,ier]=bounded_array(x,lowerx,upperx)
+% INPUT
+% x vector or matrix
+% lowerx lower limit of the entries of "x"; either a constant or
+% a matrix with the same dimension as "x"
+% upperx upper limit of the entries of "x"; either a constant or
+% a matrix with the same dimension as "x"
+% OUTPUT
+% x input matrix where entries < "lowerx" are replaced by "lowerx" and
+% entries > "upperx" are replaced by "upperx".
+% ier logical variable: "false" if no entry of "x" has been replaced and
+% "true" otherwise
+%
+% EXAMPLE
+% x0=randn(3,5);
+% lowerx=-1;
+% upperx=ones(3,5);
+% [x,ier]=bounded_array(x0,lowerx,upperx);
+% disp(x-x0)
+
+
+% Error checking
+dimsx=size(x);
+if numel(lowerx) > 1
+ dimsl=size(lowerx);
+ if ~all(dimsx == dimsl)
+ error('Incompatible dimensions of "x" and "lowerx".')
+ end
+end
+if numel(upperx) > 1
+ dimsu=size(upperx);
+ if ~all(dimsx == dimsu)
+ error('Incompatible dimensions of "x" and "upperx".')
+ end
+end
+
+ier=false;
+bool=x > upperx;
+if any(bool(:))
+ if numel(upperx) == 1
+ x(bool)=upperx;
+ else
+ x(bool)=upperx(bool);
+ end
+ ier=true;
+end
+
+bool=x < lowerx;
+if any(bool(:))
+ if numel(lowerx) == 1
+ x(bool)=lowerx;
+ else
+ x(bool)=lowerx(bool);
+ end
+ ier=true;
+end
56 boxtracking.m
@@ -0,0 +1,56 @@
+function [rboxes,rthresholds,rmasks,rmasked]=boxtracking(img)
+
+r=19;
+r1=18;
+r2=15;
+radrange=[11 22];
+
+for f=1:size(img,3)
+ clear boxes masks obw masked;
+
+ I=img(:,:,f);
+
+ [accum,circen,cirrad]=CircularHough_Grd(I,radrange);
+
+ for n=1:size(circen,1)
+ x=floor(circen(n,1));
+ y=floor(circen(n,2));
+ boxes(:,:,n)=I(y-r:y+r, x-r:x+r);
+ end
+ if (exist('boxes'))
+ ;
+ else
+ continue;
+ end
+
+ for n=1:size(boxes,3)
+ obw(:,:,n)=im2bw(boxes(:,:,n), graythresh(boxes(:,:,n)));
+ end
+
+ x0=(size(obw,1)-1)/2; y0=(size(obw,2)-1)/2;
+
+ mask=obw;
+ for n=1:size(obw,3)
+ for x=1:size(obw,1)
+ for y=1:size(obw,2)
+ d=sqrt( (x-x0)^2 + (y-y0)^2 );
+ if ( (d>r1) || (d<r2) )
+ mask(y,x,n)=0;
+ end
+ end
+ end
+ end
+
+ masked=boxes;
+ masked(mask==0)=0;
+
+ for n=1:size(masked,3)
+ ints(n)=sum(sum(masked(:,:,n)));
+ end
+
+ % save data
+ rmasks{f}=mask;
+ rboxes{f}=boxes;
+ rthresholds{f}=obw;
+ rmasked{f}=masked;
+end
155 bpass.m
@@ -0,0 +1,155 @@
+function res = bpass(image_array,lnoise,lobject,threshold)
+%
+% NAME:
+% bpass
+% PURPOSE:
+% Implements a real-space bandpass filter that suppresses
+% pixel noise and long-wavelength image variations while
+% retaining information of a characteristic size.
+%
+% CATEGORY:
+% Image Processing
+% CALLING SEQUENCE:
+% res = bpass( image_array, lnoise, lobject )
+% INPUTS:
+% image: The two-dimensional array to be filtered.
+% lnoise: Characteristic lengthscale of noise in pixels.
+% Additive noise averaged over this length should
+% vanish. May assume any positive floating value.
+% May be set to 0 or false, in which case only the
+% highpass "background subtraction" operation is
+% performed.
+% lobject: (optional) Integer length in pixels somewhat
+% larger than a typical object. Can also be set to
+% 0 or false, in which case only the lowpass
+% "blurring" operation defined by lnoise is done,
+% without the background subtraction defined by
+% lobject. Defaults to false.
+% threshold: (optional) By default, after the convolution,
+% any negative pixels are reset to 0. Threshold
+% changes the threshhold for setting pixels to
+% 0. Positive values may be useful for removing
+% stray noise or small particles. Alternatively, can
+% be set to -Inf so that no threshholding is
+% performed at all.
+%
+% OUTPUTS:
+% res: filtered image.
+% PROCEDURE:
+% simple convolution yields spatial bandpass filtering.
+% NOTES:
+% Performs a bandpass by convolving with an appropriate kernel. You can
+% think of this as a two part process. First, a lowpassed image is
+% produced by convolving the original with a gaussian. Next, a second
+% lowpassed image is produced by convolving the original with a boxcar
+% function. By subtracting the boxcar version from the gaussian version, we
+% are using the boxcar version to perform a highpass.
+%
+% original - lowpassed version of original => highpassed version of the
+% original
+%
+% Performing a lowpass and a highpass results in a bandpassed image.
+%
+% Converts input to double. Be advised that commands like 'image' display
+% double precision arrays differently from UINT8 arrays.
+
+% MODIFICATION HISTORY:
+% Written by David G. Grier, The University of Chicago, 2/93.
+%
+% Greatly revised version DGG 5/95.
+%
+% Added /field keyword JCC 12/95.
+%
+% Memory optimizations and fixed normalization, DGG 8/99.
+% Converted to Matlab by D.Blair 4/2004-ish
+%
+% Fixed some bugs with conv2 to make sure the edges are
+% removed D.B. 6/05
+%
+% Removed inadvertent image shift ERD 6/05
+%
+% Added threshold to output. Now sets all pixels with
+% negative values equal to zero. Gets rid of ringing which
+% was destroying sub-pixel accuracy, unless window size in
+% cntrd was picked perfectly. Now centrd gets sub-pixel
+% accuracy much more robustly ERD 8/24/05
+%
+% Refactored for clarity and converted all convolutions to
+% use column vector kernels for speed. Running on my
+% macbook, the old version took ~1.3 seconds to do
+% bpass(image_array,1,19) on a 1024 x 1024 image; this
+% version takes roughly half that. JWM 6/07
+%
+% This code 'bpass.pro' is copyright 1997, John C. Crocker and
+% David G. Grier. It should be considered 'freeware'- and may be
+% distributed freely in its original form when properly attributed.
+
+if nargin < 3, lobject = false; end
+if nargin < 4, threshold = 0; end
+
+normalize = @(x) x/sum(x);
+
+image_array = double(image_array);
+
+if lnoise == 0
+ gaussian_kernel = 1;
+else
+ gaussian_kernel = normalize(...
+ exp(-((-ceil(5*lnoise):ceil(5*lnoise))/(2*lnoise)).^2));
+end
+
+if lobject
+ boxcar_kernel = normalize(...
+ ones(1,length(-round(lobject):round(lobject))));
+end
+
+% JWM: Do a 2D convolution with the kernels in two steps each. It is
+% possible to do the convolution in only one step per kernel with
+%
+ % gconv = conv2(gaussian_kernel',gaussian_kernel,image_array,'same');
+ % bconv = conv2(boxcar_kernel', boxcar_kernel,image_array,'same');
+%
+% but for some reason, this is slow. The whole operation could be reduced
+% to a single step using the associative and distributive properties of
+% convolution:
+%
+ % filtered = conv2(image_array,...
+ % gaussian_kernel'*gaussian_kernel - boxcar_kernel'*boxcar_kernel,...
+ % 'same');
+%
+% But this is also comparatively slow (though inexplicably faster than the
+% above). It turns out that convolving with a column vector is faster than
+% convolving with a row vector, so instead of transposing the kernel, the
+% image is transposed twice.
+
+gconv = conv2(image_array',gaussian_kernel','same');
+gconv = conv2(gconv',gaussian_kernel','same');
+
+if lobject
+ bconv = conv2(image_array',boxcar_kernel','same');
+ bconv = conv2(bconv',boxcar_kernel','same');
+
+ filtered = gconv - bconv;
+else
+ filtered = gconv;
+end
+
+% Zero out the values on the edges to signal that they're not useful.
+lzero = max(lobject,ceil(5*lnoise));
+
+filtered(1:(round(lzero)),:) = 0;
+filtered((end - lzero + 1):end,:) = 0;
+filtered(:,1:(round(lzero))) = 0;
+filtered(:,(end - lzero + 1):end) = 0;
+
+% JWM: I question the value of zeroing out negative pixels. It's a
+% nonlinear operation which could potentially mess up our expectations
+% about statistics. Is there data on 'Now centroid gets subpixel accuracy
+% much more robustly'? To choose which approach to take, uncomment one of
+% the following two lines.
+% ERD: The negative values shift the peak if the center of the cntrd mask
+% is not centered on the particle.
+
+% res = filtered;
+filtered(filtered < threshold) = 0;
+res = filtered;
82 canny_edge.m
@@ -0,0 +1,82 @@
+function [Ix,Iy,NVI,Ibw,I_temp]=canny_edge(img)
+
+% The algorithm parameters:
+% 1. Parameters of edge detecting filters:
+% X-axis direction filter:
+ Nx1=10;Sigmax1=1;Nx2=10;Sigmax2=1;Theta1=pi/2;
+% Y-axis direction filter:
+ Ny1=10;Sigmay1=1;Ny2=10;Sigmay2=1;Theta2=0;
+% 2. The thresholding parameter alfa:
+ alfa=0.1;
+
+% % Get the initial image lena.gif
+% [x,map]=gifread('lena.gif');
+%w=ind2gray(x,map);
+%figure(1);colormap(gray);
+%subplot(3,2,1);
+%imagesc(w,200);
+%title('Image: lena.gif');
+
+% Take grayscale image as argument.
+w=img;
+figure(1);colormap(gray);
+subplot(3,2,1);imagesc(w);
+
+% X-axis direction edge detection
+subplot(3,2,2);
+filterx=d2dgauss(Nx1,Sigmax1,Nx2,Sigmax2,Theta1);
+Ix= conv2(w,filterx,'same');
+imagesc(Ix);
+title('Ix');
+
+% Y-axis direction edge detection
+subplot(3,2,3)
+filtery=d2dgauss(Ny1,Sigmay1,Ny2,Sigmay2,Theta2);
+Iy=conv2(w,filtery,'same');
+imagesc(Iy);
+title('Iy');
+
+% Norm of the gradient (Combining the X and Y directional derivatives)
+subplot(3,2,4);
+NVI=sqrt(Ix.*Ix+Iy.*Iy);
+imagesc(NVI);
+title('Norm of Gradient');
+
+% Thresholding
+I_max=max(max(NVI));
+I_min=min(min(NVI));
+level=alfa*(I_max-I_min)+I_min;
+subplot(3,2,5);
+Ibw=max(NVI,level.*ones(size(NVI)));
+imagesc(Ibw);
+title('After Thresholding');
+
+% Thinning (Using interpolation to find the pixels where the norms of
+% gradient are local maximum.)
+subplot(3,2,6);
+[n,m]=size(Ibw);
+for i=2:n-1,
+for j=2:m-1,
+ if Ibw(i,j) > level,
+ X=[-1,0,+1;-1,0,+1;-1,0,+1];
+ Y=[-1,-1,-1;0,0,0;+1,+1,+1];
+ Z=[Ibw(i-1,j-1),Ibw(i-1,j),Ibw(i-1,j+1);
+ Ibw(i,j-1),Ibw(i,j),Ibw(i,j+1);
+ Ibw(i+1,j-1),Ibw(i+1,j),Ibw(i+1,j+1)];
+ XI=[Ix(i,j)/NVI(i,j), -Ix(i,j)/NVI(i,j)];
+ YI=[Iy(i,j)/NVI(i,j), -Iy(i,j)/NVI(i,j)];
+ ZI=interp2(X,Y,Z,XI,YI);
+ if Ibw(i,j) >= ZI(1) & Ibw(i,j) >= ZI(2)
+ I_temp(i,j)=I_max;
+ else
+ I_temp(i,j)=I_min;
+ end
+ else
+ I_temp(i,j)=I_min;
+ end
+end
+end
+imagesc(I_temp);
+title('After Thinning');
+colormap(gray);
+
146 cefit_4C.m
@@ -0,0 +1,146 @@
+function cefit_4C(crop_ce)
+%CEFIT_4C Create plot of datasets and fits
+% CEFIT_4C(CROP_CE)
+% Creates a plot, similar to the plot in the main curve fitting
+% window, using the data that you provide as input. You can
+% apply this function to the same data you used with cftool
+% or with different data. You may want to edit the function to
+% customize the code and this help message.
+%
+% Number of datasets: 1
+% Number of fits: 3
+
+
+% Data from dataset "crop_ce":
+% Y = crop_ce:
+% Unweighted
+%
+% This function was automatically generated on 19-Mar-2010 14:16:24
+
+% Set up figure to receive datasets and fits
+f_ = clf;
+figure(f_);
+set(f_,'Units','Pixels','Position',[469 119 680 477]);
+legh_ = []; legt_ = {}; % handles and text for legend
+xlim_ = [Inf -Inf]; % limits of x axis
+ax_ = axes;
+set(ax_,'Units','normalized','OuterPosition',[0 0 1 1]);
+set(ax_,'Box','on');
+axes(ax_); hold on;
+
+
+% --- Plot data originally in dataset "crop_ce"
+x_1 = (1:numel(crop_ce))';
+crop_ce = crop_ce(:);
+h_ = line(x_1,crop_ce,'Parent',ax_,'Color',[0.333333 0 0.666667],...
+ 'LineStyle','none', 'LineWidth',1,...
+ 'Marker','square', 'MarkerSize',6);
+xlim_(1) = min(xlim_(1),min(x_1));
+xlim_(2) = max(xlim_(2),max(x_1));
+legh_(end+1) = h_;
+legt_{end+1} = 'crop_ce';
+
+% Nudge axis limits beyond data limits
+if all(isfinite(xlim_))
+ xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
+ set(ax_,'XLim',xlim_)
+else
+ set(ax_, 'XLim',[0.73999999999999999, 27.260000000000002]);
+end
+
+
+% --- Create fit "fit 1"
+ok_ = isfinite(x_1) & isfinite(crop_ce);
+if ~all( ok_ )
+ warning( 'GenerateMFile:IgnoringNansAndInfs', ...
+ 'Ignoring NaNs and Infs in data' );
+end
+st_ = [1.9645708521867939 -0.095342743493336671 -1.4253374742810385 -0.44632821275186407 ];
+ft_ = fittype('exp2');
+
+% Fit this model using new data
+cf_ = fit(x_1(ok_),crop_ce(ok_),ft_,'Startpoint',st_);
+
+% Or use coefficients from the original fit:
+if 0
+ cv_ = { 83655.532131371612, -0.13866236586278774, -83654.65044407727, -0.13866501234521542};
+ cf_ = cfit(ft_,cv_{:});
+end
+
+% Plot this fit
+h_ = plot(cf_,'fit',0.95);
+legend off; % turn off legend from plot method call
+set(h_(1),'Color',[1 0 0],...
+ 'LineStyle','-', 'LineWidth',2,...
+ 'Marker','none', 'MarkerSize',6);
+legh_(end+1) = h_(1);
+legt_{end+1} = 'fit 1';
+
+% --- Create fit "fit 2"
+ok_ = isfinite(x_1) & isfinite(crop_ce);
+if ~all( ok_ )
+ warning( 'GenerateMFile:IgnoringNansAndInfs', ...
+ 'Ignoring NaNs and Infs in data' );
+end
+st_ = [0.58092756599744022 -0.10000000000000001 0.9735777080550454 -0.10000000000000001 ];
+ft_ = fittype('a*x*exp(b*x)+c*x*exp(d*x)',...
+ 'dependent',{'y'},'independent',{'x'},...
+ 'coefficients',{'a', 'b', 'c', 'd'});
+
+% Fit this model using new data
+cf_ = fit(x_1(ok_),crop_ce(ok_),ft_,'Startpoint',st_);
+
+% Or use coefficients from the original fit:
+if 0
+ cv_ = { 1.8786802379116481, -1.0844933225167586, 0.41359663918692008, -0.16344125490294698};
+ cf_ = cfit(ft_,cv_{:});
+end
+
+% Plot this fit
+h_ = plot(cf_,'fit',0.95);
+legend off; % turn off legend from plot method call
+set(h_(1),'Color',[0 0 1],...
+ 'LineStyle','-', 'LineWidth',2,...
+ 'Marker','none', 'MarkerSize',6);
+legh_(end+1) = h_(1);
+legt_{end+1} = 'fit 2';
+
+% --- Create fit "fit 2 copy 1"
+ok_ = isfinite(x_1) & isfinite(crop_ce);
+if ~all( ok_ )
+ warning( 'GenerateMFile:IgnoringNansAndInfs', ...
+ 'Ignoring NaNs and Infs in data' );
+end
+st_ = [0.5411802150500159 -0.10000000000000001 ];
+ft_ = fittype('a*x*exp(b*x)',...
+ 'dependent',{'y'},'independent',{'x'},...
+ 'coefficients',{'a', 'b'});
+
+% Fit this model using new data
+cf_ = fit(x_1(ok_),crop_ce(ok_),ft_,'Startpoint',st_);
+
+% Or use coefficients from the original fit:
+if 0
+ cv_ = { 0.5583592834424036, -0.18900594409469085};
+ cf_ = cfit(ft_,cv_{:});
+end
+
+% Plot this fit
+h_ = plot(cf_,'fit',0.95);
+legend off; % turn off legend from plot method call
+set(h_(1),'Color',[0.666667 0.333333 0],...
+ 'LineStyle','-', 'LineWidth',2,...
+ 'Marker','none', 'MarkerSize',6);
+legh_(end+1) = h_(1);
+legt_{end+1} = 'fit 2 copy 1';
+
+% Done plotting data and fits. Now finish up loose ends.
+hold off;
+leginfo_ = {'Orientation', 'vertical'};
+h_ = legend(ax_,legh_,legt_,leginfo_{:}); % create and reposition legend
+set(h_,'Units','normalized');
+t_ = get(h_,'Position');
+t_(1:2) = [0.640931,0.726066];
+set(h_,'Interpreter','none','Position',t_);
+xlabel(ax_,''); % remove x label
+ylabel(ax_,''); % remove y label
20 cell2num.m
@@ -0,0 +1,20 @@
+function x=cell2num(cx)
+% Convert a cell array of numbers into a numeric array or
+% a cell array of logical variables into an array of logical variables
+%
+% Written by: E. Rietsch: March 1, 2006
+% Last updated:
+%
+% x=cell2num(cx)
+% INPUT
+% cx cell array whose entries a numbers or logical variables
+% OUTPUT
+% x numeric or logical array of the same dimension with the same
+% numbers or logical variables
+%
+% EXAMPLE
+% cell2num({1,2;3,4})
+% cell2num({true,false;false,true})
+
+x=[cx{:}];
+x=reshape(x,size(cx));
22 centerpos.m
@@ -0,0 +1,22 @@
+function [x0,y0]=centerpos(mask)
+
+sizex=size(mask,2);
+sizey=size(mask,1);
+
+sumx=0;
+sumy=0;
+numwhite=0;
+
+for i=1:sizex
+ for j=1:sizey
+ test=mask(j,i);
+ if (test>0)
+ sumx=sumx+i;
+ sumy=sumy+j;
+ numwhite=numwhite+1;
+ end
+ end
+end
+
+x0=sumx/numwhite;
+y0=sumy/numwhite;
36 check_dist.m
@@ -0,0 +1,36 @@
+function [dist,avs,stddevs] = check_dist(db,rawimg,boxsize)
+
+[Ny,Nx] = size(rawimg);
+numCircles = size(db,1);
+bs=boxsize;
+k=0;
+for i=1:numCircles
+
+ if bs==0
+ boxsize = db(i,3);
+ end
+ x0=db(i,1); y0=db(i,2);
+
+ if ((x0-boxsize)<0) || ((y0-boxsize)<0) || ((x0+boxsize)>Nx) || ((y0+boxsize)>Ny)
+ continue;
+ end
+ clear localpx; k1=0;
+ for x=round(x0-boxsize):round(x0+boxsize)
+ for y=round(y0-boxsize):round(y0+boxsize)
+ dR = sqrt( (x-x0)^2 + (y-y0)^2 );
+ if dR < boxsize
+ k=k+1;
+ k1=k1+1;
+ pxlist(k)=rawimg(y,x);
+ localpx(k1)=rawimg(y,x);
+ end
+ end
+ end
+
+ sd(i)=std(double(localpx)); m(i)=mean(double(localpx));
+
+ i
+
+end
+
+dist = pxlist; stddevs=sd; avs=m;
143 cntrd.m
@@ -0,0 +1,143 @@
+function out=cntrd(im,mx,sz,interactive)
+% out=cntrd(im,mx,sz,interactive)
+%
+% PURPOSE: calculates the centroid of bright spots to sub-pixel accuracy.
+% Inspired by Grier & Crocker's feature for IDL, but greatly simplified and optimized
+% for matlab
+%
+% INPUT:
+% im: image to process, particle should be bright spots on dark background with little noise
+% ofen an bandpass filtered brightfield image or a nice fluorescent image
+%
+% mx: locations of local maxima to pixel-level accuracy from pkfnd.m
+%
+% sz: diamter of the window over which to average to calculate the centroid.
+% should be big enough
+% to capture the whole particle but not so big that it captures others.
+% if initial guess of center (from pkfnd) is far from the centroid, the
+% window will need to be larger than the particle size. RECCOMMENDED
+% size is the long lengthscale used in bpass plus 2.
+%
+%
+% interactive: OPTIONAL INPUT set this variable to one and it will show you the image used to calculate
+% each centroid, the pixel-level peak and the centroid
+%
+% NOTE:
+% - if pkfnd, and cntrd return more then one location per particle then
+% you should try to filter your input more carefully. If you still get
+% more than one peak for particle, use the optional sz parameter in pkfnd
+% - If you want sub-pixel accuracy, you need to have a lot of pixels in your window (sz>>1).
+% To check for pixel bias, plot a histogram of the fractional parts of the resulting locations
+% - It is HIGHLY recommended to run in interactive mode to adjust the parameters before you
+% analyze a bunch of images.
+%
+% OUTPUT: a N x 4 array containing, x, y and brightness for each feature
+% out(:,1) is the x-coordinates
+% out(:,2) is the y-coordinates
+% out(:,3) is the brightnesses
+% out(:,4) is the sqare of the radius of gyration
+%
+% CREATED: Eric R. Dufresne, Yale University, Feb 4 2005
+% 5/2005 inputs diamter instead of radius
+% Modifications:
+% D.B. (6/05) Added code from imdist/dist to make this stand alone.
+% ERD (6/05) Increased frame of reject locations around edge to 1.5*sz
+% ERD 6/2005 By popular demand, 1. altered input to be formatted in x,y
+% space instead of row, column space 2. added forth column of output,
+% rg^2
+% ERD 8/05 Outputs had been shifted by [0.5,0.5] pixels. No more!
+% ERD 8/24/05 Woops! That last one was a red herring. The real problem
+% is the "ringing" from the output of bpass. I fixed bpass (see note),
+% and no longer need this kludge. Also, made it quite nice if mx=[];
+% ERD 6/06 Added size and brightness output ot interactive mode. Also
+% fixed bug in calculation of rg^2
+% JWM 6/07 Small corrections to documentation
+
+
+if nargin==3
+ interactive=0;
+end
+
+if sz/2 == floor(sz/2)
+warning('sz must be odd, like bpass');
+end
+
+if isempty(mx)
+ warning('there were no positions inputted into cntrd. check your pkfnd theshold')
+ out=[];
+ return;
+end
+
+
+r=(sz+1)/2;
+%create mask - window around trial location over which to calculate the centroid
+m = 2*r;
+x = 0:(m-1) ;
+cent = (m-1)/2;
+x2 = (x-cent).^2;
+dst=zeros(m,m);
+for i=1:m
+ dst(i,:)=sqrt((i-1-cent)^2+x2);
+end
+
+
+ind=find(dst < r);
+
+msk=zeros([2*r,2*r]);
+msk(ind)=1.0;
+%msk=circshift(msk,[-r,-r]);
+
+dst2=msk.*(dst.^2);
+ndst2=sum(sum(dst2));
+
+[nr,nc]=size(im);
+%remove all potential locations within distance sz from edges of image
+ind=find(mx(:,2) > 1.5*sz & mx(:,2) < nr-1.5*sz);
+mx=mx(ind,:);
+ind=find(mx(:,1) > 1.5*sz & mx(:,1) < nc-1.5*sz);
+mx=mx(ind,:);
+
+[nmx,crap] = size(mx);
+
+%inside of the window, assign an x and y coordinate for each pixel
+xl=zeros(2*r,2*r);
+for i=1:2*r
+ xl(i,:)=(1:2*r);
+end
+yl=xl';
+
+pts=[];
+%loop through all of the candidate positions
+for i=1:nmx
+ %create a small working array around each candidate location, and apply the window function
+ tmp=msk.*im((mx(i,2)-r+1:mx(i,2)+r),(mx(i,1)-r+1:mx(i,1)+r));
+ %calculate the total brightness
+ norm=sum(sum(tmp));
+ %calculate the weigthed average x location
+ xavg=sum(sum(tmp.*xl))./norm;
+ %calculate the weighted average y location
+ yavg=sum(sum(tmp.*yl))./norm;
+ %calculate the radius of gyration^2
+ %rg=(sum(sum(tmp.*dst2))/ndst2);
+ rg=(sum(sum(tmp.*dst2))/norm);
+
+ %concatenate it up
+ pts=[pts,[mx(i,1)+xavg-r,mx(i,2)+yavg-r,norm,rg]'];
+
+ %OPTIONAL plot things up if you're in interactive mode
+ if interactive==1
+ imagesc(tmp)
+ axis image
+ hold on;
+ plot(xavg,yavg,'x')
+ plot(xavg,yavg,'o')
+ plot(r,r,'.')
+ hold off
+ title(['brightness ',num2str(norm),' size ',num2str(sqrt(rg))])
+ pause
+ end
+
+
+end
+out=pts';
+
33 compare_labels.m
@@ -0,0 +1,33 @@
+function result = compare_labels(stack1,stack2)
+
+zmax1=size(stack1,3);
+zmax2=size(stack2,3);
+if zmax1 >= zmax2
+ zm=zmax1;
+else
+ zm=zmax2;
+end
+
+for z=1:zm
+
+ subplot(1,2,1);
+ if z<=zmax1
+ imshow(label2rgb(stack1(:,:,z),'jet','w'));
+ title(num2str(z));
+ else
+ imshow(label2rgb(stack1(:,:,z),'jet','w'));
+ title(num2str(zmax1));
+ end
+
+ subplot(1,2,2);
+ if z<=zmax2
+ imshow(label2rgb(stack2(:,:,z),'jet','w'));
+ title(num2str(z));
+ else
+ imshow(label2rgb(stack2(:,:,z),'jet','w'));
+ title(num2str(zmax2));
+ end
+
+ pause(0.05);
+end
+
34 compare_stacks.m
@@ -0,0 +1,34 @@
+function result = compare_stacks(stack1,stack2)
+
+zmax1=size(stack1,3);
+zmax2=size(stack2,3);
+if zmax1 >= zmax2
+ zm=zmax1;
+else
+ zm=zmax2;
+end
+
+for z=1:zm
+
+ subplot(1,2,1);
+ if z<=zmax1
+ imshow(stack1(:,:,z));
+ title(num2str(z));
+ else
+ imshow(stack1(:,:,zmax1));
+ title(num2str(zmax1));
+ end
+
+ subplot(1,2,2);
+ if z<=zmax2
+ imshow(stack2(:,:,z));
+ title(num2str(z));
+ else
+ imshow(stack2(:,:,zmax2));
+ title(num2str(zmax2));
+ end
+ f(z)=getframe(gcf);
+ pause(0.5);
+end
+
+result=f;
34 compare_stacks_labels.m
@@ -0,0 +1,34 @@
+function result = compare_stacks_labels(stack,labels)
+
+zmax1=size(stack,3);
+zmax2=size(labels,3);
+if zmax1 >= zmax2
+ zm=zmax1;
+else
+ zm=zmax2;
+end
+
+for z=1:zm
+
+ subplot(1,2,1);
+ if z<=zmax1
+ imshow(stack(:,:,z));
+ title(num2str(z));
+ else
+ imshow(stack(:,:,zmax1));
+ title(num2str(zmax1));
+ end
+
+ subplot(1,2,2);
+ if z<=zmax2
+ imshow(label2rgb(labels(:,:,z),'jet','w'));
+ title(num2str(z));
+ else
+ imshow(label2rgb(labels(:,:,zmax2),'jet','w'));
+ title(num2str(zmax2));
+ end
+ f(z)=getframe(gcf);
+ pause(0.1);
+end
+
+result=f;
20 connect_circles.m
@@ -0,0 +1,20 @@
+function connect_circles(dbs)
+% connect_circles.m
+% Draws bonds between halves of dumbbell particles.
+% You should issue the command imshow(image) first.
+
+hold on;
+
+numc= size(dbs,1);
+
+for i=1:numc
+
+ if dbs(i,7)==0 || dbs(i,3)==0;
+ continue;
+ end
+ xset = [dbs(i,1) dbs(i,4)];
+ yset = [dbs(i,2) dbs(i,5)];
+ plot(xset,yset,'r-');
+end
+plot(dbs(:,1), dbs(:,2), 'b+');
+hold off;
55 convert_tifs.m
@@ -0,0 +1,55 @@
+
+%This script converts files with names in the format
+%"untitled000.tif" into another form of .tif that this computer can read
+%and puts them in a folder called 'resaved using matlab' inside the same folder the images are in.
+
+filepath1 = input('Enter the file path for the folder where the images are\n i.e. C:\\Documents and Settings\\chansen3.HISPEEDCASTER\\Desktop\\01202009\\\n', 's');
+
+totalimagenumber= input('What is the highest number in the filename? For example in the file "untitled025.tif?" you would type "25"');
+
+%makes folders
+ file_path4=strcat(filepath1,'\resaved using matlab');
+ mkdir(file_path4);
+file_path5=strcat(filepath1,'\extraneous images');
+ mkdir(file_path5);
+
+
+imagenumber=0;
+
+while imagenumber<=totalimagenumber;
+ xx=0;
+ while xx==0;
+ str_imagenumber=num2str(imagenumber,'%03d');
+ file=strcat('untitled',str_imagenumber,'.tif');
+ file_path=strcat(filepath1,'\',file);
+
+ % Verify that the file exists.
+ fid = fopen(file_path, 'r');
+ if (fid == -1)
+ if ~isempty(dir(file_path))
+ error('MATLAB:imread:fileOpen', ['Can''t open file "%s" for reading;\nyou' ...
+ ' may not have read permission.'], ...
+ file_path);
+ else
+ imagenumber=imagenumber+1;
+ end
+ else
+ file_path = fopen(fid);
+ fclose(fid);
+ xx=1;
+ end
+
+ end
+ I=imread(file_path,'tif');
+ I16=immultiply(I,16);
+
+
+ str_newimage=num2str(imagenumber,'%03d');
+ file3=strcat('untitled',str_newimage,'.tif');
+ file_path3=strcat(file_path4,'\',file3);
+ imwrite(I16,file_path3);
+ 'added image written'
+ imagenumber
+ imagenumber=imagenumber+1;
+end
+
67 create_help_menu.m
@@ -0,0 +1,67 @@
+function create_help_menu(figure_handle)
+
+
+% Check if a help file exists
+path=fileparts(which('volume_browser'));
+if ~exist(path,'dir')
+ alert('Help-file directory has not been found. No help button created.')
+ return
+end
+
+callback=@help4BV;
+menuid = uimenu(figure_handle,'Label',' Need help? ','ForegroundColor','red');
+
+uimenu(menuid,'Label','General','Callback', ...
+ {callback,path,figure_handle,'help4VB_general'})
+uimenu(menuid,'Label','Explore volume','Callback', ...
+ {callback,path,figure_handle,'help4VB_explore'})
+
+uimenu(menuid,'Label',' &About ','Separator','on','enable','on', ...
+ 'CallBack',{@about_callback,figure_handle});
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function help4BV(varargin)
+
+% path=fileparts(which('volume_browser'));
+path=varargin{3};
+figure_handle=varargin{4};
+funct=varargin{5};
+
+% Check if a help file exists
+hfile=[funct,'.txt'];
+filename=fullfile(path,'HelpFiles',hfile);
+
+% Read help file
+try
+ help_info=textread(filename,'%s','delimiter','\n');
+catch
+ alert(['Help file "',hfile,'" has not been found.'])
+ return
+end
+
+handle=helpdlg(help_info,'Volume Browser Help');
+
+add_handle2delete1(handle,figure_handle)
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function about_callback(varargin)
+
+% global V3D_HANDLES
+
+createmode.Interpreter='tex';
+createmode.WindowStyle='modal';
+handle=msgbox({...
+'Matlab function to vizualize 3D Data; release 1.02';
+'Copyright \copyright 2006 2007 Eike Rietsch';
+' ';
+'This program is free software; you can redistribute it and/or modify it under the terms of the BSD License.'
+},createmode);
+
+% Store handle of message box so that it can be deleted (if it
+% still exists) upon closing of the browser
+
+add_handle2delete1(handle,varargin{3})
16 d2dgauss.m
@@ -0,0 +1,16 @@
+%%%%%%% The functions used in the main.m file %%%%%%%
+% Function "d2dgauss.m":
+% This function returns a 2D edge detector (first order derivative
+% of 2D Gaussian function) with size n1*n2; theta is the angle that
+% the detector rotated counter clockwise; and sigma1 and sigma2 are the
+% standard deviation of the gaussian functions.
+function h = d2dgauss(n1,sigma1,n2,sigma2,theta)
+r=[cos(theta) -sin(theta);
+ sin(theta) cos(theta)];
+for i = 1 : n2
+ for j = 1 : n1
+ u = r * [j-(n1+1)/2 i-(n2+1)/2]';
+ h(i,j) = gauss(u(1),sigma1)*dgauss(u(2),sigma2);
+ end
+end
+h = h / sqrt(sum(sum(abs(h).*abs(h))));
35 demo_thresh_box.m
@@ -0,0 +1,35 @@
+function intensities=demo_thresh_box(img)
+
+intensities = zeros( size(img, 3), 2 ); % mean, std
+
+for n=1:size(img,3)
+ img0=img(:,:,n);
+
+
+ [accum,circen,cirrad]=CircularHough_Grd(img(:,:,10),[11 22]);
+ clear boxint;
+
+ for c=1:size(circen,1)
+
+ r0=floor(cirrad(c));
+ x=floor(circen(c,1));
+ y=floor(circen(c,2));
+
+ clear box;
+
+ box=img0(y-r0:y+r0, x-r0:x+r0);
+ nbox=niblack(box,2,-0.2);
+ obox=im2bw(box,graythresh(box));
+ andbox=nbox & obox;
+
+ box(andbox==0) = 0;
+
+ boxint(c) = sum(sum(box));
+
+ end
+
+ intensities(n,1) = mean(double(boxint));
+ intensities(n,2) = std(double(boxint));
+
+ n
+end
12 devices_volfrac.m
@@ -0,0 +1,12 @@
+function result=devices_volfrac(pvol,volumes,numbers)
+
+for i=1:size(volumes,2)
+ for j=1:size(numbers,2)
+
+ vf=numbers(1,j)*pvol/volumes(1,i);
+ result(i,j) = vf;
+
+ end
+
+end
+
2  dgauss.m
@@ -0,0 +1,2 @@
+function y = dgauss(x,std)
+y = -x * gauss(x,std) / std^2;
141 dna_boxtracking.m
@@ -0,0 +1,141 @@
+function [pairs,average_ratio,positions,intensities,res,windowlist,masklist]=dna_boxtracking(imagestack,radrange,maxmove,boxradius,minr,maxr)
+
+numframes=size(imagestack,3);
+
+%%% DEBUG masklist windowlist
+m=0; w=0; masklist=[]; windowlist=[];
+
+
+% Get positionlist and tracklist from HoughTracking.
+[positionlist,tracklist]=HoughTracking(imagestack,radrange,1:numframes,0);
+
+% Now send tracklist to the particle tracking routine.
+param.mem=5; param.dim=2; param.good=5; param.quiet=0;
+res=track(tracklist,maxmove,param);
+
+% Now heal the tracks. When a particle drops away entirely, we assume it
+% remains stationary.
+
+numparticles=max(res(:,4));
+positions=[];
+k=0;
+for n=1:numparticles
+ k=k+1;
+ trace=res(find(res(:,4)==n),:);
+ t0=trace(1,3);
+ positions(k,1)=trace(1,1);
+ positions(k,2)=trace(1,2);
+ positions(k,3)=trace(1,3);
+ positions(k,4)=n;
+ %last=t0;
+ for f=(t0+1):numframes
+ k=k+1;
+ trindex=find(trace(:,3)==f);
+ if trindex
+ positions(k,1)=trace(trindex,1);
+ positions(k,2)=trace(trindex,2);
+ positions(k,3)=f;
+ positions(k,4)=n;
+ else
+ positions(k,1)=positions(k-1,1);
+ positions(k,2)=positions(k-1,2);
+ positions(k,3)=f;
+ positions(k,4)=n;
+ end
+ end
+end
+
+% Now that we have healed tracks we can find intensities on rings for each
+% point.
+r=boxradius;
+r1=maxr; r2=minr;
+intensities=zeros(size(positions,1),3);
+
+interval=10;
+done=size(positions,1);
+'Now calculating ring intensities.'
+
+for k=1:size(positions,1)
+
+ if mod( k, interval )==0
+ strcat(num2str(k), ' done of ', num2str(done))
+ end
+
+ x=floor(positions(k,1)); y=floor(positions(k,2));
+ if (x< (r+1) )
+ x=r+1;
+ end
+ if (y< (r+1) )
+ y=r+1;
+ end
+
+ f=positions(k,3); n=positions(k,4);
+ window = imagestack( y-r:y+r, x-r:x+r, f );
+
+ w=w+1;
+ %windowlist(:,:,w)=window;
+
+ %bw=im2bw(window, graythresh(window));
+ bw=im2bw(window, graythresh(imagestack(:,:,f)));
+
+ x0=(size(bw,1)-1)/2; y0=(size(bw,2)-1)/2;
+ mask=bw;
+ for x=1:size(bw,1)
+ for y=1:size(bw,2)
+ d=sqrt( (x-x0)^2 + (y-y0)^2 );
+ if ( (d>r1) || (d<r2) )
+ mask(y,x)=0;
+ end
+ end
+
+ end
+
+ m=m+1;
+
+ masked=window;
+ masked(mask==0)=0;
+
+ %masklist(:,:,m)=masked;
+
+ int_here=mean( mean( masked ) );
+
+ intensities(k,1)=positions(k,4);
+ intensities(k,2)=positions(k,3);
+ intensities(k,3)=int_here;
+end
+
+ % FINISH HERE
+
+%intensities=zeros(size(positions,1),3);
+%for k=1:size(positions,1)
+% intensities(k,1)=positions(k,4);
+% intensities(k,2)=positions(k,3);
+% x=positions(k,1); y=positions(k,2);
+% n=positions(k,3);
+% intensities(k,3)=mean(mean( imagestack( floor(y-(boxsize/2)):ceil(y+(boxsize/2)), floor(x-(boxsize/2)):ceil(x+(boxsize/2)),n) ));
+%end
+
+% Get user-generated pairing information.
+beep;beep;
+pairs=dna_pairing(imagestack,positions,intensities);
+
+% Finally, generate an average-ratio across time.
+for f=1:numframes
+ sumup=0;
+ num=0;
+ for k=1:size(pairs,3)
+ nowindex=find(pairs(:,1,k)==f);
+ if nowindex
+ sumup = sumup + pairs(nowindex,6,k);
+ num = num + 1;
+ end
+ end
+
+ average_ratio(f,1)=f;
+ if (num>0)
+ average_ratio(f,2)=double(sumup)/double(num);
+ end
+ average_ratio(f,3)=num;
+end
+
+
82 dna_pairing.m
@@ -0,0 +1,82 @@
+function pairs=dna_pairing(img,positions,intensities)
+
+% pairs data structure:
+% Each pair has its own matrix.
+% Column 1 is time.
+% Columns 2 and 4 are the particle labels..
+% Columns 3 and 5 are the intensities.