[Registax] Wavelet Processing from Stacked Images

close all
clear all

moviefile='jupiter1.avi';

movieinfo=aviinfo(moviefile);
numframes=movieinfo.NumFrames;
mov=aviread(moviefile,1); % open movie

ref_image=double(mov.cdata(:,:,1));% user red layer to register
ref_image_color=double(mov.cdata);% 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)

% 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); % read frame number i

    reg_image=double(mov.cdata(:,:,1));% user red layer to register
    reg_image_color=double(mov.cdata);% user red layer to register
    reg_image=reg_image-min(min(reg_image)); %subtract offset
%    imagesc(reg_image);
%    colormap(gray)
%    pause(0.01)

%register image to the reference image
    fft_regimage=(fft2(reg_image));
    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;
   
    imagesc(sum_image/(number_images_summed*256));
    pause(0.001)
   

   % pause
end

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,'c:\temp\sumimage.jpg');
   
figure;
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
imagesc(abs(sum_image_sharp))
colormap(gray)
title('sharpened sum of images');