Issue 139: Stereo PIV implementation

+piv/piv_FFTmulti.m

@@ -52,7 +52,11 @@
 max_repetitions=6; %maximum amount of repetitions of the last pass
 repetition=0;
 %repeat_last_pass=0; %set in GUI: enable repetition of last pass
-%delta_diff_min=0.025;  %set in GUI: the quality increase from one pass to the other should at least be this good. This is sort of the slope of the "quality"
+try 
+    delta_diff_min;
+catch
+    delta_diff_min=0.025;  %set in GUI: the quality increase from one pass to the other should at least be this good. This is sort of the slope of the "quality"
+end
 delta_diff=1; %initialize with bad value
 for multipass = 1:passes
 	%this while loop will run at least once. when repeat_last_pass is 0, then the while loop will break after the first execution.

+piv/piv_analysis.m

@@ -1,5 +1,5 @@
 function [x, y, u, v, typevec,corr_map] = piv_analysis(dir, filename1, filename2,...
-    preprocess_setting,piv_setting,nr_of_cores,graph)
+    preprocess_setting, piv_setting, nr_of_cores, graph, tform, output_view)
     % wrapper function to do PIV preprocess and PIV fft for a pair of image
 
     % INPUT
@@ -8,59 +8,43 @@
     % filename2: the second image 
     % preprocess_setting: cell of dimension 10 x 2
     % piv_setting: cell of dimension 13 x 2
-    % nr_of_cores: number of cores specified by user 
+    % nr_of_cores: number of cores specified by user
     % graph: bool, whether to display graphical output ( not available
     % for parallel worker)
+    % tform: tform2d for unwarping images
+    % output_view: imref2d for the output view for the transformed images
 
     % OUTPUT
     % x, y: coordinates of vectors
     % u, v: resulted components of vector field
     % typevec: type vector
     % corr_map: corellation map
+    arguments
+        dir (1,1) string 
+        filename1 (1,1) string
+        filename2 (1,1) string
+        preprocess_setting cell
+        piv_setting cell
+        nr_of_cores 
+        graph 
+        tform 
+        output_view
+    end
 
-    image1=imread(fullfile(dir, filename1)); % read images
-    image2=imread(fullfile(dir, filename2));
-    image1 = preproc.PIVlab_preproc (image1,...
-        preprocess_setting{1,2},...
-        preprocess_setting{2,2},...
-        preprocess_setting{3,2},...
-        preprocess_setting{4,2},...
-        preprocess_setting{5,2},...
-        preprocess_setting{6,2},...
-        preprocess_setting{7,2},...
-        preprocess_setting{8,2},...
-        preprocess_setting{9,2},...
-        preprocess_setting{10,2}); %preprocess images
-    image2 = preproc.PIVlab_preproc (image2,...
-        preprocess_setting{1,2},...
-        preprocess_setting{2,2},...
-        preprocess_setting{3,2},...
-        preprocess_setting{4,2},...
-        preprocess_setting{5,2},...
-        preprocess_setting{6,2},...
-        preprocess_setting{7,2},...
-        preprocess_setting{8,2},...
-        preprocess_setting{9,2},...
-        preprocess_setting{10,2});
-    [x, y, u, v, typevec,corr_map,~] = piv.piv_FFTmulti(image1,image2,...
-        piv_setting{1,2},...
-        piv_setting{2,2},...
-        piv_setting{3,2},...
-        piv_setting{4,2},...
-        piv_setting{5,2},...
-        piv_setting{6,2},...
-        piv_setting{7,2},...
-        piv_setting{8,2},...
-        piv_setting{9,2},...
-        piv_setting{10,2},...
-        piv_setting{11,2},...
-        piv_setting{12,2},...
-        piv_setting{13,2},0,...
-				piv_setting{14,2},...
-				piv_setting{15,2}); %actual PIV analysis
+    image1 = imread(fullfile(dir, filename1)); % read images
+    image2 = imread(fullfile(dir, filename2));
+    % If a tform argument exists then apply it to the images
+    if nargin == 9
+        image1 = imwarp(image1, tform, 'cubic', OutputView=output_view);
+        image2 = imwarp(image2, tform, 'cubic', OutputView=output_view);
+    end
+    image1 = preproc.PIVlab_preproc(image1, preprocess_setting{1:10,2});
+    image2 = preproc.PIVlab_preproc(image2, preprocess_setting{1:10,2});
+    [x, y, u, v, typevec,corr_map,~] = piv.piv_FFTmulti(...
+        image1, image2, piv_setting{1:15,2}...
+    ); %actual PIV analysis
 
     if graph && nr_of_cores == 1 % won't run in parallel mode
-
         imagesc(double(image1)+double(image2));colormap('gray');
         hold on
         quiver(x,y,u,v,'g','AutoScaleFactor', 1.5);
@@ -69,8 +53,7 @@
         title(['Raw result ' filename1],'interpreter','none')
         set(gca,'xtick',[],'ytick',[])
         drawnow;
-
     end
 
 
-end
+end

.gitignore

@@ -1,2 +1,33 @@
-
+# Autosave files
 *.asv
