Running out of memory when partitioning an image with ImagePartition

Question

I have an Image3D object of image type "Bit16" with the dimensions 1472 x 1472 x 1472 that I would like to split into small image cubes of dimension 64 x 64 x 64 which would give me a total of 23^3 = 12167 small Image3D objects. I want to apply an adaptive thresholding algorithm I already implemented which computes separate thresholds for each image block and then interpolates the computed threshold values (adaptive thresholding). For this purpose I use the function ImagePartition. The following code illustrates the partitioning:

croppped = ImageCrop[image, {1472, 1472, 1472}];
parts = ImagePartition[cropped, {64, 64, 64}];

I use ImageCrop to add a border to the original image which results in the mentioned Image3D object of dimensions 1472 x 1472 x 1472. If I try to run the partitioning my Mathematica kernel crashes with an out of memory error. I have 96 GB RAM installed on my machine. I'm not sure, but the memory should be a sufficient amount of memory available to do the partitioning.

What surprised me is that i can create multiple huge Image3D objects and the kernel doesn't crash, but when I partition an Image3D, using ImagePartition it does. See the following code for the memory usage when I create another Image3D object:

MemoryInUse[]
cropped2 = ImageCrop[image, {1472, 1472, 1472}];
MemoryInUse[]

14461079560

20840097904

Questions:

1. Is there a simple explanation for this problem?
2. What would be a suitable workaround in order to perform the splitting?

Edit (08/01/2014)

First of all thanks very much UDB for this nice solution you posted. I now had the time to try out the function MyImagePartition on my data and here are the results I obtained so far. First, to get a feeling for the timings and memory usage I created the same test volume and partitioned it on my machine. I then checked whether the original volume can be assembled again to give the volume again.

volume = Image3D[RandomInteger[{0, 255}, {8, 9, 10}*64], "Byte"];
First@AbsoluteTiming[ip = ImagePartition[volume, {64, 64, 64}]]

10.899000

MaxMemoryUsed[]    
First@AbsoluteTiming[mip = MyImagePartition[volume, {64, 64, 64}]]
MaxMemoryUsed[]

248907056 1.505000 542529936

Image3D[ImageAssemble@ip] == volume

True

Image3D[ImageAssemble@mip] == volume

True

Then I tested the function MyImagePartition on a real image I have. In between the partitioning of the test volume and the real data i quit the kernel to get a fresh Mathematica session. Fortunately, the function is very fast in partitioning my data and memory consumption seems to be no issue. When reassembling it again I checked for the image dimensions of the assembled image.

ImageDimensions@testImage

{458, 403, 202}

MaxMemoryUsed[]
First@AbsoluteTiming[mip = MyImagePartition[testImage, {64, 64, 64}]]
MaxMemoryUsed[]

249027112 0.278000 542635024

ImageDimensions@Image3D[ImageAssemble@mip]

{448, 384, 192}

Obviously the dimensions of the reassembled image do not match the original image dimensions. I have to admit that I didn't fully understand the usage of the second parameter 'dwdh'. Is it the option for image padding such that the image matches up with the partitioning block size? I tried several options but failed continuously ;) Besides that issue, would it be possible to add the functionality to create an overlapping partitioning of the image (e.g. give the overlap of the blocks in voxels as an additional parameter). With the built-in function ImagePartition this works quite good as follows:

parts = ImagePartition[spheroid3d, {64, 64, 64}, {60, 60, 60}, 
   Padding -> 0];

Answer 1

MyImagePartition

In the meantime, before Mathematica 10 will come out, you can enjoy my on-foot solution MyImagePartition, which both saves memory and time using the PartitionMap function from the Developer context:

