Techniques for blending images that have overlap
Show older comments
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?
Accepted Answer
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.
15 Comments
Paul Safier
on 14 Sep 2021
Paul Safier
on 14 Sep 2021
DGM
on 14 Sep 2021
That part of the code just generates test images. Any 4D image stack should work.
Paul Safier
on 14 Sep 2021
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
on 15 Sep 2021
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
on 15 Sep 2021
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
on 15 Sep 2021
Paul Safier
on 17 Sep 2021
DGM
on 17 Sep 2021
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
on 17 Sep 2021
@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:
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
on 17 Sep 2021
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
on 17 Sep 2021
DGM
on 18 Sep 2021
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
on 22 Sep 2021
More Answers (0)
Categories
Find more on Convert Image Type in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
