Registax by POC

close all
clear all

pkg load video

% load video -----------------------------
moviefile='jupiter1.avi';


movieinfo=aviinfo(moviefile);
numframes=movieinfo.NumFrames;
mov=aviread(moviefile,1); % 한프레임 얻기

% 레드채널 등 푸리에변환을 위한 더블화 ---------------
%ref_image=double(mov.cdata(:,:,1));% user red layer to register
ref_image=double(mov(:,:,1));
ref_image_color=double(mov);% user red layer to register

% 레드채널 복소공액 및 정규화 ---------------
fft_refimage=conj(fft2(ref_image)); % complex conjugate of fft of reference image
fft_refimage_norm=fft_refimage./abs(fft_refimage); %normalisation needed for "phase only" correlation

% 버퍼
xshifts=[]; %shifts to be detected in cosecutive images
yshifts=[];
imagesize=size(ref_image);
sum_image=ref_image_color; %buffer for all images summed after registration (=shift to ref image position)

% POC 피크 잡음제거를 위한 주파수도메인 로우패쓰 필터 마스크 생성-------------------
% low pass filter to reduce effect of noise on correlation peak position
lowpassfilter=ones(imagesize);
midx=imagesize(2)/2;
midy=imagesize(1)/2;
cutoff=max(imagesize)/10;
for x=1:imagesize(2)
    for y=1:imagesize(1)
        radius=sqrt((x-midx)^2+(y-midy)^2);
        if radius>cutoff
            lowpassfilter(y,x)=exp(-(radius-cutoff)/10);
        end
    end
end

%imagesc(lowpassfilter);
%pause;


% 모든 프레임에 대해서 ----------------------------------------
number_images_summed=1; % counter for averaging of the summed images

for i=2:1:numframes  
    mov=aviread(moviefile,i); % i 번쩨 프레임 로드

    reg_image=double(mov(:,:,1));% 레드 채널
    reg_image_color=double(mov);% 컬러 채널
    reg_image=reg_image-min(min(reg_image)); % 옾셋 제거
    %imagesc(reg_image);colormap(gray)
    %pause(0.01)

    %register image to the reference image
    fft_regimage=(fft2(reg_image)); % red fft2
    fft_regimage_norm=fft_regimage./abs(fft_regimage);

    fft_POC=fft_regimage_norm.*fft_refimage_norm;
   
    fft_POC=fftshift(fft_POC).*lowpassfilter;
    fft_POC=ifftshift(fft_POC);
   
    POC=ifft2(fft_POC);
    %imagesc(log(abs(POC)))
    % pause
    %colormap(gray)
    [X,xshift]=max(max(POC));
    [Y,yshift]=max(max(POC.'));
   
    if (xshift>imagesize(2)/2)
        xshift=xshift-imagesize(2);
    end
   
    if (yshift>imagesize(1)/2)
        yshift=yshift-imagesize(1);
    end
   
    xshift=xshift;
    yshift=yshift;
       
   
    xshifts=[xshifts xshift];
    yshifts=[yshifts yshift];

%   add all 3 color layers to the sum buffer
    for ypos=1:imagesize(1)
       for xpos=1:imagesize(2)
            if (((xpos+xshift)>0)&&((xpos+xshift)<imagesize(2))&&((ypos+yshift)>0)&&((ypos+yshift)<imagesize(1)))
                for layer=1:3
                    sum_image(ypos,xpos,layer)=sum_image(ypos,xpos,layer)+reg_image_color(ypos+yshift,xpos+xshift,layer);
                end
            end
       end
    end
   
    number_images_summed=number_images_summed +1;
   
    pause(0.001)
    number_images_summed

    %imagesc(sum_image/(number_images_summed*256));
    %pause(0.001)
   

   % pause
end % for all frames -----------

sum_image=sum_image/(number_images_summed*256);

sum_image=sum_image/max(max(max(sum_image))); % must be normalised to 1 before imwrite
imwrite(sum_image,'sumimage.jpg');
   
figure(1);plot(xshifts,yshifts);title('Shifts detected')

%high pass filter the sum of images to enhance details
highpassfilter=ones(imagesize);
midx=imagesize(2)/2;
midy=imagesize(1)/2;
cutoff=max(imagesize)/10;
for x=1:imagesize(2)
    for y=1:imagesize(1)
        radius=sqrt((x-midx)^2+(y-midy)^2);
        if radius>cutoff
            highpassfilter(y,x)=highpassfilter(y,x)+(radius-cutoff)/50;
        end
    end
end

% sharpen all three color layers
sum_image_sharp=sum_image;
for layer=1:3
    fft_sum_image=fftshift(fft2(sum_image(:,:,layer)));
    sum_image_sharp(:,:,layer)=ifft2(ifftshift(fft_sum_image.*highpassfilter));
end


%show enhanced sum image
sum_image_sharp=sum_image_sharp/max(max(max(sum_image_sharp)));
figure(2);imagesc(abs(sum_image_sharp));colormap(gray)
title('sharpened sum of images');