MyImagePartition[im_, wh_, dwdh_List: {0, 0, 0}] := 
 Module[{it = ImageType@im,
         cs = First@Options[im, ColorSpace],
         il = First@Options[im, Interleaving],
         is = First@Options[im, ImageSize]},
  If[! (ImageQ@im), Return@$Failed];
  If[Length@ImageDimensions@im == 2,
   Return@Developer`PartitionMap[
     Image[#, it, cs, il, is] &,
     ImageData[im, it],
     If[Length[wh] == 0, {wh, wh}, Reverse@wh],
     If[dwdh == {0, 0, 0}, Reverse@If[Length[wh] == 0, {wh, wh}, wh], Reverse@dwdh]
   ],
   Return@Developer`PartitionMap[
     Image3D[#, it, cs, il, is] &,
     ImageData[im, it],
     If[Length[wh] == 0, {wh, wh, wh}, Reverse@wh],
     If[dwdh == {0, 0, 0}, Reverse@If[Length[wh] == 0, {wh, wh, wh}, wh], Reverse@dwdh]
   ]
  ]
 ]

I think, the function above, though in part written rather semi-professionally, will work in analogy (sorry, the width and height specification only accepts fixed tile sizes so far) to the built-in ImagePartion, both for Image and Image3D data (for all image data types, such as "Byte", "Real", "Bit16" and so on), and also both for GrayLevel and RGBColor samples (including alpha channels), whereas all other colorspaces like "LUV", "CMYK" etc. work as well.

Since my 8GB RAM equipped machine (with 9.0 for Mac OS X x86 (64-bit) (January 24, 2013)) does not allow to make tests with 1472^3 volume "Bit16" data, instead I take some reduced "Byte"-valent test data set (so that no swapping will occur):

volume = Image3D[RandomInteger[{0, 255}, {8, 9, 10}*64], "Byte"];

Let us see how the built-in function works:

First@AbsoluteTiming[ip = ImagePartition[volume, {64, 64, 64}]]

24.633489

Now let's compare with my work-around:

First@AbsoluteTiming[mip = MyImagePartition[volume, {64, 64, 64}]]

3.923255

mip == ip

True

I am pretty sure your 96GB will be far enough to let you partition your 3189506048 "Bit16"-valent voxels into 64^3 blocks using MyImagePartition. Enjoy!

Addition

I have just checked the maximum required bytes for both functions using MaxMemoryUsed[] (in two separate kernel sessions):

For ImagePartition I obtained 6686483552, while for MyImagePartition I got just 1922143136...

Extension: Padding option added

As requested by g3kk0 I have extended the compatibility to the original ImagePartition function by adding a Padding option. Here is the code:

Options[MyImagePartition] = {Padding -> None};
MyImagePartition[im_, wh_, dwdh_List: {0, 0, 0}, OptionsPattern[]] := 
 Module[{it = ImageType@im, icf, ics = First@Options[im, ColorSpace], 
   iil = First@Options[im, Interleaving], 
   is = First@Options[im, ImageSize], ip = OptionValue[Padding]},
  If[is === (ImageSize -> Automatic), is = ImageSize -> All];
  If[! (ImageQ@im), Return@$Failed];
  If[Length@ImageDimensions@im == 2,
   Return@Developer`PartitionMap[
     Image[#, it, ics, iil, is] &,
     ImageData[If[ip === None,
       im,
       ImageReflect[
        ImageReflect[
         ImageCrop[ImageReflect[ImageReflect[im, Top], Left], 
          MapThread[(Ceiling[#1/#3] - 1)*#3 + #2 &, {ImageDimensions@
             im, If[Length[wh] == 0, {wh, wh}, wh], 
            If[dwdh == {0, 0, 0}, If[Length[wh] == 0, {wh, wh}, wh], 
             dwdh]}], Padding -> ip], Top], Left]
       ],
      it],
     If[Length[wh] == 0, {wh, wh}, Reverse@wh],
     If[dwdh == {0, 0, 0}, Reverse@If[Length[wh] == 0, {wh, wh}, wh], 
      Reverse@dwdh]
     ],
   icf = First@Options[im, ColorFunction];
   Return@Developer`PartitionMap[
     Image3D[#, it, icf, ics, iil, is] &,
     ImageData[If[ip === None,
       im,
       ImageReflect[
        ImageReflect[
         ImageReflect[
          ImageCrop[
           ImageReflect[ImageReflect[ImageReflect[im, Top], Left], 
            Front], 
           MapThread[(Ceiling[#1/#3] - 1)*#3 + #2 &, {ImageDimensions@
              im, If[Length[wh] == 0, {wh, wh, wh}, wh], 
             If[dwdh == {0, 0, 0}, 
              If[Length[wh] == 0, {wh, wh, wh}, wh], dwdh]}], 
           Padding -> ip], Top], Left], Front]
       ],
      it],
     If[Length[wh] == 0, {wh, wh, wh}, Reverse@wh],
     If[dwdh == {0, 0, 0}, 
      Reverse@If[Length[wh] == 0, {wh, wh, wh}, wh], Reverse@dwdh]
     ]
   ]
  ]

All known settings for Padding can be used. I did some effort to get identical results compared to ImagePartition, therefore I added some ImageReflect calls (forth and back), since otherwise due to some different rounding inside ImageCrop (used for the padding) under certain circumstances MyImagePartition vs. ImagePartition results would have slightly differed. However, this might be considered as somehow pendatic, and of course is time consuming (remove these ImageReflect calls if you want)...

What still is not implemented in MyImagePartition is any flexible tile width/height specification.

Enjoy!

Answer 2

This is not a direct answer to your question but it can help you. I see your question is about adaptive thresholding. I propose finding threshold values without exact partitioning.

img = Import["http://homepages.inf.ed.ac.uk/rbf/HIPR2/images/son1.gif"]

The simplest is GaussianFilter which is analog to $T = mean$ in your link.

GF = GaussianFilter[img, 30]

Binarize[ImageSubtract[ImageAdd[img, 0.03], GF], 0]

Here 0.03 is the adjustable parameter.

Instead of Gaussian filter you can use

MM = ImageAdd[ImageMultiply[MinFilter[img, 20], 0.5],ImageMultiply[MaxFilter[img, 20], 0.5]]

which is analog to $T=\dfrac{min+max}{2}$.

There are also MedianFilter but unfortunately its time rapidly grows with the sample size

SetAttributes[timeAvg, HoldFirst]
timeAvg[func_] := 
 Do[If[# > 1, Return[#/5^i]] & @@ AbsoluteTiming@Do[func, {5^i}], {i, 0, 15}]

Monitor[ListLogPlot[
  Transpose@
   Table[{timeAvg@GaussianFilter[img, n], timeAvg@MinFilter[img, n], 
     timeAvg@MedianFilter[img, n]}, {n, 10}], Joined -> True, 
  PlotRange -> {0.01, All}], n]

Answer 3

In Mathematica 9, ImagePartition uses ImageData and Partition so it is slow and requires a lot of memory.
This issue will be fixed in Mathematica 10.

Answer 4

Here's a way to do the adaptive thresholding using the built-in ImageFilter function. To apply this to the "Sonnet for Lena" as in ybeltukov's answer, define an auxiliary function that does the local thresholding:

adaptThresh[x_, threshRat_] := Module[{center, thisPix}, 
   center = (First[Dimensions[x]] + 1)/2;
   thisPix = x[[center, center]];
   Boole[thisPix > threshRat Mean[Flatten[x]]]];

x is the local portion of the image and the output is 1 whenever threshRat times the Mean of the portion is less than the pixel value. To apply this to the image

img = Import["http://homepages.inf.ed.ac.uk/rbf/HIPR2/images/son1.gif"];
ImageFilter[adaptThresh[#, 0.95] &, img, 2]

Since your file is so large, it would probably make sense to use ImageFileFilter which operates directly on the file (rather than reading everything into memory).