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

3 views (last 30 days)
%% 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)

Answers (2)

Bjorn Gustavsson
Bjorn Gustavsson on 5 Jul 2022
Edited: Bjorn Gustavsson on 5 Jul 2022
You should have some good use for the histogram-matching contributions found on the file exchange. For example:
HTH
  4 Comments
Bjorn Gustavsson
Bjorn Gustavsson on 8 Jul 2022
If you attach the images it might be possible for others to give more hands-on advice.

Sign in to comment.


chee gang ngui
chee gang ngui on 6 Jul 2022
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?

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by