# Running out of memory when partitioning an image with ImagePartition

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?
-
Congratulations, you posted the 10,000th question. –  Sjoerd C. de Vries Sep 13 at 14:23

## 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@DeveloperPartitionMap[ 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@DeveloperPartitionMap[ 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 full analogy 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... - Really nice answer, thanks! – g3kk0 Nov 6 at 19:01 @g3kk0: Please test your original big partitioning task and let me know what First@AbsoluteTiming@MyImagePartition as well as MaxMemoryUse[] report. – UDB Nov 6 at 19:10 I will try to do that. I have to say that I used another workaround instead, so I will have to first look for the code ;) – g3kk0 Nov 6 at 19:25 add comment 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]


-
Thanks for your interesting answer. I will definitely consider that for my problem and check what results I get. –  g3kk0 Sep 13 at 19:16
+1 I usually use a Closing operation instead of a GaussianFilter: It's robust, it removes the characters completely and it doesn't smooth the background brightness (as much). It's slower than GaussianFilter, but with a rectangular mask, processing time shouldn't depend on filter size, either. (Haven't tested that, though.) –  nikie Nov 6 at 6:39
In Mathematica 9, ImagePartition uses ImageData and Partition so it is slow and requires a lot of memory.