g_stabilize a ete revise en profondeur.

src/g_rect/g_rect.m

@@ -167,12 +167,14 @@ fprintf('\n')
 % image(imread(imgFname));
 imagesc(imread(imgFname));
 colormap(gray);
+hold on
 for i = 1:ngcp
-  text(i_gcp(i),j_gcp(i),num2str(i),'color','r','horizontalalignment','center');
+  plot(i_gcp(i),j_gcp(i),'r.');
+  text(i_gcp(i),j_gcp(i),num2str(i),'color','r','horizontalalignment','right');
 end
 title('Ground Control Points','color','r');
 
-print -dpng image1.png
+print -dpng -r300 image1.png
 
 fprintf('\n')
 ok = input('Is it ok to proceed with the rectification (y/n): ','s');
@@ -359,11 +361,12 @@ else
 end
 %%
 
+% Compute the effective resolution
+Delta = g_res(LON, LAT, field);
+
 fprintf('\n')
 fprintf('Saving rectification file in: %s\n',outputFname);
 
-errGeoFit = 0.0;
-
 save(outputFname,'imgFname','firstImgFname','lastImgFname',...
                  'LON','LAT',...
                  'LON0','LAT0',...
@@ -371,7 +374,7 @@ save(outputFname,'imgFname','firstImgFname','lastImgFname',...
                  'i_gcp','j_gcp',...
                  'hfov','lambda','phi','H','theta',...
                  'errGeoFit','errPolyFit',...
-                 'precision');
+                 'precision','Delta');
 
 clear LON LAT
 

src/g_rect/g_res.m

-%function Delta = g_res(LON,LAT,i,j,field);
-function [Delta DeltaX DeltaY] = g_res(LON,LAT,i,j,field);
-
-[m n] = size(LON);
-
-if i == 1
-  Delta1 = g_dist(LON(i+1,j),LAT(i+1,j),LON(i,j),LAT(i,j),field);
-elseif i == m
-  Delta1 = g_dist(LON(i,j),LAT(i,j),LON(i-1,j),LAT(i-1,j),field);
-else
-  Delta1 = g_dist(LON(i+1,j),LAT(i+1,j),LON(i-1,j),LAT(i-1,j),field);
-  Delta1 = Delta1/2;  
-end
-
-if j == 1
-  Delta2 = g_dist(LON(i,j+1),LAT(i,j+1),LON(i,j),LAT(i,j),field);
-elseif j == n
-  Delta2 = g_dist(LON(i,j),LAT(i,j),LON(i,j-1),LAT(i,j-1),field);
-else
-  Delta2 = g_dist(LON(i,j+1),LAT(i,j+1),LON(i,j-1),LAT(i,j-1),field);
-  Delta2 = Delta2/2;  
-end
-
-Delta = sqrt(Delta1^2 + Delta2^2);
-DeltaX = sqrt(Delta1^2);
-DeltaY = sqrt(Delta2^2);
-
+function res = g_res_field(LON, LAT, field)
+% g_res_field - Function calculating the loss of resolution with de distance
+% from the camera
+%
+% Inputs:
+% 
+%    LON   : Longitude matrix obtained from g_rect
+%    LAT   : Latitude matrix obtained from g_rect
+%    field : 'true' for field case. 'false' for lab case.
+%
+% Outputs:
+%    
+%    res : The field of resolution which is degrading with distance from
+%    the camera
+% 
+% Author: Daniel Bourgault
+%         Institut des sciences de la mer de Rimouski
+% 
+% email: daniel_bourgault@uqar.ca 
+% Website: http://demeter.uqar.ca/g_rect/
+% February 2013
+%
+% Author: Elie Dumas-Lefebvre
+%         Institut des Science de la Mer de Rimouski
+% 
+% Note: Matricial formulation of the calculation for a faster execution
+%
+% email: elie.dumas-lefebvre@uqar.ca
+% March 2019
+%
+
+% LAT and LON difference in the x-axis
+dLATi = diff(LAT, 1, 2);
+dLONi = diff(LON, 1, 2);
+
+% LAT and LON difference in the y-axis
+dLATj = diff(LAT, 1, 1);
+dLONj = diff(LON, 1, 1);
+
+% Mean LAT in both x and y axis
+LAT_meani = 0.5*(LAT(:, 2:end) + LAT(:, 1:end-1));
+LAT_meanj = 0.5*(LAT(2:end, :) + LAT(1:end-1, :));
+
+% Add a point at boundaries such that the resolution matrix 
+% has the same size of the LAT LON matrices 
+dLATi(:,end+1)     = dLATi(:,end);
+dLONi(:,end+1)     = dLONi(:,end);
+dLATj(end+1,:)     = dLATj(end,:);
+dLONj(end+1,:)     = dLONj(end,:);
+LAT_meani(:,end+1) = LAT_meani(:,end);
+LAT_meanj(end+1,:) = LAT_meani(end,:);
+
+% Conversion from degree to meters
+meterPerDegLat  = 1852*60.0;
+meterPerDegLoni = meterPerDegLat.*cosd(LAT_meani);
+meterPerDegLonj = meterPerDegLat.*cosd(LAT_meanj);
+
+if field 
+    
+    % Conversion from degrees to meters
+    dxi = dLONi.*meterPerDegLoni;
+    dxj = dLONj.*meterPerDegLonj;
+    dyi = dLATi.*meterPerDegLat;
+    dyj = dLATj.*meterPerDegLat;
+    
+else
+
+    % Lab case. Data are already in meters. No conversion required.    
+    dxi = dLONi;
+    dxj = dLONj;
+    dyi = dLATi;
+    dyj = dLATj;
+    
+end
+
+% Distances in x and y axis
+deltai = sqrt(dxi.^2 + dyi.^2);
+deltaj = sqrt(dxj.^2 + dyj.^2);
+
+% field of resolution
+res = sqrt(deltai.^2 + deltaj.^2);
+
 end
