2

I have a MxNxD volume and I need a lxhxw sliding window that goes through all the voxels of the volume. In each sliding window I need to compute the Root Mean Square Contrast. Which is the smarter way to do it?

I would like to limit the use of for loops because the volume is pretty big, 1024x1024x146.

BugsFree
  • 540
  • 1
  • 6
  • 24
  • What operations do you want to do on these voxels? Just to be clear, by sliding window you mean that you want to move the window by one row (or column/plane) for each operation such that it will overlap the previous window. – beaker Jun 10 '15 at 22:45
  • I have to compute the Root Mean Square Contrast [http://en.wikipedia.org/wiki/Contrast_%28vision%29] for each window. Yes, each time the window shifts of one voxel in one of the three directions, thus the next window overlaps the previous. – BugsFree Jun 10 '15 at 22:59
  • Does the 3D window happen to be [multiplicatively separable](http://calculus.subwiki.org/wiki/Multiplicatively_separable_function)? – Luis Mendo Jun 11 '15 at 06:59
  • I don't think that in this specific case my problem is multiplicatively separable‌​ but this can be useful for me in other cases. Only to be sure to have understood correctly, if my function was multiplicatively separable‌​ I could apply the function separately in the x,y,z directions and then sum together the result (e.g. to compute the gradient magnitude I can apply the window first along the rows and then along the columns)? – BugsFree Jun 11 '15 at 08:28
  • Did you consider using `convn`? http://www.mathworks.com/help/matlab/ref/convn.html – rayryeng Jun 11 '15 at 13:56
  • I don't think it is possible to compute the RMS contrast with a convolution, I thought about that but I didn't find the solution otherwise `convn` would be great. – BugsFree Jun 11 '15 at 21:41
  • What is the RMS contrast? I didn't know that's what you were trying to compute. Your problem statement doesn't make reference to this in any way. Please update your question to illuminate what exactly you're doing within the sliding windows. – rayryeng Jun 12 '15 at 19:27

2 Answers2

0

It sounds like you want to divide into voxels that are l x h x w. If this is the case and assuming your matrix is called q where size(q) = [M, N, D].

function N = test(q, l, h, w)
    qPrime = mat2cell(q, ...
        [l*ones(1, floor(size(q,1)/l)), mod(size(q,1), l)], ...
        [h*ones(1, floor(size(q,2)/h)), mod(size(q,2), h)], ...
        [w*ones(1, floor(size(q,3)/w)), mod(size(q,3), w)]);

    N = cellfun(@RMS, qPrime, 'uni', 0);
end

function N = RMS(M)
    var = M - mean(M(:));
    N = sqrt(sum(var(:).^2)) ./ numel(M);
end

This makes a cell array with each voxel in a cell. Apply the function to each voxel using cellfun(@foo, qPrime, 'uni', 0).

user1543042
  • 3,422
  • 1
  • 17
  • 31
  • Thank you for your answer, but this is not what I want to do. Your solution works great to split the volume, but I need a 3d convolution between the `l`x`h`x`q` window and the volume and each time the central voxel of the window gets the value computed using the root mean square contrast. – BugsFree Jun 11 '15 at 09:20
  • Isn't the window `l` x `h` x `w`? And what do you mean by a 3D convolution? That doesn't seem at all like the RMS contrast. – user1543042 Jun 11 '15 at 11:47
  • I didn't explain well. I need to create a 3D map where each voxel contains the value of the RMS contrast compute over a neighbourhood of `l`x`h`x`w`. Theorically the sliding window have to shift of a voxel for each iteration but this is what I would like to avoid because in that case `for` loops would be repeat more than 150 millions of times. – BugsFree Jun 11 '15 at 21:28
  • So you want the output to be `M` x `N` x `D`? – user1543042 Jun 11 '15 at 22:13
  • Yes, exactly, I need a `M`x`N`x`D` output matrix! – BugsFree Jun 13 '15 at 11:12
0

I have been able to solve this specific problem using the stdfilt function of MATLAB because the Root Mean Square Contrast is the standard deviation of the pixel intensities inside the sliding window. The stdfiltfunction allows you to define a

multidimensional array of zeros and ones where the nonzero elements specify the neighbors. (MATLAB Help)

BugsFree
  • 540
  • 1
  • 6
  • 24