input = imread('a1_input.jpg');
input_size = size(input);
height = floor(input_size(1)/3);
width = input_size(2);

orig_B = input(1:height,:);
orig_G = input(height+1:2*height,:);
orig_R = input(2*height+1:3*height,:);

range = 15;
margin = 20;

angle_range = 10;


fraction_factor = 2;

height = height * fraction_factor;
width = width * fraction_factor;


R = zeros(height,width);
G = zeros(height,width);
B = zeros(height,width);

for i=1:height
    for j=1:width
        R(i,j) = orig_R(ceil(i/fraction_factor),ceil(j/fraction_factor));
        G(i,j) = orig_G(ceil(i/fraction_factor),ceil(j/fraction_factor));
        B(i,j) = orig_B(ceil(i/fraction_factor),ceil(j/fraction_factor));
    end
end

template = R(margin+1:height-margin,margin+1:width-margin);

error_G = zeros(2*range+1,2*range+1,2*angle_range+1);
temp_G = zeros(height-2*margin,width-2*margin);
for theta=-angle_range:angle_range
    rotate_G = imrotate(G,theta,'bilinear','crop');
    for x=-range:range
        for y=-range:range
            temp_G(1:height-2*margin,1:width-2*margin) = rotate_G(margin+1+x:height-margin+x,margin+1+y:width-margin+y);
            error_G(range+1+x,range+1+y,angle_range+1+theta)=sum(sum((template-temp_G).^2));
        end
    end
end

error_B = zeros(2*range+1,2*range+1,2*angle_range+1);
temp_B = zeros(height-2*margin,width-2*margin);
for theta=-angle_range:angle_range
    rotate_B = imrotate(B,theta,'bilinear','crop');
    for x=-range:range
        for y=-range:range
            temp_B(1:height-2*margin,1:width-2*margin) = rotate_B(margin+1+x:height-margin+x,margin+1+y:width-margin+y);
            error_B(range+1+x,range+1+y,angle_range+1+theta)=sum(sum((template-temp_B).^2));
        end
    end
end

[error_G3 idx3] = min(error_G, [], 3);
[error_G2 idx2] = min(error_G3, [], 2);
[error_G1 idx1] = min(error_G2, [], 1);
idx1;
idx2=idx2(idx1);
idx3=idx3(idx1,idx2);
G_idx1 = idx1-range-1;
G_idx2 = idx2-range-1;
G_idx3 = idx3-angle_range-1;

[error_B3 idx3] = min(error_B, [], 3);
[error_B2 idx2] = min(error_B3, [], 2);
[error_B1 idx1] = min(error_B2, [], 1);
idx1;
idx2=idx2(idx1);
idx3=idx3(idx1,idx2);
B_idx1 = idx1-range-1;
B_idx2 = idx2-range-1;
B_idx3 = idx3-angle_range-1;

G = imrotate(G,G_idx3,'bilinear','crop');
B = imrotate(B,B_idx3,'bilinear','crop');

output = zeros(height-2*margin,width-2*margin,3);
output(:,:,1) = R(margin+1:height-margin,margin+1:width-margin);
output(:,:,2) = G(margin+1+G_idx1:height-margin+G_idx1,margin+1+G_idx2:width-margin+G_idx2);
output(:,:,3) = B(margin+1+B_idx1:height-margin+B_idx1,margin+1+B_idx2:width-margin+B_idx2);
% output=uint8(output);

process_size = size(output);
final_height = floor(process_size(1) / fraction_factor);
final_width = floor(process_size(2) / fraction_factor);

final_output = zeros(final_height,final_width,3);

for k=1:3
    for i=1:final_height
        for j=1:final_width
            final_output(i,j,k) = sum(sum(output((i-1)*fraction_factor+1:i*fraction_factor, ...
                (j-1)*fraction_factor+1:j*fraction_factor,k)))/(fraction_factor*fraction_factor);
        end
    end
end

final_output = uint8(final_output);

imshow(final_output);


    