\ No newline at end of file

src/g_stabilize/g_computemotion.m

 function[ A, T ] = g_computemotion( fx, fy, ft, roi )
+
 [ydim,xdim] = size(fx);
 [x,y] = meshgrid( [1:xdim]-xdim/2, [1:ydim]-ydim/2 );
 
 %%% TRIM EDGES
-fx = fx( 3:end-2, 3:end-2 );
-fy = fy( 3:end-2, 3:end-2 );
-ft = ft( 3:end-2, 3:end-2 );
+fx  = fx( 3:end-2, 3:end-2 );
+fy  = fy( 3:end-2, 3:end-2 );
+ft  = ft( 3:end-2, 3:end-2 );
 roi = roi( 3:end-2, 3:end-2 );
-x = x( 3:end-2, 3:end-2 );
-y = y( 3:end-2, 3:end-2 );
+x   = x( 3:end-2, 3:end-2 );
+y   = y( 3:end-2, 3:end-2 );
 
 ind = find( roi > 0 );
-x = x(ind); y = y(ind);
-fx = fx(ind); fy = fy(ind); ft = ft(ind);
-
-xfx = x.*fx; xfy = x.*fy; yfx = y.*fx; yfy = y.*fy;
-M(1,1) = sum( xfx .* xfx ); M(1,2) = sum( xfx .* yfx ); M(1,3) = sum( xfx .* xfy );
-M(1,4) = sum( xfx .* yfy ); M(1,5) = sum( xfx .* fx ); M(1,6) = sum( xfx .* fy );
-M(2,1) = M(1,2); M(2,2) = sum( yfx .* yfx ); M(2,3) = sum( yfx .* xfy );
-M(2,4) = sum( yfx .* yfy ); M(2,5) = sum( yfx .* fx ); M(2,6) = sum( yfx .* fy );
-M(3,1) = M(1,3); M(3,2) = M(2,3); M(3,3) = sum( xfy .* xfy );
-M(3,4) = sum( xfy .* yfy ); M(3,5) = sum( xfy .* fx ); M(3,6) = sum( xfy .* fy );
-M(4,1) = M(1,4); M(4,2) = M(2,4); M(4,3) = M(3,4);
-M(4,4) = sum( yfy .* yfy ); M(4,5) = sum( yfy .* fx ); M(4,6) = sum( yfy .* fy );
-M(5,1) = M(1,5); M(5,2) = M(2,5); M(5,3) = M(3,5);
-M(5,4) = M(4,5); M(5,5) = sum( fx .* fx ); M(5,6) = sum( fx .* fy );
-M(6,1) = M(1,6); M(6,2) = M(2,6); M(6,3) = M(3,6);
-M(6,4) = M(4,6); M(6,5) = M(5,6); M(6,6) = sum( fy .* fy );
+x   = x(ind); 
+y   = y(ind);
+fx  = fx(ind); 
+fy  = fy(ind); 
+ft  = ft(ind);
+
+xfx = x.*fx; 
+xfy = x.*fy; 
+yfx = y.*fx; 
+yfy = y.*fy;
+
+M(1,1) = sum( xfx .* xfx ); 
+M(1,2) = sum( xfx .* yfx ); 
+M(1,3) = sum( xfx .* xfy );
+M(1,4) = sum( xfx .* yfy ); 
+M(1,5) = sum( xfx .* fx ); 
+M(1,6) = sum( xfx .* fy );
+M(2,1) = M(1,2); 
+M(2,2) = sum( yfx .* yfx ); 
+M(2,3) = sum( yfx .* xfy );
+M(2,4) = sum( yfx .* yfy ); 
+M(2,5) = sum( yfx .* fx ); 
+M(2,6) = sum( yfx .* fy );
+M(3,1) = M(1,3); 
+M(3,2) = M(2,3); 
+M(3,3) = sum( xfy .* xfy );
+M(3,4) = sum( xfy .* yfy ); 
+M(3,5) = sum( xfy .* fx ); 
+M(3,6) = sum( xfy .* fy );
+M(4,1) = M(1,4); 
+M(4,2) = M(2,4); 
+M(4,3) = M(3,4);
+M(4,4) = sum( yfy .* yfy ); 
+M(4,5) = sum( yfy .* fx ); 
+M(4,6) = sum( yfy .* fy );
+M(5,1) = M(1,5); 
+M(5,2) = M(2,5); 
+M(5,3) = M(3,5);
+M(5,4) = M(4,5); 
+M(5,5) = sum( fx .* fx ); 
+M(5,6) = sum( fx .* fy );
+M(6,1) = M(1,6); 
+M(6,2) = M(2,6); 
+M(6,3) = M(3,6);
+M(6,4) = M(4,6); 
+M(6,5) = M(5,6); 
+M(6,6) = sum( fy .* fy );
 
 k = ft + xfx + yfy;
 
