% filename: stabilizeMovie_GCBPM.m % author: Alan Brooks % description: Given gray scale double matrix, M(width,height,frameNum), % this function stabilizes the sequence using the Grey-Coded Bit-Plane % Matching (GC-BPM) algorithm. % version history: % 08-Mar-2003 created % 09-Mar-2003 tested and optimized % 10-Mar-2003 output more vars % 11-Mar-2003 added NSS log algorithm % (improve from 0.8 FPS to 3 FPS on 1GHz PIII) function [Ms,Va,Vg,V] = stabilizeMovie_GCBPM(M) % init global variables debug_disp = 0; % display debugging plots % Gray Coded Approach based on: % [8] S. Ko, S. Lee, S. Jeon, and E. Kang. Fast digital image stabilizer % based on gray-coded bit-plane matching. IEEE Transactions on % Consumer Electronics, vol. 45, no. 3, pp. 598-603, Aug. 1999. % -------------------------------------------------------------------------- % init algorithm specific variables (uint8's for speed) bit = 5; % best grey code bit (det by paper) N = 112; %64; % matching block size (assumed NxN square) D1 = 0.95; % damping coeff (0 < D1 < 1) logSearchEnable = 1; % 1 = use logrigtmic 3-2-1 matching search % 0 = use original exhaustive matching search nSteps = 3; % number of steps for log search % (3 is pixel, 2 is 2 pixel, etc...) rotEnable = 0; % 0 = use original translational GCBPM % 1 = enable matrix-based rotation correction (Poor performance) % 2 = enable simple/efficient rot. correction [h,w,nFr] = size(M); % height, width, # of frames if ( ~rem(w,2) & ~rem(h,2) ) S = uint8(zeros(h/2,w/2,2,4)); % subimage holder (must be even) else error('video width/height must both be even # of pixels') end hw = waitbar(0,'Please wait...'); p = (h/2-N)/2; % max search window displacement bxor = uint8(zeros(N)); % intermediate value Cj = 1e9*ones(2*p+1); % correlation measure V = zeros(4,2,nFr); % motion vectors Vg = zeros(nFr,2); %[0 0]; % global motion vector Va = zeros(nFr,2); % integrated motion vectors Ms = uint8(zeros(h,w,nFr)); % init stabilized image sequence % loop over subsequent frames for fr = 1:nFr waitbar((fr-1)/nFr,hw) % display progress % get grey coded bit plane [Mg] = uint8(getGrayCodeBitPlane(M,bit,fr,debug_disp)); % break into 4 sub-images, (S1, S2, S3, S4) % SX(:,:,2,1) is S1 current image, SX(:,:,1,1) is past % could probably optimize as a reshape S(:,:,2,1) = Mg( 1:h/2, 1:w/2 ); % UL, S1 S(:,:,2,2) = Mg( 1:h/2, w/2+1:end ); % UR, S2 S(:,:,2,3) = Mg( h/2+1:end, 1:w/2 ); % LL, S3 S(:,:,2,4) = Mg( h/2+1:end, w/2+1:end ); % LR, S4 if fr > 1 % start algorithm after 1st frame for j = 1:4 % loop over sub-images if ~logSearchEnable % exhaustive search (slower, more precision) for m_pos = 1:2*p+1 % loop over possible displacements for n_pos = 1:2*p+1 % calculate correlation measure bxor = bitxor( ... % could be very fast HW S(p+1:p+N,p+1:p+N,2,j) , ... S(m_pos:m_pos+N-1,n_pos:n_pos+N-1,1,j) ); Cj(m_pos,n_pos) = sum(bxor(:)); % could mult by 1/N^2 end end % displacement loops % find min Cj arguments: m_pos_min & n_pos_min [tmp,m_pos_min] = min(Cj); [tmp,n_pos_min] = min(tmp); clear tmp; m_pos_min=m_pos_min(n_pos_min); else % log search (faster) % currently works only for 256x256 images firstJmp = 4; prev_m_pos = 9; prev_n_pos = 9; % start in center for iter = 1:nSteps % iter'th step Cj = 1e9*ones(2*p+1); % reset correlation measure curJmp = firstJmp./2.^(iter-1); for m_pos = prev_m_pos-curJmp:curJmp:prev_m_pos+curJmp for n_pos = prev_n_pos-curJmp:curJmp:prev_n_pos+curJmp % calculate correlation measure bxor = bitxor( ... % could be very fast HW S(p+1:p+N,p+1:p+N,2,j) , ... S(m_pos:m_pos+N-1,n_pos:n_pos+N-1,1,j) ); Cj(m_pos,n_pos) = sum(bxor(:)); end end % find min Cj arguments: m_pos_min & n_pos_min [tmp,m_pos_min] = min(Cj); [tmp,n_pos_min] = min(tmp); clear tmp; m_pos_min=m_pos_min(n_pos_min); % store for loop prev_m_pos = m_pos_min; prev_n_pos = n_pos_min; end % iteration loop end % log search if % store motion vector for this sub-image V(j,:,fr) = [m_pos_min n_pos_min]-p-1; % V[1] V[2] end % sub-image loop % calc current global motion vector Vg(fr,:) = median([V(:,:,fr);Vg_prev]); % apply damping to generate this frame's integrated motion vector Va(fr,:) = D1 * Va_prev + Vg(fr,:); end % fr>1 if % store current image as past image S(:,:,1,:) = S(:,:,2,:); % gray coded sub-images Vg_prev = Vg(fr,:); % global motion vector Va_prev = Va(fr,:); % integrated motion vector % stabilize the image with the current motion vector %Ms(:,:,fr) = M(:,:,fr); % add Xform !!! ??? switch rotEnable case 0 % original translational correction % (not sub-pixel for now) r = round(Va(fr,1)); % num rows moved c = round(Va(fr,2)); % num columns moved Ms(max([1 1+r]):min([h h+r]),max([1 1+c]):min([w w+c]),fr) = ... M(max([1 1-r]):min([h h-r]),max([1 1-c]):min([w w-c]),fr); case 1 % correct for rotational & translational motion: matrix method % POOR PERFORMANCE %B = IMROTATE(A,ANGLE,METHOD,'crop') cnst = 12; %128; theta = zeros(1,20); Mrs = uint8(zeros(size(M))); Mrs(:,:,1) = M(:,:,1); for fr=2:20 B=[V(:,:,fr) + cnst*[-1 1; 1 1; -1 -1; 1 -1]]';B=B(:); A=zeros(2*4,4); A(1:2:end,1)=V(:,1,fr-1) + cnst*[-1 1 -1 1]'; %1st col A(2:2:end,1)=V(:,2,fr-1) + cnst*[1 1 -1 -1]'; A(2:2:end,2)=V(:,1,fr-1) + cnst*[-1 1 -1 1]'; %2nd col A(1:2:end,2)=V(:,2,fr-1) - cnst*[1 1 -1 -1]'; A(1:2:end,3)=1; %3rd col A(2:2:end,4)=1; %4th col X=A\B theta(fr)=atan2(X(1),X(2))*180/pi-90; theta(fr) %err=V(1,:,fr)-X(3:4)' Mrs(:,:,fr) = ... imrotate(M(:,:,fr),sum(theta(1:fr)),'bilinear','crop'); figure,imshow(Mrs(:,:,fr)) end case 2 % correct for rotational: simple/efficient rot. correction end end % frame loop close(hw) % close progress bar return