From 0e74686b62112e15dc5edba60267f3a46d84ba3d Mon Sep 17 00:00:00 2001 From: "Jun, Brian H" Date: Tue, 27 Apr 2021 22:28:57 -0400 Subject: [PATCH] Add files via upload --- calculate_phase_mask.m | 140 +++ calculate_phase_mask_ellipse.m | 137 +++ calculate_phase_mask_ellipse_3D.m | 142 +++ calculate_phase_mask_gaussian.m | 140 +++ calculate_phase_mask_gaussian_3D.m | 164 +++ calculate_phase_quality.m | 57 + calculate_phase_quality_3D.m | 71 ++ calculate_phase_quality_mex.mexw64 | Bin 0 -> 100352 bytes clsm_img_gen.m | 137 +++ findGaussianWidth.m | 140 +++ gaussianWindowFilter.m | 61 + phaseOnlyFilter.m | 22 + phase_quality_mask_ensemble.m | 152 +++ piv_3d_10.m | 1552 ++++++++++++++++++++++++++ piv_3d_10_ensemble.m | 1663 ++++++++++++++++++++++++++++ subpixel.m | 386 +++++++ wrapped_phase_difference.m | 20 + 17 files changed, 4984 insertions(+) create mode 100644 calculate_phase_mask.m create mode 100644 calculate_phase_mask_ellipse.m create mode 100644 calculate_phase_mask_ellipse_3D.m create mode 100644 calculate_phase_mask_gaussian.m create mode 100644 calculate_phase_mask_gaussian_3D.m create mode 100644 calculate_phase_quality.m create mode 100644 calculate_phase_quality_3D.m create mode 100644 calculate_phase_quality_mex.mexw64 create mode 100644 clsm_img_gen.m create mode 100644 findGaussianWidth.m create mode 100644 gaussianWindowFilter.m create mode 100644 phaseOnlyFilter.m create mode 100644 phase_quality_mask_ensemble.m create mode 100644 piv_3d_10.m create mode 100644 piv_3d_10_ensemble.m create mode 100644 subpixel.m create mode 100644 wrapped_phase_difference.m 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 0000000000000000000000000000000000000000..75d85904a1a33c2fc8c7e1c05276f230d6fdfd4a GIT binary patch literal 100352 zcmeFadz@X<{r`V*L0ra6!XZL92#x!=i|d38gA9WPMHttj4O5~YC1GNuGA2ilN-8eH zph%G#b&I-9h)fVQb*Ug2Ejp2u5mZX_d|#J+)|r(=%cuQ5e!oAyeLU*5*Zb_(bFa1b z`@PQI=j=17J}#H)lS-xf;{W-4Dm4?z|5fO}|N1}mczNwje!q6=fz@BwbY`!SFKqh# zGtW9dQ+M9EKRoZ`F_}|OKIfcs$7W7BJ#*d#=VZ=0Co|&kqcUU8J?-@EOG?%llAvoU zQ>oK_-Y<2^ftwy%I9hjV`^=iXd+(ImzE>)B2%`P)|FY!&Zn@pt{ibh~6dL?rCp;Vs z{^tmd_+KJnBeau9N3T>ltf^M~XAL?0pO#*!rHmK;*M|KFuIZH;SBU)=bo;TVUpN+= z@@;W2^i7(Kqkzw(QnlNka{l>j*RP8IRP^i^qyArgD%H0Ad8eH`_GD}%eUI zs~Q$(MmAPvsv1WPtZJ+-(=JsFvop^=ozv?R7T06{p1X`}8a;4i|7m-*ni6+uHz4?f8-pdI}HPFuU&Arx{I-Qd1avK7EKg z)s0SdYlzvwU+jVrHH}(Y-MCINl!RH;gc-}0tbcm1R6e}{j4+jr%%)dGpt@;5b>mU} zVBm2ocWf(k;2$V2cjXIgVZ!F0)+=>z)0+A8$9$z>b`9p89{w-sUj@(0a>vf5Pl6GL zp9aaNuP4lhzXE~e@IKI0;gD0fl`T14mS6aI6=AZ=c^F%M=Y~{fm9AY&ssq=+L+Zz!_+!_|V$=#zP;~x5u$>IfAqJ zFOL1<@ZzyC-An(eV^Mr8ZTVQEC$+2NJmHJB8F^|l>`@r@0G#+aIGIcKFD=`)ggY!b zptS$CZ{%ed(WG5R>@|AQxVe~+Dco~-X&ENI|BP-NYDsBXK~KjiI`(!4PgsoYG4)GE zmTvo2GRqj+9>g;vIj(ZfZ!*TAh2tdS<-0M+qw2BGraceYV2vV&Zws!n8sQ$-aM$lAWnvm7RK~&MNc$xcyeT`#RC?wPE*d zR=)dtSRIR{)aqoJaKK8@GoyW{d*4ymci~C>H4S4jn1(?}Ma=xMnE9=5=eRIxOmhDj zRoKSsUt%i5ndwZDRC4X8(zElexC`pV7fKT#sv`x^&xjuqezM z(S-juHqXY34UhA`y~eGBdF2+0Ywq|}M>LI}x9Ffkq85q-XTIcE9CJid<=o`Gys3I# zK7ASv-Z=3wx|^CMjm>RpJuYb+KQ|uSJPxi|Qq8aOV^=2v*oF4zeBc6LAFd;E?YQc zXN!8GgqgDPAzy_f<1Lv3EsVN)N`I&Woi@{o zYs!>wku^&X30X7oSKlFPcv#GpCv;=#dk@XJ|1+k#AK&_hoqT|q?!E9Ny7msX(q3H#J@RMNyFU=gRY*OPY*yBzE#rD z_@AsblbtjpfKk#UYfUnjT5BeMi>w)dVW1m)CpsK?_-)^bnt5DgjpQ4%X0oz|`>cEp zAFNW=OfF1dmb>QDf4Jt~WzFPb)@ZGv@4KYUMmR7=outjtI35T7t=^g1kTz}KoHp!< z(xz?YLw@QYr44s0jJmolpMLX~-!g4z#Nb<#n-X3e;Ld>Jt95N9SL76W9^l*6~EW<*_gC>>h9>)k%h|)O(icp(61Imc?|%(_nXRVkm{x*abY!J zMCrX#xnE_km#a^G6L(*JZNP}8>bWEOWNny8dv#MqzXO~4uEJwCjnDOYoVJFKdx6ki z*=yGDhK0SVPknY|lXgBV*%3D@W!uh8Zdod4V;Hk>nekZZDr>XPK*%wf_v_*|zrA{qpoA7iYRl6z8J0l3e6f4O3Z6O5T~qNm+=jn3U%XU)hxRxVgP?$CW1~%MVFu(f&?K z2P&sWWBZ1r{9{GWq+}<pP#eBo@^b83!^`n{!{hv~@0hst?`{o8 z>B*6B%eMBherST83j2bhBb5F(yQ+N7^2UB z;JbB2-}mg>sOTFgU*kv2-RkMNu4;UZgRRzF*8$qfa(&Nr?dHarYHa0wrscXe@m5Q& zE6pwS?fEP50q}&`ZMf9?bVTVT_~@2`k8T@gm*M~YhMzS2#FPBlI9E7+{+xW+J7IAh zZzp`6E5rR-6a8wh zT+CzQ&1o^-oO?EjTsENpyi9lb*{uU;bm#%EfRmM%^qJ_s5ah!8h!u zRq#yovT&0QEqv&>54uwPSa4=lHp?AzT|eAGL&Z!V-*0%s#eI8~{$du+%!((7IgT@I znQG~If~E{+h8E7z98CSpR3$lZ6+?VM!I@eygfTi1B*+sz%w%5@z&kjI`+_d*l}3ur+&W}iW7r2)EDe*^E0W2OVYqwu41 z!yn~$*1^5kaea$j9#ZV`#eKg`pO@|ZjeWj$efGJP2U~O=-eJZ}XHVLa634_u{d3{> zjGi#2&4y$aQAHzqVR<+eV?^@P4``d;t_uD3P7mZOTQMK4JqKBO72o0Vl?(F7%*y$A zd4BQbL8Vvy2E+Bqemosh()g|Ua!es;BMjIaIe|<^9Ev$5ld7nJKYvmu}euy*B<4{f2 zsJn+6cx}M#>Ar zx4lNsANS%)I=*|yvZdqaz*+17R>>$0nZK9;s)0h^~-&muciHyjj<~7*5?xeE?E<$J2r#Qrrh^{j0p= zNi1L&}_~Q(4`KZ|Y{^nyWZVzZ;+)78k_A;~&Kpr+IIg z{q;K@7aDv{k1_jAm<4`7KfFz_i8HVV9QpKmm*Ch*@h49?EWmA1eXeWU&2@t`GuO2l zwvt6@;Ri4m_vhj~KeqCpPQ#9`)LtAsA;3? z@ev|-PfqU0Q?qtd&MnO1L8bSuJ1Km)<9)dr`gjvb?fJnDK3<$i1U_Q)$BW#dz|tj4 z_%KHon5eb+rIR%=;eB91>EyLAQJ9m?T-QM43by9vy6)aoTOGNsTezhg!r$X*JKPZR z-p==be#qzF-1LpE$9$G0*>@pNW43Uc*MhWCc5Z zRhV_?+e#L-L(!GuqQ+Y!?wISk011IpI*>yh-mn*L_y!>(Ry;o$&w=It!VTZ>cHE}n zz{%5KSQLJ?(@H0=*!5!utLud~mx1mf#b-%dJ?_i3BW^5D;~w#eD*1QETC@@FFpm81 z-brql+(_1RNPd0ce{v%kuJfA)_^y#(?)-~cT1^?Y!dp}E zj6AbpF_~}-gO6_V_cGC2Om=0JSo}4#S3x~$ym6# zJR^Bw`J2oC=vM6rL7OW$xye(ih1<);tMlmd*@=a?hi0`d;1xbp_mpgNkp9FFH5A>t zWe5KWPr7q>%imi52RFmF3AxhV?yKJV#`3uT;KmZS6)V4)%vY{>13BoDsfFbOXB;=8 zybkONnbl`d;Z|;id&}ec;*(Ll#f6*8N#AsHIqBQoTu#E&^tj`zD%@PQU7(vwQP4RY z2A|C)pVIslnc z(kCSPeD%BDRIdCU?a$lx>~L}?ScVQ`YWXPuK4kf(n@fBk+Vi2-X&|kn;|cO-`k3Tj zb^JTtTyEu#`+{;*a&!4p&UAKgc58OUS2F+CZ!TN5DUN%&n@hhZU$UM~|F?Q$`~S<$ z%S7{W#|7)x8{@?RP_~xzX4}AsTn{O^JKJUMK zb9oz%eQ@E~c7G#UxD@*k8IJ*+i2)=Ro+%HB;vs3PFr=RQhWp6@`hcSd2`{X0qq%`* zVflCS9$(Uav5wLZcVge+UXza)6h7mK?lpUSoBM5VH~EOcr(5S56- zu9n~8j&rh7$b;vK0RKwwIDdUE4!^=3XYwJ?!}x&!Zm@9EIRjh&@1~PCowp%reABs` zw{X+BEw^;jS$WP1H=Vg0rr*DJnXs7K{1wS^-!}3$NBm(<2Fl+a;q~~pN0)2$@o$g1 z3-2j9uUY)GsON8wjvs>_G*&_ z?Zu-Wj!o+k%hit`n7Uwk#%NpRV{353R3*lF|)tIfZ-Z zneCO!Xhjh%!;7`9e~Z2l!)NFRgL%-E_!e~1_Xod`uUytJe%WT@{xC9I{f0h^X@6=n z9$VA>(T{ONeIZtj$Md?_YS$xcJO;YAWzX*2a<;m6tdA*)d(b#~x$f0?SIbmMw&V>| zR-t=OyL;=ad-$RY9~HGfg;g%;-r=Ep*$aWf_hJa3<~LDs;w$9r;9a~q3FChDdn`qD z-yb+tPs(-8{u)>9je61ob6Uj?C0{gWSD6U|9!p}skdgQ{{FZ02HDNYSFMevJyyHah zfBLUABeTjLW<@7n!`Db^Gb7mM{mez~IX-X-XEgb!qmI2yz8!Cx z>Q>M1PQO36!}I1Zw}&R{Zx5T_3>YYU3b71F){eM7?C@vYh2|^E-Rn+J$D(zA876cn z=4gJdtM}T-zU<&-oNRoH-(Iy~h^|QQrHj@?q*oGYsQxT9_NMyr&!@(&vhZu{(xZNM zMPEbxq|Q`=_AH9R=aBw5(vKwcm~iAm&vb5@s#dn-1B|YaqmOy6OO>OK6%wrdsZ#a? zQ<6_lKTCPU6h#1qbTa~&qhpgf!sZxnB28EMGv2U@`T4*cFCV|%rNGFRP5#Wb~YYi&Ka8`6< z+INtC`I1ln76ag;L7WbLx^L}2BbnL+g#_YPiHhpmL_Ghu6eo2lPU^DK$>(rB3e}X< zPX(91Dp zyc4;!LRRj_j`1fL!yCTpb?JI5$;V}jRzqqq2OFxFrN*wld_vwCCFE`OqmL7-SfMP|;_^e0@KSA^_5gT>$9%t!5?`@eX z$(FPNTgcKzXX(bu(ne%y1!uqgDSU*;EX4!By?>}IRlM+d`QOh_qQ0|*iW#fa;9WGC zUMRDxoZ$@QQ1F2uS5GVevz4)B?vbQXUu+?_kM{oT9}bG=k`BktIG0GU2G;1o2Djn%0A}Pw<3m#yt_E3 zZQ3l%>5*u`!=1yN{`M)&X+Ist-?lY2tC!hT_V@m5c}~B=&OPQ7&yifP9uC+W{V)8k zZp{3Jn_=D%9Wo@G+F{9NVJ>u9HVl@|q-Yjxijku{@Q%S}Gr}iU$ByabV+Eb4@g7IV zXC=Ky3>&@llIQq=QfJ{6P`BQEN$H9F;olgL{T3G;@dWd_=)-@0!2k3ggSSJK?xI~+ z-iwvFjvX^!)^F(9Lp;L#mG*2tX6Px*6qXu1Gr_&GgMW`{I%9b0_S;<4`@G(Z(uia> z_Z)iq_s7{saGN9uTP5BDE=L}s6G}vcW-uYqp1~YMTjh?$kn0~WE&TiWhQ+n~VKqL} z-h;Uf>OM*p{N|LObMwO*Xoia4c=Pjb)Yz8Hhiv*Qbki4c!%oFk!|Ymdv5)?j{64&9 z`HK};2IkXi)6!VUOXhBj9VZl*)J^!MKYq+kXs+u^P%vKX`{60V;$RI2U|mbnMBDgiE)s_|k1LFWvaEW;nSG#Wu4YlQC*Oe%F9=J0BMdNAMAbGEdh# zC+ghH#a?`uQ>7ZgmoO7m^rzmSg3_t4j-^rMa%b2zmNm z<{2quc=4Vu$-}zka~-Ld!35{ipZZL80?Z|cXLKMA-{r#(Jh)-a)CFs!^DMcc`Sd*O zUO&D(H5P9utK5vT$4e-4p1?=&65V)M93Z8{!Tz!n(QaJpXHdE8;_@|vuaOSRZq$t9 zr$#iso=>0aKCphd#`5gIK59Pyz6iGuZaB1G)QbassBSzfTSf|Z;;i7xR`1Y@)=%E} zyLq2>KzrZoxfd3-4XLW}Tp# z=sA=8@@#2xSPW=jVK_2Jwy3;tWLyJt zk<|crw4+vjw899Nj+=vBhmL4u|CVW%w_1sVC#tHUy1O{L<@o&Mc%R6F8}P{guELbw zI~r$?CkA!I^$6*BU%Lxyc>To#Q)5f~8ve@5vxyx;>`t->Cx022ZSJG<)MIgKSggOw zf|ts4z=OD@9~=-%3tOt*b4_-ewRvI^(5puj8csSKBKdnR_%klw`Vea^e~7P+(_~<| z#n3(~G*)tmZ!WPMT=6@|zbvt{uzM)6*YU++vBZw{@&6YkwwC+Qf$98qCH9veXx#s` z68q-q`3a8n9j50H>>j3Pa5O!;yU+iN>G|8S>ch96p4l*9F8}lC!8J6>I{c68_0h+K zhm+dy!tzUZ@(oYJ__AK(3s02y#q+^fJ;pzVYRzC3!J`V}TW}pdW<$JtW_QMF&#UfF zt8T}Y$LfgTQURY6zVjgNVaI2Bk3AH-{t;na*zuBppI!0|ws7grr3aVpYf0Mhj>=1z z@pNPd|6>r2!RnZQOs~|U{n)Qdh9O>`0EYsUArJuH0<^I+0IYAo%mUzA01Dq<4Z%{o zz@Pl!mo)q}P17j+zk2B5_+iXp_+d<2@lni~<6gpLc?Yf!+B~ppzTmC$0vsiUAB`;s zI;3eNK6p&;go7dNCV@Q(%U~`!#y}h+cX;Dk_zZ~O@bgK}$HqNo>_eC@yvMe&_xdT( zNhkWE=4k%8s2a>UZ@f0(aN;@V2aW@1#vi{bc;Kj5{?y_xSIyvsczVqDnO@aH@il*j z<&U?YTvsQr6dGn%^!!XWTxTXXEtq4D8h?(f@NGc*tYHY`(=Q#R@$wCx*Os_C8r6+q z;lgs}C?y0gC-ChKQv>qTMw^5I|0F^ zZRPw#I})9nXj`H)6Kzd&dZI0fPD@nv0p6sfU2iRZzTZBv4s(`7`I%m_SzFvZwzyeS z+~fyA$-eyPBiY=exLID@99rDuy>PNGZ;z8rez1{j^7D^mGgaK=BdbZ2UkoOjd|Wcw z>@05bk;o+4QQVwc+-xgu&Ma=W7B{CCH(QFE(~6t9;^w5{CV#P*9QeZGW?gaf%;IKk zar4;XW=(Oks<>HE+}xwMSzg>6THMSOHwP9s%Zi)*ikqq8=2DbL()I4*W>;~O^^@!{ zzqr{^+?-q7Y%6ZgEN-?IH>VdjTZ)_0ikrFO=A`0geR1=`;$~fO^UUIAZE^G1;$}^8 zv#PjRQQX|4xLID@99rDW6gLMJH_M8f{fe8Z;^xw2Vfwm@n_b1t&f@0$;$}y2b8d07 zt++X}xY=6VoL<~)DQ-?HZsv-clZu=5#mx(gn{|cFY{@cw-j(b61)3MM{t!axk3K>?60w3g0DM*|0nh{5d1IU51PS$xcon{ zzrubAek>aNgY9P^_x&0FS`N7~{gC8z5O=TeXkEFkS*1EeLadfV0 zFZ(n0OYq@+@IR>=+^Gx%{{Z-Rd^XpWvbD@s2|hOu)}PtRK(M}cxUJ93(5$CymEboS z!Mc;J383^VkTCycO;-y^I*7s;aTkMzMuQQ;9%K(3VvPxf1DuaP>kQu+( zgp3dF79(ed1b@pRV%+L}Fc=tp%&1pCZZ)Ic3`t}2q=<2d84L!-lhqna*T(F}iDuN9 zA;I^Y5o2d`7$=Rvz_^%<_E}fGUv)&T>ohZJ&5)q1BgT6;7)!e+V|Or^+FRIOW4i=Z z9oet9ox#wBr6Xb_Rp}k;Hzd?d7&hP}7mUuk8#5`$e+n(^Mp4 z#&!ux8ZC@glI1T0=N;-aAL+fzytf2Z?7b!UWko)H6ZS@$!#A5-jT2MeUxLy_tC(r0 z4%W;v7&`Zs0YEhv(ETM}ZI+;b(JBHqHh{qZxRJJO$uk=V_h?@#!I=b`^{)GLxP%*7{{BzU|@6{ zqOnZET~Ds-Tr*nDkf8J<#t3v6Cy&9vm`q0dtdU(z_Y2HuF++m#j~JihU|=v97*${} z-FLG+XS)QIAKCA=oxxyVJ+Xgrr)IL=b_t3;ve(+qV6Z=Tu#Y~|_Bz`osQ1X8ww=LX zZ<0Nqrs6!-+AcwH`0B6KK&}U8YfnG ze+fz%Eo5twX=8Dlhhs4qI`@$QKus~A+yDuRIRflp0D}QAgSKqRpEh8MjJ!h=oH0X! ziZ+9H&`9?!Jwprz##v^pBSt?n%FK|U&LhTb_k+Q}*usok@e-arx?TN9nIS>#M~okv z!C+v#dyvLbJ&=rgGrIr59hxO5f3z^6$yi3A!1-G3CcfW ztZD{>f$;+{nC=JJ-eJ20l^@w3LtoJa27`T2V(({ro9z-5ePsWM?F&Mv0Z{vM+@DR;P-F%e6rH%evbFfd2b0S+j~p&Bt;wc zM!LiI0k;|_)_Z>mN*S%<*_pV$riQ`LIsAl=fVO=4b_41Rkf4|&K!pJe2EY@vWlR3L z9#drWZJME4GbE^JGx)g?(tUT&5QBknkr{tXld*vrHD*Xq=Mm#o_k+Q}*v*Xl@lvkq z)mzn%3Ns|A{fIH%3%6OP2q8>7!)c^^KIRYGE z0D}SWB5m1{ckwJxuIro@%}|RO5>&Jqt2o`OJVOiyMw1y+)*)kSGje7~Q0Ecj1NVc$ zz&OZ^*?0-p|2L~2^=3#=`w`S&?65X5)i|*mpQ^K%DJW)0 z5paqD33)&zIol<_!76V6ayw_PuPcw_Spw zkL$c0+lB`9^Y&|OKkdRjicX|ne% z-n+(oOHkR~TavBL;W>uAIa>NQ@ctFvUxHFb3)z}v+Dx2A9E-uwxt(akcR%Stl zO5lfGlRjvILnOn#OP;6$_xo=KVr0cEDQ$57H0f%P4;8Ax1m%wwCNvq# zbr>`zl)=DwcQ3@+XPtNW+lS@4o-?D<3<=6VVjN-ygMrZu2GhOC_72-6sQk#DvYo+T zKQXb7w7t!C35q_lKZxnV2pA0ZO%i*F?X9*;Q16lbY}*+O_Rc+>?o{MLi|rDWI$AiM zBwIZR-$y2U|IT~oytf3E?Y$*>G`}zRYNBsz?_clzB`9UIif89poJJgr!O*z_45-hi z7YI{&em?OaD1~3=^y$Q&cth)wN5Sjx?ip!_4o zT^KYbl)=DQvIk=Av)*`qCwwyVrWqMCBq;xgal9D}2F7o}V7lLEdztMLRDNXdZ##p* zer{qPZF|ag35q_l&%u;q1Plgyd1Bw#_U_O5TC)W8jutYia58OYFxWrY-RVw6x_8T0VQ5)^jCXf%Vtz*vopY)QG&<&G&zmliW5sBSaxhYEOMQ}5YgFfg9x zp|d3~ufnuF#f+R864ZRe*w_pP1LIONmZ={bn^A9u1nWS=nB%c97#Mq*(Odm^;~I^n z&I}2bf{1bNYMN691LKQb5o@2d_q279`kiLfnjyhb5HV`aU@$Q52ZO19i|sYGORx$= z_VsLMFxbZ?_ET)Huw8-$AhJKdsz$(IugD{bkC0dR^|=q~GJ_)4SW1vP**cMhmH*q}-oy zK3U(|CHV%k#tc3df$OE6)zxl1_{}1xpvWOb##Lr87)H7l8A#jSOxwG!R=RYWAwi{^ z!GAvysXx)P#b98}<)O1Blf?L;869RwQ2Y^Nb2AtWjLXf~R*cQfXfs2Cr66KF@3Al# z7!_vR*NbC$2alQH)U}!+!BP-09_Xt%WiT+7?SxqStUL6H&HHAwm?6PZ5HU_SgTcU< z2?kUD4%>6KORx$=_6==kFxdI`gJ}Pu?e(@xumD8%jy@UzgTcOMVjpUIo$V6Teq^6$ zJA=Xg)s9YmD(0`&b_vQJEo7l2n>|*EcSvG*&W0KrBq($nB-tz<)xq^VcWg@Bd)rlE zmjv~V7Sc9Jxi^s$I4Famm!r%uaA~=^T5g5}MIJG(F@wRt=ubwrWS@Kr>C$?Y(j{Yt z1eI>anoj*IJzES0#tS?&KHpc${n(5$GbAYfh%wX*1_Ptfj2+dFp=PAakYFi@7#$u9 zgMm?L#@8vP{(D)CrTbIv&@91HfEMO78OtN+Fy@rO!01Ir;nSNB&FC~kf~6p0oNWe! zf${hb&fU9h@338hRUoowY-cdoFG=h_w!O`E2^N6J{wjuw5il6+`z7}6Y;Uz)g4&Pl zSK7{Cu%~3tr>RK&7TYB#d$e#oiGOm4e>Y}(*UF1f7;-j9Q0O*D{CAo1>5aK#lYX$6 zPY<`N-YyC18?EAL>O@Z9pbUmy9%lw0aLcE^-cCoZGed$Rj~LU;U@$P&BO_a~pHl9@ zCZ$WQ84^^w8GJw#tNt%MTMP!qOFVS8K0h&|#taFHKVocW27`fdl^MH!)-cqLAPtC}gAwls+jGfG2Ffgt$qvA_4b~2;h z3<;Kki1CKU!eC$=Zbpgv@zE6;OPv`KECmtcY19bjl)=DQlZ^IRr=7QdO|I)RGiuF{ zU@3?gW6fYNFxt0r>OW|EjqMVw0+D@d+Zhb@Ns0ZZwpZ9L!2%H3-@$ZY1Plf{|E>tH zmv*c0+!q09ye3f%@t zHp?fA^XZ}7u_8uCcV5A;D4*G3KI1a3UEDj50C`pYZ&{ zj21H_SPCM>&&*&jFy?LPdmamef$;+~&cqK5b6ra(Xe`~`+@V>5r2s9=X)=}AF)!C+rkc6<&4MJBh{E484SHV#|#7e zn^9+m1VtV(?lFVGz}T9$Y{{^XnJ&*>s&uI}LxM^-qs*y)n`eu`!1!QOF+TW+jEQE{ zm?1&&M~vZSFc=uOnz0F9%5@DlqrwacmV$_}&|_gRFitdMu=?@!B^pb)84@f75#uG) zDo!MWfiZ}T_E}~6*0zuPkugJpr66K7n89FR{3+wqpY1M`*)G8<5ZQ;>&S0=#o7gAX zp0Ztn1t7A2geu1f7!3Af6Z?U-cjI5_V(n8<`w%2m;bhv*V6bl>J3fbjA}2d-m!Ry? z!to^8>{9$1C!aoZ6Q_O`6ow8PBq($nBw77Iete$E9h(yOadx%YB|&|oRqW+Ix;01! zLod%aqgy>#&x}?xBq;KTai19s2F4DwWlLVgOS!H;j8nR_m?1%>o6+B?Ki#v%U|=lT zSd2Rsk#VINIWr_E{)lm)84L!-on{QdOS!HC&8Rm+f~6p0eB`k(7#OFT(e)u2z0IgI zLxQCsV!VcW$BASxFg78heb#KP`u)tPHA8}>AYx21gTcUfbFfqY*^4!uHMUEz3PkpO zZD%mpZ%pi0+FoJ11Pef9{}fe@5il6+qZ9j~wwK#3LG4HO>9#W%?0h1fKMK*4>A2Rn zU4pVl3&)dWvt3HOF&o+NHz-oS%mxVx-3CcE+bP4|+_5QfpJZ3cE(z)zEu?;ua{v6B z2FYOP<%MP#xS<){_$T#PXwv))?uOU!67LxQCsVr0!=Ffcj?IraZ=p{6ruy9BF1 zWIw=m27~>U#QqE0>ur}{0f_8NvDzRN84UJQ5_^s9b+$`T`;mQy?F4eJ{s1#N&5&Rzh#0>#gTcUfe|@L^%NJ-mJ8YL= z6^QIrwlf&)cO>>Hwzt_X!2%H3zs72VGsgMFy%_#6g` zOm4AVg0e>o$CG6Bk4Wd!|FfP`{|hJ#IU6J>bQ>huY@1U5K2F90ZF@F648P}OnV}=C9A2CLm!C+uKWX8^TDc3d1j0!U(SPCM>*B%RlfpMN0 zZ>t3Rn^A6t1WQ50_y8T|QqYsJJsE|co2+j}#taFTf{1ag84QMge7LStfBt!z&NAC2 zSOp?`we1WB``wBCI@?pWORxY$_TIKL80T@utc zTE$+LxR(ruUS45_fm@r=YK8vXPK%#fhc%}6`- zA9Ft#42=8$G2VNdjGN5JnIS>(M~vglU@$QF06qG#2VTl`9cM*m3Xg(MWjGf77pSA4z&+rM)MrPESA;D4*F>W-2!NB;qzf=FM zb2Xhcwo9-IMD` zPO{@UOeiurW4i=pj}{D)Y<9Xbc|w^}zn2YVHb_wDHb}DB78z<_KKNt**zGfUF#%*SFm?1&&M~suqU@$P+O4W~j@KUbpBs1E~kYFi@7^|7V zU|@_hur}{0f_8t+0J0F|1`0mW_z9O64ZWVpM?bv znap6Y?=Cw&hk+s|Yi*aH?9qZjlFd$2CjY#|slTcXH8x04=r%~Q*_;fcU_k2IeIXn8 zajjhv)Hhm4+a%@IcP|+Xy}a5C19vf_+zbhdJYu|v<&Z_*lTl4uw&cwDOqa!H>s)2b zkf74dSU;uD_2;-B42FKJZbn*+>1LFfAwls+jML0uFfg83Oa0h#0sC>987VU)SPCM> zT4pd97!%Am>@_knW^{kZ9hxOr3eduwCS&>hEzJjmfw2!6?Xzyu<7-=*(P@SROF_iA z!wd!k;~#4}_5XU7rnAF#308r~euC``2K(<5`*hpeY?ojGi0u7sXE4|=PV8sd-fFu9 zwIA8%V78GE3S6g$h~AR^zv6`7&y#~Ix{3F@`&*YmSuFD!N53@wrt6ZuP|Lc z`Jv8Ltr-$jx)}qV`p>%`3{@09p zGbC6FBF1`VFc=t3W^Av1Y;8uJ84@f75hIVc2h0b9fpHKSg=_tCGiuF{U@3?gt!6M7 z7=6ItQ#7BPs_Cq;U4m60vY%l)gTX#~72AJjdxh;1EC7*xL)#e)_KAu8eA~-ym!S3| zd&i$N0SpHFfwJS@;ejG2Gqy`m_Gn?oBzRt5nLM?xQ$KA(nGF&Yx(yOMuP?)CFd+5q zzJd*0>)RzkeWQic-@qxirF+R>=w*u;1|DEWcNcePmY~SdDl*=FRU>6EFixN?Tk`rJ zm~vm8taH_Ah6I&v#vu3ORriCzz!+%8oR`UX*o+P{Bq;ufalRQ02F5FW)Q<^boNq>( z84@f75o1F$7z~V@8DA^)cQB*X3<;Kkh|$;c!C+vFB%|=7#9hs3F++l-AYwd-w@V}g zgMrZx45t3yYc-uY+a*{9BKwbQXE4~G?QQ$Rw%6M(!2%H3Gqy7r>{llC3vI8nU4q(= z?5|?BS^GWhhsvH$Q<0Okwo6dw7u%2iZ_#g9L?cg9Oj(%WxJ9NPW8} zvw>@UyCkS@w2-z*%5CFbG8hJVyBP*nnNe))gS%ReiEzD^$meoBU3bQ>giUSEcD zU_k2IeKi}n*0)Q7`bMjGn#$cv2176JGQ+@XGwRHcpvWUe*NYk{gMo1dZP}8^bD46h zoS=Txnjt}@o3W8o{~h;(!NAzejBS9=brI8#D)$VBq($nBzRt5 zhVx-S>f8M*HgK(Pmjv~VRDiGN(vYo+TpZ^cr=NzM%tgu~z1t7BTU^|1seqCa} z!uE38C8+(#{@!!yJcGf0ob35D6*-x)U4pVli-k_X^ZLr<-z;_NZ)HQ74H6W(4H7)B zFT;f}AocCOfel>i+a*DLqlMH@Qf^Q8lEKi+2h1?=7&E%x=MK#h6ggT&#*$|>QU(K~ zj<#&cAYJg39<6?Knjt}@o3XJ|f1&%qU|?)(#+S30a?g(vqr(gdia%mpVFrVN(e;J; zaXMbgbzNabn;8-;1rcKhGZ+kvUz@R~`f-36t!7BD6hw@Ho(~2Cqn3>JSs&iaKOT6n z87*c=uoOg$XYqE4WMD8bHUfjGzlQBO+a*{9BKxJbGZ^e|{oVHGzpt6Bw_Sn-AhPdb zJA=V~Q(~WNd!6kP)P7`Nhy@RwXE4}Ll0BcMA}4EYm!Ry?Vxd#;yuLE|x1T%px3{6j z1_=t?1__?mmti~%NPW9o*ub^ET@utcT1eX@<@R$g84SIA#0&#ZFr(ZI35q;od@)BO zWiT)1J%=)c?r+U@$OtG~?+f$mlpyj50GMDE^2s*$f5) z<8Pm-A7|sGT-RhXQf5f76hw?&%wRAueq%--^`pv+?)SJuvjj^4TA0&hEQ38C3HdBNU#({j6dM54#~h^U~CEoQ@_;q4%;PI1tR<9wlf&)fBw|=jw3Xa zZMI9W07Uj-wlf&)wY?q+y(SkvO=k=Az z_b+kk?`%WP1_=t?1__?mm*FxPkotDt&IYda?UJCr(JG#%5$+{}p_fmXVc^MT)R`ec zkw=VWPiv$M2F69S;ls^GnQ|M{s2{avNKol!WSsh+x*rS%#%^Xz664jw#i%hug5r-D zSDV3LU@ZAW{Wupd<+`pmqrwacmV$^e%nSwtAF!mjFE(TY|yVU)T0B+a*{9B739l42CZJwcGYr57SJh zY?ojGi0u2?&S0=lPwc<4z56eGtyzNFM+-?+IGK2VMCTa{_OoQib7fHEWT)*Cls#H7 zNbtNq{&ffZ`!63m_4lx$!v+Zo-3AGs*O%c67?ApQ-^B*5_3e_NzR@c7@?iIp!O+X6 z%`or`Gg{4%pvWUeFEbbnj7w?DmRviN>9TRP`q5&B1eI>arcV8(W-u5SdzrDn7;lUe zBWH#L#UC+#Wd?(R@x@2#$5_0S>-v=$^=3%06hw@D&0sJv?l336_G0aex^N2F6S?2C5&& zn9*v61WQ50*!FQ86Yuvu8Dq(4pVg&n{n2K$m?6PZ5Ha3*ObiA?Kgz*i>Thg&&UOh_ zfyjQ1?F>+AcxaqXmNm&+99bpIYeD-`|EB8zd-n8zgvMUxuq;K; zU{sp%@groscc2($W=K%{5u?TZU@$OxnQ;kT%5}AvkupPqr66KdnZaOSJl>_TY^8pj zU`F@5+@V>5r2s9=X)=}_ab<)P$zWjojEutP`lp!DX@&$#LB#ko4h9B;fw3zXO#RJm z@338hRUooYvz@_U|NCEUf3Ff~TnpGP!2%H3t8Hg6*dI>px7*%oy9Bi#*?Zg0V6a~( zdp=FY(OYbnpzP7YXeD@FUzyzgfm8n=8*(;CQ0O*D@VveZzk&g&Z}%f?puX*rpuW*6 zo~C1QHGzXN7<&1V83tZpMx7ZF6nVrbGlRjvxQe!HN!Pthx$Q=%AGKykQ0Zn2aq6$; z*JmVE^hp+ZXN+G_D0~mtX;i>_^(pV6Z=y*zdBvyOXaqOHli0A*l){ z({=`feZ1`XG!;khv|WO-M~j6{!SnjcW{Rc!v+Zo-3AGs*O%c27?ApQKfwm- z+b#*}8?9n5N8^eF2W2qy@-;IIyvU4JGbAYTh>#z}hT_^3jRoEZ`nf5f=U{a`RK)-+?1F1+qCquvY&mV$_Jq!|nb zM*E*NmaFdLSk5q`&I}2bf{3vfu8c6J3pQPL=$SWL_!O+VEW*B&>8Qt%2hh_sx^7A{{j~|)QYK8<$LByy)jbKh042(%+w9nd8x4Aztqs0sfmV$`!NvjwP z2F3_5nEE@}p0izoRUoq8X*+|#zDi>MbZ?+>EnvF@3qWK)!FC3Nea_oH`UAGt*)Bot zNA~`k8wV6PzFOUKQP0sIqrwac zmV$_JvKb5p#-Co-SSoH|KQ1t%+zbhpf{1YluJSOa3~y{?|cAocBjjSbYdT@utcT1fpQ<$j9uiGwm2df9D;f!CPPy?{G3OHkx!6&c%`!C+wA zLtD0F)ti`dhwh|)bebVSrJJ#}Q-2H37K4FtrWtKFkg=K>9cD;S{1M}6_k+Q}7-GgP zO53M*R6p9xkYFi@7(X(D!NBNzS!4P8_3X!`X0)0i!BP-0zK^Rs%qfF`aU&UpPk1Jp z(PD-KOF_g~#S8`m<5)15`Ul#cvt5E!AhJLDYjuIaVBa9IuV#C_?Gh{ik-g4#27~>L zmwfcO=SP zZIIx3eHk8s0jY2I0yc21ZqudM$iacWMYzBjY zaUX5jlJ;p#mzo{akBk`-RJs}4IQ6&oY%v%ZW6T)eLdIHVl$jwx@kflg?gxW`v6UHj z=)!C6_UcE<3<;Kkh*4(-gMsmWhsJV`-oY+6qx((n&@91HfEMO78OslFg@!q0FfeW* zqwopOm1cCBA;D4*G1f4H!N3>|22=k~+dFKRU=@h$ZMZl?7Z?oojT8G?wzt_X!2%H3 zFR-1#V1MUDAN_^xfX2VAVY>viAKACIoxxylkv*TL!ro%L1Z9sFx+}r+`pV>w|KQX= z*M^)85)`@(5)eH%W zJYwu)27`g|C~eu2^i-zYak#odKU&O?pwi9Q&Z)nPXN$qW_#ZPa`6U_anUOO?g5r-D zFS#EK2FA{2v?^_1!fhD(QE!F>OF_iA$P5Mp)S3t?ML=qY-cdo?~*;Arox`FU4pVl3*D9Ad3|N_ zKju31e`Z6O4H6W(4H7)BFT-3IkotBnVguLuc1cj*Xd(5Jl)DlausA4#p_hGmu%-lf zry1R^bBAULiX5#XV}COk42(9~@b6&#g6VS7mg+~R84^^w8QVMchk3Rb42*GRY^`^W z4bA83PMRy5-bG~<5Dvi42;j(HI|dbxWc*gu-%qc0c&G_Lh+m!S3|`!L%X4EB3v&!?%d*V-;Y*`tN- zO7OhCGWnm+IQ1{Fp~eOY3f%??p4XS*1sIUFc6YOZH+Z`wsBg58wn@tU0_PJ4Wia$| zbsnrK0p4Rqxfv1^dBix#3t-EKc<<{J)b)? zORyB6g*i>eGUivzU>m!eq`U* zb_Rp}0on6uD(o$`OHlS`p}P`1udht*lNcIo$k`x4q1zzIW>Ye}^c2UOxIbqDZ}4_W zP~T`3Ptz3S1P;nz=;hixSW^PL&x|@VBq;KTG13eM1LHZ`vL%N!GUd+7s2{avNKol! zlsokg@N6*{7(X|o{2DU0Fr&r{35q{ryz7227#RDT@uX7!-A&Yw3Ns{F3L-|M84L!- zKWAwyXY0bN*^F{CBv=X}#(C&4=9Iy}_&phg$JcH%BV&dHOF_iYKhcRxfS!z>g2B{3 z&Gs_eC0GR_`|Ic{y1-!Q!tRNE3)@q+ORxY$_N#1XFxZ#2`RE^P3^cCwU*l`d64X9g zNUFlgw4K3Ve?<0tnhJZT?GltdTIj9>&+9MahnxKp!z3FzY>=SPZIIx3eHmVRl4DNX zU$TK~eY+&6Z?uZNybc$zI4Famm+SCgO$qQJGg{4%pvWV}5oRzL7%$S6Ejj$>Oxtq? zs~;_9NKol!?3mJRZk1@_+TUTquvY& zmV$_Jl^F~MMjtbtpGd~7X4IJ>!BP-0{wK%zo62BdJoyA-g=_ukX4IM?!BP-0^iOo+ z5}+sJVlbHcXW3q3y9BF1WPckMn&<+9p$mH__N{HNuw8-$AhKU$JA=Xg_2WMJq75~P z<+e*u`;ooMb_Rp}3EA^$DsnPoy98yA7P>3J^ZLryYkj*UsBg58`bo;&bQSlX%3$c_Kpw0q0X}L*_n){!vjjzsR*`YE84L!-E3{=x zzTd!fIe(D)(P@SRm2SpPPW@`n7K4G2GvmI?$k@S*4l^Vu{)n;2{a`RKMw;;~UdnYX z+CcqiGed%gi zUSEcHAK{o2cVD|I?2@3q(L&lLDR&z#U~y0eLoWx@)|3F-2*eG&84?tE#5mpz1_R@D z+Oj1VUCMO1aD6c{W=K%!X6)?LKhm?sU|>u&W7H%vb}^&O3<-)qVsyJ73D}w8A1N~=SPCM>G&2|sj1n_0xq|(;$Bgb*xkIxAO95J#(_}0YuEb$*t*>2@YyD?u zB3AgvLw{#Rrx_9~1rbC4L?1qQ<~4ovL3*xqKl z1Pef9Z?>JmU|%h%!zw-yXx(dpuW*6o~G%@D;$);(96weYf6C65QrOkGbAYT zh*4_>gMsljZP}7FbdfNAT`_9Skf74d*u|-TjAx6%z_{Lw+mv#{%&0L#g5r-DpSvFn z2FB55ysYy1d>!?p!VC$Pf{4*<27`eyz>Mc~;dP%G{4o-Id^ZeP!}S ziQz^YI&6@j&~1?5d3_lc-On*6?zQb|vrB^dMyuG%8OTi>l)=!;@6pzj0G}feH}qym zP~;KgbTb$XjK9!^e=U7HSIW!!i_v0+1eI>au1@_EJX;I~#?5BjaS<8&nvpX@g5r-D zU%DR*2FCGb^wCArmu2cly%`cL1rg&`GZ+kv_04$aXY9vAX4IJ>!BP-0nl8s-aIN2y z@$&BwD_rY8W=5?U5-bG~L;pl4E&+NnegOtkf4uEAwo9-IMD``P&_ov)48y2N?EBhY zVY>tiKxDtub_RpJe`3#163L-`hHwBne z1_NXMy@(b5@lgJw#hB9;GbC6FBF1nt7z~UlU@-MBvpr|K1gk(~|Ds-9U@+K6CH4bs zueV)-1t7BDV>^SvzFuPQYkQsT64ZWVKiPH$gZ;JNI`yd-d#&vfls#IQF$td6S0-NU#({jC(v51_NVLGxk4^jJ7p2mhP9hL$d@+0a}>TWGquJ z!7jMg*DlGm{#*AT);{YY{R@h7%;+>jf~6p0=%47sB|uNcbzm^{uduztb_rI2$i58k zw&((bVHn3H_A1-kY?ojGi0t>-&S0=_nAlghz14OJYCp1{VLOAtzTj@BK9yx}y98yA z7P>3J^ZLre@`BYRt-i`8cj1Q@sF=S5+4}76%r;H z6p|4V#k465%n6JcW+wAUh(;yBRJcf`ty*f0cM=T@c~YBdY-5X>HpHg4)TWl&+!k%M z#H30!-lF%m+~3-3pFK0^NT9vn_rLvezVER1`t7yWUVEK=*4byDbB5E{MlRcs1gRjw z#<#gI1Zo?XvW*YCro#1L@}QFSX~mF_TbQ-BzRw6j^u2m(8&rt ztRM+abXJfA@9Rs2p>}dwXVH*T$8-F3S=Ayx&v<+Mu7)raFOZ`J>E4kD+lvc&1-l4Q*Tv}?H3Ih1wvly*LsT8Glk1L)}GT&#WR5Cl6aY%;ViN-e;cF8=aq)TQ#GpUlfo|zQMT*-`(%rs^O{{}8?GBX2`naE7PWCpiLNBbmm zl$l=1{Fa#>$^4v|ZprLprb{x9zhtgqrdl$Y%#=&!0%l4ja~3lW$yk`NOXm1fG6uPl`6Dwql6j4p zY{~qTnM}!aF=Lg?cbG|+OdB()lKCbxDUx}R86lbbm>HDMode7aNXB5MUoxf4^hxFx zW_l&_d1iVfGn1Ka$z(9oC7JV?*&~_Lndy|wU!RoSj!I^LnTTY5&rGXiUS+05G6$Gx zkj!3Y>Ljy^nOezgX2vg>N0_OWOan9Jl96wToA0E1xyB)DoXprIvzVD&$t++dM>2Dm z$(GDj%w$T&%8XSq=P;8lnNyfamCT2mWeidzbBq}wncp!pD3j`hAsII%xWhMH`H_WFGeZOI>R$u{#q!Tzo@$>#nS zEt0%I1pf~{sr90Pp|mMn+k-ZH&}J`aZ+_E~R9oXYvF6ia&HK6e+ith7J7z@#_U3y! z3!A5w7(x^@&&G^nA{u?)z6E#-@JfwOV7#!os1uG7H{8)(PXKGQ)xm{2xY7XH0NMwY zmLw``rOF1nnqY1q2_}ooNrGr1*1c$rnj0sHic2Z$()9+ql3bO&-Cg0HhJD&!bdCvEPY4(mWEE z7z`IlFpngzdFqmx6!1rgD{4FkKt6zi`WwWg`?sD(Efir8PkomL2BwmbnmcsvRBC(D zv$l7hLyyEa|MKDTM%z1(It{8tO&8=(dUi>_&HYO5y%Pn>MQ_{O2mI#?VY9qU>@nMl z11?p52h`U1_4{E4Rm=MF8#f+ezn;%#Tkb~)s6oK<*1i4m;OU}i0|8e(eX0=8P2lG9 zEH9GbH}dL_Sc)nS%x@}e-mV&W3pB6(Z3`M{fG}oXf{uIh)al%#?zkQ`ZrBGdziFc+ zjT?KIlcaIue}f}E<=2Ww$k2_jOGn4voOUh_!psBx8xXYqO8kNBhgq@@dHd&G2pK~D zXd?A>0O-sC_=*Y6OTUznP9TNNY0t@u$wke@DST|AEQzX(D`C3v>%A}m3oRO7$#@^* zwHmiG-p_cs#Q(jsLVp_P3kyE{#9L@;!`4HNKkh0mfT29%MU%jMr+ska2-@h;du4@k<#` z0gj|e=?wFL6p;QaKaxqalTQ1ljgsJwmXm;0>0`MXKv82!HI%J!%I|(Cl_RVy2^Otb z`S(|FFtB{Z%C%5Fmz61@^*_zZEvzgF7O+^kO)C?Cma~T`hgrFfl_kMC7Aw28G67J| zj48`L=diM$l_kNPgq5$<$^<~!4`rkW3fL5K`@aTdpUYWU5==^1`Q6LSJ`(_CdrUcl zlz+g=4px=~lN45dS}PNPv$_)+IEs%c91{{&BT`U*cQvkEZ_)s0k`+K|f&3x?$p-R5 z0+I=2R|4V(vM~X%18Is!zE7?ITt#5>)D-q8Z~u5Yo&f%Gnd~>N024&2MBZcsDeyBg z^_G-FtcO8(S`K4Y7(-D1lP^^C;&@@>#%%P>y!7H)p2bDv=i5wEfTh2w1Qtw*JI1Gqbi^R6E(%t?X; zFK{qN1m^a#Mqj*!?(8?XvM0W>aYHw0Hf@vylMuFZkTMmnC0MiufFYHFn5j-W*vy_p zGo5TklE#f)Y^GV8AwZc?ojZ-qL=w%kvKdKmQnQ(b+6)0?X8(qKM2v>V{tXWUqf;fn zNfw(PA)v5wBS~N+SfIS9>B83vnr=x!gf>4SW7+^2Gp0k*SpD;)Fn;aMqo2hDkj~hn zsb=9+RcNvu#srn}F-=fgcpU)dX9HQQ`>)!x`+f8mW(`)NIm zFQQ;Y!A|NwD}3OEg2n^(7mh%GTOmp?@cveND`o8fmQC;cL_mzmO9e=bh)>qPx^_S=4_WGi95uTrqX?cF*vg4q~ zZKwTaZD##F$0y!*abYE9>&M z*;CWixL-1eR7h?fwUSxfj`!oT^ogxQ?H;HVLJSkh*nyQXT6e5Ys;MxUxA)CyH?AOu ziqiX?$nWSAI^9p9t1mh(LIs%eo^J}|N`XH>;13WOs4p52p-c$8f2XoAO$ziu08{*B z$3TzC@ck7xEMVw}9In6Tuw~MN_kLCglqODkD1;K~4(HV$#jO4}I6BWQ7h>K}#%l;F z-Bk^uTaV)EA+vytMmnN#k#chu|d^ z7?$-)*Z&1_a3D=7H>XX6298J%AEv)v!j3hk{mCtHx{-(cn;QQm<9&?(RO8Px-p}}V zHNJ)M0me6KyqR%Yl!5&X8V@s0cNxHI7|;80J^p*KxDUt|G=l%iE0F??;J-CXF1z$J2yM$S}-MI-nx?^5R7i)A#>`XoY*w4N`6{A#rtA4A|9Z3hwDg9vgD zXb*xMo&q;_vnmdM1kx(NN7GUBn66nLs{y`V1jRE7S%i_Feb=*PjsK8;Z53Hf)8Q0Of{_G*MePls@;Mkw?o z)2a~)J;{7oBNTd)32B5vZ#gpcJP5{4pe0<5)IjMZ6C%*bRbX>k2KQtStmCA(9GuJu zi(upvj@LENouGA(<8+11F81u0RYn2F@i61vz>hqXCL8>juKad_(z6=>Mx4@8R(ecV zuH?#G_I|y_iy5~wey7G~G45czSmUXTmoh#_jk@jAv|)VPE32FAB)d=}#^jDJnz=P};O_+1+R&?WtkFkYeY!;D86 zzggqYGv3MgERBDQ@jZ-B(fHkrcQHOr<4YOuX8i5XnQ^$D@gBlsL-A!*&qFZ>l^AlJ z5S2sj@HsT(_8~RwLyG7I?WPn_h9UP=R?OR9|DfOv;93pJ=?idxRoaE!l`{dnf3Nl`yUvPC|2L6Epan*Cm$WjJ|^)`(^ z%6J#!U(vXiahY9bQcKNKSL(`0xippwUxDAD@%4<`8PC(Wi*X0zmuY++k`0E;f6N43r)zA0=jsJ-8TE@Sv@hyzkG5$@BH!$A7_&SZ>$#@In zZjCQwyp{1=H9nK^2;2|61c; zWxShl`q%C(?`OP+@Ysmir0Qdd5(o1^1{0IuY7Mps-XK|_&F^CKDhb}K@oz95fy$9T zohh@_ZL=;k z>$ohn5_)m{?nscZA}B`|B}zPyAd$`C3@ETUEe)N(aO*+85Qn82|LBVn?`3?h#xpn$ zeT-kO@v|6LS!%Mz-(&j&EdLQEwR{-gW_*zGH#FYGxL`kDW;`xS?a@dIL?}zOX~di- zKCBUQp158k<~)&}zQ7Y$Hs^_7(ukR*3N^y>#1B4`5c{jRIx$P7E3i3DMQ{)IP-y%Q z9N8|$`!Q|gh;%b9vs4qc)I9Yky0Q|31BWS>!}1-C=QD0+{978ogmDMs^%@_?xXMy1 zH2%A#(tbJ1muUQlj8`*$gU0`jaX;giXuOVbm8CwT@p8uNSpI#?R=Lx28E;_xHyT&{ zZDIT+ji1HxDobtC_@D2P{zq8;>l#1Ac$D$GHU19U>14c8;|Cd6S?U&z-_LfsSbny~ zn^?Y^@e4G53(NNq9?Meyg*8y>cS4qWi^0S!^>YozWvRVNy*ceP6$TzdPifo?1MnjY zPLo;c`=u&N?F8Rwle0u*DS2Z#JWK7Xzh@t^)R~Y(30Z19O2|_CVp*#9wvT10>8km- zEcKI=gnrZ7Z>--|1m#FaqQvW%%ElzJIh+R-*qrv>?eg%dEVWzXzhb}Lk!DofFcJ6R(s5Ftcp#LQCfV7`ei zL(7~?9oC4MrC!#EIhWd_5i?7*X@s-XqEi!M{{UAfW~qPzo6{C!cp$OtVK0|zoc=>L z@Gi#ZYW!lxWtMvL6xp@uy7DFrj#!pDQ{z^~?Hn$eM$@rwP7{nf82=BA|8|MAqq5YG zH2wp|%UOPh#-Cukn(;?8ejnq0#_!R18RIHT85;jQ<8>@wq;V_b4UAu-aSP)uj8E10 z8;hm?Doah&_)CmOSpH9`W*oLK9%cNO8viomos2)PaTnt%OKs73KI2_1->h*P``OKS zSmU2z`5wY!S!#)@Psmca3?^o&D>T@GaY7qL`guydIjxeLZDss#lg#ek!gvJukxjUh z!yVt*@8SNQrSV4@?`8ZyFA+ZZ2U`JEbng7HDdAJuq>alwAx%XnOt@@OOlB6RGQYQ)S^c8!=> z>ROGMSt>&#W|o?)5i?5(jc}Ig7?%+H*D-00<^N3f_@Dxt({^Hb#In?Gjkhx1#oj-z z@wJT0EOisL)I9Y*T{+~C#&WrlDviI+xSesDM)P5Mk#Ps(*J(V$xXMz~G`^1Ua+W_; z<4(q_8GrW-)Baq>{fz%g;~9*rEcHW;Phh-`<=Zv>)*{)j2F4%KcsJuMjIY-CHpW$! zD${r~;}Mo$pm87LQO2`0?qIx=@$)r)72_&P4Pi>l@jR39E|&kJ#*Z$P{&zF}GmRf$ zyoc~umfEH26SCA}3?^o&4H|60IH3*0!LQ9PXY(pcm1=wohdBay`FLpEzB&}mch681 zc>goCu<_iC!p1LVq!i7>^O(gM=>?a)STt$wFO0&b3#X?-sOb9Fg6A#ITla<-C|Gxt z9(m{ecs!wDJPCtt6Lj$+XlmY4`(^v`hTi!2^T}NE{ki?NW7I5yfM)Tmhl4S2)QZwz zUfm5MbUI1{R0@wfI>laZ8?@lhfQ9O9|4Q@D0o$A@p`Mtbc=!^)h%wn8v`u=3#V0*8 z*>?Ref_RpC43D>mnrz45Y;NLve&p22?PLW0(GjqmJ8ojiEDO#<-Yn;hpLn_@RS+WU zUWwI_ge4X7UyEUuGr$Kd4oQnauSuDBi6uvzJ^9RWrwyFiH_;!HmHHP>oS15vAxA36re&p{u%+wO@#Hx2aX0ITN_&Br ze0<%>^5`RAO%*BeTDgkYY?iwAQ1slo%p z>MtFkK2_CwmEPsU^l^VFw+j^BQ<>QQ4I|XA9ie`fs!!HGGg#fK-SxxT|D4iasql

