Update demo.m, figures/demo.png, images/Chaib2023/Case1/img1.im7,...

Update demo.m, figures/demo.png, images/Chaib2023/Case1/img1.im7, images/Chaib2023/Case1/img2.im7, images/Chaib2023/Case1/img3.im7, images/Chaib2023/Case2/img1.im7, images/Chaib2023/Case2/img2.im7, images/Chaib2023/Case2/img3.im7, images/Chaib2023/Case3/img1.im7, images/Chaib2023/Case3/img3.im7, images/Chaib2023/Case3/img2.im7, scripts/fdetect.m

demo.m

+%% Resetting the console
+
+clear
+close all
+clc
+
+%% Setup
+
+% Adding script path
+addpath("scripts/");
+
+% Computing width of morphological mask
+delta_L = [.44, .44, .47];
+spatial_res = 0.0703; % mm/px
+w = [.5, 1.5, 1.5].*delta_L./spatial_res;
+
+%% Reading and plotting data
+figure();
+
+t = tiledlayout(3,3);
+t.Padding = 'tight';
+
+for case_id = 1:3
+    for img_id = 1:3
+        % Image path
+        path = "images/Chaib2023/Case"+case_id+"/img"+img_id+".im7";
+        % Reading .im7 image
+        v = loadvec(char(path));
+        I_raw = rot90(v.w);
+        I_raw = I_raw-min(I_raw(:));
+        I_cropped = I_raw(250:699,150:599);
+        I_cropped = uint8(rescale(I_cropped).*255);
+        % Detecting flame front
+        F = fdetect(I_cropped, w(case_id));
+        % Plotting
+        nexttile();
+        BG = I_cropped;
+        imshow(BG,'InitialMagnification','fit');
+        hold on;
+        [Y,X] = find(F);
+        plot(X,Y,'.r','MarkerSize',10);
+        set(gca,'YDir','reverse','TickLabelInterpreter','latex','FontSize',30);
+        axis on;
+        caxis([min(BG(:)) max(BG(:))]);
+        xlabel('x [px]',"Interpreter","latex","FontSize",30);
+        ylabel('y [px]',"Interpreter","latex","FontSize",30);
+        box on;
+        ax = gca;
+        axis on;
+        ax.LineWidth = 4;   
+    end
+
+end
+
+
+set(gcf,'OuterPosition',[500 500 1000 1000]);
+
+
+

scripts/fdetect.m

