Techniques for blending images that have overlap

I have some images that all overlap each other by a known amount (e.g. 15% of the pixels on along a side). No registration should be needed since I know exactly the overlap amout, i.e. placement. My method of merging them into a single image is by taking the ratio of the mean values of the images' overlapping regions and then correcting one of the images by that ratio. For example:

ratioAB = mean(matA_ovlap(:))/mean(matB_ovlap(:));
matB_new = matB*ratioAB;

This is better than simply piecing them together, but still clearly leaves a lot to be desired. See result below. Any recommendations for a more sophisticated blending?

[KSSV]

Try weighted averaging.

[Paul Safier]

@KSSV what would the weights be? The distance within the overlap region?

Answer 1

[DGM]

This is a simple approach to doing edge-interpolation. This could also be done using interpolation tools in fewer lines, but right now I don't have the brain cells to spare in any effort to make this elegant. Just bear in mind that there's plenty of room for improvement. Anyone is free to post something better.

For this example, I'm going to use RGB color images. That might not be what you're using.

tiling = [3 3]; % image tiling [nup nacross]
overlap = [64 64]; % in pixels [rows cols]

% generate colored test images as a 4D RGB frame stack
% some of this uses MIMT; ignore that.
nframes = prod(tiling);
inpict = imread('cameraman.tif'); % raw single-channel source
A = zeros([imsize(inpict,2) 3 nframes]);
A(:,:,:,1) = repmat(im2double(inpict),[1 1 3]).*permute([1 0.3 0.8],[1 3 2]);
for f = 2:nframes
	A(:,:,:,f) = imtweak(A(:,:,:,1),'lchab',[1 1 (f-1)/nframes]); % rotate hue
end

% do setup
overlap = round(overlap/2)*2; % make sure it's even
s0 = [size(inpict,1) size(inpict,2)]; % source frame geometry
st = 2*(s0-overlap) + (s0-2*overlap).*(tiling-2) + overlap.*(tiling-1); % output image geometry