!H8(mPy{*uMuysQ=0c^(`aRe|3cVua8iVnV{SbDZg8W`5zuZf8LxdN5EGV zaMoDoyFErQtH>GhI4iP>XJyZvU6B>?`8*ZA+AO2SljU(&WYusjvDnL8?o~5tger-# z*>X8U&gEfm)+%QJw(xh1c1)k#L1Sgu=_#x91q|g}JTzA89&|lR8G>0=8ml!GzZ$o< z%oixDan>qNW55_Ko6%&G_&Y{BI-1^4*>YdNcsQg%%92XFIVSJvid7VsHe%k6?-y$@{eSa)sugpR!{lw)9R<89yjMn?O%*K zk{RAU*`$qbNN6+VwyXFJmuG1qEYmz9~06e#7HUa#K#w9?2!M=TPSXM6eCq#8e1>wKHgxQKEE*)|5E| z0p}Wpx!lW_BYAp)@qECyI!Nd9lAsX?X8BinoR!tSU{<~nToLm5v*rfePOtT5WbqQm zJS&~NSz$j?k>RoufV>fO`@C5}zflTgl?ol0-TF#W<2>WfCGzjcIq*9i6Ie4$YJ^=D6F|IGS7DSTU!`gVmI!|Erg z_LIfu<`MkEa;n^t#gB5?XR#-+v8rnjq6qB+NFKF$mGDQXJLMuEoXr`3~evU;*l+wl6%d>+qq?$TUi9?b^l zsN0V@pEL}}3ZLtA-IWu@8RmQu1(#7<;0;}CW3%QB_1JDv{t=U5=`H6RD!#N*#aB^P zEzJAvWcBAN!5XDEZJ2%<>hUa8B&#R;pH@%DmX80Y>7zPnJ+)6;GJQI})J9R>;=;VS zb6oBkXH~#3vKR1i`umPg5Z892){l9B$5)8~kKa9$()A!RPQWUm6ZBOn|Zz5 z9pdBZSCnfx)>QT=Jabt6awV9ozIufEvsC>h%5Gnf{ZO#H=1@6VJ&gxC-pSghu|n6a zWc76JeVRUv*H6mU4T^l-7Td>z`57R#aHLcupYfUF@zg z<}CL4ik#jxi^3s)IFuKt3fCCkP*5Za-5&6Hs}_c9DvZFQD**@mGgYwAg&uB)jzYl9GRdnd9&JCV<7e5ei^I(1>WU8Ve^IpYaBjIY=WZL2>Lv$ z3>VI)O1fzdire_cT703Pu++9JgI8w*K3^y!>B_ClH0Fho(jX_DVaRAKCOt7n>y>gv z=D!oyHFK5Bc^5e=f;8tvx5BilIEH>Cu4#XvAX#8fWZBZ|q}!ljb5nq$a4WHUtFVBA z1!MFT3oiJns1pxUYz~W=ahcx}rrSI^>T;?qfiYOahzq!D%e;P!SJPkuHU^vQPzZU% zv>1ns*&bHbhjm#ZW3Z?CwKJ`;rV2nnoY6gxfQU87)cElb%udVWzfW zwm#wZZxsH0!Y$uuwd0OIuc{>iqv@gaCbl&9cX+}pb8()~BHd{1v5`?~s6NQw!zl}Z-}!&B_8Su97%ctJ&Mt#4;M{RDm~61O#$MjDH0lu z!Qr_c&#!&-oJB)yvmdJ+9O@~X2!g=#*5WjHlsSb(&_8NM(gzgOO?pV3dCzi)xvj{Yx z?sLQ!Y|F4%@XxUH?(*RZi`&4(Yy5hJz)LrKN!R2t#v(^H_J!D&cvpCRtG!lU>Y-s|_4-2A<-V}jC2xqsfl9;b zGMp|f`K}NPeb%7g3B+S84_SR))U0yjJ_!pF)*7rxSbcuOYo)Ok6a`*+>RO#@eL~C& z_<}*DVfDDZD}-#%=?TF68mm$34&pA{VFYU2L3*AIQLh{KmM;q`^9s;xCfjW4i9 z%yqhCBUZ0m-VnB0qY^#`aF?KZVlIlHonXj{=NwiKRzsx!;mT@O_64jicR=3pVTBTR z|6z|3h0BI0a)zq4^U5sx&ykvkIMVk0i+DZcKNS7%4CaRv|HJ9bv$PO>naq#7jOo3X zGtE$Zzv4~#FJHy=QAIPdn15OEUs1d%f8{K$U#RHIvzhlR{_3llH|1Zsmigm~`g52+ z?>eThDc+PXDOoJ?-IZ>mrXpYLLo++?2jj01y8En;VQ3B;OPi)PCr~Sz>!bO zjf#qZlVl99muL`fp(CpdIO%dj6$?n9kS4ddgjX@#!UM-Zkc=>JE1jZnhA3EEh&>M< zXXYDs!4J%?7U7i+q7)bz7xc@G-3Z(zQG|>JQUVz%TO{&rMPh-ixF~Po9l~X(z3%v~@G$_OMJDpw}cBxt|R~3ZLjo|6mhiO2vu%qed{K7|9D)g@Yjq+!_iT zouX?*9vwRx8i)lAC&Uf4GJb0iz)(2gg4BGhC4elh4hMo%xZQAJ5J@P^6yYPb1m1#N zg*ODQ6f^LwV!>R8$ScZ6s=P%kDajM&na^j0MLl4lpLL|oJWouz`$0aBOgbDl`DApQ zdeMP?0+F`ia_2U)eUqB?!+8@QPEGst+laJ@T(*}xpQ*`fY8p^-CYy};@eEt6U4)+!LmAI+#UbqFDC2plnC-xOqA25;srU}|T9onJ zR7|~9i2G5-vr{pIo&_kJQ z(H3uy_;{Pd$J-@7nQf}0Exrv#PX0X^fs+w98G(}#I2nPH5jYuvlMy%>fs+yV|0x1i z<&31MU_Q_H2EzO}O4Xkv|AY}J$Ne`Q?nr7RJ(3g2jg&^pBejvbNNXe#*%RrC^hWw3 zgAuVQeUo)l&ZgW=rJKq()o!ZW)Ve9MY0svvO}(4?HVtkPo6|R2H|K25-CVl4d~@yQ zy3LL)ekGy*<0;f9b2on`nNW0ZP^;#+PSrRYtPpHtpi(Awxw>% z+?Kt~zRj_%dYgY+%eK~Sk!@Yuy0`Ui>)R%_r)*E%p1nP1d+v7o_UO*eo!vWocJ}Wa z*qO2`byxZ>-!8<-ZxK^mGgN_SN6 zXxY)aBeDYqPj3Ho5vZ7d5e*fQ*<%ru3QkeG1AI01I_$J1d6r&QW8mi^jM~-LUXAsn z$yf`+8X>J2()u8+UD4XtgRt=w_I=p+`;{GZMB_z zY^hg!pR+YCLeB??M?pJiJ zqOFQPrRet*eO1wZMJIfh?OGL`qo~=&|F-?&VeJj)zohC{Df+mgPb=D^=-Y~pdxraS zp`x=CU8?9^iaxC9Zbc6&dR)=T&$69aiY`?2i;AvQ^j{TyLecLjdPvbX6g{r!dGwQx z*ybqeP_#nPfT9hGZdUZiioU9-d3;XBhvS(nKJTb_jl(%a+j)v!si+yp6P25*^okYr zD7seBM-|cg@BhWb zMCXK>+7c!Av?}9$J$#LTFA}7sUuf-;P#|0x!s|!$=?66;_Q%)Ztrt8m(uis+PJFwr zEA&<2Jrqw|-Lu^Ge7E6o&0PasykWs@Y|&cz@(a5X`I01_=Fs&bFJz-jQz~Cm=HUZi z{=k^N*_u@$?qwW-n2Q*#m`A@~G~eT_%5&kpkYG?rEl}DrjK%nhYH=Z+SY4}7kAb+- zhsE>o?pa8wh$^ndXIe&}I8+jHR<2N0L9Uta#``ez>Z7UA!FZwJTxBF`Ea4jbs7ff{ zTVp!V%&6UP`il+qF~4j_wux(}GQ38$=L;0#$C_@FuZ_u$s`^b_UjR?5j6i&iOSQWO zZzQ@a=hCZn!T9)BNpqc{M{%3EAp60XR@ zTXj|563ids>mTF#0!%LS&?eQ{PMSH%IRq<;70lV_E7RNYGcZN`D|!Vk)sjU?d)!V; zGfJ3K(UzwI#*(xro4?I`SypLJiK|nu(W;tpadsWzY37RLN9dT{2=M?T6el@=oDBbIQ4gVi68rV{!+%TRcV!X8s|GNUx&eU-AcdlP(qkfn>#Tce&(=Dsq@t z3CQqxIGkes_?mG1*;;~Ewq_5*=HaJCvWH=L^r^il7>1V6hta}k!D(FR!`V+~l@N!C zr{YqdSAH`|Z&@O1;aM=%l!RP%r`JVanPq_PZ#=4w{Mf$DI`?UB$ew&I1h!r8O9 G{(k`bvRI!0 literal 0 HcmV?d00001 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