Home > Matching_SYM_LSM > src > 2D_LSM_62 > 2D_LSM_62_sym_warp > LSM_62_sym_warp_main.m

LSM_62_sym_warp_main

PURPOSE ^

% LS-matching 6+2 of two images g and h, symmetric

SYNOPSIS ^

function [est_x,est_sigma_0,NN,Mf,Red,N_iter,max_ratio] =LSM_62_sym_warp_main(g,h,vg,vh,A_a,R_a,sigma_smooth,max_iter,Nf_MIN,pt)

DESCRIPTION ^

% LS-matching 6+2 of two images g and h, symmetric

 model: observed images g and h referring to unknown reference f

 Basic model:

 g(y) -> (A,R) -> h(z)

 geometric model:    z   = z(y,A), A = A(p1, ..p6)
 radiometric model:  h   = h(g,R), R = R(p7,p8)
 unknown parameters: p   = p(1:8), to be estimated



 Relation to unknown signal f:

 g(y) -> (B,S) -> f(x) -> (B,S) -> h(z)

 geometric transformation:    B(a) = A^(-1/2)(u)
 radiometric transformation;  S(a) = R^(-1/2)(u)
 unknown parameters:          u = u(1:8), estimated internally

 affine geometric and radiometric transformations
 g(y) = a_7^-1 * (f(B    *   y + a) - a_8) + n_g(y)
 h(z) = a_7    * (f(B^-1 *  (z - a) + a_8) + n_h(z)

 B_a  = [a1 a3 a5;    S = [a7 a8
         a2 a4 a6;          0  1]
          0  0 1]

 all images are square:
     size      (M x M), M odd
     radius     N=(M-1)/2,
     centre at (N+1,N+1)

 call:
 [est_x,est_sigma_0,NN,Mf,Red,N_iter] = ...
    LSM_62_sym_warp_main(g,h,vg,vh,A_a,R_a,sigma_smooth,max_iter,pt)

 g            = Ng x Ng matrix of first image, type = double
 h            = Nh x Nh matrix of second image, type = double
 vg           = 256-vector of noise variances [0:1:255] of g
 vh           = 256-vector of noise variances [0:1:255] of h
 A_a          = 3x3 matrix, approximate transformation geometric parameters
                   A_a  = [a b c;
                           d e f;
                           0 0 1]
                   zh = A_a * yh
 R_a          = 1x2 vector of approximate contrast and intensity deviation
                   h = R_a(1) * g + R_a(2)
 sigma_smooth = smoothing parameter for true image at first iteration
 max_iter     = maximum number of iterations
 Nf_MIN       = minimum size of overlap (radius of window f)
 pt           in (0,1,2) plot output control

 est_x       = estimated transformation parameters p with CovM
               est_x.x  = parameters p of A and R
               est_x.C  = covariance matrix
               est_x.Ah = 3x3 homogeneous geometric affinity
               est_x.R  = 1x2 radiometric parameters
               1..6: geometric parameters
               7..8: intenisty parameters
 est_sigma_0   estimated factor sigma_0 (should be close to 1)
 NN            number of observing pixels
 Mf            size of estimated image f
 Red           redundancy
 N_iter        number of iterations needed
 max_ratio     maximum ratio of corrections of parameters/sigma

 wf 7/2020

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %% LS-matching 6+2 of two images g and h, symmetric
0002 %
0003 % model: observed images g and h referring to unknown reference f
0004 %
0005 % Basic model:
0006 %
0007 % g(y) -> (A,R) -> h(z)
0008 %
0009 % geometric model:    z   = z(y,A), A = A(p1, ..p6)
0010 % radiometric model:  h   = h(g,R), R = R(p7,p8)
0011 % unknown parameters: p   = p(1:8), to be estimated
0012 %
0013 %
0014 %
0015 % Relation to unknown signal f:
0016 %
0017 % g(y) -> (B,S) -> f(x) -> (B,S) -> h(z)
0018 %
0019 % geometric transformation:    B(a) = A^(-1/2)(u)
0020 % radiometric transformation;  S(a) = R^(-1/2)(u)
0021 % unknown parameters:          u = u(1:8), estimated internally
0022 %
0023 % affine geometric and radiometric transformations
0024 % g(y) = a_7^-1 * (f(B    *   y + a) - a_8) + n_g(y)
0025 % h(z) = a_7    * (f(B^-1 *  (z - a) + a_8) + n_h(z)
0026 %
0027 % B_a  = [a1 a3 a5;    S = [a7 a8
0028 %         a2 a4 a6;          0  1]
0029 %          0  0 1]
0030 %
0031 % all images are square:
0032 %     size      (M x M), M odd
0033 %     radius     N=(M-1)/2,
0034 %     centre at (N+1,N+1)
0035 %
0036 % call:
0037 % [est_x,est_sigma_0,NN,Mf,Red,N_iter] = ...
0038 %    LSM_62_sym_warp_main(g,h,vg,vh,A_a,R_a,sigma_smooth,max_iter,pt)
0039 %
0040 % g            = Ng x Ng matrix of first image, type = double
0041 % h            = Nh x Nh matrix of second image, type = double
0042 % vg           = 256-vector of noise variances [0:1:255] of g
0043 % vh           = 256-vector of noise variances [0:1:255] of h
0044 % A_a          = 3x3 matrix, approximate transformation geometric parameters
0045 %                   A_a  = [a b c;
0046 %                           d e f;
0047 %                           0 0 1]
0048 %                   zh = A_a * yh
0049 % R_a          = 1x2 vector of approximate contrast and intensity deviation
0050 %                   h = R_a(1) * g + R_a(2)
0051 % sigma_smooth = smoothing parameter for true image at first iteration
0052 % max_iter     = maximum number of iterations
0053 % Nf_MIN       = minimum size of overlap (radius of window f)
0054 % pt           in (0,1,2) plot output control
0055 %
0056 % est_x       = estimated transformation parameters p with CovM
0057 %               est_x.x  = parameters p of A and R
0058 %               est_x.C  = covariance matrix
0059 %               est_x.Ah = 3x3 homogeneous geometric affinity
0060 %               est_x.R  = 1x2 radiometric parameters
0061 %               1..6: geometric parameters
0062 %               7..8: intenisty parameters
0063 % est_sigma_0   estimated factor sigma_0 (should be close to 1)
0064 % NN            number of observing pixels
0065 % Mf            size of estimated image f
0066 % Red           redundancy
0067 % N_iter        number of iterations needed
0068 % max_ratio     maximum ratio of corrections of parameters/sigma
0069 %
0070 % wf 7/2020
0071 
0072 function [est_x,est_sigma_0,NN,Mf,Red,N_iter,max_ratio] = ...
0073     LSM_62_sym_warp_main(g,h,vg,vh,A_a,R_a,sigma_smooth,max_iter,Nf_MIN,pt)
0074 
0075 %% ensure g and h are double
0076 g=double(g);
0077 h=double(h);
0078 
0079 % check for postiveness
0080 if det(A_a) < 0
0081     disp('Geometric affinity is not positive definite')
0082     est_x = 0;
0083     est_sigma_0 = 0;
0084     Red = -1;
0085     NN = 0;
0086     Mf = 0;
0087     N_iter = 0;
0088     return
0089 end
0090 
0091 %% test CPU-time
0092 NO_TEST_CPU = 1;
0093 sum_cpu=0;
0094 
0095 %% determine approximate values: take squareroot of given affinities
0096 
0097 % take squareroot from given affinities
0098 B_a = sqrtm(A_a);
0099 a7  = sqrt(R_a(1));
0100 a8  = R_a(2)/(1+a7);
0101 S_a = [a7,a8]';
0102 
0103 % set approximate values
0104 Be    = B_a(1:2,:);
0105 xa    = [Be(:);S_a];
0106 
0107 if pt > 0
0108     [im_lr,weite,border] = fill_left_right_images_0(g,h,max_iter);  
0109 end
0110 
0111 if pt > 1
0112     fig_conv = show_main_convergence_0(xa);
0113 end
0114 
0115 %% start iterations =============================================
0116 est_sigma_0 = 1;
0117 max_ratio   = 10;
0118 
0119 for iter=1:max_iter
0120     
0121     if pt > 0
0122         disp(['Iteration No.: ',num2str(iter), ' -----------------------------'])
0123     end
0124     %% find observation positions
0125     if pt > 1
0126         tic;
0127     end
0128     if max_ratio > 0.5
0129         [Nf,yi,zi,w] = LSM_62_sym_warp_find_observation_positions(...
0130             g,h,vg,vh,xa,Nf_MIN,NO_TEST_CPU*pt);
0131     end
0132     if Nf < 5
0133         disp('Overlap is too small')
0134         est_x = 0;
0135         est_sigma_0 = 0;
0136         Red = -1;
0137         NN = 0;
0138         Mf = Nf;
0139         N_iter = iter;
0140         max_ratio = 100;
0141         return
0142     end
0143     
0144     %% determine estimated fe and smoothed fes
0145     
0146     [fe,fes,fg,fh] = LSM_62_sym_warp_estimate_f(...
0147         g,h,vg,vh,Nf,xa,sigma_smooth,NO_TEST_CPU*pt);
0148     if pt > 1
0149         cpu=toc;
0150         sum_cpu=sum_cpu+cpu;
0151     end
0152     % adapt smoothing parameter
0153     sigma_smooth=sigma_smooth*0.5;
0154     
0155     if pt > 0
0156         im_lr = fill_left_right_images(...
0157             im_lr,fg,fh,weite,iter);
0158     end
0159     
0160     %% estimate parameters
0161     old_est_sigma_0=est_sigma_0;
0162     if pt > 1
0163         old_cpu=cpu;
0164         tic;
0165     end
0166     
0167     % --------------  main routine  -------------------
0168     [est_xa,d_xa,est_sigma_0,NN,Red,Xm,dg]  = ...
0169         LSM_62_sym_warp_estimate_parameters(...
0170         fe,fes,g,h,yi,zi,w,xa,NO_TEST_CPU*pt);
0171     % --------------  main routine  -------------------
0172     
0173     if pt > 1
0174         cpu=cpu+toc;
0175         old_cpu=cpu;
0176         sum_cpu=sum_cpu+cpu;
0177     end
0178     
0179     if pt > 1
0180         show_main_convergence(...
0181             fig_conv,xa,est_xa,old_est_sigma_0,est_sigma_0,iter,old_cpu,cpu,sum_cpu);
0182     end
0183     xa = est_xa.x;
0184     
0185     %
0186     %% close iteration
0187     
0188     if pt > 1
0189         check=(Xm'*diag(w)*dg);
0190         check_AtWv=check'
0191     end
0192     max_ratio = max(abs(d_xa./sqrt(diag(est_xa.C))));
0193     %% plot/print
0194     if pt > 0
0195         display(['N = ',num2str(NN),', sigma_0 = ', num2str(est_sigma_0),', ratio = ',num2str(max_ratio)])
0196     end
0197    
0198     if  max_ratio < 0.1
0199         break
0200     end
0201     
0202 end
0203 N_iter = iter;
0204 Mf = size(fe,1);
0205 
0206 
0207 if pt > 0
0208     figure('name','undistorted images')
0209     imshow(im_lr)
0210     
0211     for iter = 1:max_iter
0212         text(border+(iter+0.5)*weite,1.75*border+weite,num2str(iter))
0213     end
0214 end
0215 %% Finale
0216 est_x = finalize_sym_lsm(...
0217     iter,...
0218     xa,...
0219     est_xa,...
0220     pt);
0221 
0222 return

Generated on Sun 19-Jul-2020 23:00:25 by m2html © 2005