+*.m~
+*.autosave
+*.slx.r*
+*.mdl.r*
+
+# Derived content-obscured files
+*.p
+
+# Compiled MEX files
+*.mex*
+
+# Packaged app and toolbox files
+#*.mlappinstall
+#*.mltbx
+
+# Deployable archives
+*.ctf
+
+# Generated helpsearch folders
+helpsearch*/
+
+# Code generation folders
+slprj/
+sccprj/
+codegen/
+
+# Cache files
+*.slxc
+
+# Cloud based storage dotfile
+.MATLABDriveTag

Example_scripts/+stereopivfunctions/compute_output_view.m

@@ -0,0 +1,19 @@
+function [output_view] = compute_output_view(image_size, tform1, tform2)
+%COMPUTE_OUTPUT_VIEW Summary of this function goes here
+%   Detailed explanation goes here
+    [xl1, yl1] = outputLimits(tform1, [1 image_size(2)], [1 image_size(1)]);
+    [xl2, yl2] = outputLimits(tform2, [1 image_size(2)], [1 image_size(1)]);
+    xLimits = [xl1; xl2];
+    yLimits = [yl1; yl2];
+    xMin = min([1; xLimits(:,1)]);
+    xMax = max([image_size(2); xLimits(:,2)]);
+    yMin = min([1; yLimits(:,1)]);
+    yMax = max([image_size(1); yLimits(:,2)]);
+    width  = round(xMax - xMin);
+    height = round(yMax - yMin);
+    output_view = imref2d(...
+        [height width], ...
+        [min(xLimits(:,1)), max(xLimits(:,2))], ...
+        [min(yLimits(:,1)), max(yLimits(:,2))]...
+    );
+end

Example_scripts/+stereopivfunctions/get_tan_angles.m

@@ -0,0 +1,25 @@
+function [alpha,beta] = get_tan_angles(coords, camera_position, lightplane_normal, e_X, e_Y)
+%get_tang_angles get the alpha and beta angles
+% alpha and beta are the resolved angles between the vector connecting 
+% the camera to the point in the laser light field where a velocity 
+% vector has been calculated and the normal of the lightplane.
+% Arguments
+%   coords: the 3-D coordinates of the sampling points 
+%   camera_position: the 3-D coordinates of the camera
+%   lightplane_normal: the normal vector for the laser light plane
+
+    % calculate the normalised vectors connecting the camera to the
+    % sampling points
+    camera_to_coords = coords - camera_position;
+    camera_to_coords = camera_to_coords / vecnorm(camera_to_coords);
+
+    % Resolve the vectors to vertical and horizontal components in the
+    % plane of the light field using the rows and columns of the sampling
+    % points as the Y and X directions in the plane
+
+    % Use the dot and cross products to calculate the angles to the plane
+    % tan = sin/con = a X b / a . b
+
+    alpha = abs(coords(:,1)) ./ distance;
+    beta = abs(coords(:,2)) ./ distance;
+end

Example_scripts/+stereopivfunctions/image_feature_extraction.m

@@ -0,0 +1,36 @@
+function [img_p, features, points_obj] = image_feature_extraction(img)
+%IMAGE_FEATURE_EXTRACTION preprocess an image and detect MSER features
+%   This function has been designed to aid feature detection on images of
+%   vapour. 
+% Arguments:
+%   img: An image object process
+% Ouputs:
+%   img_p: The processed image used in the feature detection
+%   features: The detected features
+%   points_obj: THe points associated with the features
+
+    test = false;
+
+    % Try imhistmatch
+    img_p = imgaussfilt(imflatfield(img,300), 21, FilterSize=51);
+    img_p = imadjust(img_p);
+    regions = detectMSERFeatures(...
+        img_p, ...
+        ThresholdDelta=1, ...
+        RegionAreaRange=[500, 1e6] ...
+    );
+    [features, points_obj] = extractFeatures(img_p, regions);
+    if test
+        show_regions(img_p, regions) % Show the detected regions
+    end
+end
+
+
+function show_regions(img, regions)
+    figure; 
+    imshow(img); 
+    hold on;
+    plot(regions,'showPixelList',true,'showEllipses',false);
+    hold off
+    drawnow
+end

Example_scripts/+stereopivfunctions/image_intersection.m

@@ -0,0 +1,19 @@
+function [img1m, img2m, piv_mask] = image_intersection(img1, img2)
+%IMAGE_INTERSECTION The mask for a pair of images in a common view where
+%both images contain non-zero data
+%   Detailed explanation goes here    
+    mask = make_intersction_mask(img1, img2);
+
+    img1m = zeros(size(img1), 'uint8');
+    img2m = zeros(size(img2), 'uint8');
+    img1m(mask) = img1(mask);
+    img2m(mask) = img2(mask);
+    piv_mask = not(mask);
+end
+
+
+function [img_mask] = make_intersction_mask(img1, img2)
+    img_a = img1 > 0;
+    img_b = img2 > 0;
+    img_mask = img_a & img_b;
+end

Example_scripts/+stereopivfunctions/recover_3d_velocities.m

