How to homogeneous two different images into same color after image stitching

%% The above image shows different colors once it stitches. I did use adobe lightroom and the color of the image is homogeneous.

%% The below code was how I stitched my two images into one image.

%% read images

imgPath = 'C:\Users\nguic\Documents\MATLAB\test_png\';

imgDir = dir([imgPath,'*.png']);

% select reference

reference_order = 1;

ref = imread([imgPath imgDir(reference_order).name]);

% convert ref into gray scale

ref_g = im2gray(ref);

% compute the all possible feature points

ref_points = detectSIFTFeatures(ref_g);

% select souce image

source_order = [2];

panel=zeros(10000,10000,3); % adjustable

% mannually translate, ensuring the all locations are positive

final_H_for_ref = [1 0 3000;0 1 4000;0 0 1];

%create a location list for ref

for ii = 1:size(ref_g,2)

if ii == 1

loc_list_ref_to = transpose([ones(size(ref_g,1),1) (1:size(ref_g,1))']);

else

loc_list_ref = transpose([ii.*ones(size(ref_g,1),1) (1:size(ref_g,1))']);

loc_list_ref_to = [loc_list_ref_to loc_list_ref];

end

end

loc_list_ref_hc = [loc_list_ref_to ; ones(1,size(loc_list_ref_to,2))];

% get the transformed location of ref

loc_list_ref_transed_hc = final_H_for_ref*loc_list_ref_hc;

ref_col_right_lim = max(loc_list_ref_transed_hc(1,:));

ref_col_left_lim= min(loc_list_ref_transed_hc(1,:));

ref_row_right_lim= max(loc_list_ref_transed_hc(2,:));

ref_row_left_lim= min(loc_list_ref_transed_hc(2,:));

panel(ref_row_left_lim:ref_row_right_lim,ref_col_left_lim:ref_col_right_lim,:) = ref(:,:,:);

for aa = 1:length(source_order)

source = imread([imgPath imgDir(source_order(aa)).name]);

% convert source into gray scale

source_g = im2gray(source);

% compute the all possible feature points

source_points = detectSIFTFeatures(source_g);

%source_points = detectHarrisFeatures(source_g);

% create feature discroptor

[features1,valid_points1] = extractFeatures(ref_g,ref_points);

[features2,valid_points2] = extractFeatures(source_g,source_points);

% match feature points

indexPairs = matchFeatures(features1,features2);

% RANSACE to select inliers

iteration_time = 100; %% you can also adjust this itereation time

for ii1=1:iteration_time

r = randi([1 length(indexPairs)],1,4);

for ii2=1:length(r)

pair(ii2,:) = indexPairs(r(ii2),:);

% extract the pixel location (x,y) in image space

point_ref_location(ii2,:) = valid_points1.Location(pair(ii2,1),:);

point_sou_location(ii2,:) = valid_points2.Location(pair(ii2,2),:);

% convert (x,y) into (x,y,1)

point_ref_location_hc(ii2,:) = [point_ref_location(ii2,:) 1];

point_sou_location_hc(ii2,:) = [point_sou_location(ii2,:) 1];

% normalize pixel location

Norm_matrix = [2/size(ref,2) 0 -1; 0 2/size(ref,1) -1; 0 0 1];

%Norm_matrix = [1 0 0; 0 1 0; 0 0 1];

point_ref_location_hc_nor(ii2,:)=Norm_matrix*point_ref_location_hc(ii2,:)';

point_sou_location_hc_nor(ii2,:)=Norm_matrix*point_sou_location_hc(ii2,:)';

% construct Ai(a) matrix

a(:,:,ii2) = [zeros(1,3) -point_sou_location_hc_nor(ii2,:) point_ref_location_hc_nor(ii2,2)*point_sou_location_hc_nor(ii2,:);

point_sou_location_hc_nor(ii2,:) zeros(1,3) -point_ref_location_hc_nor(ii2,1)*point_sou_location_hc_nor(ii2,:);

-point_ref_location_hc_nor(ii2,2)*point_sou_location_hc_nor(ii2,:) point_ref_location_hc_nor(ii2,1)*point_sou_location_hc_nor(ii2,:) zeros(1,3)];

end

% construct A matrix and SVD to A

A = [a(:,:,1);a(:,:,2);a(:,:,3);a(:,:,4)];

[U,S,V] = svd(A);

% extract the last column of the V and resize to H matrix

h_vector = V(:,9);

H = [h_vector(1:3)';

h_vector(4:6)';

h_vector(7:9)'];

% H_final = H * Norm_matrix;

% set up the threshold to identify inliers

% transfer all original pixel location of the correspondance point

matchedPoints1 = valid_points1.Location(indexPairs(:,1),:);

matchedPoints2 = valid_points2.Location(indexPairs(:,2),:);

matchedPoints1_hc_nor= Norm_matrix*transpose([matchedPoints1 ones(size(matchedPoints1,1),1)]);

matchedPoints2_hc_nor= Norm_matrix*transpose([matchedPoints2 ones(size(matchedPoints2,1),1)]);

%matchedPoints2_hc_nor_transed=H*matchedPoints2_hc_nor;

matchedPoints2_nor_transed=H*matchedPoints2_hc_nor;

hc = matchedPoints2_nor_transed(3,:);

matchedPoints2_hc_nor_transed=matchedPoints2_nor_transed./(hc);

v1 = matchedPoints2_hc_nor_transed-matchedPoints1_hc_nor;

v2 = v1.*v1;

v3 = sum(v2,1);

error_dis = sqrt(v3);

% use error_dis to compute how much inliers based on threshold

Treshold = 0.2; %% the only thing you can adjust

inlier_vector = find(error_dis<=Treshold);

number_inlier = length(inlier_vector);

number_inlier_list(ii1) = number_inlier;

if ii1>1

if number_inlier_list(ii1) == max(number_inlier_list)

saved_index = inlier_vector;

end

else

saved_index = inlier_vector;

end

end

max(number_inlier_list)

% extract the inliers from the image

for ii = 1:length(saved_index)

idx = saved_index(ii);

inlier_ref(ii,:) = matchedPoints1_hc_nor(:,idx);

inlier_sou(ii,:) = matchedPoints2_hc_nor(:,idx);

% construct a matrix by inliers

a_m(:,:,ii) = [zeros(1,3) -inlier_ref(ii,3)*inlier_sou(ii,:) inlier_ref(ii,2)*inlier_sou(ii,:);

inlier_ref(ii,3)*inlier_sou(ii,:) zeros(1,3) -inlier_ref(ii,1)*inlier_sou(ii,:) ;

-inlier_ref(ii,2)*inlier_sou(ii,:) inlier_ref(ii,1)*inlier_sou(ii,:) zeros(1,3)];

% compute the A matrix by a

if ii == 1

A_final = a_m(:,:,ii);

else

A_final = [A_final;a_m(:,:,ii)];

end

end

% SVD to A and get final H

[U1,S1,V1] = svd(A_final);

h_vector_final= V1(:,9);

H_after_ransac = [h_vector_final(1:3)';

h_vector_final(4:6)';

h_vector_final(7:9)']

New_H = inv(Norm_matrix)*H_after_ransac*Norm_matrix;

% bonus part, compute the mean error after transformation

matchedPoints1_err = valid_points1.Location(indexPairs(:,1),:);

matchedPoints2_err = valid_points2.Location(indexPairs(:,2),:);

matchedPoints1_hc_err= transpose([matchedPoints1_err ones(size(matchedPoints1_err,1),1)]);

matchedPoints2_hc_err= transpose([matchedPoints2_err ones(size(matchedPoints2_err,1),1)]);

matchedPoints2_nor_transed_err=New_H*matchedPoints2_hc_err;

hc_err = matchedPoints2_nor_transed_err(3,:);

matchedPoints2_hc_transed_err=matchedPoints2_nor_transed_err./(hc_err);

v1_err = matchedPoints2_hc_transed_err-matchedPoints1_hc_err;

v2_err = v1_err.*v1_err;

v3_err = sum(v2_err,1);

error_dis_err = sqrt(v3_err);

total_err_bonus(aa)= mean(error_dis_err);

% compute the invert of H for source

H_inv = inv(final_H_for_ref*New_H);

for ii1 = 1:size(panel,1)

for ii2 = 1:size(panel,2)

% create pixel location(pl)

pl_panel = [ii2;ii1;1];

% convert panel pixel location back to source image

pl_sou = H_inv*pl_panel;

pl_sou_hc = round(pl_sou/pl_sou(3));

if pl_sou_hc(1)>0 && pl_sou_hc(2)>0 && pl_sou_hc(2)<size(source,1) && pl_sou_hc(1)<size(source,2)

if panel(ii1,ii2,:) == zeros(1,1,3)

panel(ii1,ii2,:) = source(pl_sou_hc(2),pl_sou_hc(1),:);

else

panel(ii1,ii2,:) = source(pl_sou_hc(2),pl_sou_hc(1),:);

panel(ii1,ii2,1) = 0.5*source(pl_sou_hc(2),pl_sou_hc(1),1)+0.5*panel(ii1,ii2,1);

panel(ii1,ii2,2) = 0.5*source(pl_sou_hc(2),pl_sou_hc(1),2)+0.5*panel(ii1,ii2,2);

panel(ii1,ii2,3) = 0.5*source(pl_sou_hc(2),pl_sou_hc(1),3)+0.5*panel(ii1,ii2,3);

end

end

end

end

end

imshow(uint8(panel))

mean(total_err_bonus)

Answer 1

[Bjorn Gustavsson]

You should have some good use for the histogram-matching contributions found on the file exchange. For example:

https://se.mathworks.com/matlabcentral/fileexchange/61928-rgb-histogram-matching?s_tid=ta_fx_results, matchhistograms, or perhaps the imhistmatch function out of the image processing toolbox.

HTH

[Bjorn Gustavsson]

As far as I understand the tools I linked to they could be used on the separate images before the stitching (modify one to fit the other). Just try and see what you get. If the result is what I anticipate the modified image should have colours similar to the "target-image" which then would make the stitching simpler.

[chee gang ngui]

%Could you show me how you run the program? So, I have two images, the first image named 1 and the second image named 2. Both images are store in this directory C:\Users\nguic\Documents\MATLAB\test_png\.

%demo

clc

close all

%load image

filename='image.jpg';

img=imread(filename);

ref=imread('reference.jpg');

%rgb_hist matching

out_image = RGBHistMatch(img,ref);

filename_=strcat(filename(1:end-4),'_rgb_hist_match.jpg');

imwrite(out_image, filename_);

%matchin the histogram of a given colored image i to a reference image r for

%generate the result image

%%

% Author: Mahmoud Afifi, York University

close all

I_=I; R_=R; %take a backup

Result=zeros(size(I));

for Channel=1:3

if size(R,2)>1

%an image, replace the image with its histogram equalization transformation

if size(R,3)>1

R_=R(:,:,Channel);

end

%Calc CDF (G) of the reference image

[countR,binsR]=imhist(R_);

%pR=countR/(size(R_,1)*size(R_,2));

pR=countR/numel(R_);

else

%else, the given R is the reference PDF

pR=R_;

end

G=255*cumsum(pR);

G=round(G);

if size(I,3)>1

I_=I(:,:,Channel);

end

%Calc the given image's histogram equalization transformation

[S] = HistEquTrans(I_);

%Map S to G

F=zeros(size(G)); %the map (F) where F represents S->G

min=1000000;minIndex=-1;

for i=1:size(S,1)

for j=1:size(G,1)

T=abs(S(i,1)-G(j,1));

if T==0

minIndex=j;

break

elseif T<min

minIndex=j;

min=T;

end

end

F(i,1)=minIndex;

minIndex=-1;

min=1000000;

end

%Map the value of mapped S (F) to corresponding pixels in the original image

result=zeros(size(I_(:)));

for k=1:size(F,1)

indecies=find(I_==k-1);

result(indecies)=F(k,1);

end

Result(:,:,Channel)=reshape(uint8(result),size(I_));

end

Result=Result/255;

end

function [s,result]=HistEquTrans(i)

% return the histogram equalization transformation

L=256;

if(size(i,3)>1)

i=rgb2gray(i);

end

%reshape the image into single dimension

temp=i(:);

[counts,bins]=imhist(temp); %calc histogram of the given image

p=counts/size(temp,1); %calc PDF of the image

s=(L-1)*cumsum(p); %calc s=(L-1)*CDF(p)

s=round(s); %round the float into the nearest integare

end

[chee gang ngui]

I have figured out the how to use the link you provided to me. Ty, but it seems like the result is not that great.

[Bjorn Gustavsson]

If you attach the images it might be possible for others to give more hands-on advice.

Answer 2

[chee gang ngui]

I am new with MATLAB. Do you think the code that you provided will able to intergrate with my code or I should do it seperately? Should I do it before image stitching or after?