Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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];
share|improve this question
9  
Congratulations, you posted the 10,000th question. –  Sjoerd C. de Vries Sep 13 '13 at 14:23

4 Answers 4

up vote 14 down vote accepted

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!

share|improve this answer
    
Really nice answer, thanks! –  g3kk0 Nov 6 '13 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 '13 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 '13 at 19:25
    
I added some information on the application of the function on my data in the question above. –  g3kk0 Jan 8 at 15:23
    
@g3kk0: Hope the extended version fits your needs. Generally, I think the images should be cropped prior to partitioning in a way so that no padding is required... –  UDB Jan 15 at 22:26

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"]

enter image description here

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

GF = GaussianFilter[img, 30]

enter image description here

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

enter image description here

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]

enter image description here

share|improve this answer
    
Thanks for your interesting answer. I will definitely consider that for my problem and check what results I get. –  g3kk0 Sep 13 '13 at 19:16
1  
+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 '13 at 6:39

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.

share|improve this answer

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]

enter image description here

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).

share|improve this answer
    
Does the time of this method grows fast with the radius of the filter? –  ybeltukov Jan 15 at 23:24
    
@ybeltukov -- Yes, because the filter is n by n (where n is the parameter given to the filter). On the other hand, as long as the underlying changes are smooth, it shouldn't be necessary to use a large n (I used n=2 for the above example). –  bill s Jan 15 at 23:44
    
It was my idea half a year ago, but I reject it because of bad scalability. It will be quite slow for OP's 1472^3 data. Anyway, +1 for implementing this! –  ybeltukov Jan 15 at 23:53
    
@bills: Thanks for the nice answer. Since my data is quite large, I will try it with the ImageFileFilter idea you proposed and see what comes out. I will add the results to the question. –  g3kk0 Jan 16 at 10:52

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.