How to Stitch 25 images with different dynamic range

Hi all,

Actually the question about how to take care of already stitched picture. I don't have any questions about the stitched process itself, I know how to do it. However after I do stitch I see the transition between the pictures. I would like to make a smooth transition between the pictures (see red arrows).

Every picture has 768X768 pixels. Every picture has a different dynamic rangle. So first of all I'd like to bring all the pictures to the same dynamic range, then to make slope compenstation on each picture and maybe then to stitch all of them. See the picture above one of the 25 pictures. I did mesh (X,Y,Z) and turned around this picture to show you the slope. I'd like to remove that slope.

I saw in one of the topics (https://www.mathworks.com/matlabcentral/answers/86498-how-to-blend-image) Image Analyst suggested to "You can use conv2() to blur the edge zone of the stitched pair, or you can stitch with a weighted average so that each edge has a ramp where it goes from full strength to zero." But how exactly to do this? Can you give me some example how to do it?

Here I attached 2 matrices for example to show you the slope and different dynamic ranges between them. So I'd like to bring all the matrices I have (25) to the same dynamic range and make slope compensation of each matrix I have.

Thank you very much.

P.S. If you want to stitch these two matrices you should put them as following:

Answer 1

[Bjorn Gustavsson]

You could do something like this to level the regions:

dI = mean(Im(row_above,:)-Im(row_above+1,:)); % Perhaps median works better, depends on intensity variance
Im(1:row_above,:) = Im(1:row_above,:) - dI/2;
Im(row_below:end,:) = Im(row_below:end,:) + dI/2;

Then you simply have to step through all the seams to level the entire image. You should also keep track of all the dIs to know where the absolute intensity-correction goes.

HTH

[Dimani4]

Hi Bjorn,

Thank you for your answer. I dont understand something.

The first line:dI = mean(Im(row_above,:)-Im(row_above+1,:));

Here you mean I should create for loop for calculating the mean values of the differences between the two adjacent rows?

dI(1)=mean(Im(1,:)-Im(2,:));

dI(2)=mean(Im(2,:)-Im(3,:));

dI(3)=mean(Im(3,:)-Im(4,:))..... ?

Sorry, but the rest of the lines I dont understand?

Im(1:row_above,:) = Im(1:row_above,:) - dI/2;

Im(row_below:end,:) = Im(row_below:end,:) + dI/2;

Should I do loop as well?

The first one should be like this?

Im(1:2,:) = Im(1:2,:) - dI(1)/2;

Second one?

Im(2:3,:) = Im(2:3,:) - dI(2)/2;

Im(1:end,:) = Im(1:end,:) + dI(1)/2;

Third one?

Im(3:4,:) = Im(3:4,:) - dI(3)/2;

Im(2:end,:) = Im(2:end,:) + dI(3)/2;....

Thank you.

[Bjorn Gustavsson]

No, not exactly. In the image you show above it seems as if you have stitched together 5 rows of images above each other, right?

That would make the boundaries be on row-numbers:

Rows_above = 768*[1:4];
for iRow = 1:numel(Rows_above)
  row_above = Rows_above(iRow);
  row_below = row_above + 1;
  dI = mean(Im(row_above,:)-Im(row_above+1,:)); % Perhaps median works better, depends on intensity variance
  Im(1:row_above,:) = Im(1:row_above,:) - dI/2;
  Im(row_below:end,:) = Im(row_below:end,:) + dI/2;
end

This would take the mean (or median if that works better...) of the difference between intensities above and below the stitching, and then shift the two parts half of that. This works if you have a constand zero-level offset for each different image. Then you have to repeat the procedure for each stitching. You can also use more than one line above and below - just check that you get a good scalar off-set to work with

[Dimani4]

Thank you Bjorn.

Yes, exactly 5 rows of images (actually 5X5 5 rows, 5 columns).Each picture is 768X768pixels

See the results I did according your idea.

Take a look at the results.

Here is without your suggestion:

Here zoomed region

Here is with your suggestion (with median):

And here is the full stitching:

[Bjorn Gustavsson]

OK, so your spiky higher-intensity-regions vary much between the sub-images, while the smoother low-intensity-bands are reasonably similar. After looking at the 2 matrices in your files it is not exactly clear to me what you want to achieve...

Here's the 2 stitched images with columns 190 and 300 to the right and rows 768 and 769 below.

and here's a zoom in on a part of the intersection-region:

If you want to level the low-intensity regions you have to adapt my snipped to also select only the low-intensity-regions (here columns ~140-218 etc and calculate the average difference between the pixels in those regions...)

If you want something else you have to give a more detailed explanation.

Answer 2

[Dimani4]

Bjorn,

All I want is to do kind of plane removal of each matrix to get the intensity distribution more homogenical and then to get all the matrices to the same dynamic range to eliminate the differences in the intensities between them as you pointed out in the Figure1 and then to get all the pictures stitched.

Actually I did some kind of very primitive slope removal of every line in each matrix and then got them stitched.

The algorithm is like that:

for sc=1:size(pl_laser_matr,2)
    
pl_laser_matr{sc}=minus_slope(pl_laser_matr{sc}); 

end

function matri_compens2=minus_slope(matrix)

matri_compens=matrix;

for i=1:size(matri_compens,2)  %minus vertical slope
    
    lineB=(matri_compens(:,i))';
    x1=1:size(matri_compens,1);
    p1=polyfit(x1,lineB,1);
    p1(1);
    for j=1:size(matri_compens,1)
       
        new_line1(j,i)=lineB(j)-p1(1)*x1(j);
        
    end
   
end

matri_compens2=new_line1;

And here is the result of the first matrix:

From the first sight it seems OK however you got some artifacts in the edges, look at the zoom I did in the picture above:

It happens because I apply this algorithm to the straight vertical lines (columns) but my structure is little bit inclined so it creates these artifacts.

Here the whole stitch: (25 frames)

And here is zoom:

You see all these artificial lines (red arrows). Not so good for me.

Thank you for your help.

[Bjorn Gustavsson]

OK, now I understand that part of the problem. Since your columns of bright-spots are not aligned with the columns this becomes a more labour-intensive task. I'd go about it something like this:

1, rotate the 768x768 images such that the spike-columns align with the image-columns, you can use imrotate, see the help and documentation for that function. I'd use 'bicubic' for method, and 'loose' for the bbox.

2, do your adjustments column-by-column to your hearts content. (I havent fully understood what you want to do to what purpose, but I'm happy to assume you have a plan.)

3, rotate back to the original orientation. It seems to me that you'll have to do that to get the stitching right.

4, crop that rotated-back image to your initial 768x768 size.

5, stitch.

6, celebrate?

[Dimani4]

Thank you Bjorn,

Actually I tried something similar, but maybe I found out more elegant way to do it. What do you think?

I try to fit it to some surface, find out the formula for that surface, find out the points in each (X,Y) and then just subtract them from each point of my Z value.

So what I do:

1. Find out surface fit to my points:

X=1:size(MyMatrix,2);
Y=1:size(MyMatrix,1);
[x,y]=meshgrid(X,Y);
f2=fit( [x(:), y(:)], MyMatrix(:), 'poly11' ); %linear surface fit
formula(f2)

ans =

p00 + p10*x + p01*y;

coeffs = coeffvalues(f2)

coeffs =

    0.1121    0.0000   -0.0001  
    
figure, sd2=plot(f2,[x(:),y(:)],MyMatrix(:));%plotting my matrix with the fit

Now lets plot only the fit.

tyr2=sd2.Parent
mesh(tyr2.Children(2).XData,tyr2.Children(2).YData,tyr2.Children(2).ZData)

Since I know the formula and the coefficients I can make a new Z value corresponding to my surface fit:

Z_new=0.1121+0*x-0.0001*y;

Then I'll subtract MyMatrix-Z_new.

That's the way I do kind of plane compensation.

also I can do quadratic polynomial if I'll do:

f3 = fit( [x(:), y(:)], handles.matrices{1}(:), 'poly22' );

And do the same procedure (i.e. find out new Z_new and substract from MyMatrix).. this how I do quadratic plane compensation... What do you think?

The picture below is quadratic plane after fit (f3).

[Bjorn Gustavsson]

If that gives you a correction you're happy with, then all's well that ends well?

Answer 3

[Dimani4]

Hey Bjorn,

Not yet. Now I have a problem of different color scales. Do you know how to bring all the pictures to the same range?

Look at the stitch:

This is the typical surface plane (only first matrix). I guess the other surface planes of other matrices look approx. the same:

Also I saw when I subtract my matrix from the surface plane I got negative numbers, what I did

func_f=str2func(['@(x,y)' fformulaD]);
func_f_fitted_surface_matrix=func_f(x,y);
diff_matr=matrix-func_f_fitted_surface_matrix;

comp_matrix=diff_matr+abs(min(diff_matr(:)));

Thank you.

[Bjorn Gustavsson]

If you want to normalize max and min of each image to 0 to 1, you can do this:

ImNorm = (Im-min(Im(:)))/(max(Im(:))-min(Im(:)));

[Dimani4]

Thank you for your suggestion.I use it.

However I still have a problem:

I did histograms graph of first 5 frames and found out they are little bit different.

I try to bring all the frames to the same reference of histogram:

pl_laser_matr{sc} = (pl_laser_matr{sc}-min(pl_laser_matr{sc}(:)))/(max(pl_laser_matr{sc}(:))-min(pl_laser_matr{sc}(:)));

if sc>1
       
pl_laser_matr{sc} = imhistmatch(pl_laser_matr{sc},pl_laser_matr{1});   
    
end

%----------------------------------------------------------------------------------------------

if handles.checkbox7.Value==1
    
 if sc<=str2double(handles.edit8.String)   
    
    figure('Name',strcat('Frame No:',num2str(sc)));imagesc(pl_laser_matr{sc});colorbar;
    figure('Name',strcat('Histigram of Frame No:',num2str(sc)));imhist(pl_laser_matr{sc});

 end
 
end

[Dimani4]

However it didnt solve the problem, more than that the histograms I got are different from the reference I give.

As you can see it didnt bring all the frames to the same histogram disribution as reference (my reference is first frame).

What can you suggest? Maybe there is another way to get rid of this transitions between the frames?

Thank you.

[Image Analyst]

How were the images captured in the first place? Did you take a "blank shot" so that you can do a background correction to flatten it? All cameras have shading where it's darker at the endge than at the center and you need to take a photo of a uniform field to know the percentage of light at each pixel and to correct for that (in most cases division but for some cases, such as x-ray or fluorescence, subtraction).

[Dimani4]

Thank you. I expected you to see my post, because as I mentioned in the first post you have answered the similar issue here (https://www.mathworks.com/matlabcentral/answers/86498-how-to-blend-image).

I didnt take the photos. I will ask the guy who took it.

[Bjorn Gustavsson]

Stay off the histogram-matching. That is not the right way to solve this problem.

Image Analysts suggestion to make a proper background correction is important.

However, you now seem to be back at the initial problem-description. Meaning that you now have intensity-transitions at the stitchings. I suggested that you can remove these transitions by leveling the dark bands. Something like this:

Rows_above = 768*[1:4];
for iRow = 1:numel(Rows_above)
  row_above = Rows_above(iRow);
  row_below = row_above + 1;
  dI = mean(Im(row_above,idx_darkBands{iRow}))-...
            Im(row_above+1,idx_darkBands{iRow}))); % Perhaps median works better, depends on intensity variance
  Im(1:row_above,:) = Im(1:row_above,:) - dI/2;
  Im(row_below:end,:) = Im(row_below:end,:) + dI/2;
end

I suggest you do this on each individual column of images before stitching the columns together. You also have to create the cell-array of indices, idx_darkBands,to the dark bands for each image-pair.

This should get rid of the horizontal intensity-steps at the stitches.

[Dimani4]

Thank you Bjorn. Here is my another try. I did fit to polynom (x=1,y=0) (formula=0.087352+1.8237e-05*x)

Here is the full picture:

and a zoom:

You are absolutely right, I got back to the initial problem and I need to get rid off the transitions in the stitching. Can you please explain what do you mean by creation of idx_darkBands. I dont really understand it. One more thing, as you can see I got regions with brighter background in some of the stiches. What do you think, how can I get rid of that?

Thank you.

About the proper background correction I'm still waiting for the answer.

[Bjorn Gustavsson]

In my comment from "27 Jan 2021 at 16:06" I have a horizontal line-plot in the embedded figure. In that one there are 2 rows plotted, one from just above the stitch, one from just below. In those lines you have low-intensity-regions at (approximately) 150-215, 370-430 and 590-660. Removing the intensity-difference between regions would be my first step. Since your bands have a slight slope you will have to save these indices for each stitching. So in order to run my attempt you would have to make a cell-array with those indices for each stitching. Perhaps a 2-D cell-array would be best, organized to match your image-stitching. Something like this:

iRow = 1;
iCol = 1;
idx_darkBands{iRow,iCol} = [150:215,370:430,590:660];

Then you'd have to add the indices for the dark bands at each stitching in the vertical direction between the first column of images by incremental iRow and then start over with the second column of images. If you do this you'll have to modify my snippet to account for the cell-array being 2D.

[Dimani4]

Hey Bjorn,

Thank you for your suggestion.

I want to take these parts automatically, so I thought to draw profile and just exclude these regions but as it turned out to exclude them its not so easy thing because they are very noisy.

xMatr=1:size(init_matr1,2);
yMatr(1:size(init_matr1,2))=768;
[cx,cy,c]=improfile(init_matr1,xMatr,yMatr);   
figure;imagesc(init_matr1);
% colormap(hot(255));
colorbar;
hold on;
plot(xMatr,yMatr,'LineStyle','--','LineWidth',2,'Color','green');
figure;
plot(cx,c);

I took row 768. init_matr1 is the stitched matrix of (768*5X768*5). Here if you want you can download it (https://wetransfer.com/downloads/c0ad07dcadd336b16cd700fd4e6c47c220210204120137/a0983e6fc4cf07d5e59103ef84449f1120210204120205/07b050).

And zoomed region:

What do you think, how can I take these low-intensity regions automatically?

Thank you.