@@ -0,0 +1,49 @@
+function [u,v,w] = recover_3d_velocities(u1, v1, a1, b1, u2 ,v2 ,a2 , b2)
+%RECOVER_3D_VELOCITIES Recover 3-D velocities for a pair of 2-D velocities
+%and the tangent of the angles between the 2-D plane and the respective
+%cameras
+%   Arguments
+%       u1, v1, u2, v2: The horizontal and vertical components of the
+%       velocity for cameras 1 and 2
+%       a1, b1, a2, b2: The angles between a horizontal and vertical plane
+%       that are perpendicular to both u and v
+%   Returns:
+%       u, v, w: The 3-D velocities
+
+    % Calculate the new 3-d velocities using 
+    % u = (u1 * tan(alpha2) + u2 * tan(alpha1)) / (tan(alpha1) + tan(alpha2))
+    % v = (v1 * tan(beta2) + v2 * tan(beta2)) / (tan(beta1) + tan(beta2))
+    % w = (u1 - u2) / (tan(alpha1) + tan(alpha2))
+    % or
+    % w = (v1 + v2) / (tan(beta1) + tan(beta2))
+    % if one of the denominators approaches zero (eg )  then use
+    % v = 0.5 * (v1+v2) + 0.5 * w * (tan(beta1) - tan(beta2))
+    % to avoid the singularity
+    % See equations 7.3 -> 7.8 (pages 213, 214) in Particle Image
+    % Velocimetry: Apractical Guide (2007), M. Raffel, C. E. Willert, S. T.
+    % Wereley, J. Kompenhans, Springer
+    da = a1 + a2;
+    db = b1 + b2;
+    alpha_mask = abs(atan(a1)) < 0.002 | abs(atan(a2)) < 0.002;
+    beta_mask = abs(atan(b1)) < 0.002 | abs(atan(b2)) < 0.002;
+
+    u = (u1 .* a2 + u2 .* a1) ./ da;
+    v = (v1 .* b2 + v2 .* b1) ./ db;
+
+    w = (u1 - u2) ./ da;
+
+    if any(alpha_mask,'all') || any(beta_mask, 'all')
+        u_mean = 0.5 * (u1 + u2);
+        u_diff = 0.5 * (u1 - u2);
+        v_mean = 0.5 * (v1 + v2);
+        v_diff = 0.5 * (v1 - v2);
+        u_coef = (b1 - b2) ./ da;
+        v_coef = (a1 - a2) ./ db;
+
+        u(alpha_mask) = u_mean(alpha_mask) + v_diff(alpha_mask) .* v_coef(alpha_mask);
+        v(beta_mask) = v_mean(beta_mask) + u_diff(beta_mask) .* u_coef(beta_mask);
+
+        w(alpha_mask) = (u(alpha_mask) - u(alpha_mask)) ./ da(alpha_mask);
+        w(beta_mask) = (v1(beta_mask) - v2(beta_mask)) ./ db(beta_mask);
+    end
+end

Example_scripts/+stereopivfunctions/run_correlation.m

@@ -0,0 +1,55 @@
+function [x,y,u,v,typevector, correlation_map] = run_correlation(image_dir, filenames1, filenames2, p, s, nr_of_cores, tform, output_view)
+%RUN_CORRELATION Taken from the "Example_PIVlab_commandline" script with
+%small addition for passing the required transform data
+%   Detailed explanation goes here
+    arguments
+        image_dir
+        filenames1
+        filenames2
+        p cell
+        s cell
+        nr_of_cores
+        tform 
+        output_view
+    end
+
+    n = length(filenames1);
+    x = cell(n,1);
+    y = x;
+    u = x;
+    v = x;
+    typevector = x; % typevector will be 1 for regular vectors, 0 for masked areas
+    correlation_map = x; % correlation coefficient
+    if nr_of_cores > 1
+        try
+		    local_cluster = parcluster('local'); % single node
+		    corenum = local_cluster.NumWorkers; % fix : get the number of cores available
+	    catch
+		    warning('on');
+		    warning('parallel local cluster can not be created, assigning number of cores to 1');
+		    nr_of_cores = 1;
+        end
+    end
+    if nr_of_cores > 1
+	    if misc.pivparpool('size') < nr_of_cores
+		    misc.pivparpool('open', nr_of_cores);
+	    end
+
+	    parfor i = 1:numel(x)  % index must increment by 1
+		    [x{i}, y{i}, u{i}, v{i}, typevector{i},correlation_map{i}] = ...
+			    piv.piv_analysis(image_dir, filenames1{i}, filenames2{i},p,s,nr_of_cores,false);
+        end
+    else % sequential loop
+	    for i = 1:numel(x)  % index must increment by 1
+		    [x{i}, y{i}, u{i}, v{i}, typevector{i}, correlation_map{i}] = ...
+			    piv.piv_analysis(image_dir, filenames1{i}, filenames2{i},p,s,nr_of_cores,false, tform, output_view);
+            if i == 1
+                disp("Processing: ")
+            else
+                fprintf(repmat('\b', 1, length(progress_string) + 1));
+            end
+            progress_string = [int2str(i / numel(x) * 100), ' %'];
+            disp(progress_string);
+        end
+    end
+end