myn = repmat(linspace(1,0,overlap(1)).',[1 s0(2)]); % transition masks
mys = repmat(linspace(0,1,overlap(1)).',[1 s0(2)]);
mxw = repmat(linspace(1,0,overlap(2)),[st(1) 1]);
mxe = repmat(linspace(0,1,overlap(2)),[st(1) 1]);
colpict = zeros([st(1) s0(2) size(A,3) tiling(2)]); % preallocate
outpict = zeros([st size(A,3)]);

% assemble the frames into columns
f = 1;
for n = 1:tiling(2)
	for m = 1:tiling(1)
		if m == 1
			sr = s0(1)-overlap(1);
			ry = 1:sr; 
			colpict(ry,:,:,n) = A(1:sr,:,:,f);
								
		elseif m == tiling(1)
			sr = s0(1)-overlap(1);
			ry = ry(end)+1:ry(end)+sr;
			colpict(ry,:,:,n) = A(overlap(1)+1:end,:,:,f);
			
		else
			sr = s0(1)-2*overlap(1);
			ry = ry(end)+1:ry(end)+sr; 
			colpict(ry,:,:,n) = A(overlap(1)+1:overlap(1)+sr,:,:,f);

		end
		
		% insert transition
		if m ~= tiling(1)
			ry = ry(end)+1:ry(end)+overlap(1);
			colpict(ry,:,:,n) = A(end-overlap(1)+1:end,:,:,f).*myn ...
				              + A(1:overlap(1),:,:,f+1).*mys;
		end

		f = f+1;
	end
end

% assemble the columns into an image
for n = 1:tiling(2)
	if n == 1
		sr = s0(2)-overlap(2);
		rx = 1:sr; 
		outpict(:,rx,:) = colpict(:,1:sr,:,n);

	elseif n == tiling(2)
		sr = s0(2)-overlap(2);
		rx = rx(end)+1:rx(end)+sr;
		outpict(:,rx,:) = colpict(:,overlap(2)+1:end,:,n);

	else
		sr = s0(2)-2*overlap(2);
		rx = rx(end)+1:rx(end)+sr; 
		outpict(:,rx,:) = colpict(:,overlap(2)+1:overlap(2)+sr,:,n);

	end

	% insert transition
	if n ~= tiling(2)
		rx = rx(end)+1:rx(end)+overlap(2);
		outpict(:,rx,:) = colpict(:,end-overlap(2)+1:end,:,n).*mxw ...
					    + colpict(:,1:overlap(2),:,n+1).*mxe;
	end
end
clear colpict

In this example, each supercolumn of the image is assembled from single-source blocks and transition blocks betwen them. This takes up a lot of space with indexing, but the interpolation needs to be applied to the minimum number of pixels. The interpolation itself is a simple linear transition implemented with multiplicative composition using precalculated masks.

The supercolumns are assembled into the final image in a similar fashion.

[Paul Safier]

Thanks. Let me study this.

In this line, do you mean to use size?

A = zeros([imsize(inpict,2) 3 nframes]);

[Paul Safier]

Nevermind, I see imsize and imtweak are from your package. I will get them online...

[DGM]

That part of the code just generates test images. Any 4D image stack should work.

[Paul Safier]

Thanks! I will take a look.

Would this code work if the adjacent images have no overlap, i.e. they just show seams? For example in the image I show, after the overlap has been chopped off, if it was again separated into distinct images that belong together, can the method treat this?

[DGM]

Ah, no. I didn't bother handling the case where overlap=0. It breaks if there's zero overlap or if the overlap exceeds the image width. It's just a demo.

If the images all have the same size, then the zero-overlap case looks just a simple tiling operation, and the above code would be a slow option.

tiling = [3 3]; % image tiling [nup nacross]

% generate test images as a 4D RGB frame stack
% some of this uses MIMT; ignore that.
nframes = prod(tiling);
inpict = imread('cameraman.tif'); % raw single-channel source
A = zeros([imsize(inpict,2) 3 nframes]);
A(:,:,:,1) = repmat(im2double(inpict),[1 1 3]).*permute([1 0.3 0.8],[1 3 2]);
for f = 2:nframes
	A(:,:,:,f) = imtweak(A(:,:,:,1),'lchab',[1 1 (f-1)/nframes]);
end

% tiling can be done using cell array tools
% this might simplify if the images are already in a cell array
% instead of being in a 4D stack
outpict = cell2mat(reshape(squeeze(num2cell(A,1:3)),tiling));

% or it can be done using MIMT imtile()
outpict = imtile(A,tiling,'direction','col'); % this is MIMT imtile(), not IPT imtile()
% the output is the same

Regarding the hypothetical "... if it was again separated into distinct images ...", if this actually needs to be done in practice mat2cell() or MIMT imdetile() may be used. Since its geometry should be integer-divisible by the tiling vector, either should work.

[Paul Safier]

I see. In this hypothetical case, is there a way to blend the seams even though there is no more overlap? For example, do you know how to 'fix' the seam I show in the image below (taken from one of Matlab's examples)?

[DGM]

Oh I see what you're saying. Hmm. I imagine there are some canonical approaches. I don't know what they would be or how well they would work, though.

[Paul Safier]

I might investigate Matlab's roifilt2 function...

Thanks for the advice on all this!

[Paul Safier]

Hi @DGM , if I converted a matrix to RGB, before running your example code, does your toolbox have a routine that would be able to convert back to a matrix of doubles?

[DGM]

I'm not sure what you mean.

If an image is RGB, that just means it has 3 channels (i.e. it's an MxNx3 array), as opposed to a monochrome (I) array (MxNx1) or an RGBA array (MxNx4). It's a matter of shape, and should be independent of the numeric class.

If you're just asking how to convert an image (of any arbitrary shape) from one numeric class to another, MATLAB already has tools like im2double(), im2uint8(), etc that handle the casting and rescaling. MIMT has imcast() which is sort of a generalized version of those, but other than that convenience, it does the same thing.

If you're wondering if you need to use an RGB workflow to use the above code and you want to know how to switch back and forth because you're processing monochrome images, there's a couple things to point out.

If instead of generating an RGB test image, I just used a monochrome test image like so:

% ...
nframes = prod(tiling);
inpict = imread('cameraman.tif'); % raw single-channel source
% instead of making the test image RGB for viewing purposes
% A = zeros([imsize(inpict,2) 3 nframes]);
% A(:,:,:,1) = repmat(im2double(inpict),[1 1 3]).*permute([1 0.3 0.8],[1 3 2]);
% for f = 2:nframes
% 	A(:,:,:,f) = imtweak(A(:,:,:,1),'lchab',[1 1 (f-1)/nframes]);
% end

% just replicate the same single-channel (I) image
A = repmat(im2double(inpict),[1 1 1 nframes]);
% ...

... the rest of the code works the same, but the output is now has only one channel, same as A. So the core of the code can handle RGB images, but it doesn't need RGB images.

While in this case, you shouldn't need to convert from RGB to mono, other cases elsewhere you might. In those cases, MATLAB has things like rgb2gray() which reduces the image by extracting the luma. MIMT has mono(), which again is more generalized and can extract luma, value, lightness, or other things.

[Paul Safier]

@DGM I'm trying to understand that if I start with a floating point image, where the average/min/max values are important. Say, they represent a surface temperature, for example. Then I convert the matrices to an image format suitable for a blending routine (for example RGB or grayscale), I would want to convert the final (stitched image) back to the original scale so it would be a meaningful temperature.

For example, the attached matrices have an overlap of 20. It does not seem to work properly to deal with them directly. Not sure why the code you posted won't work with them, as is. The attached image shows the results of both ways. I had to either convert them to an RGB image, or convert them to a single channel image as such:

A(:,:,1,k) = uint8(255*mat2gray(B.hgtDep));

Not sure why I couldn't just use:

A(:,:,1,k) = B.hgtDep;

but, this is why I'm asking about converting back to a floating point matrix...

[DGM]

Ah okay. There are some things that can be done. All I changed was the tiling (product of the tiling vector needs to match the number of images), and I changed the calculation of s0 to use A instead of inpict.

tiling = [2 2]; % image tiling [nup nacross] 
overlap = [10 10]; % in pixels [rows cols]

% read test frames from files
nframes = prod(tiling);
aa = load('1_Dep.mat');
bb = load('2_Dep.mat');
cc = load('3_Dep.mat');
dd = load('4_Dep.mat');
A = cat(4,aa.hgtDep,bb.hgtDep,cc.hgtDep,dd.hgtDep);

overlap = round(overlap/2)*2; % make sure it's even
s0 = [size(A,1) size(A,2)]; % source frame geometry 
% ... the rest of the code is the same

This gives a 2x2 tiling, and note that no conversion or normalization was done. The output image should be scaled the same as the input images.

Of course, if you try to view it, imshow() will just display a solid white image since it expects floating-point images to be normalized. If you want to view it, you can just use imshow(outpict,[]), which displays the image as normalized to its extrema.

Trying to normalize each input image is problematic because you're likely losing information. Unless you specify the normaiization limits, mat2gray uses the image extrema, which will tend to vary. Once they've been normalized, denormalizing requires going back to the source images and figuring out what the global extrema were. Picking normalizing limits explicitly would avoid all that.

That said, I don't know that normalizing (or uint8 conversion) is strictly necessary. Just like with imshow(), it might depend on what other image tools you need to use.

[Paul Safier]

@DGM Sorry for poorly explaining my concern in the previous post. I don't have an issue setting up to run these 4 matrices. My issue is with the results of the blend. It works well when the matrices are scalled via mat2gray, but works poorly when run as is. If you comment/uncomment the different definitions of matrix A, it will run with and without the scaling. It will display the result at the end and you can see what I'm concerned with.

tiling = [2 2]; % image tiling [nup nacross]
overlap = [20 20]; % in pixels [rows cols]


nframes = prod(tiling);
aa = load('1_Dep.mat');
bb = load('2_Dep.mat');
cc = load('3_Dep.mat');
dd = load('4_Dep.mat');

% HERE IS WHERE THE DIFFERENCE OCCURS
% This doesn't seem to blend well
A = cat(4,aa.hgtDep,cc.hgtDep,bb.hgtDep,dd.hgtDep);
% This does blend well
%A = cat(4,mat2gray(aa.hgtDep),mat2gray(cc.hgtDep),mat2gray(bb.hgtDep),mat2gray(dd.hgtDep));

overlap = round(overlap/2)*2; % make sure it's even
s0 = [size(A,1) size(A,2)]; % source frame geometry
st = 2*(s0-overlap) + (s0-2*overlap).*(tiling-2) + overlap.*(tiling-1); % output image geometry

myn = repmat(linspace(1,0,overlap(1)).',[1 s0(2)]); % transition masks
mys = repmat(linspace(0,1,overlap(1)).',[1 s0(2)]);
mxw = repmat(linspace(1,0,overlap(2)),[st(1) 1]);
mxe = repmat(linspace(0,1,overlap(2)),[st(1) 1]);
colpict = zeros([st(1) s0(2) size(A,3) tiling(2)]); % preallocate
outpict = zeros([st size(A,3)]);

% assemble the frames into columns
f = 1;
for n = 1:tiling(2)
	for m = 1:tiling(1)
		if m == 1
			sr = s0(1)-overlap(1);
			ry = 1:sr; 
			colpict(ry,:,:,n) = A(1:sr,:,:,f);
								
		elseif m == tiling(1)
			sr = s0(1)-overlap(1);
			ry = ry(end)+1:ry(end)+sr;
			colpict(ry,:,:,n) = A(overlap(1)+1:end,:,:,f);
			
		else
			sr = s0(1)-2*overlap(1);
			ry = ry(end)+1:ry(end)+sr; 
			colpict(ry,:,:,n) = A(overlap(1)+1:overlap(1)+sr,:,:,f);

		end
		
		% insert transition
		if m ~= tiling(1)
			ry = ry(end)+1:ry(end)+overlap(1);
			colpict(ry,:,:,n) = A(end-overlap(1)+1:end,:,:,f).*myn ...
				              + A(1:overlap(1),:,:,f+1).*mys;
		end

		f = f+1;
	end
end

% assemble the columns into an image
for n = 1:tiling(2)
	if n == 1
		sr = s0(2)-overlap(2);
		rx = 1:sr; 
		outpict(:,rx,:) = colpict(:,1:sr,:,n);

	elseif n == tiling(2)
		sr = s0(2)-overlap(2);
		rx = rx(end)+1:rx(end)+sr;
		outpict(:,rx,:) = colpict(:,overlap(2)+1:end,:,n);

	else
		sr = s0(2)-2*overlap(2);
		rx = rx(end)+1:rx(end)+sr; 
		outpict(:,rx,:) = colpict(:,overlap(2)+1:overlap(2)+sr,:,n);

	end

	% insert transition
	if n ~= tiling(2)
		rx = rx(end)+1:rx(end)+overlap(2);
		outpict(:,rx,:) = colpict(:,end-overlap(2)+1:end,:,n).*mxw ...
					    + colpict(:,1:overlap(2),:,n+1).*mxe;
	end
end
clear colpict

% check how well the blending did
figure; imagesc(outpict)
colormap(jet(256))
colorbar
%figure; imshow(outpict,[])

[DGM]

I don't think the issue is with the interpolation. It looks like the data in each image has some offset. Normalizing the input images gets rid of the offsets. If that's the case, I don't know where the offset came from or what it means in this context, so I don't know if it's unwanted information.

A = cat(4,aa.hgtDep,cc.hgtDep,bb.hgtDep,dd.hgtDep);
outpict = mat2gray(imtile(A,[2 2],'direction','col'));

If the inputs are normalized, the question becomes which values should be used as references for denormalizing. If the inputs instead have their offsets removed somehow, the question becomes which tile is the calibrated one? Without figuring out where the offsets are coming from, I don't think I can answer either.

[Paul Safier]

@DGM Understood. Thanks for your help with this.