+function F = fdetect(I,w)
+
+%FDETECT Detect the position of the flame front in a grayscale images using
+%the Filtered Canny algorithm
+%
+%   FDETECT takes a raw OH-PLIF image I and a mask width w as its input, and
+%   returns a binary image F of the same size as I with 1's where the function
+%   finds edges in I and 0's elsewhere.
+%
+%   FDETECT can be broken down to four steps :
+%       1) Pre-processing (using two different schemes destined to
+%       segmentation and edge detection respectively).
+%       2) Enhanced Otsu segmentation
+%       3) Mask construction
+%       4) Filtered Canny edge detection
+%
+
+
+%% Pre-processing
+
+% Two different pre-processing schemes are used depending on which flame
+% front detection approach (segmentation vs. edge detection) is used.
+% Segmentation makes use of very mild filtering due to its high sensitivity
+% to the blurring effects induced by filtering.
+
+%%% Segmentation pre-processing scheme
+    % Filling new image array
+    Ifs=I;
+    % Applying constrast-limited adaptive histogram equalization (CLAHE)
+    Ifs=adapthisteq(Ifs,'NumTiles',[8,8],'ClipLimit',0.01);
+    % Applying a median filter iteratively for background subtraction
+    for ifilt = 1:20
+        Ifs = medfilt2(Ifs,[3 3]);
+    end
+    % Applying the non-linear diffusion (NLD) filter for smoothing
+    Ifs = imdiffusefilt(Ifs,'ConductionMethod','quadratic','NumberOfIterations',20);
+
+%%% Edge detection pre-processing scheme
+    % Filling new image array
+    Ife=I;
+    % Applying a median filter iteratively for background subtraction
+    for ifilt = 1:20
+        Ife = medfilt2(Ife,[3 3]);
+    end
+    % Applying the bilateral filter for smoothing
+    Ife = imbilatfilt(Ife,500,10);
+    % Applying the non-linear diffusion (NLD) filter for smoothing
+    Ife = imdiffusefilt(Ife,'ConductionMethod','quadratic','NumberOfIterations',50);
+
+%% Enhanced Otsu segmentation
+
+%%% Stage 1 : Detecting a preliminary Otsu contour
+
+    % Applying Otsu thresholding to the filtered image
+    Ibin = imbinarize(Ifs);        
+    % Removing small stray white pixels
+    Ibin = bwareaopen(Ibin,20);       
+    % Selecting a preliminary contour containing 
+    otsu_prelim = logical(bwperim(Ibin));       
+    % Flame object classification
+        % Setting up temporary variables
+        Ibin_temp = Ibin;
+        % Finding connected components in the temporary binary image
+        Z = bwconncomp(Ibin_temp);
+        % Counting the number of elements in each connected components
+        numPixels = cellfun(@numel,Z.PixelIdxList);
+        % Finding the index of the largest object within the binary image
+        [~,idx] = max(numPixels);
+        % Discarding all objects but the largest one in the binary image
+        for ibw = 1:Z.NumObjects
+            if ibw ~= idx
+                Ibin_temp(Z.PixelIdxList{ibw}) = 0;
+            end
+        end
+        % Classifying objects into one of three categories : main front (mf),
+        % product pockets (pp), and reactant pockets (rp)
+        mf = bwperim(imfill(Ibin_temp,'holes'));
+        pp = bwperim(Ibin_temp-Ibin);
+        rp = logical(otsu_prelim-pp-mf);
+    % Cleaning the main front (mf)
+        % Detecting the position of bottom tips of OH fluorescence
+        [Y,X] = find(mf==1);
+        ylim1 = max(Y(1:floor(length(Y)/2)));
+        xlim1=X(find(Y==ylim1,1,'first'));
+        ylim2 = max(Y(floor(length(Y)/2):end));
+        xlim2=X(find(Y==ylim2,1,'last'));
+        % Discarding unwanted contours
+        Ibin_temp2 = imfill(Ibin_temp,'holes');
+        Ibin_temp2(ylim1-5:ylim1,1:xlim1) = 1;
+        Ibin_temp2(ylim2-5:ylim1,xlim2:end) = 1;
+        Ibin_temp2(1,:) = 1;
+        Ibin_temp2(1:ylim1-5,1)=1;
+        Ibin_temp2(1:ylim2-5,end)=1;
+        mf = bwperim(imfill(Ibin_temp2,'holes'));
+        mf(ylim1-5:ylim1,1:xlim1) = 0;
+        mf(ylim2-5:ylim1,xlim2:end) = 0;
+    % Cleaning product pockets (rp) - Keeping only those upstream the main 
+    % flame front
+        pp = pp.*(1-imfill(Ibin_temp2,'holes'));
+
+        
+%%% Stage 2 : Pocket identification
+
+    % Removing the smallest objects in the burnt gas region
+    rp = bwareaopen(rp,20);
+    % Computing gradient magnitudes of potential reactant pockets
+    g_otsu = double(rp).*rescale(imgradient(Ife));
+    % Extracting the gradients of each potential reactant pocket
+    E = bwconncomp(rp);
+    % Finding and discarding false low-gradient pockets
+    th_pockets = .2;
+    rp_true = zeros(size(rp));    
+    for ipockets = 1:E.NumObjects
+        mean_gradient = mean(g_otsu(E.PixelIdxList{1,ipockets}));
+        if(mean_gradient<th_pockets)
+            E.PixelIdxList{1,ipockets} = [];
+        else
+            rp_true(E.PixelIdxList{1,ipockets}) = 1;
+        end
+    end
+
+%%% Constructing final Otsu contour
+% Combining main front, product pockets, and true reactant pockets
+F_otsu = logical(mf+pp+rp_true);
+F_otsu(1,:) = 0;
+F_otsu(:,1) = 0;
+F_otsu(end,:) = 0;
+F_otsu(:,end) = 0;
+
+%% Mask construction
+
+% The width provided by the user is used to construct a 2D binary
+% morphological mask. Recommended values of w are one to twice the
+% unstretched laminar flame thickness (conversion should be made to pixels)
+
+mask = bwdist(F_otsu)<=w;
+
+
+%% Filtered Canny algorithm
+
+%%%% Gaussian filter standard deviation (static value of 2 recommended, higher values increase smoothing)
+sigma = 2;
+
+%%% Canny algorithm (MATLAB implementation - obtained from edge.m)
+% Filtering and computing 2D gradient map
+    % Determine filter length
+    filterExtent = ceil(4*sigma);
+    x = -filterExtent:filterExtent;  
+    % Create 1-D Gaussian Kernel
+    c = 1/(sqrt(2*pi)*sigma);
+    gaussKernel = c * exp(-(x.^2)/(2*sigma^2));  
+    % Normalize to ensure kernel sums to one
+    gaussKernel = gaussKernel/sum(gaussKernel);
+    % Create 1-D Derivative of Gaussian Kernel
+    derivGaussKernel = gradient(gaussKernel);  
+    % Normalize to ensure kernel sums to zero
+    negVals = derivGaussKernel < 0;
+    posVals = derivGaussKernel > 0;
+    derivGaussKernel(posVals) = derivGaussKernel(posVals)/sum(derivGaussKernel(posVals));
+    derivGaussKernel(negVals) = derivGaussKernel(negVals)/abs(sum(derivGaussKernel(negVals)));
+    % Compute smoothed numerical gradient of image I along x (horizontal)
+    % direction. GX corresponds to dG/dx, where G is the Gaussian Smoothed
+    % version of image I.
+    GX = imfilter(double(Ife), gaussKernel', 'conv', 'replicate');
+    GX = imfilter(GX, derivGaussKernel, 'conv', 'replicate');    
+    % Compute smoothed numerical gradient of image I along y (vertical)
+    % direction. GY corresponds to dG/dy, where G is the Gaussian Smoothed
+    % version of image I.
+    GY = imfilter(double(Ife), gaussKernel, 'conv', 'replicate');
+    GY  = imfilter(GY, derivGaussKernel', 'conv', 'replicate');    
+    % Gradient magnitude
+    G = hypot(GX,GY);
+    G = G./max(G(:));
+% Filtering gradient map using the binary morphological mask
+GX = ((mask).*GX);
+GY = ((mask).*GY);
+Gf = hypot(GX,GY);
+Gf = Gf./max(Gf(:));
+% Non-maxima suppression (NMS)
+lowThresh = 0; % (Should be tuned between 0 and 0.1, lower values result in fewer discontinuities across the flame front but may introduce spurs)
+highThresh = 0;     
+E = images.internal.builtins.cannyFindLocalMaxima(GX,GY,Gf,lowThresh);    
+if ~isempty(E)
+    [rstrong,cstrong] = find(G>highThresh & E);  
+    if ~isempty(rstrong) % result is all zeros if idxStrong is empty
+        F = bwselect(E, cstrong, rstrong, 8);
+    else
+        F = false(size(E));
+    end
+else
+    F = false(size(E));
+end
+%%% Cleaning final flame edge
+F = bwareaopen(bwskel(F),40);
+
+
+end
+
+
+
+
+
+
+
+