close all
clear all

cropfactor=0.70;

fname='C:\temp\Jupiter Animation'
[frames.im frames.map]=imread([fname '.gif'],'frames','all');
width_orig=size(frames.im,1);
height_orig=size(frames.im,2);



%reference image
if (true)
    ref_ima_nr=23;
    imaRGB=ind2rgb(frames.im(:,:,1,ref_ima_nr),frames.map);
    image=imresize(imaRGB,2,'bicubic');
    ima=squeeze(image);
    ima=ima/max(max(max(ima)));
    
    width=size(ima,2);
    height=size(ima,1);
    
    %clear everything outside planet
    if (false)
        for x=1:width
            for y=1:height
                xx=x-width/2;
                yy=y-height/2;
                
                if (y<0.10*height) | (y>0.90*height) 
%                if sqrt(xx*xx+yy*yy)>0.8*width/2 
                    ima(y,x,:)=0;
                end
            end
        end
    end
    
    
    widthrange=1:width;
    heightrange=(1:height).';
    
    
    imaBW=squeeze(sum(ima,3));
    imaBW_thresholded=(imaBW>(max(max(max(imaBW)))/8)); 

    imagesc(imaBW_thresholded);
    colormap(gray)
    drawnow
    sumsum=sum(sum(imaBW_thresholded));

    centreofmass_ver0=sum(heightrange.*sum(imaBW_thresholded,2))/sumsum
    centreofmass_hor0=sum(widthrange.*sum(imaBW_thresholded,1))/sumsum

    %clear out everything out of com
    if (false)
        for x=1:width
            for y=1:height
                xx=x-width/2;
                yy=y-height/2;
                if sqrt(xx*xx+yy*yy)>0.8*width/2 
                    ima(x,y,:)=0;
                end
            end
        end
    end

    ima0=ima;
end
    

for i=1:size(frames.im,4)-1

  
    %convert lookup 8bit to RGB
    
    %compressed gif? zeros are equal to pixel(x-1)
    if (false)
        for x=1:width_orig
            for y=1:height_orig
                if frames.im(x,y,1,i)==0
%                    frames.im(x,y,1,i)=frames.im(max(1,x-1),y,1,i);
                    frames.im(x,y,1,i)=frames.im(x,max(1,y-1),1,i);
                end
            end
        end
    end

    imaRGB=ind2rgb(frames.im(:,:,1,i),frames.map);
    image=imresize(imaRGB,2,'bicubic');
    ima=squeeze(image);
    ima=ima/max(max(max(ima)));
    
    width=size(ima,2);
    height=size(ima,1);
    
    %clear everything outside planet
    if (true)
        for x=1:width
            for y=1:height
                xx=x-width/2;
                yy=y-height/2;
                
                if (y<0.10*height) | (y>0.90*height) 
%                if sqrt(xx*xx+yy*yy)>0.8*width/2 
                    ima(y,x,:)=0;
                end
            end
        end
    end
    
    
    widthrange=1:width;
    heightrange=(1:height).';
    
  
        
    if (false)
        for x=1:width
            for y=1:height
                if sqrt((x-centreofmass_hor0)^2+(y-centreofmass_ver0)^2)>(height/2)*cropfactor
                    ima(y,x,:)=0;
                end
            end
        end
    end

    imaBW=squeeze(sum(ima,3));
    imaBW_thresholded=(imaBW>(max(max(max(imaBW)))/8)); 

    sumsum=sum(sum(imaBW_thresholded));
    imagesc(imaBW_thresholded);
    colormap(gray)
    drawnow
    centreofmass_ver=sum(heightrange.*sum(imaBW_thresholded,2))/sumsum;
    centreofmass_hor=sum(widthrange.*sum(imaBW_thresholded,1))/sumsum;

    [-round(centreofmass_ver0-centreofmass_ver) -round(centreofmass_hor0-centreofmass_hor)]  
    for layer=1:3
        %shift back image
        a=round(centreofmass_ver0-centreofmass_ver);
        b=round(centreofmass_hor0-centreofmass_hor);
        ima(:,:,layer)=circshift(ima(:,:,layer),[a b]);

        if (true)
            %correct colors
            params=polyfit(ima(:,:,layer).*imaBW_thresholded,ima0(:,:,layer).*imaBW_thresholded,1);
            ima(:,:,layer)=polyval(params,ima(:,:,layer));
        end
    end

    ima=max(ima,0);
    ima=min(ima,1);
    ima=imresize(ima,0.5,'bilinear');
    f=imagesc(ima);
    axis equal
    axis tight
    drawnow
    
    imagedata=get(f,'CData');
    if (i==1) [im map]= rgb2ind(imagedata,256,'dither'); end
    framesout(:,:,1,i) = rgb2ind(imagedata,map);
end


imwrite(framesout,map,[fname '_improved.gif'],'DelayTime',0.10,'LoopCount',inf) %g443800
