Matlab code
% 1. make sure the figure drawing of deformation vector field is right. % 2. regenerate data. function motion_fusion() path(path,'Users/cirrus_mdn/Desktop/Desktop - Madans Macbook Pro/Research/Motion_Fusion10102009'); %== Parameter set ================================================= m=513; n=513; nbin = 256; idx_i = 25; dataname = 'R'; % Multiresolution starting level Multi_SL = 4; % level 1: knot space hx, hy == 2^0 == 1 ite_max_h = 10; load T; load R; T = double(T); R = double(R); % Normalize T, R, making SSD of T, R = 2000 norm_factor = sqrt(2000/sum((R(:)-T(:)).^2)); T = T * norm_factor; R = R * norm_factor; % idx_i: deformation index for generating groundtruth SNR_dB = 100; interp_f = 1; deform_case = 1; % Using deformation method case 1, 2 % Solver method: FCD: 0, LSFEM: 1 solver_method = 0; bd_cond = 0; % 0: dirichlet, 1: neumann % LSFEM grid size*2+1: 0, LSFEM full size*2+1: 1, 17x17: 2 solver_sz = 1; regrid = 0; % Fix Jacobian fix_f1 = 0; % Euler:0, Mid: 1, RK4: 2 RK_method = 0; RK_steps = 1; % f1, f2 dt, seperately: 1 simultaneously: 0 dt_f1f2 = 0; groundT = 0; tstep = 10^(0); % load([filename(1:end-2),'.mat']); load('./results/case1_2d_R_interpf1_bdcond0_SNRdB100_regrid0_fixJ0_solver_method0_RKmethod0_RKsteps1_ite_400_h8-h4.mat') %====================================================================== % Motion Fusion %====================================================================== ssd_old = sum((F(:) - R_resampled(:)).^2); alpha = ones(size(F)); %alpha = double(F > 0.005); tstep = 10^-1; succeed = 1; [Id_X, Id_Y] = ndgrid(1:m, 1:n); Phi_X = DISP_X + Id_X; Phi_Y = DISP_Y + Id_Y; % T->R, R->F while tstep > 10^(-5) if succeed % R_resampled R_resampled = interp2c(Phi_X, Phi_Y, R, ... uint32([m,n]),uint32([m,n])); end [R_x2, R_x1] = gradient(R_resampled); da = 2*(R_resampled-F).*(R_x1.*DISP_X + R_x2.*DISP_Y); alpha_tmp = alpha + tstep*da; DISP_X_tmp = DISP_X.*alpha_tmp; DISP_Y_tmp = DISP_Y.*alpha_tmp; Phi_X_tmp = DISP_X + Id_X; Phi_Y_tmp = DISP_Y + Id_Y; R_resampled_fusion_tmp ... = interp2c(Phi_X_tmp, Phi_Y_tmp, R, uint32([m,n]),uint32([m,n])); ssd_tmp = sum((F(:)-R_resampled_fusion_tmp(:)).^2); if ssd_tmp <= ssd_old Phi_Y = Phi_Y_tmp; Phi_X = Phi_X_tmp; DISP_X = DISP_X_tmp; DISP_Y = DISP_Y_tmp; alpha = alpha_tmp; ssd_old = ssd_tmp; tstep = tstep*1.1; R_resampled = R_resampled_fusion_tmp; fprintf('ssd:%5.4f tstep:%4.4e\n',... ssd_old, tstep); succeed = 1; else tstep = tstep*0.85; succeed = 0; end end show_figures(R_resampled, Phi_X, Phi_Y, Id_X, Id_Y); end function show_figures(R_resampled, Phi_X, Phi_Y, Id_X, Id_Y) close all; load truth_xy; figure(1); imshow((R_resampled),[]); hold on; deformplot2_griddata(Phi_X, Phi_Y, ... Id_X, Id_Y, 5, 1); deformplot2(true_x+Id_X, true_y+Id_Y, ... Id_X, Id_Y, 5, 1); hold off end