Regionprops handling of big images

Hello. I am processing a quite big image (120 GB) and need to find connected regions and calcularte some properties.

The "preprocessing" of the image works pretty fine. Here is the code:

bim = blockedImage("my_big_image.tif");
tmp = apply(bim,@(bs)processImgBlock(bs,filter_value));

And the function is like this

function bw = processImgBlock(bs,m)
         bw = bs.Data == m; % masked image
         bw = bwareaopen(bw,250);
         bw = bwconncomp(bw);
         bw = labelmatrix(bw);
end

So far so good. I get my label matrix ready to be processed. It is a 80k x 80k lines, single band image. According to documentation I will get as return a blocked image that I can retrieve using:

img = gather(tmp);

The issue is that if I run regionprops on img, Matlab does not process it due to memory. Note that I have 64GB in my computer.

So my next try that (quite obviously) failed was:

my_props = apply(tmp,@(bs)get_stats(bs,land_code));

% % % getstats is like follows

function s = get_stats(bs,land_code)
    
     s = regionprops('table',bs.Data,'Area','Perimeter','Eccentricity','Orientation','Centroid');
     
end

because the function should return a blocked image and not a table.

Here comes the question(s):

1. Is there any function that works similarly to apply that can perform the regionprops operations?
2. If not, any suggestion on how to handle it? Splitting image in small parts would be the obvious choice but it could split a connected region in 2 or more subimages and I'd like to avoid it (if possible).

Thank you in advance and please reach me if you need more information or if you find my explanation unclear.

Answer 1

[Walter Roberson]

Splitting image in small parts would be the obvious choice but it could split a connected region in 2 or more subimages

Any method that involves splitting the image into sub-images and processing the sub-images, runs the risk that you might miss connectivity that extends over sub-image boundaries.

blockedImage and apply and the older blockproc pass in partial images and assume the processing on the partial images can be handled independently. They do not transform the called function to be able to handle large arrays without running out of memory, and they have no idea how to synthesize the results of processing the sub-images together to form a complete result.

The sort of thing you would need would be to use regionprops with a tall array. Unfortunately, regionprops does not support tall arrays.

You are going to have to split the image up into sub-images, process that, get the results out, and then patch up the boundaries -- such as checking to see whether any connected components extend to the sub-image boundaries in adjacent sub-images, and if so then adjusting the information, such as re-labeling.

bw = bwconncomp(bw);

I am morally certain that is not going to return "correct" information when you use blocked images.

because the function should return a blocked image and not a table.

If you use the older blockproc() then it is not required that the output be the same size as the input block -- but if you use overlapping blocks, then make sure you use options to not trim the output. You have to return something numeric, but as long as you can make the output a consistent size then you could potentially pack data into a numeric array. This would not work very well for returning per-region information (unless you set a fixed maximum number to return), but it would work, for example, to return information about the largest region.

Answer 2

[Image Analyst]

How big is your bs.data binary image? I would start by trying to delete any variables the in workspace that you don't need anymore. Use the clear function. If you still don't have enough memory, try writing out the binary image to a disk file and then, in a separate program read in the binary file and get your properties. That should work since you won't have any other large variables hanging around in memory anymore.

[Marcio Teixeira]

Hello @Image Analyst!

First of all, thanks for your prompt reply. I will provide more details below.

1. My tifs are typically 80 to 120 MB. They are loaded in Matlab as 80000x80000 matrices I can load it, display (imshow), no problems there.
2. The label matrix (tmp) consists of 121 unique labels (or connected regions). Dimensions are 6400000000 bytes, class uint8. Again no issues displaying (I have 64GB of RAM on desktop and 32 MB in my notebook)
3. When you said (...) try writing out the binary image to a disk file and then, in a separate program read in the binary file and get your properties (...) you mean not Matlab?

Clearing vars helped. Bad news are it crashed (running in my notebook).

Following the error message:

----------------------------------------------------

Requested 80000x80000 (47.7GB) array exceeds maximum array size preference (31.7GB).

This might cause MATLAB to become unresponsive.

Error in regionprops>ComputePerimeterWithLabelMatrix (line 1193)

B = images.internal.builtins.regionboundaries(double(L),8);

Error in regionprops (line 358)