-b(1) = sum( k .* xfx ); b(2) = sum( k .* yfx );
-b(3) = sum( k .* xfy ); b(4) = sum( k .* yfy );
-b(5) = sum( k .* fx ); b(6) = sum( k .* fy );
+b(1) = sum( k .* xfx ); 
+b(2) = sum( k .* yfx );
+b(3) = sum( k .* xfy ); 
+b(4) = sum( k .* yfy );
+b(5) = sum( k .* fx ); 
+b(6) = sum( k .* fy );
 
 %v = inv(M) * b';
 v = M\b';
 
-%A = [1 0 ; 0 1];
-
+%A = [1 0 ; 0 1]; % For translation only 
 A = [v(1) v(2) ; v(3) v(4)];
-T = [v(5) ; v(6)];
\ No newline at end of file
+T = [v(5) ; v(6)];

src/g_stabilize/g_down.m

 function[ f ] = g_down( f, L );
+
 blur = [1 2 1]/4;
+
 for k = 1 : L
-f = conv2( conv2( f, blur, 'same' ), blur', 'same' );
-f = f(1:2:end,1:2:end);
+    f = conv2( conv2( f, blur, 'same' ), blur', 'same' );
+    f = f(1:2:end,1:2:end);
 end
\ No newline at end of file

src/g_stabilize/g_opticalflow.m

 function[ Acum, Tcum ] = g_opticalflow( f1, f2, roi, L )
+
 f2orig = f2;
 Acum = [1 0 ; 0 1];
 Tcum = [0 ; 0];
 
 for k = L : -1 : 0
-%for k = L : -1 : L
     
     %%% DOWN-SAMPLE
     f1d = g_down( f1, k );
     f2d = g_down( f2, k );
     ROI = g_down( roi, k );
+    
     %%% COMPUTE MOTION
     [Fx,Fy,Ft] = g_spacetimederiv( f1d, f2d );
     [A,T]      = g_computemotion( Fx, Fy, Ft, ROI );
     T = (2^k) * T;
     [Acum,Tcum] = g_accumulatewarp( Acum, Tcum, A, T );
+    
     %%% WARP ACCORDING TO ESTIMATED MOTION
     f2 = g_warp( f2orig, Acum, Tcum );
     

src/g_stabilize/g_roi.m

@@ -58,4 +58,4 @@ roi(ij) = 1;
 figure(2);
 imagesc(roi);
 
-save roi.mat roi imgFname
\ No newline at end of file
+save([imgFname(1:end-4),'_roi.mat'],'roi','imgFname');
\ No newline at end of file

src/g_stabilize/g_stabilize.m

@@ -9,66 +9,100 @@ function g_stabilize(img_ref_fname,varargin)
 %                  necessary), of the reference image in the current 
 %                  directory.
 %
-%   Parameters, Name,       Value
-%               roiFile,    [Defaults: roi.mat] Name of the file with roi
-%                           information.
-%               auto,       [Defaults: 0 (false)] If true, overrides any
-%                           interaction with the user (for use in scripts)
-%               suffix,     [Defaults: '_stable'] Suffix to add to the
-%                           stabilized images.
+%   Optional parameters:
+%       L_level_Gaussian:   Size of the L-level Gaussian pyramid [Default: 4]
+%       roi_file:           Name of the file with roi information [Default: 'roi.mat']
+%       fname_suffix:       Suffix to add to the stabilized images [Default: '_stable']
+%       output_folder:      Name of the output folder [Default: 'stable']
+%       equalize:           Use image equalization over the roi region [Default: 0 (false)]
+%       use_img_grad:       Perform the stabilization using the gradient of the image [Default: 0 (false)]
+%       auto:               If true, overrides any interaction with the user 
+%                           (for use in scripts) [Default: 0 (false)] 
 %
 % Output:
 %   All stabilized images are rewritten as image files in the subfolder 
-%   /stable.
+%   /stable (by default or other name is specified).
 %
-% Example:
-%   g_stabilize('IMG_4018.jpg');
+% Examples:
+%   Example 1: Using all default parameters
+%   >> g_stabilize('IMG_4018.jpg');
 %
+%   Example 2: Using a different L level Gaussian than the default value (4)
+%   >> g_stabilize('IMG_4018.jpg','L_level_Gaussian',6);
+%
+%   Example 3: Using a different L level Gaussian than the default value (4)
+%              and a different suffix than the default (_stable) and
+%              without any image equalization
+%   >> g_stabilize('IMG_4018.jpg','L_level_Gaussian',6,'fname_suffix','_very_stable','equalize',false);
+%
+
+%%
+p = inputParser;
+
+% Default values
+default_L_level_Gaussian = 4;
+default_roi_file         = 'roi.mat';
+default_fname_suffix     = '_stable';
+default_output_folder    = 'stable';
+default_equalize         = false;
+default_use_img_grad     = false;
+default_auto             = false;
+
+addRequired(p,'img_ref_fname',@ischar);
+addParameter(p,'L_level_Gaussian',default_L_level_Gaussian,@isnumeric);
+addParameter(p,'roi_file',default_roi_file,@ischar);
+addParameter(p,'fname_suffix',default_fname_suffix,@ischar);
+addParameter(p,'output_folder',default_output_folder,@ischar);
+addParameter(p,'equalize',default_equalize,@islogical);
+addParameter(p,'use_img_grad',default_use_img_grad,@islogical);
+addParameter(p,'auto',default_auto,@islogical);
+
+parse(p,img_ref_fname,varargin{:})
 
-%% Some parameters
-% Pyramid level for the stabilization.
-% See documentation in the g_stabilize folder.
-% L = 4 worked well on initial application but may require adjustments for
-% particular application.
-L = 4;
+img_ref_fname = p.Results.img_ref_fname;
+L             = p.Results.L_level_Gaussian;
+roi_file      = p.Results.roi_file;
+fname_suffix  = p.Results.fname_suffix;
+output_folder = p.Results.output_folder;
+equalize      = p.Results.equalize;
+use_img_grad  = p.Results.use_img_grad;
+auto          = p.Results.auto;
 
-% Number of iterative improvement. Probably ok just 1 iteration.
-niter = 1;
+disp(['   ']);
+disp(['   -----------------------------------------------------------------------']);
+disp(['   IMAGE STABILIZATION PARAMETERS']);
+disp(['     Filename of the reference image:              ',img_ref_fname]);
+disp(['     L-level Gaussian:                             ',num2str(L)]);
+disp(['     Name of the ROI filename:                     ',roi_file]);
+disp(['     Suffix of stabilized filenames:               ',fname_suffix]);
+disp(['     Output folder:                                ',output_folder]);
+disp(['     Equalize (1) or not (0) the roi region:       ',num2str(equalize)]);
+disp(['     Use (1) or not (0) the gradient of the image: ',num2str(use_img_grad)]);
+disp(['   -----------------------------------------------------------------------']);
 
+%% Equalization parameters
 % Parameters for equalization. See help clahs for details.
-% This may not be needed fpr most application but may help other cases.
+% This may not be needed for most applications but may help in some cases 
+% where the lighting conditions changes a lot.
 % Need to play with it if problems. 
-equalize = false;
-nry = 4;
-nrx = 4;
+if equalize
+    nry = 8;
+    nrx = 8;
+end
+
+%%
 
 % Extract the path and extension of the reference image
 [pathstr,name,ext] = fileparts(img_ref_fname);
 
-% Check for optional options
-fname_suffix = '_stable';
-auto = 0;
-roifile = 'roi.mat';
-if length(varargin) > 1
-    for i=1:2:length(varargin)
-        if strcmpi(varargin{i},'roiFile')
-            roifile = varargin{i+1};
-        elseif strcmpi(varargin{i},'auto')
-            auto = varargin{i+1};
-        elseif strcmpi(varargin{i},'suffix')
-            fname_suffix = varargin{i+1};
-        end
-    end
-end
-
-hw = waitbar(0,'Preparing images');
 
 %% Load the Region of interest (roi) file
-display(' Loading the roi.mat file containing the region of interest (roi)');
+disp(['   ']);
+display('   Loading the file containing the region of interest (roi)...');
 if length(pathstr) == 0  % Images are in the current folder
-    load(roifile);
+    load(roi_file);
 else
-    load([pathstr,'/',roifile]);
+    load([pathstr,'/',roi_file]);
 end
 
 % Read the reference image
@@ -76,71 +110,76 @@ im2 = imread(img_ref_fname);
 
 % Convert into gray scale if not already a B/W image
 if size(im2,3) ~= 1
-    im2 = g_rgb2gray(im2);
+    im2_BW = g_rgb2gray(im2);
 end
 
+im2_BW_cooked = im2_BW;
+clear im2;
+
 % May equalize with this equalizer. Careful this doesn't always work well
 % and may compromise the stabilization. Use with care. 
 if equalize
-    im2 = clahs(im2,nry,nrx);
+    im2_BW_cooked = clahs(im2_BW_cooked,nry,nrx);
 end
 
-im2 = double(im2);
+im2_BW        = double(im2_BW);
+im2_BW_cooked = double(im2_BW_cooked);
 
 % Remove mean and divide by the standard deviation.
-% This generally help the algorithm to compare images
-inorm    = find(roi > 0);
-im2_norm = (im2 - nanmean(im2(inorm)))./nanstd(im2(inorm));
+% This generally help the algorithm to compare images that have different
+% light conditions
+%inorm    = find(roi > 0);
+%im2_BW_cooked = (im2_BW_cooked - nanmean(im2_BW_cooked(inorm)))./nanstd(im2_BW_cooked(inorm));
 
-% Take the norm of the gradient of the image. This helps the
-% stabilization algorithm, nut again you may want to explore
-% without it
-[im2x, im2y] = gradient(im2_norm);
-grad_im2     = sqrt(im2x.^2 + im2y.^2);
+if use_img_grad
+    % Take the norm of the gradient of the image. This may help the
+    % stabilization algorithm but not always. Use this option with care. 
+    [im2x, im2y]  = gradient(im2_BW_cooked);
+    im2_BW_cooked = sqrt(im2x.^2 + im2y.^2);    
+end
+
+frames(2).im = im2_BW_cooked;
 
-% The second frame is the reference frame for the stabilization algorithm.
-frames(2).im = grad_im2;
 
-%% Some figures to make sure everything seems ok
+%% Some figures to make sure everything is ok
 figure(1);
-imagesc(im2);
+imagesc(im2_BW);
 colormap(gray);
-title('Reference image. Intensity normalized over ROI');
+title('Reference image in B/W');
 
 figure(2);
-imagesc(log10(grad_im2+eps));
+imagesc(im2_BW_cooked);
 colormap(gray);
-title('Gradient of the reference image on a log scale');
+title('Reference image manipulated (if so)');
 
 figure(3);
-imagesc(roi);
-title('roi');
+% Produce a mask over non-ROI region and plot it
+iNaN = find(roi == 0);
+roi_with_NaN = roi;
+roi_with_NaN(iNaN) = NaN;
+imagesc(im2_BW_cooked.*roi_with_NaN);
+colormap(gray);
+title('Reference image manipulated with mask over non-ROI');
 
-if ~auto
-    answer = input('Happy with roi (y/n)? ','s');
-    if answer == 'n'
-        return
-    end
-end
+%if ~auto
+%    disp([' '])
+%    answer = input('   Happy with roi and parameters (y/n)? ','s');
+%    if answer == 'n'
+%        return
+%    end
+%end
+
+display(['    ']);
 
 %% Create the folder where to write the stable images
 if length(pathstr) == 0  % Images are in the current folder
-    [status,message,messageid] = mkdir('stable');
+    [status,message,messageid] = mkdir(output_folder);
+    save([output_folder,'/g_stabilize_info.mat'],'p','roi');
 else
-    [status,message,messageid] = mkdir([pathstr,'/stable']);
+    [status,message,messageid] = mkdir([pathstr,'/',output_folder]);
+    save([pathstr,'/',output_folder,'/g_stabilize_info.mat'],'p','roi');
 end
 
-% Check for images already stablilized that would be found in the 
-% stable folder. This is not yet fully implemeted. The goal would
-% be not to perform the stabilization to images already stabilized
-%if strcmp(message,'Directory already exists.')
-%    if length(pathstr) == 0  
-%        img_done = dir(['stable/*',ext]);
-%    else
-%        img_done = dir([pathstr,'/stable/','*',ext]);
-%    end
-%end
-
 %% Read all filemanes in the current directory with the same extension as the ref img
 if length(pathstr) == 0  % Images are in the current folder
     all_img = dir(['*',ext]);
@@ -155,15 +194,8 @@ time_per_img = 0;
 
 for i = 1:N
     
-    
     time_rem = (N-i+1)*time_per_img;
-    if i > 1
-        waitbar((i-1)/N,hw,sprintf('Stabilizing image #%d. %.0fh%02.2fmn remaining',...
-            i,fix(time_rem/3600),rem(time_rem,3600)/60));
-    else
-        waitbar((i-1)/N,hw,sprintf('Stabilizing image #%d.',...
-            i,fix(time_rem/3600),rem(time_rem,3600)/60));
-    end
+    
     total = tic;
     
     if length(pathstr) == 0
@@ -174,32 +206,62 @@ for i = 1:N
         im_info = imfinfo([pathstr,'/',all_img(i).name]);
     end
     
+    if i == 1
+        display(['   Processing image: ',all_img(i).name]);
+        display(['   Time remaining:   Yet undetermined. Please wait until next image.']);
+        display(['  ']);
+    else
+        display(['   Processing image: ',all_img(i).name]);
+        display(['   Time remaining:   ',datestr(time_rem/86400,'HH:MM:SS')]);
+        display(['  ']);
+    end
+
     % Convert into gray scale if not already a B/W image
     if size(im1,3) ~= 1
-        im1 = g_rgb2gray(im1);
+        im1_BW = g_rgb2gray(im1);
     end
     
+    im1_BW_cooked = im1_BW;