7941943TPS-RPM.zip_TPS matlab_TPS RPM_plate_rpm_tps

% Robust Point Matching (RPM) Demo (version 20000427): % ---------------------------------------------------- % Copyright (C) 2000 Haili Chui, Anand Rangarajan % % Authors: Haili Chui and Anand Rangarajan % Date: 04/27/2000 % % Contact Information: % % Haili Chui: chui@noodle.med.yale.edu % Anand Rangarajan: anand@noodle.med.yale.edu % % Terms: % % The source code (M-files) are provided under the % terms of the GNU General Public License with an explicit % clause permitting the execution of the M-files from within % a MATLAB environment. See the LICENSE file for details. % % % ---------------------------------------------------------------------- % cMIX.m % ---------------------------------------------------------------------- % Gaussian Mixture Robust Point Matching. % % Usage: % % [c,d,m] = cMIX (x, y, frac, Tinit, Tfinalfac); % % Optional ones: % % [c,d] = cMIX (x, y, frac, Tinit, Tfinalfac); % [c,d] = cMIX (x, y, frac, Tinit, Tfinalfac, disp_flag); % [c,d] = cMIX (x, y, frac, Tinit, Tfinalfac, disp_flag, m_method); % % [w] = cMIX (x, y, z, sigma, frac, Tinit, Tfinalfac); % [w] = cMIX (x, y, z, sigma, frac, Tinit, Tfinalfac, disp_flag); % [w] = cMIX (x, y, z, sigma, frac, Tinit, Tfinalfac, disp_flag, m_method, % icp_sigma); % % Notes: the program will set transformation type automatically. % ['tps', 'rbf']. % % 'icp' -- 'icp0', 'icp3', 'icp5'. --> to set k_sigma. % % 04/27/00 function [o1,o2,o3] = cMIX (in1,in2,in3,in4,in5,in6,in7,in8,in9); fig(1); clf; whitebg('k'); set(gcf,'color',[0 0 0]); % set(gcf,'DoubleBuffer','on') % Init control parameters: % ------------------------ perT_maxit = 5; relax_maxit = 1; anneal_rate = 0.93; lamda1_init = 1; lamda2_init = 0.01; % default: % -------- disp_flag = 1; m_method = 'mixture'; debug_flag = 0; % check input: % ------------ % --- [c,d] = cMIX (x, y, frac, Tinit, Tfinalfac) --------------------------- if (nargin == 5) x = in1; y = in2; frac = in3; T_init = in4; T_finalfac = in5; trans_type = 'tps'; z = x; sigma = 1; % --- [c,d] = cMIX (x, y, frac, Tinit, Tfinalfac, disp_flag) -------------- elseif (nargin == 6) & (length(in3) == 1) x = in1; y = in2; frac = in3; T_init = in4; T_finalfac = in5; disp_flag = in6; trans_type = 'tps'; z = x; sigma = 1; % --- [c,d] = cMIX (x, y, frac, Tinit, Tfinalfac, disp_flag, m_method) ------ elseif (nargin == 7) & (length(in3) == 1) x = in1; y = in2; frac = in3; T_init = in4; T_finalfac = in5; disp_flag = in6; m_method = in7; trans_type = 'tps'; z = x; sigma = 1; % --- [w] = cMIX (x, y, z, sigma, frac, Tinit, Tfinalfac) ----------------- elseif (nargin == 7) & (length(in3) > 1) x = in1; y = in2; z = in3; sigma = in4; frac = in5; T_init = in6; T_finalfac = in7; trans_type = 'rbf'; disp_flag = 1; % --- [w] = cMIX (x, y, z, sigma, frac, Tinit, Tfinalfac, disp_flag) ----- elseif (nargin == 8) & (length(in3) > 1) x = in1; y = in2; z = in3; sigma = in4; frac = in5; T_init = in6; T_finalfac = in7; disp_flag = in8; trans_type = 'rbf'; % --- [w] = cMIX (x,y,z,sigma,frac,Tinit,Tfinalfac,disp_flag,m_method); --- elseif (nargin == 9) & (length(in3) > 1) x = in1; y = in2; z = in3; sigma = in4; frac = in5; T_init = in6; T_finalfac = in7; disp_flag = in8; m_method = in9; trans_type = 'rbf'; % --- [theta,tx,ty] = cMIX (x, y, frac, Tinit, Tfinalfac) ---------------------- elseif (nargin == 5) & (nargout == 3) x = in1; y = in2; frac = in3; T_init = in4; T_finalfac = in5; trans_type = 'r+t'; theta = 0; t = zeros (2,1); s = 1; disp_flag = 1; z = x; sigma = 1; % --- [theta,tx,ty] = cMIX (x, y, frac, Tinit, Tfinalfac, disp_flag) --------- elseif (nargin == 6) & (nargout == 3) x = in1; y = in2; frac = in3; T_init = in4; T_finalfac = in5; disp_flag = in6; trans_type = 'r+t'; theta = 0; t = zeros (2,1); s = 1; z = x; sigma = 1; % --- [theta,tx,ty] = cMIX (x, y, frac, Tinit, Tfinalfac,disp_flag, m_method); --- elseif (nargin == 7) & (nargout == 3) x = in1; y = in2; frac = in3; T_init = in4; T_finalfac = in5; disp_flag = in6; m_method = in7; trans_type = 'r+t'; theta = 0; t = zeros (2,1); s = 1; z = x; sigma = 1; else disp ('# ERROR #: cMIX -- wrong input!'); help cMIX; return; end; % take care of 'icp' k_sigma stuff. if (strcmp(m_method(1:3), 'icp')) if length(m_method) == 3 k_sigma = 0; m_method = 'icp'; else k_sigma = str2num(m_method(4)); m_method = 'icp'; end; else k_sigma = 0; end; % Init: % ----- % init x,y,z: [xmax, dim] = size(x); x = x (1:frac:xmax, :); [xmax, dim] = size(x); [ymax, tmp] = size(y); y = y (1:frac:ymax, :); [ymax, tmp] = size(y); [zmax, tmp] = size(z); if strcmp(trans_type, 'tps') z = x; end; % init m: m = ones (xmax, ymax) ./ (xmax * ymax); T0 = max(x(:,1))^2; moutlier = 1/sqrt(T0)*exp(-1); % /xmax *0.001; % moutlier = 1/xmax*0.01; m_outliers_row = ones (1,ymax) * moutlier; m_outliers_col = ones (xmax,1) * moutlier; % init transformation parameters: theta = 0; t = zeros (2,1); s = 1; c_tps = zeros (xmax,dim+1); d_tps = eye (dim+1, dim+1); w = zeros (xmax+dim+1, dim); % icp: perT_maxit = 1: if strcmp (m_method(1:3),'icp') perT_maxit = 1; end; % ------------------------------------------------------------------- % Annealing procedure: % ------------------------------------------------------------------- T = T_init; T_final = T_init / T_finalfac; vx = x; vy = y; it_total = 1; flag_stop = 0; while (flag_stop ~= 1) for i=1:perT_maxit % repeat at each termperature. % Given vx, y, Update m: if debug_flag; disp ('calc m ...'); end; m = cMIX_calc_m (vx, y, T, m_method, ... m_outliers_row, m_outliers_col,it_total,k_sigma); % Given m, update transformation: vy = m * y ./ ( (sum(m'))' * ones(1,dim)); lamda1 = lamda1_init*length(x)*T; lamda2 = lamda2_init*length(x)*T; %lamda1 = lamda1_init; %lamda2 = lamda2_init; if debug_flag; disp ('calc c,d ...'); end; [c_tps, d_tps, w] = cMIX_calc_transformation (trans_type, ... lamda1, lamda2, sigma, x, vy, z); % w(1:length(z),:) = w(1:length(z),:)*0 % d_tps = d_tps * 0 + eye(dim+1,dim+1); % c_tps = c_tps *0; if debug_flag; disp ('calc new x ...'); end; [vx] = cMIX_warp_pts (trans_type, x, z, c_tps, d_tps, w, sigma); end % end of iteration/perT T = T * anneal_rate; % Determine if it's time to stop: % ------------------------------- % when T <= cluster variance, stop. if T < T_final; flag_stop = 1; end; str = sprintf ('%s, T = %.4f:\t lamda1: %.4f lamda2: %.4f', ... m_method,T, lamda1, lamda2); disp (str); % Display: % -------- it_total = it_total + 1; if (disp_flag == 1) | (flag_stop == 1) | ... (disp_flag ==2 & mod(it_total,10)==0) | ... (it_total == 1) figure(1); clf; cMIX_plot_simple (2, x, y, z, vx, m, 1/xmax, T, ... trans_type, c_tps, d_tps, w, sigma, m_method); pause(0.5); end; end % end of annealing. % return outputs: % --------------- if strcmp (trans_type, 'tps') o1 = c_tps; o2 = d_tps; o3 = []; elseif strcmp (trans_type, 'rbf') o1 = w; o2 = []; o3 = []; elseif strcmp (trans_type, 'r+t') o1 = theta; o2 = tx; o3 = ty; end; o3 = m; %%%%% % 1 % %%% cMIX_calc_m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% % % Update m (correspondence). % % Usage: % [m] = cMIX_cal