ComputePerimeterWithLabelMatrix(L,stats, ...

Related documentation

----------------------------------

Any workaround for this? Do you believe that running in my desktop (with 64 GB) I can avoid this error?

Thanks again!

[Marcio Teixeira]

Updating: running on desktop didn't work either. I got an "Out of memory"

[Image Analyst]

Those are big images. May I ask what they are pictures of, and what you are trying to measure in the image? Also, are you absolutely sure you can't get usable answers by using a downsampled image?

Anyway, looks like you got a solution from Walter that works for you.

[Marcio Teixeira]

Hello @Image Analyst !

Of course! I have images of land use and land cover. The map consists of land uses and each one has a code e.g. forest = 10, pasture = 20 etc. When I apply the maks in my code I am measuring the desired region. I am interested in perimeter, area, the centroids and orientation of these areas.

The solution from @Walter Roberson was the one I was trying to avoid: split the image. If it is the only possible path, I will go for it. Now there is a suggestion from by @Ashish Uthama. I will try it.

About downsampling I am not so sure. My idea is to collect statistics on thse areas, so it would be nice to keep them in the original 30 x 30 m per pixel.

cheers!

Márcio

[Ashish Uthama]

Depending on how large these regions are, downsampling might be worth trying. i.e if they occupy say ~>5-10% each of the full image, the resolution lost in area/perimeter might be small when compared to the full image. (its likely also a function of the details in the mask edges, if its mostly uniform/straight then the loss might be even less).

If you only have a few well defined regions, you could consider using blockedImage's zero-copy crop like this example shows: crop and mask geoTiff Files. Adapt the idea here like so: 1)create a low res image of the mask (retain the WorldStart/WorldEnd of the orginal!) Use regionprops on this to get bounding boxes for each object, then iterate on that to crop the finest level one by one to recompute stats on the finest level individually. This will work as long as your objects are well localized, so a crop is actually much smaller than the full image.

(For questions like these, it really helps to share sample images so some 'non standard' ideas can be formed :D)

Answer 3

[Ashish Uthama]

@Marcio Teixeira - what is the maximum size of an individual 'region' that you are trying to analyze? If its small enough w.r.t to the image the approach outlined in this example: Detect Nuclei in Large Whole Slide Images Using Cellpose ought to work. It shows how to use regionprops on a blockedImage using apply.

It shows how to use blockedImage to process 100kx100k+ images while ensuring that objects on the border are correctly captured only once.

In short:

1. Determine the maxium size of your individual object. This is the crucial part.
2. Process large image in blocks with additional border size of max size/2 at least (use apply())
3. In the processing function, reject objects whose centroids fall on the border region and only retain stats for objects whose centroids fall within the core block region

You ought to be able to modify the helpers in the linked example for your use case.

[Marcio Teixeira]

Hello @Ashish Uthama.

Definetely I can determine the maximum size of my object. This step is OK. Unfortunatelly I don't have the medical image toolbox in my Matlab distribution (I use one of my University).

Nevertheless I found a workaround that I detailed in my answer below.

Thank you for the suggestion and the very interesting example you provided.

Cheers,

Márcio

[Ashish Uthama]

Glad you have something that works for you!

Point of clarification: You only need Medical if you want to use cellpose for segmentation. For just blockedImage+regionprops, you only need Image processing Toolbox. So you could just replace the cellpose call in the above example with whatever you currently use for segmentation to adapt it. Or you could just study the code and modify it for your use.

Answer 4

[Marcio Teixeira]

Guys!

First of all thanks for your messages. I tried all your suggestions and I found one that is working nicely for me. After loading my image, I filtered it for my pixel of interest (let's say 20). With this binary image (bw) I do the following:

1. Instead of removing the areas with less than 250 pixels I selected the top 100 using a1 = bwareafilt(bw,100). I came to 100 images by trial and error. This was the limit that I could use regionprops without crashing.
2. I extracted the properties if interest and kept them in a table
3. I can use the a1 as a filter to clear the top 100 areas from my original image bw. Then I apply step one again and get the next patches (101 to 200).

Of course I need to clear memory every step and read again the new filtered image. It is not the most elegant solution, but it is working for me.

Again thank you all for valuable inputs.