diff --git a/calculate_phase_mask.m b/calculate_phase_mask.m new file mode 100644 index 0000000..c804955 --- /dev/null +++ b/calculate_phase_mask.m @@ -0,0 +1,140 @@ +function [PHASE_MASK] = calculate_phase_mask(PHASE_QUALITY, KERNEL_RADIUS) + % This function computes a quality-based mask of the phase angle plane of a cross correlation. + + % Rename the kernel radius + rad = KERNEL_RADIUS; + + % Extract the central portion of the phase quality + phase_quality_interior = PHASE_QUALITY(rad + 1 : end - rad - 1, rad + 1 : end - rad - 1); + + phase_quality_sub = (phase_quality_interior - ... + min(phase_quality_interior(:))); + + % Rescale the phase quality so that its range is [0, 1] + phase_quality_scaled = phase_quality_sub ./ max(phase_quality_sub(:)); + + % Measure size + [region_height, region_width] = size(phase_quality_scaled); + + % Make coordinates (Cartesian) + [x, y] = meshgrid(1 : region_width, 1 : region_height); + + % Centroid + xc = region_width / 2; + yc = region_height / 2; + + % Make coordinates (polar), + % with origin (r = 0) at geometric centroid of the region + [~, r] = cart2pol(x - xc, y - yc); + + % Threshold the phase quality map using histogram equalization +% phase_quality_bw = (histeq(phase_quality_scaled, 2)); + phase_quality_bw = im2bw(phase_quality_scaled, 0.5); + +% Set border pixels to 1 to prevent connecting regions via the border + phase_quality_bw(1, :) = 1; + phase_quality_bw(end, :) = 1; + phase_quality_bw(:, 1) = 1; + phase_quality_bw(:, end) = 1; + + % Find properties of all the connected regions in the thresholded image. + phase_quality_region_props = regionprops(~phase_quality_bw, 'Centroid', 'PixelIdxList'); + + % Count the number of regions + num_regions = length(phase_quality_region_props); + + % Allocate a vector that will contain the median values + % of the radial coordinates of each region. + region_weighted_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the centroids of the detected regions + region_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the mean values of the radial coordinates + % of the pixels comprising the detected regions. + region_radius_median = zeros(num_regions, 1); + + % Measure the median radial coordinate of each region + for k = 1 : num_regions + + % Find the radial coordinate of the region's centroid. + % whose value is equal to the k'th iteration of this loop + region_centroid_subs = phase_quality_region_props(k).Centroid; + + % Radial coordinates of the pixels in the region + region_radius_median(k) = median(r(phase_quality_region_props(k).PixelIdxList)); + + % Convert to radial + region_centroid_radial(k) = sqrt((region_centroid_subs(1) - xc)^2 + (region_centroid_subs(2) - yc)^2); + + % Centroid coordinates weighted by the median radial coordinate of the pixels in the region + region_weighted_centroid_radial(k) = region_radius_median(k) * region_centroid_radial(k); + + end + + % Find the minimum + [~, region_weighted_centroid_idx] = min(region_weighted_centroid_radial); + + % Mask indices + mask_idx = phase_quality_region_props(region_weighted_centroid_idx).PixelIdxList; + + % Allocate a mask matrix + phase_mask_holder = zeros(region_height, region_width); + + % Apply the mask + phase_mask_holder(mask_idx) = 1; + + % Zero the border pixels just in cae + phase_mask_holder(:, 1) = 0 ; + phase_mask_holder(:, end) = 0; + phase_mask_holder(1, :) = 0; + phase_mask_holder(end, :) = 0; + + % Fill the holes + phase_mask_holder = imfill(phase_mask_holder); + + % Get the region properties again + phase_mask_region_props = regionprops(phase_mask_holder,... + 'PixelIdxList', 'MajorAxisLength', 'MinorAxisLength', ... + 'Orientation'); + + % Read the measured major axis of the ellipse + % that best fits the binary quality mask + ax_maj = phase_mask_region_props.MajorAxisLength; + + % Same for the minor axis + ax_min = phase_mask_region_props.MinorAxisLength; + + % Measure the orientation angle of + % the best fit ellipse and convert it to radians. + ax_angle = 1 * deg2rad(phase_mask_region_props.Orientation); + + % Gaussian standard deviations as fractions + % of the major and minor axes of the ellipse fit + std_maj = ax_maj / 1.8; + std_min = ax_min / 1.8; + + % Rotate the coordinates for the elliptical Gaussian + x2 = (x - xc) * cos(ax_angle) - (y - yc) * sin(ax_angle); + y2 = (x - xc) * sin(ax_angle) + (y - yc) * cos(ax_angle); + + % Calculate the Gaussian function + gaussian_mask = exp(-(x2.^2)/(std_maj^2) - (y2.^2) / (std_min^2)); + + % Insert the cropped mask into the phase + % quality matrix, so that it's the same + % size as the original phase matrix. + PHASE_MASK = zeros(size(PHASE_QUALITY)); + PHASE_MASK(rad + 1 : end - rad - 1, rad + 1 : end - rad - 1) = gaussian_mask; + +end + + + + + + + + + + diff --git a/calculate_phase_mask_ellipse.m b/calculate_phase_mask_ellipse.m new file mode 100644 index 0000000..34648be --- /dev/null +++ b/calculate_phase_mask_ellipse.m @@ -0,0 +1,137 @@ +function [PHASE_MASK] = calculate_phase_mask_ellipse(PHASE_QUALITY, KERNEL_RADIUS) + % This function computes a quality-based mask of the phase angle plane of a cross correlation. + + % Rename the kernel radius + rad = KERNEL_RADIUS; + + % Extract the central portion of the phase quality + phase_quality_interior = PHASE_QUALITY(rad + 1 : end - rad - 1, rad + 1 : end - rad - 1); + + phase_quality_sub = (phase_quality_interior - ... + min(phase_quality_interior(:))); + + % Rescale the phase quality so that its range is [0, 1] + phase_quality_scaled = phase_quality_sub ./ max(phase_quality_sub(:)); + + % Measure size + [region_height, region_width] = size(phase_quality_scaled); + + % Make coordinates (Cartesian) + [x, y] = meshgrid(1 : region_width, 1 : region_height); + + % Centroid + xc = region_width / 2; + yc = region_height / 2; + + % Make coordinates (polar), + % with origin (r = 0) at geometric centroid of the region + [~, r] = cart2pol(x - xc, y - yc); + + % Threshold the phase quality map using histogram equalization +% phase_quality_bw = (histeq(phase_quality_scaled, 2)); + phase_quality_bw = im2bw(phase_quality_scaled, 0.5); + +% Set border pixels to 1 to prevent connecting regions via the border + phase_quality_bw(1, :) = 1; + phase_quality_bw(end, :) = 1; + phase_quality_bw(:, 1) = 1; + phase_quality_bw(:, end) = 1; + + % Find properties of all the connected regions in the thresholded image. + phase_quality_region_props = regionprops(~phase_quality_bw, 'Centroid', 'PixelIdxList'); + + % Count the number of regions + num_regions = length(phase_quality_region_props); + + % Allocate a vector that will contain the median values + % of the radial coordinates of each region. + region_weighted_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the centroids of the detected regions + region_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the mean values of the radial coordinates + % of the pixels comprising the detected regions. + region_radius_median = zeros(num_regions, 1); + + % Measure the median radial coordinate of each region + for k = 1 : num_regions + + % Find the radial coordinate of the region's centroid. + % whose value is equal to the k'th iteration of this loop + region_centroid_subs = phase_quality_region_props(k).Centroid; + + % Radial coordinates of the pixels in the region + region_radius_median(k) = median(r(phase_quality_region_props(k).PixelIdxList)); + + % Convert to radial + region_centroid_radial(k) = sqrt((region_centroid_subs(1) - xc)^2 + (region_centroid_subs(2) - yc)^2); + + % Centroid coordinates weighted by the median radial coordinate of the pixels in the region + region_weighted_centroid_radial(k) = region_radius_median(k) * region_centroid_radial(k); + + end + + % Find the minimum + [~, region_weighted_centroid_idx] = min(region_weighted_centroid_radial); + + % Mask indices + mask_idx = phase_quality_region_props(region_weighted_centroid_idx).PixelIdxList; + + % Allocate a mask matrix + phase_mask_holder = zeros(region_height, region_width); + + % Apply the mask + phase_mask_holder(mask_idx) = 1; + + % Zero the border pixels just in cae + phase_mask_holder(:, 1) = 0 ; + phase_mask_holder(:, end) = 0; + phase_mask_holder(1, :) = 0; + phase_mask_holder(end, :) = 0; + + % Fill the holes + phase_mask_holder = imfill(phase_mask_holder); + + % Get the region properties again + phase_mask_region_props = regionprops(phase_mask_holder,... + 'PixelIdxList', 'MajorAxisLength', 'MinorAxisLength', ... + 'Orientation'); + + % Read the measured major axis of the ellipse + % that best fits the binary quality mask + ax_maj = phase_mask_region_props.MajorAxisLength/2; + + % Same for the minor axis + ax_min = phase_mask_region_props.MinorAxisLength/2; + + % Measure the orientation angle of + % the best fit ellipse and convert it to radians. + ax_angle = 1 * deg2rad(phase_mask_region_props.Orientation); + + % Rotate the coordinates for the elliptical Gaussian + x2 = (x - xc) * cos(ax_angle) - (y - yc) * sin(ax_angle); + y2 = (x - xc) * sin(ax_angle) + (y - yc) * cos(ax_angle); + + % Calculate the Gaussian function + elliptical_mask = zeros(size(phase_mask_holder)); + elliptical_mask((x2.^2 / ax_maj^2 + y2.^2 / ax_min^2) < 1) = 1; + + % Insert the cropped mask into the phase + % quality matrix, so that it's the same + % size as the original phase matrix. + PHASE_MASK = zeros(size(PHASE_QUALITY)); + PHASE_MASK(rad + 1 : end - rad - 1,... + rad + 1 : end - rad - 1) = elliptical_mask; + +end + + + + + + + + + + diff --git a/calculate_phase_mask_ellipse_3D.m b/calculate_phase_mask_ellipse_3D.m new file mode 100644 index 0000000..188c417 --- /dev/null +++ b/calculate_phase_mask_ellipse_3D.m @@ -0,0 +1,142 @@ +function [PHASE_MASK] = calculate_phase_mask_ellipse_3D(PHASE_QUALITY, KERNEL_RADIUS) + % This function computes a quality-based mask of the phase angle plane of a cross correlation. + + % Rename the kernel radius + rad = KERNEL_RADIUS; + + % Extract the central portion of the phase quality + phase_quality_interior = PHASE_QUALITY(rad + 1 : end - rad - 1, rad + 1 : end - rad - 1, rad + 1 : end - rad - 1); + + phase_quality_sub = (phase_quality_interior - ... + min(phase_quality_interior(:))); + + % Rescale the phase quality so that its range is [0, 1] + phase_quality_scaled = phase_quality_sub ./ max(phase_quality_sub(:)); + + % Measure size + [region_height, region_width, region_slice] = size(phase_quality_scaled); + + % Make coordinates (Cartesian) + [x, y, z] = meshgrid(1 : region_width, 1 : region_height, 1 : region_slice); + + % Centroid + xc = region_width / 2; + yc = region_height / 2; + zc = region_slice / 2; + + % Make coordinates (polar), + % with origin (r = 0) at geometric centroid of the region + [~, r] = cart2pol(x - xc, y - yc, z - zc); + + % Threshold the phase quality map using histogram equalization +% phase_quality_bw = (histeq(phase_quality_scaled, 2)); + phase_quality_bw = im2bw(phase_quality_scaled, 0.5); + +% Set border pixels to 1 to prevent connecting regions via the border + phase_quality_bw(1, :, :) = 1; + phase_quality_bw(end, :, :) = 1; + phase_quality_bw(:, 1, :) = 1; + phase_quality_bw(:, end, :) = 1; + phase_quality_bw(:, :, 1) = 1; + phase_quality_bw(:, :, end) = 1; + + % Find properties of all the connected regions in the thresholded image. + phase_quality_region_props = regionprops(~phase_quality_bw, 'Centroid', 'PixelIdxList'); + + % Count the number of regions + num_regions = length(phase_quality_region_props); + + % Allocate a vector that will contain the median values + % of the radial coordinates of each region. + region_weighted_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the centroids of the detected regions + region_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the mean values of the radial coordinates + % of the pixels comprising the detected regions. + region_radius_median = zeros(num_regions, 1); + + % Measure the median radial coordinate of each region + for k = 1 : num_regions + + % Find the radial coordinate of the region's centroid. + % whose value is equal to the k'th iteration of this loop + region_centroid_subs = phase_quality_region_props(k).Centroid; + + % Radial coordinates of the pixels in the region + region_radius_median(k) = median(r(phase_quality_region_props(k).PixelIdxList)); + + % Convert to radial + region_centroid_radial(k) = sqrt((region_centroid_subs(1) - xc)^2 + (region_centroid_subs(2) - yc)^2 + (region_centroid_subs(3) - zc)^2); + + % Centroid coordinates weighted by the median radial coordinate of the pixels in the region + region_weighted_centroid_radial(k) = region_radius_median(k) * region_centroid_radial(k); + + end + + % Find the minimum + [~, region_weighted_centroid_idx] = min(region_weighted_centroid_radial); + + % Mask indices + mask_idx = phase_quality_region_props(region_weighted_centroid_idx).PixelIdxList; + + % Allocate a mask matrix + phase_mask_holder = zeros(region_height, region_width, region_slice); + + % Apply the mask + phase_mask_holder(mask_idx) = 1; + + % Zero the border pixels just in cae + phase_mask_holder(:, 1, :) = 0 ; + phase_mask_holder(:, end, :) = 0; + phase_mask_holder(1, :, :) = 0; + phase_mask_holder(end, :, :) = 0; + phase_mask_holder(:, :, 1) = 0; + phase_mask_holder(:, :, end) = 0; + + % Fill the holes + phase_mask_holder = imfill(phase_mask_holder); + + % Get the region properties again + phase_mask_region_props = regionprops(phase_mask_holder,... + 'PixelIdxList', 'MajorAxisLength', 'MinorAxisLength', ... + 'Orientation'); + + % Read the measured major axis of the ellipse + % that best fits the binary quality mask + ax_maj = phase_mask_region_props.MajorAxisLength/2; + + % Same for the minor axis + ax_min = phase_mask_region_props.MinorAxisLength/2; + + % Measure the orientation angle of + % the best fit ellipse and convert it to radians. + ax_angle = 1 * deg2rad(phase_mask_region_props.Orientation); + + % Rotate the coordinates for the elliptical Gaussian + x2 = (x - xc) * cos(ax_angle) - (y - yc) * sin(ax_angle); + y2 = (x - xc) * sin(ax_angle) + (y - yc) * cos(ax_angle); + + % Calculate the Gaussian function + elliptical_mask = zeros(size(phase_mask_holder)); + elliptical_mask((x2.^2 / ax_maj^2 + y2.^2 / ax_min^2) < 1) = 1; + + % Insert the cropped mask into the phase + % quality matrix, so that it's the same + % size as the original phase matrix. + PHASE_MASK = zeros(size(PHASE_QUALITY)); + PHASE_MASK(rad + 1 : end - rad - 1,... + rad + 1 : end - rad - 1) = elliptical_mask; + +end + + + + + + + + + + diff --git a/calculate_phase_mask_gaussian.m b/calculate_phase_mask_gaussian.m new file mode 100644 index 0000000..8313881 --- /dev/null +++ b/calculate_phase_mask_gaussian.m @@ -0,0 +1,140 @@ +function [PHASE_MASK] = calculate_phase_mask_gaussian(PHASE_QUALITY, KERNEL_RADIUS) + % This function computes a quality-based mask of the phase angle plane of a cross correlation. + + % Rename the kernel radius + rad = KERNEL_RADIUS; + + % Extract the central portion of the phase quality + phase_quality_interior = PHASE_QUALITY(rad + 1 : end - rad - 1, rad + 1 : end - rad - 1); + + phase_quality_sub = (phase_quality_interior - ... + min(phase_quality_interior(:))); + + % Rescale the phase quality so that its range is [0, 1] + phase_quality_scaled = phase_quality_sub ./ max(phase_quality_sub(:)); + + % Measure size + [region_height, region_width] = size(phase_quality_scaled); + + % Make coordinates (Cartesian) + [x, y] = meshgrid(1 : region_width, 1 : region_height); + + % Centroid + xc = region_width / 2; + yc = region_height / 2; + + % Make coordinates (polar), + % with origin (r = 0) at geometric centroid of the region + [~, r] = cart2pol(x - xc, y - yc); + + % Threshold the phase quality map using histogram equalization +% phase_quality_bw = (histeq(phase_quality_scaled, 2)); + phase_quality_bw = im2bw(phase_quality_scaled, 0.5); + +% Set border pixels to 1 to prevent connecting regions via the border + phase_quality_bw(1, :) = 1; + phase_quality_bw(end, :) = 1; + phase_quality_bw(:, 1) = 1; + phase_quality_bw(:, end) = 1; + + % Find properties of all the connected regions in the thresholded image. + phase_quality_region_props = regionprops(~phase_quality_bw, 'Centroid', 'PixelIdxList'); + + % Count the number of regions + num_regions = length(phase_quality_region_props); + + % Allocate a vector that will contain the median values + % of the radial coordinates of each region. + region_weighted_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the centroids of the detected regions + region_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the mean values of the radial coordinates + % of the pixels comprising the detected regions. + region_radius_median = zeros(num_regions, 1); + + % Measure the median radial coordinate of each region + for k = 1 : num_regions + + % Find the radial coordinate of the region's centroid. + % whose value is equal to the k'th iteration of this loop + region_centroid_subs = phase_quality_region_props(k).Centroid; + + % Radial coordinates of the pixels in the region + region_radius_median(k) = median(r(phase_quality_region_props(k).PixelIdxList)); + + % Convert to radial + region_centroid_radial(k) = sqrt((region_centroid_subs(1) - xc)^2 + (region_centroid_subs(2) - yc)^2); + + % Centroid coordinates weighted by the median radial coordinate of the pixels in the region + region_weighted_centroid_radial(k) = region_radius_median(k) * region_centroid_radial(k); + + end + + % Find the minimum + [~, region_weighted_centroid_idx] = min(region_weighted_centroid_radial); + + % Mask indices + mask_idx = phase_quality_region_props(region_weighted_centroid_idx).PixelIdxList; + + % Allocate a mask matrix + phase_mask_holder = zeros(region_height, region_width); + + % Apply the mask + phase_mask_holder(mask_idx) = 1; + + % Zero the border pixels just in cae + phase_mask_holder(:, 1) = 0 ; + phase_mask_holder(:, end) = 0; + phase_mask_holder(1, :) = 0; + phase_mask_holder(end, :) = 0; + + % Fill the holes + phase_mask_holder = imfill(phase_mask_holder); + + % Get the region properties again + phase_mask_region_props = regionprops(phase_mask_holder,... + 'PixelIdxList', 'MajorAxisLength', 'MinorAxisLength', ... + 'Orientation'); + + % Read the measured major axis of the ellipse + % that best fits the binary quality mask + ax_maj = phase_mask_region_props.MajorAxisLength/2; + + % Same for the minor axis + ax_min = phase_mask_region_props.MinorAxisLength/2; + + % Measure the orientation angle of + % the best fit ellipse and convert it to radians. + ax_angle = 1 * deg2rad(phase_mask_region_props.Orientation); + + % Gaussian standard deviations as fractions + % of the major and minor axes of the ellipse fit + std_maj = ax_maj / 1.00; + std_min = ax_min / 1.00; + + % Rotate the coordinates for the elliptical Gaussian + x2 = (x - xc) * cos(ax_angle) - (y - yc) * sin(ax_angle); + y2 = (x - xc) * sin(ax_angle) + (y - yc) * cos(ax_angle); + + % Calculate the Gaussian function + gaussian_mask = exp(-(x2.^2)/(std_maj^2) - (y2.^2) / (std_min^2)); + + % Insert the cropped mask into the phase + % quality matrix, so that it's the same + % size as the original phase matrix. + PHASE_MASK = zeros(size(PHASE_QUALITY)); + PHASE_MASK(rad + 1 : end - rad - 1, rad + 1 : end - rad - 1) = gaussian_mask; + +end + + + + + + + + + + diff --git a/calculate_phase_mask_gaussian_3D.m b/calculate_phase_mask_gaussian_3D.m new file mode 100644 index 0000000..7776bea --- /dev/null +++ b/calculate_phase_mask_gaussian_3D.m @@ -0,0 +1,164 @@ +function [PHASE_MASK] = calculate_phase_mask_gaussian_3D(PHASE_QUALITY, KERNEL_RADIUS) + % This function computes a quality-based mask of the phase angle plane of a cross correlation. + + % Rename the kernel radius + rad = KERNEL_RADIUS; + + % Extract the central portion of the phase quality + phase_quality_interior = PHASE_QUALITY(rad + 1 : end - rad - 1, rad + 1 : end - rad - 1, rad + 1 : end - rad - 1); + + phase_quality_sub = (phase_quality_interior - ... + min(phase_quality_interior(:))); + + % Rescale the phase quality so that its range is [0, 1] + phase_quality_scaled = phase_quality_sub ./ max(phase_quality_sub(:)); + + % Measure size + [region_height, region_width, region_slice] = size(phase_quality_scaled); + + % Make coordinates (Cartesian) + [x, y, z] = meshgrid(1 : region_width, 1 : region_height, 1 : region_slice); + + % Centroid + xc = region_width / 2; + yc = region_height / 2; + zc = region_slice / 2; + + % Make coordinates (polar), + % with origin (r = 0) at geometric centroid of the region + [~, r] = cart2pol(x - xc, y - yc, z - zc); + + % Threshold the phase quality map using histogram equalization +% phase_quality_bw = (histeq(phase_quality_scaled, 2)); + + + %%%%%%% addition for thresholding in 3D %%%%%%%% + + +% level = graythresh(phase_quality_scaled); + phase_quality_bw = zeros(region_height, region_width, region_slice); + for slice_k = 1:region_slice + phase_quality_bw(1:region_height, 1:region_width, slice_k) = im2bw(phase_quality_scaled(:,:,slice_k), 0.5); %ORIGINALLY 0.5 + end +% Set border pixels to 1 to prevent connecting regions via the border + phase_quality_bw(1, :, :) = 1; + phase_quality_bw(end, :, :) = 1; + phase_quality_bw(:, 1, :) = 1; + phase_quality_bw(:, end, :) = 1; + phase_quality_bw(:, :, 1) = 1; + phase_quality_bw(:, :, end) = 1; + + % Find properties of all the connected regions in the thresholded image. + phase_quality_region_props = regionprops(~phase_quality_bw, 'Centroid', 'PixelIdxList'); + + % Count the number of regions + num_regions = length(phase_quality_region_props); + + % Allocate a vector that will contain the median values + % of the radial coordinates of each region. + region_weighted_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the centroids of the detected regions + region_centroid_radial = zeros(num_regions, 1); + + % Allocate a vector to hold the mean values of the radial coordinates + % of the pixels comprising the detected regions. + region_radius_median = zeros(num_regions, 1); + + % Measure the median radial coordinate of each region + for k = 1 : num_regions + + % Find the radial coordinate of the region's centroid. + % whose value is equal to the k'th iteration of this loop + region_centroid_subs = phase_quality_region_props(k).Centroid; + + % Radial coordinates of the pixels in the region + region_radius_median(k) = median(r(phase_quality_region_props(k).PixelIdxList)); + + % Convert to radial + region_centroid_radial(k) = sqrt((region_centroid_subs(1) - xc)^2 + (region_centroid_subs(2) - yc)^2 + (region_centroid_subs(3) - zc)^2); + + % Centroid coordinates weighted by the median radial coordinate of the pixels in the region + region_weighted_centroid_radial(k) = region_radius_median(k) * region_centroid_radial(k); + + end + + % Find the minimum + [~, region_weighted_centroid_idx] = min(region_weighted_centroid_radial); + + % Mask indices + mask_idx = phase_quality_region_props(region_weighted_centroid_idx).PixelIdxList; + + % Allocate a mask matrix + phase_mask_holder = zeros(region_height, region_width, region_slice); + + % Apply the mask + phase_mask_holder(mask_idx) = 1; + + % Zero the border pixels just in cae + phase_mask_holder(:, 1, :) = 0 ; + phase_mask_holder(:, end, :) = 0; + phase_mask_holder(1, :, :) = 0; + phase_mask_holder(end, :, :) = 0; + phase_mask_holder(:, :, 1) = 0; + phase_mask_holder(:, :, end) = 0; + + % Fill the holes + phase_mask_holder = imfill(phase_mask_holder); + + %%% test for 3D by Brian %%% + stacked_phase_mask_holder = zeros(size(phase_mask_holder(:,:,1))); + + for stack_slice = 1:region_slice; + + stacked_phase_mask_holder = stacked_phase_mask_holder + phase_mask_holder(:,:,stack_slice); + end + + stacked_phase_mask_holder = im2bw(stacked_phase_mask_holder, 0.5); + + % Get the region properties again + phase_mask_region_props = regionprops(stacked_phase_mask_holder,... + 'PixelIdxList', 'MajorAxisLength', 'MinorAxisLength', ... + 'Orientation'); + + % Read the measured major axis of the ellipse + % that best fits the binary quality mask + ax_maj = phase_mask_region_props.MajorAxisLength/2; + + % Same for the minor axis + ax_min = phase_mask_region_props.MinorAxisLength/2; + + % Measure the orientation angle of + % the best fit ellipse and convert it to radians. + ax_angle = 1 * deg2rad(phase_mask_region_props.Orientation); + + % Gaussian standard deviations as fractions + % of the major and minor axes of the ellipse fit + std_maj = ax_maj / 1.00; + std_min = ax_min / 1.00; + + % Rotate the coordinates for the elliptical Gaussian + x2 = (x - xc) * cos(ax_angle) - (y - yc) * sin(ax_angle); + y2 = (x - xc) * sin(ax_angle) + (y - yc) * cos(ax_angle); + z2 = (z - zc); + % Calculate the Gaussian function + gaussian_mask = exp(-(x2.^2)/(std_maj^2) - (y2.^2) / (std_min^2) - (z2.^2) / (std_min^2)); + + % Insert the cropped mask into the phase + % quality matrix, so that it's the same + % size as the original phase matrix. + PHASE_MASK = zeros(size(PHASE_QUALITY)); + PHASE_MASK(rad + 1 : end - rad - 1,... + rad + 1 : end - rad - 1, rad + 1 : end - rad - 1) = gaussian_mask; + +end + + + + + + + + + + diff --git a/calculate_phase_quality.m b/calculate_phase_quality.m new file mode 100644 index 0000000..21e2713 --- /dev/null +++ b/calculate_phase_quality.m @@ -0,0 +1,57 @@ +function phase_quality_array = calculate_phase_quality(wrapped_phase_angle, kernel_radius) + + +% Calculate the horizontal and vertical wrapped phase differences +% phase_diff_rows = phase_diff_kernel_wrapped(wrapped_phase_angle, 1, 'sobel'); +% phase_diff_cols = phase_diff_kernel_wrapped(wrapped_phase_angle, 2, 'sobel'); + +phase_diff_rows = wrapped_phase_difference(wrapped_phase_angle, 1); +phase_diff_cols = wrapped_phase_difference(wrapped_phase_angle, 2); + +% Measure the array size +% and subtract 2 from each dimensions due to the differencing +[M, N] = size(wrapped_phase_angle); + +% Array height and width +array_height = M; +array_width = N; + +% Min and maxes of rows and columns +col_min = 1 + kernel_radius; +col_max = array_width - kernel_radius; + +row_min = 1 + kernel_radius; +row_max = array_height - kernel_radius; + +% Allocate the phase quality array +phase_quality_array = zeros(array_height, array_width); + +% Loop over the rows +for m = row_min : row_max + for n = col_min : col_max + + % Extract the sub-region of phase differences + row_diffs = phase_diff_rows(m - kernel_radius : m + kernel_radius,... + n - kernel_radius : n + kernel_radius); + + col_diffs = phase_diff_cols(m - kernel_radius : m + kernel_radius,... + n - kernel_radius : n + kernel_radius); + + % Standard deviations + % The 1 indicates normalization by N, not by N-1 + row_diff_std = std(row_diffs(:), 1); + col_diff_std = std(col_diffs(:), 1); + + % Populate the phase quality array + phase_quality_array(m, n) = row_diff_std + col_diff_std; + + end +end + +end + + + + + + diff --git a/calculate_phase_quality_3D.m b/calculate_phase_quality_3D.m new file mode 100644 index 0000000..9a4a6a0 --- /dev/null +++ b/calculate_phase_quality_3D.m @@ -0,0 +1,71 @@ +function phase_quality_array = calculate_phase_quality_3D(wrapped_phase_angle, kernel_radius) + + +% Calculate the horizontal and vertical wrapped phase differences +% phase_diff_rows = phase_diff_kernel_wrapped(wrapped_phase_angle, 1, 'sobel'); +% phase_diff_cols = phase_diff_kernel_wrapped(wrapped_phase_angle, 2, 'sobel'); + +phase_diff_rows = wrapped_phase_difference(wrapped_phase_angle, 1); +phase_diff_cols = wrapped_phase_difference(wrapped_phase_angle, 2); +phase_diff_slice = wrapped_phase_difference(wrapped_phase_angle, 3); + +% Measure the array size +% and subtract 2 from each dimensions due to the differencing +[M, N, O] = size(wrapped_phase_angle); + +% Array height and width +array_height = M; +array_width = N; +array_slice = O; + +% Min and maxes of rows and columns +col_min = 1 + kernel_radius; +col_max = array_width - kernel_radius; + +row_min = 1 + kernel_radius; +row_max = array_height - kernel_radius; + +slice_min = 1 + kernel_radius; +slice_max = array_slice - kernel_radius; + +% Allocate the phase quality array +phase_quality_array = zeros(array_height, array_width, array_slice); + +% Loop over the rows +tic + +for m = row_min : row_max + for n = col_min : col_max + for o = slice_min : slice_max + + % Extract the sub-region of phase differences + row_diffs = phase_diff_rows(m - kernel_radius : m + kernel_radius,... + n - kernel_radius : n + kernel_radius, o - kernel_radius : o + kernel_radius); + + col_diffs = phase_diff_cols(m - kernel_radius : m + kernel_radius,... + n - kernel_radius : n + kernel_radius, o - kernel_radius : o + kernel_radius); + + slice_diffs = phase_diff_slice(m - kernel_radius : m + kernel_radius,... + n - kernel_radius : n + kernel_radius, o - kernel_radius : o + kernel_radius); + % Standard deviations + % The 1 indicates normalization by N, not by N-1 + row_diff_std = std(row_diffs(:), 1); + col_diff_std = std(col_diffs(:), 1); + slice_diff_std = std(slice_diffs(:), 1); + + % Populate the phase quality array + phase_quality_array(m, n, o) = row_diff_std + col_diff_std + slice_diff_std; + end + + end +end + +toc + +end + + + + + + diff --git a/calculate_phase_quality_mex.mexw64 b/calculate_phase_quality_mex.mexw64 new file mode 100644 index 0000000..75d8590 Binary files /dev/null and b/calculate_phase_quality_mex.mexw64 differ diff --git a/clsm_img_gen.m b/clsm_img_gen.m new file mode 100644 index 0000000..7c120c2 --- /dev/null +++ b/clsm_img_gen.m @@ -0,0 +1,137 @@ +function [Img1,Img2]=clsm_img_gen(U,V,W) +%clsm_generator function generates synthetic confocal laser scanning microscopy images (1D/2D/3D) +%Developed/updated by Brian Jun (2016) based on synthetic PIV image gnerator version previously updated by Sam Raben (2011) +% +% +%Inputs +%Assign uniform velocity U,V,W (pixels/frame) +%%%%%%%%%%%%%%%%%%%%%%%%%% +% physical parameters +%%%%%%%%%%%%%%%%%%%%%%%%%% + d = 2.2; %particle diameter (pix) + Lx = 32; %image size width x (pix) + Ly = 32; %image size height y (pix) + Lz = 32; %image depth z (pix) + C = 0.001; %density (particles/pix^3) + Vout = 24; %distance outside image to generate particles (pix) + distr = 'uniform'; % normal or uniform + + pix_res = 0.31; % pixel resolution (um) + D = 5.01*10^-6; %Diffusion coefficient of particles (um^2/us) + T_s = 0.01; %pixel dwell time of the scanner (us) +%%%%%%%%%%%%%%%%%%%%%%%%%% + + if numel(d) ~= 1 + d_ave = mean(d); + N=round(C*(Ly+2*d_ave)*(Lx+2*d_ave+2*Vout)*(Lz+2*d_ave)); + switch lower(distr) + case 'normal' + d = d_ave + 3.*randn(N,1); + case 'uniform' + d = rand(N,1).*(d(2)-d(1)) + d(1); + end + elseif nargin < 11 + N=round(C*(Ly+2*d)*(Lx+2*d+2*Vout)*(Lz+2*d)); + d = repmat(d,N,1); + else + N = numel(X); + d = repmat(d,N,1); + end + + if nargin < 11 + X = (Lx-16)*rand(N,1)+8; + Y = (Ly-16)*rand(N,1)+8; + Z = (Lz-16)*rand(N,1)+8; + end + + N = numel(X); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %Generate Image 1 % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + D_p = D*T_s/(pix_res^2); %Diffusity per pixel scan (pix^2) + sigma = sqrt(2*D_p); + + Img1 = zeros(Lx+1,Ly+1,Lz+1); + I1 = zeros(Lx+1,Ly+1,Lz+1); + + d_one = d(1); + + for t_3 = 1:1:Lz + for t_2 = 1:1:Ly + for t_1 = 1:1:Lx + x0 = U/(Lx*Ly*Lz) + sigma*randn(size(X,1),size(X,2)); + y0 = V/(Lx*Ly*Lz) + sigma*randn(size(Y,1),size(Y,2)); + z0 = W/(Lx*Ly*Lz) + sigma*randn(size(Z,1),size(Z,2)); + + X = X + x0; + Y = Y + y0; + Z = Z + z0; + + for p=1:N + if X(p)<=t_1+4 && X(p)>=t_1-4 + for i=floor(X(p)-0.75*d_one):ceil(X(p)+0.75*d_one) + for j=floor(Y(p)-0.75*d_one):ceil(Y(p)+0.75*d_one) + for k=floor(Z(p)-0.75*d_one):ceil(Z(p)+0.75*d_one) + + if i>=0 && i<=Lx && j>=0 && j<=Ly && k>=0 && k<=Lz && sqrt((i-X(p))^2+(j-Y(p))^2+(k-Z(p))^2)<0.75*d_one + I1(i+1,j+1,k+1) = I1(i+1,j+1,k+1) + 0.4*exp(-8*((i-X(p))^2+(j-Y(p))^2+(k-Z(p))^2)/d_one ^2); + end + end + end + end + end + end + + Img1(t_1, t_2, t_3) = I1(t_1, t_2, t_3); + I1 = zeros(Lx+1,Ly+1,Lz+1); + + end + end + end + + Iimg = 1; + Img1 = permute(Img1,[2 1 3])/Iimg; + + +%%%%%%%%%%GENERATE THE SECOND IMAGE %%%%%%%%%%%%%%%%% + + Img2 = zeros(Lx+1,Ly+1,Lz+1); + I1 = zeros(Lx+1,Ly+1,Lz+1); + parlist = zeros(N,1); + + for t_3 = 1:1:Lz + for t_2 = 1:1:Ly + for t_1 = 1:1:Lx + + x0 = U/(Lx*Ly*Lz) + sigma*randn(size(X,1),size(X,2)); + y0 = V/(Lx*Ly*Lz) + sigma*randn(size(Y,1),size(Y,2)); + z0 = W/(Lx*Ly*Lz) + sigma*randn(size(Z,1),size(Z,2)); + + X = X + x0; + Y = Y + y0; + Z = Z + z0; + + for p=1:N + if X(p)<=t_1+4 && X(p)>=t_1-4 + for i=floor(X(p)-0.75*d_one):ceil(X(p)+0.75*d_one) + for j=floor(Y(p)-0.75*d_one):ceil(Y(p)+0.75*d_one) + for k=floor(Z(p)-0.75*d_one):ceil(Z(p)+0.75*d_one) + if i>=0 && i<=Lx && j>=0 && j<=Ly && k>=0 && k<=Lz && sqrt((i-X(p))^2+(j-Y(p))^2+(k-Z(p))^2)<0.75*d_one + parlist(p) = p; + I1(i+1,j+1,k+1) = I1(i+1,j+1,k+1) + 0.4*exp(-8*((i-X(p))^2+(j-Y(p))^2+(k-Z(p))^2)/d_one ^2); + end + end + end + end + end + end + Img2(t_1, t_2, t_3) = I1(t_1, t_2, t_3); + I1 = zeros(Lx+1,Ly+1,Lz+1); + end + end + end + + Iimg = 1; + Img2 = permute(Img2,[2 1 3])/Iimg; diff --git a/findGaussianWidth.m b/findGaussianWidth.m new file mode 100644 index 0000000..8ef5f46 --- /dev/null +++ b/findGaussianWidth.m @@ -0,0 +1,140 @@ +function [STDY, STDX] = findGaussianWidth(IMAGESIZEY, IMAGESIZEX, WINDOWSIZEY, WINDOWSIZEX) +% FINDGAUSSIANWIDTH determines the standard deviation of a normalized Gaussian function whose +% area is approximately equal to that of a top-hat function of the desired +% effective window resolution. +% +% John says that the volume under a this Gaussian window will not be +% exactly the volume under the square window. Check into this. +% +% INPUTS +% xregion = Width of each interrogation region (pixels) +% yregion = Height of each interrogation region (pixels) +% xwin = Effective window resolution in the x-direction (pixels) +% ywin = Effective window resolution in the y-direction (pixels) +% +% OUTPUTS +% sx = standard deviation of a normalized Gaussian for the x-dimension of the window (pixels) +% sy = standard deviation of a normalized Gaussian for the y-dimension of the window (pixels) +% +% EXAMPLE +% xregion = 32; +% yregion = 32; +% xwin = 16; +% ywin = 16; +% [sx sy] = findGaussianWidth(xregion, yregion, xwin, ywin); +% +% SEE ALSO +% + +% Initial guess for standard deviaitons are half the respective window sizes +% sx = xwin/2; +% sy = ywin/2; + +STDX = 50 * WINDOWSIZEX; +STDY = 50 * WINDOWSIZEY; + +% Generate x and y domains +x = - IMAGESIZEX/2 : IMAGESIZEX/2; +y = - IMAGESIZEY/2 : IMAGESIZEY/2; + +% Generate normalized zero-mean Gaussian windows +xgauss = exp(-(x).^2/(2 * STDX^2)); +ygauss = exp(-(y).^2/(2 * STDY^2)); + +% Calculate areas under gaussian curves +xarea = trapz(x, xgauss); +yarea = trapz(y, ygauss); + +if WINDOWSIZEX < xarea + + % Calculate initial errors of Gaussian windows with respect to desired + % effective window resolution + xerr = abs(1 - xarea / WINDOWSIZEX); + + % Initialize max and min values of standard deviation for Gaussian x-window + sxmax = 100 * IMAGESIZEX; + sxmin = 0; + + % Iteratively determine the standard deviation that gives the desired + % effective Gaussian x-window resolution + + % Loop while the error of area under curve is above the specified error tolerance + while xerr > 1E-5 + + % If the area under the Gaussian curve is less than that of the top-hat window + if xarea < WINDOWSIZEX + % Increase the lower bound on the standard deviation + sxmin = sxmin + (sxmax - sxmin) / 2; + else + % Otherwise, increase the upper bound on the standard deviation + sxmax = sxmin + (sxmax - sxmin) / 2; + end + + % Set the standard deviation to halfway between its lower and upper bounds + STDX = sxmin + (sxmax - sxmin) / 2; + + % Generate a Gaussian curve with the specified standard deviation + xgauss = exp(-(x).^2/(2 * STDX^2)); + + % Calculate the area under this Gaussian curve via numerical + % integration using the Trapezoidal rule + xarea = trapz(x,xgauss); + + % Calculate the error of the area under the Gaussian curve with + % respect to the desired area + xerr = abs(1 - xarea / WINDOWSIZEX); + + end + +end + +if WINDOWSIZEY < yarea + + yerr = abs(1 - yarea / WINDOWSIZEY); + + % Initialize max and min values of standard deviation for Gaussian y-window + symax = 100 * IMAGESIZEY; + symin = 0; + + % Iteratively determine the standard deviation that gives the desired + % effective Gaussian y-window resolution + + % Loop while the error of area under curve is above the specified error tolerance + while yerr > 1E-5 + + % If the area under the Gaussian curve is less than that of the top-hat window + if yarea < WINDOWSIZEY + % Increase the lower bound on the standard deviation + symin = symin + (symax - symin) / 2; + else + % Otherwise, increase the upper bound on the standard deviation + symax = symin + (symax - symin) / 2; + end + + % Set the standard deviation to halfway between its lower and upper bounds + STDY = symin + (symax - symin) / 2; + + % Generate a Gaussian curve with the specified standard deviation + ygauss = exp(-(y).^2/(2 * STDY^2)); + + % Calculate the area under this Gaussian curve via + % numerical integration using the Trapezoidal rule + yarea = trapz(y, ygauss); + + % Calculate the error of the area under the Gaussian curve with + % respect to the desired area + yerr = abs(1 - yarea / WINDOWSIZEY); + + end + +end + +end + + + + + + + + diff --git a/gaussianWindowFilter.m b/gaussianWindowFilter.m new file mode 100644 index 0000000..8e26786 --- /dev/null +++ b/gaussianWindowFilter.m @@ -0,0 +1,61 @@ +function WINDOW = gaussianWindowFilter(DIMENSIONS, WINDOWSIZE, WINDOWTYPE) +% gaussianWindowFilter(DIMENSIONS, WINDOWSIZE, WINDOWTYPE) creates a 2-D +% gaussian window +% +% INPUTS +% DIMENSIONS = 2 x 1 Vector specifying the Dimensions (in rows and columns) of the gaussian window filter +% WINDOWSIZE = 2 x 1 vector specifying the effective window resolution of +% the gaussian window. WINDOWSIZE can either specify the resolution in +% pixels or as a fraction of the filter dimensions. This option is +% controlled by the input WINDOWTYPE. +% WINDOWTYPE = String specifying whether WINDOWSIZE specifies a +% resolution in pixels ('pixels') or as a fraction of the window +% dimensions ('fraction'). +% +% OUTPUTS +% WINDOW = 2-D gaussian window +% +% SEE ALSO +% findwidth + +% Default to an absolute size window type +if nargin < 3 + WINDOWTYPE = 'fraction'; +end + +% Signal height and width +height = DIMENSIONS(1); +width = DIMENSIONS(2); + +% Determine whether window size is an absolute size or a fraction of the +% window dimensions +if strcmp(WINDOWTYPE, 'fraction') + windowSizeX = width .* WINDOWSIZE(2); + windowSizeY = height .* WINDOWSIZE(1); +elseif strcmp(WINDOWTYPE, 'pixels'); + windowSizeX = WINDOWSIZE(2); + windowSizeY = WINDOWSIZE(1); +else + error('Invalid window type "%s"\n', WINDOWTYPE); +end + +% Standard deviations +[sy, sx] = findGaussianWidth(height, width, windowSizeY, windowSizeX); + +% Calculate center of signal +xo = round(width/2); +yo = round(height/2); + +% Create grid of x,y positions to hold gaussian filter data +[x,y] = meshgrid(1:width, 1:height); + +% Calculate gaussian distribution (X) +WindowX = exp( - ((x - xo).^2 / (2 * sx^2))); + +% Calculate gaussian distribution (Y) +WindowY = exp( - ((y - yo).^2 / (2 * sy^2))); + +% 2-D Gaussian Distribution +WINDOW = WindowX .* WindowY; + +end \ No newline at end of file diff --git a/phaseOnlyFilter.m b/phaseOnlyFilter.m new file mode 100644 index 0000000..8d67938 --- /dev/null +++ b/phaseOnlyFilter.m @@ -0,0 +1,22 @@ +function PHASE_ONLY_CORRELATION_PLANE = phaseOnlyFilter(COMPLEX_CROSS_CORRELATION_PLANE) +% PHASE_ONLY_CORRELATION_PLANE = phaseOnlyFilter(COMPLEX_CROSS_CORRELATION_PLANE) +% performs phase-only filtering of a spectral-domain cross correlation signal + +% Calculate the spectral magnitude of the complex cross correlation +% in the frequency domain +spectral_magnitude = sqrt(COMPLEX_CROSS_CORRELATION_PLANE .* ... + conj(COMPLEX_CROSS_CORRELATION_PLANE)); + +% If the value of the cross correlation is zero, set the magnitude to 1. +% This avoids division by zero when calculating the phase correlation +% below (i.e., when dividing the cross correlation by its magnitude). +spectral_magnitude(COMPLEX_CROSS_CORRELATION_PLANE == 0 ) = 1; + +% Divide cross correlation by its nonzero magnitude to extract the phase information +PHASE_ONLY_CORRELATION_PLANE = COMPLEX_CROSS_CORRELATION_PLANE ./ ... + spectral_magnitude; + +end + + + diff --git a/phase_quality_mask_ensemble.m b/phase_quality_mask_ensemble.m new file mode 100644 index 0000000..1eac8fc --- /dev/null +++ b/phase_quality_mask_ensemble.m @@ -0,0 +1,152 @@ +%Phase quality mask code developed by Matt Giarra (2016) + + +% Input data directory +input_dir = 'Enter here'; + +% Input data base name +input_base_name = 'Enter here'; + +% Input data number format +num_format = '%04d'; + +% Input data extension +input_extension = '.tif'; + +% Start image +start_image = 1; + +% End image +end_image = 1000; + +% Frame step +frame_step = 1; + +% Correlation step +correlation_step = 1; + +% Region size +region_height = 512; +region_width = 512; + +% Window fractions +window_fraction = [0.5, 0.5]; + +% Correlation method +correlation_method = 'spc'; + +% Cutoff radius +rc_outer = 8; +rc_inner = 3; + +% Coordinates +[x, y] = meshgrid(1 : region_width, 1 : region_height); + +% Centroid +xc = region_height / 2 + 1; +yc = region_width / 2 + 1; + +% Angular coordinates +[~, r] = cart2pol(x - xc, y - yc); + +% Weighting filter +plane_fit_weights = ones(size(x)); +plane_fit_weights(r > rc_outer) = 0; +plane_fit_weights(r < rc_inner) = 0; + +% First image numbers +image_numbers_01 = start_image : frame_step : end_image; + +% Second image numbers +image_numbers_02 = image_numbers_01 + correlation_step; + +% Number of images +num_images = length(image_numbers_01); + +% Create a region window +region_window = gaussianWindowFilter([region_height, region_width], ... + window_fraction, 'fraction'); + +% Allocate a correlation plane +cross_correlation = zeros(region_height, region_width, 'double'); + +% Make file paths +for k = 1 : num_images + + % File name of the first images + file_name_01 = [input_base_name num2str(image_numbers_01(k), num_format) input_extension]; + + % File name of the second images + file_name_02 = [input_base_name num2str(image_numbers_02(k), num_format) input_extension]; + + % File paths + file_path_01{k} = fullfile(input_dir, file_name_01); + file_path_02{k} = fullfile(input_dir, file_name_02); + +end + +% Do the correlations +for k = 1 : num_images + + % Check existence of both filepaths + if exist(file_path_01{k}, 'file') && exist(file_path_02{k}, 'file') + + % Inform the user + fprintf('Processing image %d of %d...\n', k, num_images); + + % Read the first image + img_01 = double(imread(file_path_01{k})); + + % Read the second image + img_02 = double(imread(file_path_02{k})); + + % Extract the parts of the image used in the correlation + region_01 = img_01(1:region_height,1:region_width) .* region_window; %img_01(region_rows, region_cols) .* region_window; + region_02 = img_02(1:region_height,1:region_width) .* region_window;%img_02(region_rows, region_cols) .* region_window; + + % Take the FFT of both regions + ft_01 = fftn(region_01, [region_height, region_width]); + ft_02 = fftn(region_02, [region_height, region_width]); + + % Conjugate-multiply the regions' Fourier Transforms + % to produce the complex cross correlation + cross_correlation = cross_correlation + ft_02 .* conj(ft_01); + + end + +end + +Sx = region_width; +Sy = region_height; +D = [1,1]*2.8; +Peaklocator = 1; +Peakswitch = 0; +cnorm = ones(Sy,Sx); +fftindy = [ceil(Sy/2)+1:Sy 1:ceil(Sy/2)]; +fftindx = [ceil(Sx/2)+1:Sx 1:ceil(Sx/2)]; + +% Extract the phase from the ensemble cross correlation plane +spectral_phase_plane = phaseOnlyFilter(cross_correlation); + +% % Zero the parts of the phase plane that are outsize the cutoff radius +% spectral_phase_plane(r > rc_outer) = 0; + +% Extract the phase angle from the phase plane +%phase_angle_plane = angle(spectral_phase_plane); + +phase_angle = fftshift(angle(spectral_phase_plane)); +kernel_radius = 3; +phase_quality = calculate_phase_quality_mex(phase_angle, kernel_radius); +phase_mask = calculate_phase_mask_ellipse... + (phase_quality, kernel_radius); + +filtered_spectral_phase_plane = fftshift(phase_mask).*spectral_phase_plane; + + +G = ifftn(filtered_spectral_phase_plane,'symmetric'); +G = G(fftindy,fftindx); +G = abs(G); + +%subpixel estimation +[U,V,dX,dY]=subpixel(G,Sx,Sy,cnorm,Peaklocator,Peakswitch,D); + diff --git a/piv_3d_10.m b/piv_3d_10.m new file mode 100644 index 0000000..832e231 --- /dev/null +++ b/piv_3d_10.m @@ -0,0 +1,1552 @@ +function piv_3d_10(piv_parameters); +% This function is designed to process two or three dimensional PIV data using +% the parameters specified by the 'piv_parameters' data structure. +% +% This code is based upon the code 'basic_3d_rpc_processing_05.m'. +% +% Updates on previous versions: +% +% Version 08 +% +% Saves the PIV processing parameters data structure to the PIV +% vector field write directory for later reference. +% +% Version 09 +% +% Added the ability to use several different outlier vector +% replacement methods including a local mean interpolation (default +% in previous versions), a Laplacian interpolation, and a Delaunay +% interpolation. +% +% Version 10 +% +% This updated the order of the spatial and temporal dimensions to be +% consistent with other PIV software versions and to minimize the +% number of overhead operations (since MATLAB prefers some orders +% over others). +% +% The ability to import other data formats was also included. The +% code can now import 2D images in one of the following formats: +% 'bmp', 'jpg', 'jp2', 'png', 'ppm', 'tif', 'mat', 'dcm', 'cdf', or +% 'h5' while 3D images may be one of the following formats: 'mat', +% 'dcm', 'cdf', or 'h5'. +% +% Any NaN velocity vector values that are produced during the +% correlations (likely due to simulated data not being seeded within +% a window) are considered outliers and replaced accordingly. +% +% Things to work on, modify, fix, et cetera: +% +% Change the code that is optimized for C to code that is optimized for +% matlab. +% +% Check on the speed of the if-then statements added to make the code run +% for 2D images as well. +% +% Add peak fitting algorithms besides the 3-point Gaussian. +% +% Add the ability to record multiple correlation peaks, specifically to +% be used during validation. Possibly also record the peak ratio. +% +% Add the ability to do multiple iterations of validation and replacement +% for each pass. +% +% Add in the ability to perform deform correlations. +% +% Authors: Rod La Foy +% First Written On: 31 October 2014 +% Last Modified On: 21 November 2014 + +% This extracts the vector field saving directory from the data processing structure +vector_field_write_directory=piv_parameters.general_parameters.vector_field_write_directory; +% This saves a copy of the processing parameters to the vector field write +% directory for future reference +save([vector_field_write_directory,'piv_processing_parameters.mat'],'piv_parameters'); + +% This extracts the number of passes to perform +pass_number=piv_parameters.general_parameters.pass_number; + +% This is the initial image frame to load +initial_image_frame=piv_parameters.general_parameters.initial_image_frame; +% This is the final image frame to load +final_image_frame=piv_parameters.general_parameters.final_image_frame; +% This is the number of frames to step each iteration +image_frame_step=piv_parameters.general_parameters.image_frame_step; +% This is the number of frames to step between correlation frame pairs +image_correlation_step=piv_parameters.general_parameters.image_correlation_step; + +% This is the Boolean value stating whether to perform window deformation +perform_window_deformation=piv_parameters.general_parameters.perform_window_deformation; +% If window deformation is to be performed, this loads the window +% deformation parameters +if perform_window_deformation; + % This is the minimum number of window deformation operations to + % perform + window_deformation_iteration_min=piv_parameters.general_parameters.window_deformation_iteration_min; + % This is the maximum number of window deformation operations to + % perform + window_deformation_iteration_max=piv_parameters.general_parameters.window_deformation_iteration_max; + % This is the convergence threshhold of the window deformation process + window_deformation_threshhold=piv_parameters.general_parameters.window_deformation_threshhold; +end; + +% This is a Boolean value stating whether to perform pyramid correlations +perform_pyramid_correlations=piv_parameters.general_parameters.perform_pyramid_correlations; + +% This is the image filename extension +image_filename_extension=piv_parameters.general_parameters.image_filename_extension; +% This is the image variable name (if the variable is saved within a data +% object that contain multiple variables, ie a 'mat' file) +image_variable_name=piv_parameters.general_parameters.image_variable_name; + +% This is a list of the images within the image read directory +image_filename_list=dir([piv_parameters.general_parameters.image_read_directory,piv_parameters.general_parameters.image_filename_prefix,'*',image_filename_extension]); + +% This is the filename of the first image file (for the purpose of determing the image +% file size) +first_image_filename=[piv_parameters.general_parameters.image_read_directory,image_filename_list(1).name]; + +% This finds the size of the images being loaded which can be 'bmp', 'jpg', +% 'jp2', 'png', 'ppm', or 'tif' images in 2D or 'mat' or 'dcm' images in 3D +if any(strcmpi(image_filename_extension,{'bmp','jpg','jp2','png','ppm','tif'})); + + % This loads in the first image + I_First_Image=imread(first_image_filename); + % This measures the image size + image_size=size(I_First_Image); + % This clears the image from memory + clear('I_First_Image'); + % This adds the third dimension as 1 if it does not exist + if length(image_size)==2; + % This adds the third dimension of the image size as 1 + image_size(3)=1; + end; + % This clears the first image from memory + clear('I_First_Image'); + +elseif strcmpi(image_filename_extension,'mat'); + + % This creates a matlab mat object containing information about the image + first_image_object=matfile(first_image_filename); + % This extracts the image size information + image_size=size(first_image_object,image_variable_name); + % This adds the third dimension as 1 if it does not exist + if length(image_size)==2; + % This adds the third dimension of the image size as 1 + image_size(3)=1; + end; + +elseif strcmpi(image_filename_extension,'dcm'); + + % This reads int the DICOM image + I_First_Image=dicomread(first_image_filename); + % This extracts the image size information + image_size=size(I_First_Image); + % This adds the third dimension as 1 if it does not exist and removes + % the third dimension if the image is 3D since the 3rd dimension will + % correspond to the colormap + if length(image_size)==2; + % This adds the third dimension of the image size as 1 + image_size(3)=1; + elseif length(image_size)==4; + % This removes the 3rd dimension since this corresponds to the + % colormap fod DICOM files which is not used for PIV data + image_size(3)=[]; + end; + % This clears the first image from memory + clear('I_First_Image'); + +elseif strcmpi(image_filename_extension,'cdf'); + + % This reads in the Common Data Format image into a cell array + I_First_Image=cdfread(first_image_filename,'Variables',image_variable_name); + % This extracts the image size information from the cell array + image_size=size(I_First_Image{1}); + % This adds the third dimension as 1 if it does not exist and removes + % the third dimension if the image is 3D since the 3rd dimension will + % correspond to the colormap + if length(image_size)==2; + % This adds the third dimension of the image size as 1 + image_size(3)=1; + end; + % This clears the first image from memory + clear('I_First_Image'); + +elseif strcmpi(image_filename_extension,'h5'); + + % This reads in the HDF5 image + I_First_Image=hdf5read(first_image_filename,image_variable_name); + % This extracts the image size information + image_size=size(I_First_Image); + % This adds the third dimension as 1 if it does not exist and removes + % the third dimension if the image is 3D since the 3rd dimension will + % correspond to the colormap + if length(image_size)==2; + % This adds the third dimension of the image size as 1 + image_size(3)=1; + end; + % This clears the first image from memory + clear('I_First_Image'); + +end; + +% These are the coordinates of the estimated velocities (initially nothing is assummed +% about the velocity field so a zero displacement is used and since the extrapolation value +% is set to 0, only one arbitrary coordinate must be specified) +if image_size(3)==1; + % This calculates the initial interpolation coordinates with only one + % kk coordinate value + [jj_position,ii_position,kk_position]=meshgrid([1,image_size(1)],[1,image_size(2)],[1]); + % This is the estimated (initially zero) velocity for the 2D array + ii_velocity=zeros(2,2,1); + jj_velocity=zeros(2,2,1); + kk_velocity=zeros(2,2,1); + + +elseif image_size(3)>1; + % This calculates the initial interpolation coordinates with the full + % 3D array of coordinates + [jj_position,ii_position,kk_position]=meshgrid([1,image_size(1)],[1,image_size(2)],[1,image_size(3)]); + % This is the estimated (initially zero) velocity for the 3D array + ii_velocity=zeros(1,1,1); % originally 2,2,2 + jj_velocity=zeros(1,1,1); % originally 2,2,2 + kk_velocity=zeros(1,1,1); % originally 2,2,2 + + +end; + +% This iterates through the passes +for pass_index=1:pass_number; + + % This displays that the current pass is being processed + fprintf('\n\nCompleting pass %d of %d passes.\n',pass_index,pass_number); + + % This extracts the effective current window resolution + window_resolution=piv_parameters.pass_parameters(pass_index).window_resolution; + % This extracts the full window size + window_size=piv_parameters.pass_parameters(pass_index).window_size; + % This extracts the window gridding method for the current pass + window_gridding_method=piv_parameters.pass_parameters(pass_index).window_gridding_method; + + % This is the Boolean value stating whether to perform validation on the vector field + validate_vector_field=piv_parameters.pass_parameters(pass_index).validate_vector_field; + + % This is the Boolean value stating whether to perform smoothing of the velocity vector field + smooth_vector_field=piv_parameters.pass_parameters(pass_index).smooth_vector_field; + + % If specified by the user to perform smoothing of the velocity field, this loads the + % smoothing specific parameters + if smooth_vector_field; + + % This loads the Gaussian smoothing standard deviation + gaussian_smoothing_kernal_std=piv_parameters.pass_parameters(pass_index).gaussian_smoothing_kernal_std; + % This loads the Gaussian smoothing kernal size + gaussian_smoothing_kernal_size=piv_parameters.pass_parameters(pass_index).gaussian_smoothing_kernal_size; + + end; + + % If specified by the user to perform validation, this loads the loads the validation + % specific parameters + if validate_vector_field; + + % This loads the vector outlier threshhold limits + minimum_outlier_threshhold=piv_parameters.pass_parameters(pass_index).minimum_outlier_threshhold; + maximum_outlier_threshhold=piv_parameters.pass_parameters(pass_index).maximum_outlier_threshhold; + + % This loads the UOD kernal size + uod_kernal_size=piv_parameters.pass_parameters(pass_index).uod_kernal_size; + % This loads the UOD expected error + uod_epsilon=piv_parameters.pass_parameters(pass_index).uod_epsilon; + % This loads the UOD residual threshhold value + uod_residual_threshhold=piv_parameters.pass_parameters(pass_index).uod_residual_threshhold; + + % This loads the outlier vector replacement method + vector_replacement_method=piv_parameters.pass_parameters(pass_index).vector_replacement_method; + % This loads the relevant vector replacement parameters based upon + % the interpolation method + if strcmp(vector_replacement_method,'local_mean'); + % This loads the minimum number of local vectors used to + % calculate the local mean + local_mean_minimum_valid_vector_number=piv_parameters.pass_parameters(pass_index).local_mean_minimum_valid_vector_number; + elseif strcmp(vector_replacement_method,'laplacian_interpolation'); + % This loads the number of adjacent points to use in + % calculating the interpolated values (this can be either 4 or + % 8 in 2D and 6, 18, or 26 in 3D) + laplacian_interpolation_adjacent_connectivity=piv_parameters.pass_parameters(pass_index).laplacian_interpolation_adjacent_connectivity; + elseif strcmp(vector_replacement_method,'delaunay_interpolation'); + % This loads the Delaunay interpolation weighting method to be + % used in interpolating values + delaunay_interpolation_weighting_method=piv_parameters.pass_parameters(pass_index).delaunay_interpolation_weighting_method; + end; + + end; + + % This calculates the grid spacing for the current pass using either + % the amount of window overlap or the window spacing + if strcmp(window_gridding_method,'window_overlap'); + % This is the ratio of the window overlap distance + window_overlap=piv_parameters.pass_parameters(pass_index).window_overlap; + % This is the spacing between the centers of the windows + grid_spacing=round(window_resolution.*(1-window_overlap)); %originally grid_spacing=round(window_resolution.*(1-window_overlap)); + % This checks if the overlap percentage is so high that grid spacing is zero (although + % this should actually never be done) + if grid_spacing==0; + % This sets the grid spacing to 1 voxel + grid_spacing=ones(size(grid_spacing)); + end; + elseif strcmp(window_gridding_method,'window_spacing'); + % This extracts the window grid spacing from the parameters structure + grid_spacing=piv_parameters.pass_parameters(pass_index).window_spacing; + end; + + + + + % This calculates the window domains + [window_min,window_max,window_center,window_number]=calculate_window_domains(image_size,window_resolution,window_size,grid_spacing); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % This iterates through the frames calculating the velocity fields. % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Z_ensemble = zeros(window_size); + + % This iterates through the image frames calculating the vector fields + for frame_index=initial_image_frame:image_frame_step:final_image_frame; + + % This displays that the current frame is being processed + fprintf('\nProcessing frame %d of %d frames.\n',frame_index,final_image_frame); + + + + % If this is the second or higher pass, this loads the + % previously calculated vector field which is used to deform + % the images in the pyramid correlations + if pass_index>1; + + % This is the filename to load the data from + data_filename_read=[vector_field_write_directory,'pass_',sprintf('%02.0f',pass_index-1),'_frame_',sprintf('%04.0f',frame_index),'_x_',sprintf('%04.0f',frame_index+image_correlation_step),'.mat']; + % This loads in the data file + load(data_filename_read); + + % This renames the data loaded from memory to be conistent + % with the variables names within this code (ie Prana uses X, + % Y, Z, U, V, W for the position and velocity values and this + % code uses the index dimensions ii, jj, kk) + ii_position=Y; + jj_position=X; + kk_position=Z; + ii_velocity=V; + jj_velocity=U; + kk_velocity=W; + + end; + + + + % This checks whether standard correlations are supposed to be + % performed + if not(perform_pyramid_correlations); + + % This performs the standard correlations of the current frame + [Z,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=standard_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity); + + %ensemble correlation addition + Z_ensemble = Z_ensemble + Z; + C=ifftn(Z_ensemble,'symmetric'); + + displacement_vector=subpixel_displacement_vector(C,window_size); + + + + elseif perform_pyramid_correlations; + + % This performs the pyramid correlations of the current frame + [ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=pyramid_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity); + + end; + + + + + + % This performs validation of the vector field if specified by the user + if validate_vector_field; + + % This displays that the vector field is being validated + disp('Validating the vector field . . . '); + + % This determines the outliers of the vector field based upon several metrics + outlier_vector_array=locate_vector_outliers(ii_velocity,jj_velocity,kk_velocity,minimum_outlier_threshhold,maximum_outlier_threshhold,uod_kernal_size,uod_epsilon,uod_residual_threshhold); + + % This replaces the outlier vectors using the specified method + if strcmp(vector_replacement_method,'local_mean'); + % This replaces the outlier vectors using a local mean of valid vectors + [ii_velocity,jj_velocity,kk_velocity]=local_mean_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,local_mean_minimum_valid_vector_number); + elseif strcmp(vector_replacement_method,'laplacian_interpolation'); + % This replaces the outlier vectors using Laplacian interpolation + [ii_velocity,jj_velocity,kk_velocity]=laplacian_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,laplacian_interpolation_adjacent_connectivity); + elseif strcmp(vector_replacement_method,'delaunay_interpolation'); + % This replaces the outlier vectors using a weighted sum based upon the + % Delaunay triangulation of the valid vectors + [ii_velocity,jj_velocity,kk_velocity]=delaunay_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,delaunay_interpolation_weighting_method); + end; + + end; + + % This performs smoothing of the vector field if specified by the users + if smooth_vector_field; + + % This displays that the vector field is being smoothed + disp('Smoothing the vector field . . . '); + + % This smooths the vector fields + [ii_velocity,jj_velocity,kk_velocity]=smooth_velocity_field(ii_velocity,jj_velocity,kk_velocity,gaussian_smoothing_kernal_size,gaussian_smoothing_kernal_std*ones(1,3)); + + end; + + % This is the filename to save the data as + data_filename_write=[vector_field_write_directory,'pass_',sprintf('%02.0f',pass_index),'_frame_',sprintf('%04.0f',frame_index),'_x_',sprintf('%04.0f',frame_index+image_correlation_step),'.mat']; + + % This renames the data from the current vector field to be + % conistent with the variables names used by other vector field + % code (ie Prana uses X, Y, Z, U, V, W for the position and + % velocity values and this code uses the index dimensions ii, jj, + % kk) + Y=ii_position; + X=jj_position; + Z=kk_position; + V=ii_velocity; + U=jj_velocity; + W=kk_velocity; + + % This renames the validation array to be consistent with the + % variables used by other vector field code (ie Prana uses Valid as + % the name for the valid vectors and this code uses + % outlier_vector_array) + if validate_vector_field; + % This creates the binary array variable 'Valid' which gives + % the locations of the vectors that were considered valid + Valid=not(outlier_vector_array); + else; + % This creates an array of NaNs for the variable 'Valid' since + % validation was not performed on the vector field + Valid=NaN(size(X)); + end; + + + % This saves the vector field data to memory + save(data_filename_write,'X','Y','Z','U','V','W','Valid','-v7.3'); + + end; + +end; + + + +function [Z,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=standard_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity); +% This function calculates standard pair correlations of the data set listed +% within the parameters variable 'piv_parameters'. The standard correlation +% consists of taking pair-wise correlations from a time series of images. +% The input parameters of this function are given by +% +% piv_parameters A data structure containing all the necessary +% parameters required to perform the PIV processing +% +% pass_index A scalar integer giving the current pass of the PIV +% processing +% +% frame_index A scalar integer giving the index of the current +% frame being processed from the directory of images. +% For pyramid correlations, the velocity field will +% always be calculated for the frame_index and +% frame_index + 1 froms, ie the changing the +% correlation step does nothing. +% +% window_number A 1 x 3 vector giving the number of PIV windows to +% process in the first, second, and third dimensions +% respectively. If the images are 2D then +% window_number(3) should always equal 1. +% +% window_min A L x 3 vector giving the minimum indices of the +% window domains in each respective dimension where L +% is the total number of windows to process. +% +% window_max A L x 3 vector giving the maximum indices of the +% window domains in each respective dimension where L +% is the total number of windows to process. +% +% window_center A L x 3 vector giving the coordinates of the center +% of the window domains in each respective dimension +% where L is the total number of windows to process. +% +% ii_position A M x N x P array giving the 1st dimension +% coordinates of the previously measured vector field +% if this is at least the second pass. If this is +% the first pass then the coordinates correspond to +% an all-zero vector field. +% +% jj_position A M x N x P array giving the 2nd dimension +% coordinates of the previously measured vector field +% if this is at least the second pass. If this is +% the first pass then the coordinates correspond to +% an all-zero vector field. +% +% kk_position A M x N x P array giving the 3rd dimension +% coordinates of the previously measured vector field +% if this is at least the second pass. If this is +% the first pass then the coordinates correspond to +% an all-zero vector field. +% +% ii_velocity A M x N x P array giving the 1st dimension velocity +% of the previously measured vector field if this is +% at least the second pass. If this is the first +% pass then the velocity will be set to all zeros. +% +% jj_velocity A M x N x P array giving the 2nd dimension velocity +% of the previously measured vector field if this is +% at least the second pass. If this is the first +% pass then the velocity will be set to all zeros. +% +% kk_velocity A M x N x P array giving the 3rd dimension velocity +% of the previously measured vector field if this is +% at least the second pass. If this is the first +% pass then the velocity will be set to all zeros. +% +% The output parameters give the measured (un-validated, un-smoothed) +% velocity field of the current frame pair and are described by +% +% ii_position A Q x R x S array giving the 1st dimension +% coordinates of the measured vector field, ie the +% center of the correlation window in the 1st +% dimension. +% +% jj_position A Q x R x S array giving the 2nd dimension +% coordinates of the measured vector field, ie the +% center of the correlation window in the 1st +% dimension. +% +% kk_position A Q x R x S array giving the 3rd dimension +% coordinates of the measured vector field, ie the +% center of the correlation window in the 1st +% dimension. +% +% ii_velocity A Q x R x S array giving the 1st dimension velocity +% of the current vector field measured at each +% correlation window. +% +% jj_velocity A Q x R x S array giving the 1st dimension velocity +% of the current vector field measured at each +% correlation window. +% +% kk_velocity A Q x R x S array giving the 1st dimension velocity +% of the current vector field measured at each +% correlation window. +% +% Authors: Rod La Foy +% First Created On: 14 March 2014 +% Last Modified On: 14 March 2014 + +% This is the image filename extension +image_filename_extension=piv_parameters.general_parameters.image_filename_extension; +% This is the image variable name (if the variable is saved within a data +% object that contain multiple variables, ie a 'mat' file) +image_variable_name=piv_parameters.general_parameters.image_variable_name; + +% This is a list of the images within the image read directory +image_filename_list=dir([piv_parameters.general_parameters.image_read_directory,piv_parameters.general_parameters.image_filename_prefix,'*',image_filename_extension]); + +% This extracts the effective current window resolution +window_resolution=piv_parameters.pass_parameters(pass_index).window_resolution; +% This extracts the full window size +window_size=piv_parameters.pass_parameters(pass_index).window_size; + +% This is the number of frames to step between correlation frame pairs +image_correlation_step=piv_parameters.general_parameters.image_correlation_step; + +% This is the bulk window offset distance +bulk_window_offset=piv_parameters.pass_parameters(pass_index).bulk_window_offset; + +% This is the correlation method to use for the current pass +correlation_method=piv_parameters.pass_parameters(pass_index).correlation_method; +% If the correlation method is the RPC method, this creates the RPC filter +if strcmp(correlation_method,'RPC'); + % This extracts the size of the particles in the images + particle_diameter=piv_parameters.general_parameters.particle_diameter; + % This creates the spectral energy filter + spectral_filter=single(fftshift(create_spectral_filter(window_size,particle_diameter))); +elseif strcmp(correlation_method,'SCC'); + % This sets the spectral filter to a null value + spectral_filter=[]; +end; + +% This is a Boolean value stating whether to zero-mean the correlation +% windows +zero_mean_windows=piv_parameters.pass_parameters(pass_index).zero_mean_windows; + +% This sets the fft library for optimal speed calculation +fftw('planner','measure'); + +% This determines the size of the Gaussian function to use for the window masking +gaussian_width=determine_gaussian_size_2(window_size,window_resolution); + +% This creates a Gaussian windowing function +gaussian_filter=single(create_gaussian_filter(window_size,gaussian_width)); + +% These are the coordinate vectors +ii_position_current=reshape(window_center(:,1),window_number); +jj_position_current=reshape(window_center(:,2),window_number); +kk_position_current=reshape(window_center(:,3),window_number); + + +% This initializes the velocity vectors +ii_velocity_current=zeros(size(ii_position_current)); +jj_velocity_current=zeros(size(ii_position_current)); +kk_velocity_current=zeros(size(ii_position_current)); + +% This determines the domain over which the known velocity field should be +% extrapolated +% ii_extrap_min=min(ii_position_current(:))-ceil(window_size(1)/2); +% ii_extrap_max=max(ii_position_current(:))+ceil(window_size(1)/2); +% jj_extrap_min=min(jj_position_current(:))-ceil(window_size(2)/2); +% jj_extrap_max=max(jj_position_current(:))+ceil(window_size(2)/2); +% kk_extrap_min=min(kk_position_current(:))-ceil(window_size(3)/2); +% kk_extrap_max=max(kk_position_current(:))+ceil(window_size(3)/2); +% % This adds extrapolated points to the previously measured vector field so +% % that if the vector field is interpolated at points near the edge of the +% % image, the interpolated points will have approximately correct values (as +% % opposed to zero or NaNs) +% [~,~,~,ii_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,ii_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max); +% [~,~,~,jj_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,jj_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max); +% [ii_position,jj_position,kk_position,kk_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,kk_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max); + +% If the bulk window offset is non-zero, this adds it to the previous pass +% velocity field estimate +if any(bulk_window_offset); + % This adds the bulk window offset to the previously measured velocity + % field + ii_velocity=ii_velocity+bulk_window_offset(1); + jj_velocity=jj_velocity+bulk_window_offset(2); + kk_velocity=kk_velocity+bulk_window_offset(3); + % This creates a Boolean value stating to apply the bulk window offset + perform_bulk_window_offset=true; +else; + % This creates a Boolean value stating to apply the bulk window offset + perform_bulk_window_offset=false; +end; + +% If this is the second pass or higher, this calculates the discrete window +% offset, otherwise the offset is just set to zero +if (pass_index==1)&&(not(perform_bulk_window_offset)); + % This sets the first frame minimum and maximum values to the + % non-offset values + frame_1_window_min=window_min; + frame_1_window_max=window_max; + % This sets the second frame minimum and maximum values to the + % non-offset values + frame_2_window_min=window_min; + frame_2_window_max=window_max; + % This sets the displacement offset array to all zeros + displacement_offset=zeros(size(window_min)); +else; + % This function adds the known velocity field to the window domains, ie + % accounts for the discrete window offset + [frame_1_window_min,frame_1_window_max,frame_2_window_min,frame_2_window_max,displacement_offset]=dwo_window_domains(window_min,window_max,window_center,image_correlation_step,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity); +end; + +% This is the filename of the first image to load +image_1_filename=[piv_parameters.general_parameters.image_read_directory,image_filename_list(frame_index).name]; +% This is the filename of the first image to load + + +image_2_filename=[piv_parameters.general_parameters.image_read_directory,image_filename_list(frame_index+image_correlation_step).name]; + +% This finds the size of the images being loaded which can be 'bmp', 'jpg', +% 'jp2', 'png', 'ppm', or 'tif' images in 2D or 'mat images in 3D +if any(strcmpi(image_filename_extension,{'bmp','jpg','jp2','png','ppm','tif'})); + + % This loads the first image + I1=double(imread(image_1_filename)); + % This loads the second image + I2=double(imread(image_2_filename)); + +elseif strcmpi(image_filename_extension,'mat'); + + % This creates a matlab file object of the first image of the correlation pairs + image_1_object=matfile(image_1_filename); + % This loads the first image + eval(['I1=double(image_1_object.',image_variable_name,');']); + % This creates a matlab file object of the second image of the correlation pairs + image_2_object=matfile(image_2_filename); + % This loads the second image + eval(['I2=double(image_2_object.',image_variable_name,');']); + +elseif strcmpi(image_filename_extension,'dcm'); + + % This loads the first image + I1=double(squeeze(dicomread(image_1_filename))); + % This loads the second image + I2=double(squeeze(dicomread(image_2_filename))); + +elseif strcmpi(image_filename_extension,'cdf'); + + % This loads the first image cell array + I1_Temp=double(cdfread(image_1_filename,'Variables',image_variable_name)); + % This loads the image from the cell array + I1=I1_Temp{1}; + % This clears the cell array from memory + clear('I1_Temp'); + % This loads the second image + I2_Temp=double(cdfread(image_2_filename,'Variables',image_variable_name)); + % This loads the image from the cell array + I2=I2_Temp{1}; + % This clears the cell array from memory + clear('I2_Temp'); + +elseif strcmpi(image_filename_extension,'h5'); + + % This loads the first image + I1=double(hdf5read(image_1_filename,image_variable_name)); + % This loads the second image + I2=double(hdf5read(image_2_filename,image_variable_name)); + +end; + +% This determines the image size +image_size=zeros(1,3); +image_size(1)=size(I1,1); +image_size(2)=size(I1,2); +image_size(3)=size(I1,3); + +% This iterates through the windows performing the cross-correlations + + + + + + + +for window_index=1:prod(window_number); + + + % This displays the calculation progress + display_calculation_progress(window_index,1:size(window_min,1)); + + % This extracts the correlation window from the first frame + I1_Window=extract_correlation_window(I1,window_size,frame_1_window_min(window_index,:),frame_1_window_max(window_index,:),image_size); + + % This extracts the correlation window from the second frame + I2_Window=extract_correlation_window(I2,window_size,frame_2_window_min(window_index,:),frame_2_window_max(window_index,:),image_size); + + % This subtracts the mean value from the windows if specified by the + % user + if zero_mean_windows; + % This subtracts the mean value from the first image window + I1_Window=I1_Window-mean(I1_Window(:)); + % This subtracts the mean value from the second image window + I2_Window=I2_Window-mean(I2_Window(:)); + end; + + + + % This performs the cross-correlation measurement + [Z,displacement_vector]=calculate_window_correlation(I1_Window,I2_Window,gaussian_filter,spectral_filter,window_size,correlation_method); + + % These are the velocities of the current correlation volume + ii_velocity_current(window_index)=displacement_vector(1)+displacement_offset(window_index,1); + jj_velocity_current(window_index)=displacement_vector(2)+displacement_offset(window_index,2); + kk_velocity_current(window_index)=displacement_vector(3)+displacement_offset(window_index,3); + +end; + +% This overwrites the input velocity field coordinates with the current +% pass's coordinates +ii_position=ii_position_current; +jj_position=jj_position_current; +kk_position=kk_position_current; +% This overwrites the input velocity field vectors with the current +% pass's vectors +ii_velocity=ii_velocity_current/image_correlation_step; + +size(ii_velocity_current) + +jj_velocity=jj_velocity_current/image_correlation_step; +kk_velocity=kk_velocity_current/image_correlation_step; + +%pyramid correlation location - originally + +function [Z,displacement_vector]=calculate_window_correlation(I1,I2,gaussian_window_filter,spectral_filter,window_size,method); +% This performs the volumetric correlation between windows I1 and I2. The +% gaussian filter is applied to the windows, the correlation is performed, +% and a sub-pixel fit is made to the correlation volume. + +% This calculates the correlation volumes using either the SCC or RPC +% correlations +if strcmp(method,'SCC'); + + + % This performs the SCC correlation on the two windows + [Z,C]=scc_correlation(I1,I2,gaussian_window_filter); + + +elseif strcmp(method,'RPC'); + + % This performs the RPC correlation on the two windows + C=rpc_correlation(I1,I2,gaussian_window_filter,spectral_filter); + +end; + + + +% This calculates the displacement vector from the correlation volumes +displacement_vector=subpixel_displacement_vector(C,window_size); + + + + +function [Z,C]=scc_correlation(I1,I2,gaussian_window_filter); +% This function performs the standard cross-correlation on the two windows 'I1' and 'I2' +% after masking the windows with a Gaussian function given by 'gaussian_window_filter'. + +% This applies the Gaussian filter to the first window + +I1=I1.*gaussian_window_filter; + + +% This applies the Gaussian filter to the second window +I2=I2.*gaussian_window_filter; + +% This calculates the FFT of the first window +Y1=fftn(I1); + +% This calculates the FFT of the second window +Y2=fftn(I2); + +% This performs the cross-correlation in the spectral domain +Z= Y2.*conj(Y1); % originally Y2.*conj(Y1) + + +% This performs the inverse FFT to return the correlation volume +C=ifftn(Z,'symmetric'); + + + + +function C=rpc_correlation(I1,I2,gaussian_window_filter,spectral_filter); +% This function performs the robust phase correlation on the two windows 'I1' and 'I2' +% after masking the windows with a Gaussian function given by 'gaussian_window_filter' and +% using the RPC spectral filter defined by 'spectral_filter'. + +% This applies the Gaussian filter to the first window +I1=I1.*gaussian_window_filter; + +% This applies the Gaussian filter to the second window +I2=I2.*gaussian_window_filter; + +% This calculates the FFT of the first window +Y1=fftn(I1); + +% This calculates the FFT of the second window +Y2=fftn(I2); + +% This performs the cross-correlation in the spectral domain +YC=Y2.*conj(Y1); + +% This is the magnitude of the cross-correlation in the spectral domain +YC_Magnitude=sqrt(YC.*conj(YC)); + +% These are the indices of the non-zero magnitude components of the +% cross-correlation +non_zero_indices=(YC_Magnitude(:)~=0); + +% This initializes the phase correlation array +R=YC; + +% THis is the phase correlation in the spectral domain +R(non_zero_indices)=YC(non_zero_indices)./YC_Magnitude(non_zero_indices); + +% This applies the spectral filter +R=R.*spectral_filter; + +% This performs the inverse FFT to return the correlation volume +C=ifftn(R,'symmetric'); + +% This takes the absolute value of the correlation volume to ensure positive values +C=abs(C); + + + +function displacement_vector=subpixel_displacement_vector(C,window_size); +% This function finds the subpixel fit to the correlation volume given by C. + +% This is the linear index of the correlation maximum (this only takes the first maximum +% value if multiple are found - this is unlikely, but should be addressed) +max_linear_index=find(C==max(C(:)),1); + +% This converts the linear index to a set of subscripted indices +[ii_max,jj_max,kk_max]=ind2sub(window_size,max_linear_index); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% This performs the three-point Gaussian fitting in the ii-direction. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% These are the three points to fit in the ii-direction +if (ii_max>1)&&(ii_max1)&&(jj_max1)&&(kk_max(window_size(1)/2); + ii_fit=ii_fit-window_size(1); +end; + +% This accounts for the periodicity of the correlation volume for the jj displacement +if jj_fit>(window_size(2)/2); + jj_fit=jj_fit-window_size(2); +end; + +% This accounts for the periodicity of the correlation volume for the kk displacement +if kk_fit>(window_size(3)/2); + kk_fit=kk_fit-window_size(3); +end; + +% This is the displacement vector of the current window +displacement_vector=[ii_fit,jj_fit,kk_fit]; + + + +function I_Window=extract_correlation_window(I,window_size,window_min,window_max,image_size); +% This function extracts the current windows from the frame 'I' where +% 'window_size' is the size of the correlation window, 'window_min' and +% 'window_max' are the 1 x 3 vectors of the window minimum and maximum +% indices, and 'image_size' is the the size of the frame 'I'. + +% These are the coordinate domains to extract +ii_min=window_min(1); +ii_max=window_max(1); +jj_min=window_min(2); +jj_max=window_max(2); +kk_min=window_min(3); +kk_max=window_max(3); + +% This checks if the ii indices are outside the domain of the image and sets them to the +% domain of the image if true +if ii_min<1; + ii_min=1; +end; +if ii_max>image_size(1); + ii_max=image_size(1); +end; +% This checks if the jj indices are outside the domain of the image and sets them to the +% domain of the image if true +if jj_min<1; + jj_min=1; +end; +if jj_max>image_size(2); + jj_max=image_size(2); +end; +% This checks if the kk indices are outside the domain of the image and sets them to the +% domain of the image if true +if kk_min<1; + kk_min=1; +end; +if kk_max>image_size(3); + kk_max=image_size(3); +end; + +% This extracts the current window from the frame +I_Window=I(ii_min:ii_max,jj_min:jj_max,kk_min:kk_max); + +% This initilizes the image size vector +I_Window_Size=zeros(1,3); +% This is the size of the extracted image +I_Window_Size(1)=size(I_Window,1); +I_Window_Size(2)=size(I_Window,2); +I_Window_Size(3)=size(I_Window,3); + +% If the extracted image is smaller than the full window size (due to the window being +% on the edge of the full image) this pads the image with zeros +if any(I_Window_Size~=window_size); + % This initializes an array of zeros the size of the full window + I_Window_Temp=zeros(window_size); + % These are the indices to insert the image into the zero padded image (the + % indices are centered in the window so that information is not lost during the + % masking process) + window_insert_min=floor((window_size-I_Window_Size)/2)+1; + window_insert_max=window_insert_min+I_Window_Size-1; + % This copies the current image from the frame into the window + I_Window_Temp(window_insert_min(1):window_insert_max(1),window_insert_min(2):window_insert_max(2),window_insert_min(3):window_insert_max(3))=I_Window; + % This renames the variable so that the zero-padded image is used + I_Window=I_Window_Temp; +end; + + + + +function p=determine_gaussian_size_2(window_size,resolution_size); +% This function is used to determine the diameter of a Gaussian function +% used to mask PIV windows to prevent aliasing during the Fourier +% transform. +% +% The input argument 'window_size' refers to the full size of +% the PIV window prior to performing the masking operation and may be a +% scaler or a vector. The input argument 'resolution_size' refers to the +% effective resolution that the user wants to attain from the masked PIV +% window, ie the window may be 64 x 64 in size (window_size=[64,64]) but +% the user wants the effective resolution after masking to be equal to 32 x +% 32 (resolution_size=[32,32]). +% +% The output argument 'p' is used to contruct the Gaussian masking function +% of the following form +% +% G(x) = exp(-((p*x)^2)/2) (1) +% +% where the independent variable lies in the domain +% +% -(window_size-1)/2 <= x <= (window_size-1)/2 (2) +% +% The variable 'p' is calculated by solving the integral equation +% +% resolution_size = int(G(x),-(window_size-1)/2,(window_size-1)/2) (3) +% +% or the equivalent expression +% +% resolution_size = sqrt(2*pi)*erf(window_size*p/(2*sqrt(2))/p (4) +% +% given by simplifying the integral. These equations are taken from +% "Assessment of advanced windowing techniques for digital particle image +% velocimetry (DPIV)" Eckstein 2009. +% +% This function is written to replace the 'findwidth' function in prana. +% The equation is solved by re-writing equation (4) in the form +% +% 0 = f(x) = erf(b*x)/x-a (5) +% +% and solving for the independent variable using Halley's method. +% +% Written By: Rod La Foy +% Written On: 6 June 2013 + +% These are the ratios of the resolution size to window size +size_ratio=resolution_size./window_size; + +% This initializes the output vector p +p=zeros(size(window_size)); + +% This pre-computes the square root of pi for later use +pi_sqrt=sqrt(pi); + +% This is the tolerance with which to find the root +tolerance=1e-4; + +% This iterates through the number of dimensions to solve for (probably +% either 2 or 3) +for dimension_number=1:length(window_size); + + % This checks if the size ratio is equal to one or greater + if size_ratio(dimension_number)>=1; + + % This sets the Gaussian width variable to 0 + p(dimension_number)=0; + + else; + + % This is the first constant term in the equation + a=resolution_size(dimension_number)/sqrt(2*pi); + % This is the second constant term in the equation + b=window_size(dimension_number)/(2*sqrt(2)); + + % This is the initial guess of value of the Gaussian width term + p_temp=1; + + % This iterates to the root + while true; + + % This is the error function of the b*p term + bp_erf=erf(b*p_temp); + % This is the exponential term of the (b*p)^2 term + bp_exp=exp((b*p_temp)^2); + + % This is the numerator of the difference term + dp_numerator=a*p_temp-bp_erf; + + % This initializes the denominator term + dp_denominator=0; + % This adds the first term to the denominator + dp_denominator=dp_denominator+a; + % This adds the second term to the denominator + dp_denominator=dp_denominator+2*b/(pi_sqrt*bp_exp); + % This adds the third term to the denominator + dp_denominator=dp_denominator+(2*b^3*p_temp^2*dp_numerator)/(2*b*p_temp-bp_exp*pi_sqrt*bp_erf); + + % This is the expected difference to the root + dp=dp_numerator/dp_denominator; + + % This is the new root approximation + p_temp=p_temp-dp; + + % This breaks if the error is less then the tolerance + if abs(erf(b*p_temp)/p_temp-a)2)&&(size(ii_coordinates,2)>2); + + % This is the estimated ii-displacement at the center of the windows + estimated_window_dii=interp2(ii_coordinates,jj_coordinates,ii_displacement,window_center(:,1),window_center(:,2),'cubic',0); + % This is the estimated jj-displacement at the center of the windows + estimated_window_djj=interp2(ii_coordinates,jj_coordinates,jj_displacement,window_center(:,1),window_center(:,2),'cubic',0); + % This is the estimated kk-displacement at the center of the windows + estimated_window_dkk=interp2(ii_coordinates,jj_coordinates,kk_displacement,window_center(:,1),window_center(:,2),'cubic',0); + + else; + + % This is the estimated ii-displacement at the center of the windows + estimated_window_dii=interp2(ii_coordinates,jj_coordinates,ii_displacement,window_center(:,1),window_center(:,2),'linear',0); + % This is the estimated jj-displacement at the center of the windows + estimated_window_djj=interp2(ii_coordinates,jj_coordinates,jj_displacement,window_center(:,1),window_center(:,2),'linear',0); + % This is the estimated kk-displacement at the center of the windows + estimated_window_dkk=interp2(ii_coordinates,jj_coordinates,kk_displacement,window_center(:,1),window_center(:,2),'linear',0); + + end; + +else; + + % If the vector field is greater then 3 x 3, then a cubic interpolation + % is performed, otherwise a linear interpolation is performed + if (size(ii_coordinates,1)>2)&&(size(ii_coordinates,2)>2)&&(size(ii_coordinates,3)>2); + + % This is the estimated ii-displacement at the center of the windows + estimated_window_dii=interp3(ii_coordinates,jj_coordinates,kk_coordinates,ii_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0); + % This is the estimated jj-displacement at the center of the windows + estimated_window_djj=interp3(ii_coordinates,jj_coordinates,kk_coordinates,jj_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0); + % This is the estimated kk-displacement at the center of the windows + estimated_window_dkk=interp3(ii_coordinates,jj_coordinates,kk_coordinates,kk_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0); + + else; + + % This is the estimated ii-displacement at the center of the windows + estimated_window_dii=interp3(ii_coordinates,jj_coordinates,kk_coordinates,ii_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0); + % This is the estimated jj-displacement at the center of the windows + estimated_window_djj=interp3(ii_coordinates,jj_coordinates,kk_coordinates,jj_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0); + % This is the estimated kk-displacement at the center of the windows + estimated_window_dkk=interp3(ii_coordinates,jj_coordinates,kk_coordinates,kk_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0); + + end; + +end; + +% This is half of the measured velocity at the window coordinates +rounded_half_estimated_displacement=round([estimated_window_dii,estimated_window_djj,estimated_window_dkk]/2); + +% This calculates the first frame minimum domain values +frame_1_window_min=window_min-floor(rounded_half_estimated_displacement); +% This calculates the first frame maximum domain values +frame_1_window_max=window_max-floor(rounded_half_estimated_displacement); + +% This calculates the second frame minimum domain values +frame_2_window_min=window_min+ceil(rounded_half_estimated_displacement); +% This calculates the second frame maximum domain values +frame_2_window_max=window_max+ceil(rounded_half_estimated_displacement); + +% This is the offset that must be added back onto the measured velocity +% field of the current pass to account for the dicrete window offset +displacement_offset=ceil(rounded_half_estimated_displacement)+floor(rounded_half_estimated_displacement); + + + +function [window_min,window_max,window_center,window_number]=calculate_window_domains(image_size,window_resolution,window_size,grid_spacing); +% This function calculates the minimum and maximum values of the window +% indices to extract from the image. + +% These are the number of windows in each dimension +window_number=ceil((image_size-window_resolution)./grid_spacing); %ORIGINALLY window_number=ceil((image_size-window_resolution)./grid_spacing)+1 + +% This initializes the window indices +window_min=zeros(prod(window_number),3); +window_max=zeros(prod(window_number),3); + +% This initializes the window centers +window_center=zeros(prod(window_number),3); + +% This is the full resolution covered by the PIV windows (which will be equal to or larger +% than the resolution of the full image) +full_resolution=grid_spacing.*(window_number-1)+window_resolution; + + + + + +% This is the distance in indices below 1 to start the windows +negative_index_distance=floor((full_resolution-image_size)/2); + +% This is a counting variable for indexing the window domains +count=0; + +% This calculates the window domains +for kk=1:window_number(3); + for jj=1:window_number(2); + for ii=1:window_number(1); + + % This increments the counting variable + count=count+1; + + % These are the minimum window indices + window_min(count,:)=([ii,jj,kk]-1).*grid_spacing-ceil((window_size-window_resolution)/2)-negative_index_distance+1; + % These are the maximum window indices + window_max(count,:)=window_min(count,:)+window_size-1; + + % These are the window center indices + window_center(count,:)=(window_min(count,:)+window_max(count,:))/2 - 0.5; + + window_min(count,:) + window_max(count,:) + + end; + end; +end; + + + +function spectral_filter=create_spectral_filter(window_size,particle_diameter); +% This function creates the RPC spectral energy filter with a size of 'window_size' using +% a particle diameter given by 'particle_diameter' where 'bit_depth' is the bit depth of +% the three-dimensional image to be processed. + +% These are the wavenumber vectors +k_ii_vector=-pi:2*pi/window_size(1):pi-2*pi/window_size(1); +k_jj_vector=-pi:2*pi/window_size(2):pi-2*pi/window_size(2); +k_kk_vector=-pi:2*pi/window_size(3):pi-2*pi/window_size(3); + +% These are the wavenumber coordinates +[k_ii,k_jj,k_kk]=ndgrid(k_ii_vector,k_jj_vector,k_kk_vector); + +% This is the squared radius of the wavenumber +rho=k_ii.^2+k_jj.^2+k_kk.^2; + +% None of the coefficients actually matter . . . so RPC filter is just a low pass filter . . . +spectral_filter=exp(-(rho*particle_diameter^2)/16); + + + +function display_calculation_progress(current_value,value_vector); +% This function displays a progress bar showing the percent complete of the +% currently running calculation where the calculation is being iterated +% over the vector 'value_vector' and the current iteration's value is equal +% to 'current_value'. For example in the for loop +% +% for ii=1:10; +% commands . . . +% end; +% +% the arguments would be equal to +% +% current_value=ii; +% value_vector=1:10; +% +% This convention holds for non-integer or non-monotonic vectors, although +% to work correctly all values of value_vector must be unique. + +% This is the number of characters to display in the pogress bar (denoted +% by either '#' or '-' characters) +progress_bar_character_number=50; + +% This is the index of the value_vector that corresponds to the +% current_value +[null_variable,value_index]=min(abs(current_value-value_vector)); + +% This is the percentage of the calculation that is completed +current_progress_decimal=value_index/length(value_vector); +% This creates the character string showing the numerical and graphical +% progress for the current iteration +current_text_string=generate_progress_string(current_progress_decimal,progress_bar_character_number); + +% If this is the first iteration, then a new line is added, the progress +% bar is displayed, and the function is exited +if value_index==1; + % This displays the portion of the string showing the percentage of + % the calculation that is complete that is new to this iteration + fprintf(current_text_string); + % This ends the function + return; +end; + +% This is the percentage of the calculation that was completed during the +% last iteration +previous_progress_decimal=(value_index-1)/length(value_vector); +% This creates the character string showing the numerical and graphical +% progress for the previous iteration +previous_text_string=generate_progress_string(previous_progress_decimal,progress_bar_character_number); + +% This compares the current progress string with the previous string and if +% they are the same, then the function exits. If they are different, then +% only the text after the difference is displayed. +if strcmp(current_text_string,previous_text_string); + + % If this is the last time that the progress bar will be displayed, this + % prints a new line + if value_index==length(value_vector); + % This prints a new line after the progress bar + fprintf('\n'); + end; + + % This exits the function without changing the progress bar + return; + +else; + % This is the total number of charcters to be displayed + string_character_number=length(current_text_string); + + % This is the index into the string where the strings first differ + first_difference_index=find(not(current_text_string==previous_text_string),1,'first'); + + % These are the locations of the double percent signs '%%' + double_percent_indices=strfind(current_text_string,'%%'); + % This removes the double percent indices that are before the first + % difference index + double_percent_indices(double_percent_indices1; + % This calculates the initial interpolation coordinates with the full + % 3D array of coordinates + [jj_position,ii_position,kk_position]=meshgrid([1,image_size(1)],[1,image_size(2)],[1,image_size(3)]); + % This is the estimated (initially zero) velocity for the 3D array + ii_velocity=zeros(1,1,1); % originally 2,2,2 + jj_velocity=zeros(1,1,1); % originally 2,2,2 + kk_velocity=zeros(1,1,1); % originally 2,2,2 + + +end; + +% This iterates through the passes +for pass_index=1:pass_number; + + % This displays that the current pass is being processed + fprintf('\n\nCompleting pass %d of %d passes.\n',pass_index,pass_number); + + % This extracts the effective current window resolution + window_resolution=piv_parameters.pass_parameters(pass_index).window_resolution; + % This extracts the full window size + window_size=piv_parameters.pass_parameters(pass_index).window_size; + % This extracts the window gridding method for the current pass + window_gridding_method=piv_parameters.pass_parameters(pass_index).window_gridding_method; + + % This is the Boolean value stating whether to perform validation on the vector field + validate_vector_field=piv_parameters.pass_parameters(pass_index).validate_vector_field; + + % This is the Boolean value stating whether to perform smoothing of the velocity vector field + smooth_vector_field=piv_parameters.pass_parameters(pass_index).smooth_vector_field; + + % If specified by the user to perform smoothing of the velocity field, this loads the + % smoothing specific parameters + if smooth_vector_field; + + % This loads the Gaussian smoothing standard deviation + gaussian_smoothing_kernal_std=piv_parameters.pass_parameters(pass_index).gaussian_smoothing_kernal_std; + % This loads the Gaussian smoothing kernal size + gaussian_smoothing_kernal_size=piv_parameters.pass_parameters(pass_index).gaussian_smoothing_kernal_size; + + end; + + % If specified by the user to perform validation, this loads the loads the validation + % specific parameters + if validate_vector_field; + + % This loads the vector outlier threshhold limits + minimum_outlier_threshhold=piv_parameters.pass_parameters(pass_index).minimum_outlier_threshhold; + maximum_outlier_threshhold=piv_parameters.pass_parameters(pass_index).maximum_outlier_threshhold; + + % This loads the UOD kernal size + uod_kernal_size=piv_parameters.pass_parameters(pass_index).uod_kernal_size; + % This loads the UOD expected error + uod_epsilon=piv_parameters.pass_parameters(pass_index).uod_epsilon; + % This loads the UOD residual threshhold value + uod_residual_threshhold=piv_parameters.pass_parameters(pass_index).uod_residual_threshhold; + + % This loads the outlier vector replacement method + vector_replacement_method=piv_parameters.pass_parameters(pass_index).vector_replacement_method; + % This loads the relevant vector replacement parameters based upon + % the interpolation method + if strcmp(vector_replacement_method,'local_mean'); + % This loads the minimum number of local vectors used to + % calculate the local mean + local_mean_minimum_valid_vector_number=piv_parameters.pass_parameters(pass_index).local_mean_minimum_valid_vector_number; + elseif strcmp(vector_replacement_method,'laplacian_interpolation'); + % This loads the number of adjacent points to use in + % calculating the interpolated values (this can be either 4 or + % 8 in 2D and 6, 18, or 26 in 3D) + laplacian_interpolation_adjacent_connectivity=piv_parameters.pass_parameters(pass_index).laplacian_interpolation_adjacent_connectivity; + elseif strcmp(vector_replacement_method,'delaunay_interpolation'); + % This loads the Delaunay interpolation weighting method to be + % used in interpolating values + delaunay_interpolation_weighting_method=piv_parameters.pass_parameters(pass_index).delaunay_interpolation_weighting_method; + end; + + end; + + % This calculates the grid spacing for the current pass using either + % the amount of window overlap or the window spacing + if strcmp(window_gridding_method,'window_overlap'); + % This is the ratio of the window overlap distance + window_overlap=piv_parameters.pass_parameters(pass_index).window_overlap; + % This is the spacing between the centers of the windows + grid_spacing=round(window_resolution.*(1-window_overlap)); %originally grid_spacing=round(window_resolution.*(1-window_overlap)); + % This checks if the overlap percentage is so high that grid spacing is zero (although + % this should actually never be done) + if grid_spacing==0; + % This sets the grid spacing to 1 voxel + grid_spacing=ones(size(grid_spacing)); + end; + elseif strcmp(window_gridding_method,'window_spacing'); + % This extracts the window grid spacing from the parameters structure + grid_spacing=piv_parameters.pass_parameters(pass_index).window_spacing; + end; + + + + + % This calculates the window domains + [window_min,window_max,window_center,window_number]=calculate_window_domains(image_size,window_resolution,window_size,grid_spacing); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % This iterates through the frames calculating the velocity fields. % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Z_ensemble = zeros(window_size); + tick = 1; + + % This iterates through the image frames calculating the vector fields + for frame_index=initial_image_frame:image_frame_step:final_image_frame; + + %comment + + % This displays that the current frame is being processed + fprintf('\nProcessing frame %d of %d frames.\n',frame_index,final_image_frame); + + + + % If this is the second or higher pass, this loads the + % previously calculated vector field which is used to deform + % the images in the pyramid correlations + if pass_index>1; + + % This is the filename to load the data from + data_filename_read=[vector_field_write_directory,'pass_',sprintf('%02.0f',pass_index-1),'_frame_',sprintf('%04.0f',frame_index),'_x_',sprintf('%04.0f',frame_index+image_correlation_step),'.mat']; + % This loads in the data file + load(data_filename_read); + + % This renames the data loaded from memory to be conistent + % with the variables names within this code (ie Prana uses X, + % Y, Z, U, V, W for the position and velocity values and this + % code uses the index dimensions ii, jj, kk) + ii_position=Y; + jj_position=X; + kk_position=Z;d + ii_velocity=V; + jj_velocity=U; + kk_velocity=W; + + end; + + + + % This checks whether standard correlations are supposed to be + % performed + if not(perform_pyramid_correlations); + + % This performs the standard correlations of the current frame + [Z,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=standard_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity); + + %ensemble correlation addition + Z_ensemble = Z_ensemble + Z; + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% if frame_index == initial_image_frame +% +% C=ifftn(Z_ensemble,'symmetric'); %for SCC +% +% C = abs(C); %for SCC +% end + + C=ifftn(Z_ensemble,'symmetric'); %for SCC + + C = abs(C); %for SCC + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + if frame_index == final_image_frame - 1; +% + Z_ensemble_Magnitude = sqrt(Z_ensemble.*conj(Z_ensemble)); + + non_zero_indices=(Z_ensemble_Magnitude(:)~=0); + + R=Z_ensemble; + + R(non_zero_indices)=Z_ensemble(non_zero_indices)./Z_ensemble_Magnitude(non_zero_indices); + + particle_diameter=piv_parameters.general_parameters.particle_diameter; + + spectral_filter=single(fftshift(create_spectral_filter(window_size,particle_diameter))); + + phase_angle = fftshift(angle(R)); +% % +% % %%%%below specifically APC%%%%%%%%%%%%% +% % + kernel_radius = 2; + + phase_quality = calculate_phase_quality_3D(phase_angle, kernel_radius); + + phase_mask = calculate_phase_mask_ellipse_3D... + (phase_quality, kernel_radius); + % + % + % + filtered_spectral_phase_plane = fftshift(phase_mask).*R; + + C = ifftn(filtered_spectral_phase_plane, 'symmetric'); + + C = abs(C); + end +% + + %%%%%%below is RPC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% R_final = R.*spectral_filter; +% +% +% C = ifftn(R_final,'symmetric'); %for RPC +% +% C = abs(C); %FOR RPC + + + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + + + displacement_vector=subpixel_displacement_vector(C,window_size); + + ii_velocity_current(1)=displacement_vector(1); + jj_velocity_current(1)=displacement_vector(2); + kk_velocity_current(1)=displacement_vector(3); + + ii_velocity=ii_velocity_current/image_correlation_step; + jj_velocity=jj_velocity_current/image_correlation_step; + kk_velocity=kk_velocity_current/image_correlation_step; + + + + elseif perform_pyramid_correlations; + + % This performs the pyramid correlations of the current frame + [ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=pyramid_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity); + + end; + + + + + + % This performs validation of the vector field if specified by the user + if validate_vector_field; + + % This displays that the vector field is being validated + disp('Validating the vector field . . . '); + + % This determines the outliers of the vector field based upon several metrics + outlier_vector_array=locate_vector_outliers(ii_velocity,jj_velocity,kk_velocity,minimum_outlier_threshhold,maximum_outlier_threshhold,uod_kernal_size,uod_epsilon,uod_residual_threshhold); + + % This replaces the outlier vectors using the specified method + if strcmp(vector_replacement_method,'local_mean'); + % This replaces the outlier vectors using a local mean of valid vectors + [ii_velocity,jj_velocity,kk_velocity]=local_mean_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,local_mean_minimum_valid_vector_number); + elseif strcmp(vector_replacement_method,'laplacian_interpolation'); + % This replaces the outlier vectors using Laplacian interpolation + [ii_velocity,jj_velocity,kk_velocity]=laplacian_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,laplacian_interpolation_adjacent_connectivity); + elseif strcmp(vector_replacement_method,'delaunay_interpolation'); + % This replaces the outlier vectors using a weighted sum based upon the + % Delaunay triangulation of the valid vectors + [ii_velocity,jj_velocity,kk_velocity]=delaunay_replacement(ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity,outlier_vector_array,delaunay_interpolation_weighting_method); + end; + + end; + + % This performs smoothing of the vector field if specified by the users + if smooth_vector_field; + + % This displays that the vector field is being smoothed + disp('Smoothing the vector field . . . '); + + % This smooths the vector fields + [ii_velocity,jj_velocity,kk_velocity]=smooth_velocity_field(ii_velocity,jj_velocity,kk_velocity,gaussian_smoothing_kernal_size,gaussian_smoothing_kernal_std*ones(1,3)); + + end; + + % This is the filename to save the data as + data_filename_write=[vector_field_write_directory,'pass_',sprintf('%02.0f',pass_index),'_frame_',sprintf('%04.0f',frame_index),'_x_',sprintf('%04.0f',frame_index+image_correlation_step),'.mat']; + + % This renames the data from the current vector field to be + % conistent with the variables names used by other vector field + % code (ie Prana uses X, Y, Z, U, V, W for the position and + % velocity values and this code uses the index dimensions ii, jj, + % kk) + Y=ii_position; + X=jj_position; + Z=kk_position; + V=ii_velocity; + U=jj_velocity; + W=kk_velocity; + + % This renames the validation array to be consistent with the + % variables used by other vector field code (ie Prana uses Valid as + % the name for the valid vectors and this code uses + % outlier_vector_array) + if validate_vector_field; + % This creates the binary array variable 'Valid' which gives + % the locations of the vectors that were considered valid + Valid=not(outlier_vector_array); + else; + % This creates an array of NaNs for the variable 'Valid' since + % validation was not performed on the vector field + Valid=NaN(size(X)); + end; + + + if frame_index == final_image_frame - 1; + + + % This saves the vector field data to memory + save(data_filename_write,'X','Y','Z','U','V','W','Valid','-v7.3'); + %save('spectral_phase_3D_v2.mat','R') + end + + end; + +end; + + + +function [Z,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity]=standard_correlations(piv_parameters,pass_index,frame_index,window_number,window_min,window_max,window_center,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity); + +basedir = 'D:\3DPIV\matlab_3d_piv_code-master\matlab_3d_piv_code-master\test_image\synthetic_confocal\100nmD\parfor_scan01\scan0.1'; +im_base = 'im_'; %default value was B + + +% This function calculates standard pair correlations of the data set listed +% within the parameters variable 'piv_parameters'. The standard correlation +% consists of taking pair-wise correlations from a time series of images. +% The input parameters of this function are given by +% +% piv_parameters A data structure containing all the necessary +% parameters required to perform the PIV processing +% +% pass_index A scalar integer giving the current pass of the PIV +% processing +% +% frame_index A scalar integer giving the index of the current +% frame being processed from the directory of images. +% For pyramid correlations, the velocity field will +% always be calculated for the frame_index and +% frame_index + 1 froms, ie the changing the +% correlation step does nothing. +% +% window_number A 1 x 3 vector giving the number of PIV windows to +% process in the first, second, and third dimensions +% respectively. If the images are 2D then +% window_number(3) should always equal 1. +% +% window_min A L x 3 vector giving the minimum indices of the +% window domains in each respective dimension where L +% is the total number of windows to process. +% +% window_max A L x 3 vector giving the maximum indices of the +% window domains in each respective dimension where L +% is the total number of windows to process. +% +% window_center A L x 3 vector giving the coordinates of the center +% of the window domains in each respective dimension +% where L is the total number of windows to process. +% +% ii_position A M x N x P array giving the 1st dimension +% coordinates of the previously measured vector field +% if this is at least the second pass. If this is +% the first pass then the coordinates correspond to +% an all-zero vector field. +% +% jj_position A M x N x P array giving the 2nd dimension +% coordinates of the previously measured vector field +% if this is at least the second pass. If this is +% the first pass then the coordinates correspond to +% an all-zero vector field. +% +% kk_position A M x N x P array giving the 3rd dimension +% coordinates of the previously measured vector field +% if this is at least the second pass. If this is +% the first pass then the coordinates correspond to +% an all-zero vector field. +% +% ii_velocity A M x N x P array giving the 1st dimension velocity +% of the previously measured vector field if this is +% at least the second pass. If this is the first +% pass then the velocity will be set to all zeros. +% +% jj_velocity A M x N x P array giving the 2nd dimension velocity +% of the previously measured vector field if this is +% at least the second pass. If this is the first +% pass then the velocity will be set to all zeros. +% +% kk_velocity A M x N x P array giving the 3rd dimension velocity +% of the previously measured vector field if this is +% at least the second pass. If this is the first +% pass then the velocity will be set to all zeros. +% +% The output parameters give the measured (un-validated, un-smoothed) +% velocity field of the current frame pair and are described by +% +% ii_position A Q x R x S array giving the 1st dimension +% coordinates of the measured vector field, ie the +% center of the correlation window in the 1st +% dimension. +% +% jj_position A Q x R x S array giving the 2nd dimension +% coordinates of the measured vector field, ie the +% center of the correlation window in the 1st +% dimension. +% +% kk_position A Q x R x S array giving the 3rd dimension +% coordinates of the measured vector field, ie the +% center of the correlation window in the 1st +% dimension. +% +% ii_velocity A Q x R x S array giving the 1st dimension velocity +% of the current vector field measured at each +% correlation window. +% +% jj_velocity A Q x R x S array giving the 1st dimension velocity +% of the current vector field measured at each +% correlation window. +% +% kk_velocity A Q x R x S array giving the 1st dimension velocity +% of the current vector field measured at each +% correlation window. +% +% Authors: Rod La Foy +% First Created On: 14 March 2014 +% Last Modified On: 14 March 2014 + +% This is the image filename extension +image_filename_extension=piv_parameters.general_parameters.image_filename_extension; +% This is the image variable name (if the variable is saved within a data +% object that contain multiple variables, ie a 'mat' file) +image_variable_name=piv_parameters.general_parameters.image_variable_name; +image_variable_name_2=piv_parameters.general_parameters.image_variable_name_2; + +% This is a list of the images within the image read directory +image_filename_list=dir([piv_parameters.general_parameters.image_read_directory,piv_parameters.general_parameters.image_filename_prefix,'*',image_filename_extension]); + +% This extracts the effective current window resolution +window_resolution=piv_parameters.pass_parameters(pass_index).window_resolution; +% This extracts the full window size +window_size=piv_parameters.pass_parameters(pass_index).window_size; + +% This is the number of frames to step between correlation frame pairs +image_correlation_step=piv_parameters.general_parameters.image_correlation_step; + +% This is the bulk window offset distance +bulk_window_offset=piv_parameters.pass_parameters(pass_index).bulk_window_offset; + +% This is the correlation method to use for the current pass +correlation_method=piv_parameters.pass_parameters(pass_index).correlation_method; +% If the correlation method is the RPC method, this creates the RPC filter +if strcmp(correlation_method,'RPC'); + % This extracts the size of the particles in the images + particle_diameter=piv_parameters.general_parameters.particle_diameter; + % This creates the spectral energy filter + spectral_filter=single(fftshift(create_spectral_filter(window_size,particle_diameter))); + +elseif strcmp(correlation_method,'SCC'); + % This sets the spectral filter to a null value + spectral_filter=[]; +end; + +% This is a Boolean value stating whether to zero-mean the correlation +% windows +zero_mean_windows=piv_parameters.pass_parameters(pass_index).zero_mean_windows; + +% This sets the fft library for optimal speed calculation +fftw('planner','measure'); + +% This determines the size of the Gaussian function to use for the window masking +gaussian_width=determine_gaussian_size_2(window_size,window_resolution); + +% This creates a Gaussian windowing function +gaussian_filter=single(create_gaussian_filter(window_size,gaussian_width)); + +% These are the coordinate vectors +ii_position_current=reshape(window_center(:,1),window_number); +jj_position_current=reshape(window_center(:,2),window_number); +kk_position_current=reshape(window_center(:,3),window_number); + + +% This initializes the velocity vectors +ii_velocity_current=zeros(size(ii_position_current)); +jj_velocity_current=zeros(size(ii_position_current)); +kk_velocity_current=zeros(size(ii_position_current)); + +% This determines the domain over which the known velocity field should be +% extrapolated +% ii_extrap_min=min(ii_position_current(:))-ceil(window_size(1)/2); +% ii_extrap_max=max(ii_position_current(:))+ceil(window_size(1)/2); +% jj_extrap_min=min(jj_position_current(:))-ceil(window_size(2)/2); +% jj_extrap_max=max(jj_position_current(:))+ceil(window_size(2)/2); +% kk_extrap_min=min(kk_position_current(:))-ceil(window_size(3)/2); +% kk_extrap_max=max(kk_position_current(:))+ceil(window_size(3)/2); +% % This adds extrapolated points to the previously measured vector field so +% % that if the vector field is interpolated at points near the edge of the +% % image, the interpolated points will have approximately correct values (as +% % opposed to zero or NaNs) +% [~,~,~,ii_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,ii_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max); +% [~,~,~,jj_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,jj_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max); +% [ii_position,jj_position,kk_position,kk_velocity]=extrapolate_velocity_field(ii_position,jj_position,kk_position,kk_velocity,ii_extrap_min,ii_extrap_max,jj_extrap_min,jj_extrap_max,kk_extrap_min,kk_extrap_max); + +% If the bulk window offset is non-zero, this adds it to the previous pass +% velocity field estimate +if any(bulk_window_offset); + % This adds the bulk window offset to the previously measured velocity + % field + ii_velocity=ii_velocity+bulk_window_offset(1); + jj_velocity=jj_velocity+bulk_window_offset(2); + kk_velocity=kk_velocity+bulk_window_offset(3); + % This creates a Boolean value stating to apply the bulk window offset + perform_bulk_window_offset=true; +else; + % This creates a Boolean value stating to apply the bulk window offset + perform_bulk_window_offset=false; +end; + +% If this is the second pass or higher, this calculates the discrete window +% offset, otherwise the offset is just set to zero +if (pass_index==1)&&(not(perform_bulk_window_offset)); + % This sets the first frame minimum and maximum values to the + % non-offset values + frame_1_window_min=window_min; + frame_1_window_max=window_max; + % This sets the second frame minimum and maximum values to the + % non-offset values + frame_2_window_min=window_min; + frame_2_window_max=window_max; + % This sets the displacement offset array to all zeros + displacement_offset=zeros(size(window_min)); +else; + % This function adds the known velocity field to the window domains, ie + % accounts for the discrete window offset + [frame_1_window_min,frame_1_window_max,frame_2_window_min,frame_2_window_max,displacement_offset]=dwo_window_domains(window_min,window_max,window_center,image_correlation_step,ii_position,jj_position,kk_position,ii_velocity,jj_velocity,kk_velocity); +end; + +% This is the filename of the first image to load + +% image_1_filename=[piv_parameters.general_parameters.image_read_directory,image_filename_list(frame_index).name]; + +% This is the filename of the first image to load + + + + + + +image_1_filename = fullfile(basedir,[im_base,num2str(frame_index,'%01d'),'.mat']); + + + +% image_1_filename % comment + +image_2_filename = fullfile(basedir,[im_base,num2str(frame_index+image_correlation_step,'%01d'),'.mat']); + + + +% image_2_filename=[piv_parameters.general_parameters.image_read_directory,image_filename_list(frame_index+image_correlation_step).name]; +% +% image_2_filename + +% This finds the size of the images being loaded which can be 'bmp', 'jpg', +% 'jp2', 'png', 'ppm', or 'tif' images in 2D or 'mat images in 3D +if any(strcmpi(image_filename_extension,{'bmp','jpg','jp2','png','ppm','tif'})); + + % This loads the first image + I1=double(imread(image_1_filename)); + % This loads the second image + I2=double(imread(image_2_filename)); + +elseif strcmpi(image_filename_extension,'mat'); + + % This creates a matlab file object of the first image of the correlation pairs + image_1_object=matfile(image_1_filename); + % This loads the first image + eval(['I1=double(image_1_object.',image_variable_name,');']); + % This creates a matlab file object of the second image of the correlation pairs + image_2_object=matfile(image_2_filename); + % This loads the second image + eval(['I2=double(image_2_object.',image_variable_name_2,');']); + +elseif strcmpi(image_filename_extension,'dcm'); + + % This loads the first image + I1=double(squeeze(dicomread(image_1_filename))); + % This loads the second image + I2=double(squeeze(dicomread(image_2_filename))); + +elseif strcmpi(image_filename_extension,'cdf'); + + % This loads the first image cell array + I1_Temp=double(cdfread(image_1_filename,'Variables',image_variable_name)); + % This loads the image from the cell array + I1=I1_Temp{1}; + % This clears the cell array from memory + clear('I1_Temp'); + % This loads the second image + I2_Temp=double(cdfread(image_2_filename,'Variables',image_variable_name)); + % This loads the image from the cell array + I2=I2_Temp{1}; + % This clears the cell array from memory + clear('I2_Temp'); + +elseif strcmpi(image_filename_extension,'h5'); + + % This loads the first image + I1=double(hdf5read(image_1_filename,image_variable_name)); + % This loads the second image + I2=double(hdf5read(image_2_filename,image_variable_name)); + +end; + +% This determines the image size +image_size=zeros(1,3); +image_size(1)=size(I1,1); +image_size(2)=size(I1,2); +image_size(3)=size(I1,3); + +% This iterates through the windows performing the cross-correlations + + + + + + + +for window_index=1:prod(window_number); + + + % This displays the calculation progress + display_calculation_progress(window_index,1:size(window_min,1)); + + % This extracts the correlation window from the first frame + I1_Window=extract_correlation_window(I1,window_size,frame_1_window_min(window_index,:),frame_1_window_max(window_index,:),image_size); + + % This extracts the correlation window from the second frame + I2_Window=extract_correlation_window(I2,window_size,frame_2_window_min(window_index,:),frame_2_window_max(window_index,:),image_size); + + % This subtracts the mean value from the windows if specified by the + % user + if zero_mean_windows; + % This subtracts the mean value from the first image window + I1_Window=I1_Window-mean(I1_Window(:)); + % This subtracts the mean value from the second image window + I2_Window=I2_Window-mean(I2_Window(:)); + end; + + + + % This performs the cross-correlation measurement + [Z,displacement_vector]=calculate_window_correlation(I1_Window,I2_Window,gaussian_filter,spectral_filter,window_size,correlation_method); + + % These are the velocities of the current correlation volume + ii_velocity_current(window_index)=displacement_vector(1)+displacement_offset(window_index,1); + jj_velocity_current(window_index)=displacement_vector(2)+displacement_offset(window_index,2); + kk_velocity_current(window_index)=displacement_vector(3)+displacement_offset(window_index,3); + + +end; + +% This overwrites the input velocity field coordinates with the current +% pass's coordinates +ii_position=ii_position_current; +jj_position=jj_position_current; +kk_position=kk_position_current; +% This overwrites the input velocity field vectors with the current +% pass's vectors +ii_velocity=ii_velocity_current/image_correlation_step; +jj_velocity=jj_velocity_current/image_correlation_step; +kk_velocity=kk_velocity_current/image_correlation_step; + +%pyramid correlation location - originally + +function [Z,displacement_vector]=calculate_window_correlation(I1,I2,gaussian_window_filter,spectral_filter,window_size,method); +% This performs the volumetric correlation between windows I1 and I2. The +% gaussian filter is applied to the windows, the correlation is performed, +% and a sub-pixel fit is made to the correlation volume. + +% This calculates the correlation volumes using either the SCC or RPC +% correlations +if strcmp(method,'SCC'); + + + % This performs the SCC correlation on the two windows + [Z,C]=scc_correlation(I1,I2,gaussian_window_filter); + + +elseif strcmp(method,'RPC'); + + % This performs the RPC correlation on the two windows + [C]=rpc_correlation(I1,I2,gaussian_window_filter,spectral_filter); + +end; + + + +% This calculates the displacement vector from the correlation volumes +displacement_vector=subpixel_displacement_vector(C,window_size); + + + + +function [Z,C]=scc_correlation(I1,I2,gaussian_window_filter); +% This function performs the standard cross-correlation on the two windows 'I1' and 'I2' +% after masking the windows with a Gaussian function given by 'gaussian_window_filter'. + +% This applies the Gaussian filter to the first window + +I1=I1.*gaussian_window_filter; + + +% This applies the Gaussian filter to the second window +I2=I2.*gaussian_window_filter; + +% This calculates the FFT of the first window +Y1=fftn(I1); + +% This calculates the FFT of the second window +Y2=fftn(I2); + +% This performs the cross-correlation in the spectral domain +Z= Y2.*conj(Y1); % originally Y2.*conj(Y1) + + +% This performs the inverse FFT to return the correlation volume +C=ifftn(Z,'symmetric'); + + + + +function [C]=rpc_correlation(I1,I2,gaussian_window_filter,spectral_filter); +% This function performs the robust phase correlation on the two windows 'I1' and 'I2' +% after masking the windows with a Gaussian function given by 'gaussian_window_filter' and +% using the RPC spectral filter defined by 'spectral_filter'. + +% This applies the Gaussian filter to the first window +I1=I1.*gaussian_window_filter; + +% This applies the Gaussian filter to the second window +I2=I2.*gaussian_window_filter; + +% This calculates the FFT of the first window +Y1=fftn(I1); + +% This calculates the FFT of the second window +Y2=fftn(I2); + +% This performs the cross-correlation in the spectral domain +YC=Y2.*conj(Y1); + +% This is the magnitude of the cross-correlation in the spectral domain +YC_Magnitude=sqrt(YC.*conj(YC)); + +% These are the indices of the non-zero magnitude components of the +% cross-correlation +non_zero_indices=(YC_Magnitude(:)~=0); + +% This initializes the phase correlation array +R=YC; + +% THis is the phase correlation in the spectral domain +R(non_zero_indices)=YC(non_zero_indices)./YC_Magnitude(non_zero_indices); + +% This applies the spectral filter +R=R.*spectral_filter; + +% This performs the inverse FFT to return the correlation volume +C=ifftn(R,'symmetric'); + +% This takes the absolute value of the correlation volume to ensure positive values +C=abs(C); + + + +function displacement_vector=subpixel_displacement_vector(C,window_size); +% This function finds the subpixel fit to the correlation volume given by C. + +% This is the linear index of the correlation maximum (this only takes the first maximum +% value if multiple are found - this is unlikely, but should be addressed) +max_linear_index=find(C==max(C(:)),1); + +% This converts the linear index to a set of subscripted indices +[ii_max,jj_max,kk_max]=ind2sub(window_size,max_linear_index); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% This performs the three-point Gaussian fitting in the ii-direction. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% These are the three points to fit in the ii-direction +if (ii_max>1)&&(ii_max1)&&(jj_max1)&&(kk_max(window_size(1)/2); + ii_fit=ii_fit-window_size(1); +end; + +% This accounts for the periodicity of the correlation volume for the jj displacement +if jj_fit>(window_size(2)/2); + jj_fit=jj_fit-window_size(2); +end; + +% This accounts for the periodicity of the correlation volume for the kk displacement +if kk_fit>(window_size(3)/2); + kk_fit=kk_fit-window_size(3); +end; + +% This is the displacement vector of the current window +displacement_vector=[ii_fit,jj_fit,kk_fit]; + + + +function I_Window=extract_correlation_window(I,window_size,window_min,window_max,image_size); +% This function extracts the current windows from the frame 'I' where +% 'window_size' is the size of the correlation window, 'window_min' and +% 'window_max' are the 1 x 3 vectors of the window minimum and maximum +% indices, and 'image_size' is the the size of the frame 'I'. + +% These are the coordinate domains to extract +ii_min=window_min(1); +ii_max=window_max(1); +jj_min=window_min(2); +jj_max=window_max(2); +kk_min=window_min(3); +kk_max=window_max(3); + +% This checks if the ii indices are outside the domain of the image and sets them to the +% domain of the image if true +if ii_min<1; + ii_min=1; +end; +if ii_max>image_size(1); + ii_max=image_size(1); +end; +% This checks if the jj indices are outside the domain of the image and sets them to the +% domain of the image if true +if jj_min<1; + jj_min=1; +end; +if jj_max>image_size(2); + jj_max=image_size(2); +end; +% This checks if the kk indices are outside the domain of the image and sets them to the +% domain of the image if true +if kk_min<1; + kk_min=1; +end; +if kk_max>image_size(3); + kk_max=image_size(3); +end; + +% This extracts the current window from the frame +I_Window=I(ii_min:ii_max,jj_min:jj_max,kk_min:kk_max); + +% This initilizes the image size vector +I_Window_Size=zeros(1,3); +% This is the size of the extracted image +I_Window_Size(1)=size(I_Window,1); +I_Window_Size(2)=size(I_Window,2); +I_Window_Size(3)=size(I_Window,3); + +% If the extracted image is smaller than the full window size (due to the window being +% on the edge of the full image) this pads the image with zeros +if any(I_Window_Size~=window_size); + % This initializes an array of zeros the size of the full window + I_Window_Temp=zeros(window_size); + % These are the indices to insert the image into the zero padded image (the + % indices are centered in the window so that information is not lost during the + % masking process) + window_insert_min=floor((window_size-I_Window_Size)/2)+1; + window_insert_max=window_insert_min+I_Window_Size-1; + % This copies the current image from the frame into the window + I_Window_Temp(window_insert_min(1):window_insert_max(1),window_insert_min(2):window_insert_max(2),window_insert_min(3):window_insert_max(3))=I_Window; + % This renames the variable so that the zero-padded image is used + I_Window=I_Window_Temp; +end; + + + + +function p=determine_gaussian_size_2(window_size,resolution_size); +% This function is used to determine the diameter of a Gaussian function +% used to mask PIV windows to prevent aliasing during the Fourier +% transform. +% +% The input argument 'window_size' refers to the full size of +% the PIV window prior to performing the masking operation and may be a +% scaler or a vector. The input argument 'resolution_size' refers to the +% effective resolution that the user wants to attain from the masked PIV +% window, ie the window may be 64 x 64 in size (window_size=[64,64]) but +% the user wants the effective resolution after masking to be equal to 32 x +% 32 (resolution_size=[32,32]). +% +% The output argument 'p' is used to contruct the Gaussian masking function +% of the following form +% +% G(x) = exp(-((p*x)^2)/2) (1) +% +% where the independent variable lies in the domain +% +% -(window_size-1)/2 <= x <= (window_size-1)/2 (2) +% +% The variable 'p' is calculated by solving the integral equation +% +% resolution_size = int(G(x),-(window_size-1)/2,(window_size-1)/2) (3) +% +% or the equivalent expression +% +% resolution_size = sqrt(2*pi)*erf(window_size*p/(2*sqrt(2))/p (4) +% +% given by simplifying the integral. These equations are taken from +% "Assessment of advanced windowing techniques for digital particle image +% velocimetry (DPIV)" Eckstein 2009. +% +% This function is written to replace the 'findwidth' function in prana. +% The equation is solved by re-writing equation (4) in the form +% +% 0 = f(x) = erf(b*x)/x-a (5) +% +% and solving for the independent variable using Halley's method. +% +% Written By: Rod La Foy +% Written On: 6 June 2013 + +% These are the ratios of the resolution size to window size +size_ratio=resolution_size./window_size; + +% This initializes the output vector p +p=zeros(size(window_size)); + +% This pre-computes the square root of pi for later use +pi_sqrt=sqrt(pi); + +% This is the tolerance with which to find the root +tolerance=1e-4; + +% This iterates through the number of dimensions to solve for (probably +% either 2 or 3) +for dimension_number=1:length(window_size); + + % This checks if the size ratio is equal to one or greater + if size_ratio(dimension_number)>=1; + + % This sets the Gaussian width variable to 0 + p(dimension_number)=0; + + else; + + % This is the first constant term in the equation + a=resolution_size(dimension_number)/sqrt(2*pi); + % This is the second constant term in the equation + b=window_size(dimension_number)/(2*sqrt(2)); + + % This is the initial guess of value of the Gaussian width term + p_temp=1; + + % This iterates to the root + while true; + + % This is the error function of the b*p term + bp_erf=erf(b*p_temp); + % This is the exponential term of the (b*p)^2 term + bp_exp=exp((b*p_temp)^2); + + % This is the numerator of the difference term + dp_numerator=a*p_temp-bp_erf; + + % This initializes the denominator term + dp_denominator=0; + % This adds the first term to the denominator + dp_denominator=dp_denominator+a; + % This adds the second term to the denominator + dp_denominator=dp_denominator+2*b/(pi_sqrt*bp_exp); + % This adds the third term to the denominator + dp_denominator=dp_denominator+(2*b^3*p_temp^2*dp_numerator)/(2*b*p_temp-bp_exp*pi_sqrt*bp_erf); + + % This is the expected difference to the root + dp=dp_numerator/dp_denominator; + + % This is the new root approximation + p_temp=p_temp-dp; + + % This breaks if the error is less then the tolerance + if abs(erf(b*p_temp)/p_temp-a)2)&&(size(ii_coordinates,2)>2); + + % This is the estimated ii-displacement at the center of the windows + estimated_window_dii=interp2(ii_coordinates,jj_coordinates,ii_displacement,window_center(:,1),window_center(:,2),'cubic',0); + % This is the estimated jj-displacement at the center of the windows + estimated_window_djj=interp2(ii_coordinates,jj_coordinates,jj_displacement,window_center(:,1),window_center(:,2),'cubic',0); + % This is the estimated kk-displacement at the center of the windows + estimated_window_dkk=interp2(ii_coordinates,jj_coordinates,kk_displacement,window_center(:,1),window_center(:,2),'cubic',0); + + else; + + % This is the estimated ii-displacement at the center of the windows + estimated_window_dii=interp2(ii_coordinates,jj_coordinates,ii_displacement,window_center(:,1),window_center(:,2),'linear',0); + % This is the estimated jj-displacement at the center of the windows + estimated_window_djj=interp2(ii_coordinates,jj_coordinates,jj_displacement,window_center(:,1),window_center(:,2),'linear',0); + % This is the estimated kk-displacement at the center of the windows + estimated_window_dkk=interp2(ii_coordinates,jj_coordinates,kk_displacement,window_center(:,1),window_center(:,2),'linear',0); + + end; + +else; + + % If the vector field is greater then 3 x 3, then a cubic interpolation + % is performed, otherwise a linear interpolation is performed + if (size(ii_coordinates,1)>2)&&(size(ii_coordinates,2)>2)&&(size(ii_coordinates,3)>2); + + % This is the estimated ii-displacement at the center of the windows + estimated_window_dii=interp3(ii_coordinates,jj_coordinates,kk_coordinates,ii_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0); + % This is the estimated jj-displacement at the center of the windows + estimated_window_djj=interp3(ii_coordinates,jj_coordinates,kk_coordinates,jj_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0); + % This is the estimated kk-displacement at the center of the windows + estimated_window_dkk=interp3(ii_coordinates,jj_coordinates,kk_coordinates,kk_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'cubic',0); + + else; + + % This is the estimated ii-displacement at the center of the windows + estimated_window_dii=interp3(ii_coordinates,jj_coordinates,kk_coordinates,ii_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0); + % This is the estimated jj-displacement at the center of the windows + estimated_window_djj=interp3(ii_coordinates,jj_coordinates,kk_coordinates,jj_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0); + % This is the estimated kk-displacement at the center of the windows + estimated_window_dkk=interp3(ii_coordinates,jj_coordinates,kk_coordinates,kk_displacement,window_center(:,1),window_center(:,2),window_center(:,3),'linear',0); + + end; + +end; + +% This is half of the measured velocity at the window coordinates +rounded_half_estimated_displacement=round([estimated_window_dii,estimated_window_djj,estimated_window_dkk]/2); + +% This calculates the first frame minimum domain values +frame_1_window_min=window_min-floor(rounded_half_estimated_displacement); +% This calculates the first frame maximum domain values +frame_1_window_max=window_max-floor(rounded_half_estimated_displacement); + +% This calculates the second frame minimum domain values +frame_2_window_min=window_min+ceil(rounded_half_estimated_displacement); +% This calculates the second frame maximum domain values +frame_2_window_max=window_max+ceil(rounded_half_estimated_displacement); + +% This is the offset that must be added back onto the measured velocity +% field of the current pass to account for the dicrete window offset +displacement_offset=ceil(rounded_half_estimated_displacement)+floor(rounded_half_estimated_displacement); + + + +function [window_min,window_max,window_center,window_number]=calculate_window_domains(image_size,window_resolution,window_size,grid_spacing); +% This function calculates the minimum and maximum values of the window +% indices to extract from the image. + +% These are the number of windows in each dimension +window_number=ceil((image_size-window_resolution)./grid_spacing); %ORIGINALLY window_number=ceil((image_size-window_resolution)./grid_spacing)+1 + +% This initializes the window indices +window_min=zeros(prod(window_number),3); +window_max=zeros(prod(window_number),3); + +% This initializes the window centers +window_center=zeros(prod(window_number),3); + +% This is the full resolution covered by the PIV windows (which will be equal to or larger +% than the resolution of the full image) +full_resolution=grid_spacing.*(window_number-1)+window_resolution; + + + + + +% This is the distance in indices below 1 to start the windows +negative_index_distance=floor((full_resolution-image_size)/2); + +% This is a counting variable for indexing the window domains +count=0; + +% This calculates the window domains +for kk=1:window_number(3); + for jj=1:window_number(2); + for ii=1:window_number(1); + + % This increments the counting variable + count=count+1; + + % These are the minimum window indices + window_min(count,:)=([ii,jj,kk]-1).*grid_spacing-ceil((window_size-window_resolution)/2)-negative_index_distance+1; + % These are the maximum window indices + window_max(count,:)=window_min(count,:)+window_size-1; + + % These are the window center indices + window_center(count,:)=(window_min(count,:)+window_max(count,:))/2 - 0.5; + + window_min(count,:) + window_max(count,:) + + end; + end; +end; + + + +function spectral_filter=create_spectral_filter(window_size,particle_diameter); +% This function creates the RPC spectral energy filter with a size of 'window_size' using +% a particle diameter given by 'particle_diameter' where 'bit_depth' is the bit depth of +% the three-dimensional image to be processed. + +% These are the wavenumber vectors +k_ii_vector=-pi:2*pi/window_size(1):pi-2*pi/window_size(1); +k_jj_vector=-pi:2*pi/window_size(2):pi-2*pi/window_size(2); +k_kk_vector=-pi:2*pi/window_size(3):pi-2*pi/window_size(3); + +% These are the wavenumber coordinates +[k_ii,k_jj,k_kk]=ndgrid(k_ii_vector,k_jj_vector,k_kk_vector); + +% This is the squared radius of the wavenumber +rho=k_ii.^2+k_jj.^2+k_kk.^2; + +% None of the coefficients actually matter . . . so RPC filter is just a low pass filter . . . +spectral_filter=exp(-(rho*particle_diameter^2)/16); + + + +function display_calculation_progress(current_value,value_vector); +% This function displays a progress bar showing the percent complete of the +% currently running calculation where the calculation is being iterated +% over the vector 'value_vector' and the current iteration's value is equal +% to 'current_value'. For example in the for loop +% +% for ii=1:10; +% commands . . . +% end; +% +% the arguments would be equal to +% +% current_value=ii; +% value_vector=1:10; +% +% This convention holds for non-integer or non-monotonic vectors, although +% to work correctly all values of value_vector must be unique. + +% This is the number of characters to display in the pogress bar (denoted +% by either '#' or '-' characters) +progress_bar_character_number=50; + +% This is the index of the value_vector that corresponds to the +% current_value +[null_variable,value_index]=min(abs(current_value-value_vector)); + +% This is the percentage of the calculation that is completed +current_progress_decimal=value_index/length(value_vector); +% This creates the character string showing the numerical and graphical +% progress for the current iteration +current_text_string=generate_progress_string(current_progress_decimal,progress_bar_character_number); + +% If this is the first iteration, then a new line is added, the progress +% bar is displayed, and the function is exited +if value_index==1; + % This displays the portion of the string showing the percentage of + % the calculation that is complete that is new to this iteration + fprintf(current_text_string); + % This ends the function + return; +end; + +% This is the percentage of the calculation that was completed during the +% last iteration +previous_progress_decimal=(value_index-1)/length(value_vector); +% This creates the character string showing the numerical and graphical +% progress for the previous iteration +previous_text_string=generate_progress_string(previous_progress_decimal,progress_bar_character_number); + +% This compares the current progress string with the previous string and if +% they are the same, then the function exits. If they are different, then +% only the text after the difference is displayed. +if strcmp(current_text_string,previous_text_string); + + % If this is the last time that the progress bar will be displayed, this + % prints a new line + if value_index==length(value_vector); + % This prints a new line after the progress bar + fprintf('\n'); + end; + + % This exits the function without changing the progress bar + return; + +else; + % This is the total number of charcters to be displayed + string_character_number=length(current_text_string); + + % This is the index into the string where the strings first differ + first_difference_index=find(not(current_text_string==previous_text_string),1,'first'); + + % These are the locations of the double percent signs '%%' + double_percent_indices=strfind(current_text_string,'%%'); + % This removes the double percent indices that are before the first + % difference index + double_percent_indices(double_percent_indicesCORRELATION_WIDTH + x_max=CORRELATION_WIDTH; + end + if y_min<1 + y_min=1; + end + if y_max>CORRELATION_HEIGHT + y_max=CORRELATION_HEIGHT; + end + points=double(SPATIAL_CORRELATION_PLANE(y_min:y_max,x_min:x_max).*WEIGHTING_MATRIX(y_min:y_max,x_min:x_max)); + + %Options for the lsqnonlin solver using Levenberg-Marquardt solver + options=optimset('MaxIter',1200,'MaxFunEvals',5000,'TolX',1e-6,'TolFun',1e-6,... + 'Display','off','DiffMinChange',1e-7,'DiffMaxChange',1,... + 'Algorithm','levenberg-marquardt'); + + %xvars is [M,betaX,betaY,CX,CY,alpha] + + % Set empty lower bounds (LB) and upper bounds (UB) + % for the least squares solver LSQNONLIN. + LB = []; + UB = []; + + %Initial values for the solver (have to convert D into Beta) + %x0 must be double for lsqnonlin, so convert M + x0 = [double(M(i)) 0.5*(sigma/D1)^2 0.5*(sigma/D2)^2 shift_locx shift_locy 0]; + + [xloc, yloc]=meshgrid(x_min:x_max,y_min:y_max); + + %Run solver; default to 3-point gauss if it fails + try + %[xvars resnorm resid exitflag output]=lsqnonlin(@leastsquares2D,x0,LB,UB,options,points(:),[yloc(:),xloc(:)],method); + [xvars]=lsqnonlin(@leastsquares2D,x0, LB, UB,options,points(:),[yloc(:),xloc(:)], method); + shift_errx=xvars(4)-shift_locx; + shift_erry=xvars(5)-shift_locy; + + %convert beta to diameter, diameter = 4*std.dev. + dA = sigma/sqrt(2*abs(xvars(2))); %diameter of axis 1 + dB = sigma/sqrt(2*abs(xvars(3))); %diameter of axis 2 + + % Find the equivalent diameter for a circle with + % equal area and return that value + D(i) = sqrt(dA*dB); + peak_angle = mod(xvars(6),2*pi); + + % Calculate the eccentricity of the elliptical Gaussian + % peak. + % dA is the length of the major axis of the best-fit + % ellipse. + % dB is the length of the minoir axis of the best-fit + % ellipse. + PEAK_ECCENTRICITY = sqrt(1 - dB^2 / dA^2); + + % These are the lengths of the projections of the + % elliptical Gaussian peak onto the horizontal and vertical + % axes. + dX = max( abs(dA*cos(peak_angle)), abs(dB*sin(peak_angle)) ); + dY = max( abs(dA*sin(peak_angle)), abs(dB*cos(peak_angle)) ); + + DX(i) = dX; + DY(i) = dY; + PEAK_ANGLE(i) = peak_angle; + + %LSqF didn't fail... + goodSize = 1; + + + %if D1 or D2 are already too big, just quit - it's + %the best we're going to do. + %Have to check in loop, if check at while statement, + %might never size it at all. + if 2*D1