Stitching sub images to reconstruct full image

I have subdivided an image into 4x4 tiles and produced a subimage and obtained a binary imaged. I add each of these binary images to a binary stack using a loop.

After the loop, how do I reconstruct/stitch the sub binary images back to the full image size, so I can create a binary image representing the while raw image so to use as a mask:

My binary sub images are expressed as :

BI(:,:,i)   %where i is 1:16 as I'm using 4x4 tiles

This is my approach that isn't working:

%Now combine binary images so to create regions to act as mask on original
%image
size(BI)   %Confirm there are 16 planes of images in the binarystack
Binary=[]; %Create empty Binary Image that will hold reconstructed sub images 
ct=0;  %counter
for jj=1:tiles
  for ii=1:tiles  
      ii
      jj
      ct=ct+1
      startingCols(jj)
      startingRows(ii)

    Binary(startingCols(ii):endingCols(ii), startingRows(jj):endingRows(jj))=BI(:,:,ct);
    figure(2)
    subplot(4,4,ct)
    imshow(BI(:,:,ct),[0,1])
    Binary=BI; 
    end
  end

Answer 1

[Mohammad Abouali]

Again, why don't you use blockproc()? :D

what you are asking, i.e. stitching the sub images into one big one, can be done in a single command using blockproc() as follows:

Binary=blockproc( reshape(1:16,4,4)', [1,1], @(x) BI(:,:,x.data) );

where Binary contains all the sub images placed in the order you show in the figure. 16 sub images are not much but if there were more you can even do this in parallel as easily as:

Binary =blockproc( reshape(1:16,4,4)', [1,1], @(x) BI(:,:,x.data), 'UseParallel',true);

[Jason]

Yes this worked, brilliant - thankyou the only thing that put me off with this approach is that I don't udnerstand it well enough, for example if you wanted to perform two operations how would the syntax look?

[Jason]

I just noticed that this also rounds to integer tiles. My raw image (ROI) is 101x101 pixels, I divide each dimension by 4, which gives 25 pixels for eacg subregion. This then results in the"combined binary" image being 100x100.

My next step then falls over as I am trying to use this binary image via drawing regions round bright objects to superimpose on the raw ROI image. As they are different sizes, its falling over:

%DrawRegions from Binary to raw image:
    labeledImage = bwlabel(Binary, 8);     % Label each blob so we can make measurements of it
    blobMeasurements = regionprops(labeledImage, ROI, 'all');     
    numberOfBlobs = size(blobMeasurements, 1);
    hold on;
    boundaries = bwboundaries(Binary);  
    numberOfBoundaries = size(boundaries);
    for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 1);
    end

It seems trying to handle the additional pixels is a real headache!

[Mohammad Abouali]

If your sub images, i.e. BI, is of size 25x25x16 then of course when you put them back together in a 4x4 tiles you would get 100x100 pixels and not 101x101. If your original image before processing was 101x101 and after tiling it and processing it you ended up with 16 tiles of 25x25 then pretty much you have dropped one row and column.

I can suggest two approaches:

Approach 1 based on what you have and the codes that you posted I think you are dropping the last raw and column of the original image. Therefore, that row and column is not being processed as part of your tile. then you have two options:

Approach1, Option 1 Let say OrigImage is your original image with size 101x101. Then you tile it into BI with size 25x25x16. You can stitch BI back using the command I gave you above:

Binary=blockproc( reshape(1:16,4,4)', [1,1], @(x) BI(:,:,x.data) );

Binary would be 100x100. Then to compare it with original image resize your original image as follow

OrigImage=OrigImage(1:end-1,1:end=1);

The above command drops the last row and column from your original image and makes it 100x100, now both Binary and OrigImage have the same size.

Approach1, Option 2 This is the same as previous option, but instead you resize the binary to make it 101x101. Here is what you should do:

% stitch the tiles
Binary=blockproc( reshape(1:16,4,4)', [1,1], @(x) BI(:,:,x.data) );

% resize Binary
Binary=imresize(Binary,[101,101],'nearest');

Now both Binary and OrigImage are 101x101.

Approach 2 The second approach is to go back and fix the tiling; so that you don't drop a row and a column from the image. I really suggest that you use blockproc and then you won't have any of these problems.

Let's say you have combined all the processing that you need to do on a tile in a function called myOperation.m. Therefore you have

BI(:,:,1)=myOperation(Tile1);
BI(:,:,2)=myOperation(Tile2);
... 
BI(:,:,16)=myOperation(Tile16); of course you do this in a loop.

This means that you need to divide the original image into Tiles (tile 1 to tile 16. then you need to stitch BI back into one image.

If you have such function then you can use blockprock as follow:

Binary=blockproc(OrigImage,[tileRowSize, tileColSize],@(x) myOperation(x.data));

Note that the above command does not requires the origimage to be divided in tile. and note that you won't get the BIs but you would get directly the stitched image.

In this approach you just have to make sure that myOperation should be able to work on variable sized tiles.

[Mohammad Abouali]

If you are just thresholding using graythresh you can write your blockproc like this

Binary=blockproc(OrigImage,[25,25],@(x) (im2bw(x.data,graythresh(x.data(:)))) )

That would threshold your image using Otsu's method and it is not global, the threshold is decided on each 25x25 tile separately. You don't need to divide the image into tiles, you don't need to spend time stitching them back.

Another approach is to use the function that ImageAnalyst told you.

[Jason]

Thankyou very much

Answer 2

[Iain]

depending how you've done it....

 BInew = permute(BI,[1 3 2]);
 BInew = reshape(BInew,16*size(BI,1),[]);

ought to work.

[Jason]

Thanks for identifying this method, I don't understand the [1 3 2] part.

I require a 4x4 array Thanks

[Iain]

I wanted to swap the 2nd and 3rd dimensions of your 3D array around, so that the pixels would be stored down a the first column of the first segment, then down the first column of the 2nd segment... etc.

You might want:

 BInew = reshape(BI,[size(BI,[1 2]) 4 4]);
 BInew = permute(BInew,[1 3 2 4]);
 BInew = reshape(BInew,size(BI,1)*4,[]);

[Jason]

Thanks for your time, Im getting this error:

Error using size
Dimension argument must be a positive integer scalar within indexing range.

if I do

size(BI) I get:
ans =

    25    25    16

[Iain]

BInew = reshape(BI,[size(BI,[1 2]) 4 4]);

needs to be:

 BInew = reshape(BI,[25 25 4 4]);

[Jason]

Im using this:

BInew = reshape(BI,[25 25 4 4]);
BInew = permute(BInew,[1 3 2 4]);
size(BInew)
figure(3)
imshow(BInew,[0,1])

and get the error:

Error using images.internal.imageDisplayValidateParams>validateCData (line 115)
Multi-plane image inputs must be RGB images of size MxNx3.

Error in images.internal.imageDisplayValidateParams (line 27)
common_args.CData = validateCData(common_args.CData,image_type);

If I do size(BInew) I get:

ans =

      25     4    25     4

[Iain]

You forgot the last reshape

 BInew = reshape(BInew,size(BI,1)*4,[]);

 BInew = reshape(BInew,100,[]);

[Jason]

Ah OK. How would I change this line to cater for an arbitrary number of tiles (so far I had tiles=4)?

BInew = permute(BInew,[1 3 2 4]);

[Iain]

You'd need to change the reshapes, not the permute. You've got 16 tiles though.... (4 x 4)

 BInew = reshape(BI,[25 25 tileshigh tileswide]);
 BInew = permute(BInew,[1 3 2 4]);
 BInew = reshape(BInew, 25*tileshigh, [])

[Jason]

Thankyou

Answer 3

[Image Analyst]

Why bother saying this:

Binary(startingCols(ii):endingCols(ii), startingRows(jj):endingRows(jj))=BI(:,:,ct);

if you're just going to say this

Binary=BI;

three lines later? Not only that, but you switched rows and columns. The first index should be rows, not columns as you have it, and the second index should be columns, not rows as you have it.

[Jason]

The Binary=BI was a typo,: I have this now:

Binary = zeros(rows, cols);
ct=0;
for jj=1:tiles
  for ii=1:tiles  
      ii;
      jj;
      ct=ct+1;
      p=[ct,startingRows(jj),endingRows(jj),startingCols(ii),endingCols(ii)]   %print out info for debugging
      Binary(startingRows(jj):endingRows(jj),startingCols(ii):endingCols(ii))=BI(:,:,ct);
  end
end

The error Is:

Subscripted assignment dimension mismatch.

Error in MultiImageanalysis>pushbutton7_Callback (line 1091)
      Binary(startingRows(jj):endingRows(jj),startingCols(ii):endingCols(ii))=BI(:,:,ct);

[Image Analyst]

Well, what are startingRows(jj),endingRows(jj),startingCols(ii), endingCols(ii), and size(BI)? Chances are they aren't the full size of a plane of BI. For example if BI is 100 by 200 then the distance between rows must be 100 inclusive, and between columns 200 inclusive. It any one of them is not, then something is wrong. You can't stuff a 100 by 200 image into a space that's 99 by 199 or whatever.

[Jason]

BI is 25x25 which are the sub images. The full image is 101x101 and consists of 4 BI's in each dimension.

OK, so 4x25=100, so there must be 1 pixel missing??

[Jason]

Just to add, its falling over on ii=4 in the first loop:

for jj=1:4
  for ii=1:4  
......

To help debug, I print out this:

p=[ct,startingRows(jj),endingRows(jj),startingCols(ii),endingCols(ii)]   %print out info for debugging

to get:

p =

             4           1          25          76         101

Subscripted assignment dimension mismatch.

I predefine:

Binary = zeros(rows, cols); = zeros(101,101)

[Image Analyst]

Jason's reply moved from an Answer to here:

I think I have identified the problem (not sur ehow to fix it): the definition of the start and end points for the grid to split up the image (in this case into 4x4 subimages , as per I.A help is:

[rows,cols]=size(ROI) % ROI is my Image
tiles = 4  % Number of tiles or chunks you want it divide up into.
% Find out where the dividing lines are.  This may be a fractional value.
startingRows = linspace(1, rows, tiles+1);
startingCols = linspace(1, cols, tiles+1);
% Now get the starting integer rows / cols.
startingRows = int32(startingRows(1:end-1))
startingCols = int32(startingCols(1:end-1))
% Get the ending rows.
endingRows = [startingRows(2:end)-1, rows]
endingCols = [startingCols(2:end)-1, cols]

%Determine number of pixels in each tile;
nx=endingRows(1)-startingRows(1);
ny=endingCols(1)-startingCols(1);

If I look at the output of the end points, the last entry should be 100 and not 101? (My raw image size is 101 x 101.

endingRows =

            25          50          75         101

endingCols =

            25          50          75         101

[Image Analyst]

If your image is 101-by-101, then why do you want to process only out to 100-by-100? Why don't you want to include the last row and column?

[Jason]

I do, but due to the error above when trying to reconstruct, thought it was easier. My wish is to include all pixels.

Binary(startingCols(ii):endingCols(ii), startingRows(jj):endingRows(jj))=BI(:,:,ct);

[Jason]

If I manually enter the values for the 1st four subimages:

for jj=1:tiles
  for ii=1:tiles  
      ct=ct+1;
      p=[ct,startingRows(jj),endingRows(jj),startingCols(ii),endingCols(ii)]   %print out info for debugging
      %Binary(startingCols(ii):endingCols(ii),startingRows(jj):endingRows(jj))=BI(:,:,ct);
      Binary(1:25,1:25)=BI(:,:,1);
      Binary(1:25,26:50)=BI(:,:,2);
      Binary(1:25,51:75)=BI(:,:,3);
      Binary(1:25,76:101)=BI(:,:,4);

    end
  end

it fails over on the 4th one.

[Jason]

This is my code for stacking the subimages into the 3rd dimension (i.e. an image stack)

tiles=4;    %image chunks
nimages=tiles*tiles
%Determine number of pixels in each tile;
nx=endingRows(1)-startingRows(1);
ny=endingCols(1)-startingCols(1);

 %Initialise Image and binary stack (IM)

IM=[]; %Initialise Image stack array 
BI=[]; %Initialise Binary array
subplot(2,4,[1,2,5,6])
hold on
ct=0;    %counter
for j=1:tiles
   for i=1:tiles 
    ct=ct+1;
    %Plot tile boundary lines on raw image
    plot([startingRows(i),startingRows(i)],[0,cols],'y-');  %plot tile lines
    plot([0,rows],[startingCols(i),startingCols(i)],'y-');  %plot tile lines

      %Add region number
      xr=double(floor(endingCols(i)-ny/2));
      yr=double(floor(endingRows(j)-nx/2));
      text(xr,yr,num2str(ct),'FontSize', 20, 'Color', [0 1 0]);

      %Divide Image into tiles
      IM1=imcrop(ROI,[startingRows(i) startingCols(j) nx ny ]);
      [thresh1, thresh2] = Threshold(hObject, eventdata, handles,IM1);

      %Create Binary Image
      ix=IM1;
      level=thresh2;
      ix(ix<level)=0;
      ix(ix>level+1)=1;
      BI1=ix;

      %Concatenate to image stack
      IM=cat(3,IM,IM1);   %add along 3rd dimension
      BI=cat(3,BI,BI1);   %create binary image stack

     end; 
  end
  hold off;

[Image Analyst]

Your 4th line has 26 columns instead of 25 like the others.

What does this say:

size(BI)

[Jason]

I will have to try tomorrow. I did notice that tnd the width and height i use in imcrop to get the subregions are defined by

nx=endingRows(1)-startingRows(1); ny=endingCols(1)-startingCols(1);

Which will be different for the 4th line of the loop

[Image Analyst]

This will happen if the number of rows or columns is not an integer multiple of the number of tiles you want to divide it up into. So for 4 tiles, you can have 100 or 104 or 108, but not 101 or 102 or 103, etc. You can still chop it up and stitch them together you just can't do it like you're doing. You should do variations of this

tallMatrix = [m1; m2];  % m1 and m2 have the same number of columns.
wideMatrix = [m1, m2];  % m1 and m2 have the same number of rows.

[Jason]

Thanks for your help, it seems like the Blockproc approach is easier and caters for all posiblities, although I understand the above apporach better, just can't get it to work.

You asked what

size(BI) was:

ans =

      25    25    16

[Image Analyst]

See attached blockproc demos. And don't forget to look at Mohammad's last comment under his answer.

[Jason]

Thankyou

Answer 4

[Jason]

Why does this approach not work:

ct=0;
for ii=1:4
 for jj=1:4
   ct=ct+1
   Binary(startingCols(jj):endingCols(jj), startingRows(ii):endingRows(ii))=BI(:,:,ct);
 end
end

imshow(Binary,[0,1])

[Iain]

That should work provided you calculate startingcols etc. correctly.

[Image Analyst]

The first index should be rows, not columns as you have it, and the second index should be columns, not rows as you have it.

Answer 5

[Alessandro Masullo]

What about using mat2cell and cell2mat? It should be much easier.

[Jason]

Hi, I don't see how this can reconstruct a full image from the subimages

Answer 6

[Jason]

I've managed to get a solution where I first have to crop the image to ensure there number of pixels corresponding to integer tiling: (100 in my case). its not ideal, but it seems a big struggle to get the non integer tiles as well.

[Image Analyst]

Why did you even split up your image in the first place?

[Jason]

local thresholding worked a lot better than global!

[Image Analyst]

OK but you don't need to split apart your image to do that. You can get a better local thresholding using adapthisteq() to flatten your image and then use a global threshold, or use blockproc like Mohammad suggested. adapthisteq is like splitting your image apart into tons of tiles that are only a pixel apart and will be much better and more accurate than splitting your image into 4 tiles and then computing the threshold for the whole tile. You should really look into these methods.

Answer 7

[KHOR WEI KOK]

Hi, I would like to ask, how you manage to show the line segmentation and numbering on your image axes?