2

Data: Say I have a 2000 rows by 500 column matrix (image)

What I need: Compute the FFT of 64 rows by 10 column chunks of above data. In other words, I want to compute the FFT of 64X10 window that is run across the entire data matrix. The FFT result is used to compute a scalar value (say peak amplitude frequency) which is used to create a new "FFT value" image.

Now, I need the final FFT image to be the same size as the original data (2000 X 500).

What is the fastest way to accomplish this in MATLAB? I am currently using for loops which is relatively slow. Also I use interpolation to size up the final image to the original data size.

Eitan T
  • 32,660
  • 14
  • 72
  • 109
hkf
  • 663
  • 2
  • 11
  • 23

2 Answers2

2

As @EitanT pointed out, you can use blockproc for batch block processing of an image J. However you should define your function handle as

fun = @(block_struct) fft2(block_struct.data);
B = blockproc(J, [64 10], fun);

For a [2000 x 500] matrix this will give you a [2000 x 500] output of complex Fourier values, evaluated at sub-sampled pixel locations with a local support (size of the input to FFT) of [64 x 10]. Now, to replace those values with a single, e.g. with the peak log-magnitude, you can further specify 

fun = @(block_struct) max(max(log(abs(fft2(block_struct.data)))));
B = blockproc(J, [64 10], fun);

The output then is a [2000/64 x 500/10] output of block-patch values, which you can resize by nearest-neighbor interpolation (or something else for smoother versions) to the desired [2000 x 500] original size 

C = imresize(B, [2000 500], 'nearest');

I can include a real image example if it will further help.

Update: To get overlapping blocks you can use the 'Bordersize' option of blockproc by setting the overlap [V H] such that the final windows of size [M + 2*V, N + 2*H] will still be [64, 10] in size. Example: 

fun = @(block_struct) log(abs(fft2(block_struct.data)));
V = 16; H = 3; % overlap values
overlap = [V H]; 
M = 32; N = 4; % non-overlapping values
B1 = blockproc(J, [M N], fun, 'BorderSize', overlap); % final windows are 64 x 10

However, this will work with keeping the full Fourier response, not the single-value version with max(max()) above.

See also this post for filtering using blockproc:Dealing with “Really Big” Images: Block Processing.

gevang
  • 4,994
  • 25
  • 33
  • 1
    What is the point of `block_struct.data`? Why do you assume it is stored in this way? – Eitan T Sep 11 '12 at 00:09
  • 2
    quoting from http://www.mathworks.com/help/toolbox/images/ref/blockproc.html "A block struct is a MATLAB structure that contains the block data as well as other information about the block." You can check the examples also. In R2012a a handle such as @fft will not work, input should be struct. The function handle to `blockproc` should be "a function that accepts a block struct as input and returns a matrix, vector, or scalar." – gevang Sep 11 '12 at 00:25
  • Interesting. I thought `blockproc` behaved just like the (now obsolete) `blkproc`... thanks for pointing that out then. I'd better correct my answer too. – Eitan T Sep 11 '12 at 00:29
  • no problem, I was under the impression that it worked by default through @fun also. – gevang Sep 11 '12 at 00:31
  • Thank you for the response gevang and EitanT. Is there a way for the blocks to overlap? I don't understand how filter2 can be used. Thanks! – hkf Sep 12 '12 at 17:57
  • See the Update in my answer for a simple way to get overlapping windows, thus smoother boundaries between the FFT patches. – gevang Sep 12 '12 at 18:41
  • Unfortunately, `blockproc` doesn't "slide" the window, it just `inflates` the blocks' borderes so that they overlap, _but their amount remains the same_. I'm not sure that was the OP's original intention. – Eitan T Sep 12 '12 at 18:48
  • @EitanT true. What do you mean by "their amount" though? Note that you can change the sliding parameter through the final window length, i.e. you have 2 parameters final length and overlap amount. These can be translated to window center and window length and this way emulate a sliding window effect. – gevang Sep 12 '12 at 19:15
  • For this purpose I believe `blockproc` would be significantly slow. I suggest using `colfilt` then, which does not perform unnecessary padding. Please see my answer for reference. – Eitan T Sep 12 '12 at 19:21
1

If you want to apply the same function (in your case, the 2-D Fourier transform) on individual distinct blocks in a larger matrix, you can do that with the blkproc function, which is replaced in newer MATLAB releases by blockproc.

However, I infer that you wish to apply apply fft2 on overlapping blocks in a "sliding window" fashion. For this purpose you can use colfilt with the 'sliding' option. Note that the function that we're applying on each block is the fft:

block_size = [64, 10];
temp_size = 5 * block_size;
col_func = @(x)cellfun(@(y)max(max(abs(fft2(y)))), num2cell(x, 1), 'Un', 0);
B = colfilt(A, block_size, 10 * block_size, 'sliding', col_func);

How does this work? colfilt processes the matrix A by rearranging each "sliding" block into a separate column of a new temporary matrix, and then applying the col_func to this new matrix. col_func in turn restores each column into the original block and applies fft2 on it, returning the largest amplitude value for each column.

Important things to note:

  1. Since this mentioned temporary matrix includes all possible "sliding" blocks, memory could be a limitation. Therefore, in order to use less memory in calculations, colfilt breaks up the original matrix A into sub-matrices of temp_size, and performs calculations separately on each. The resulting matrix B is still the same, of course.

  2. Each element in the resulting matrix B is computed from the corresponding block neighborhood. The larger your image is, the more blocks you will need to process, so the computation time will increase geometrically. I believe that you'll have to wait quite a bit until MATLAB finishes processing all sliding windows on your 2000-by-500 matrix.

Eitan T
  • 32,660
  • 14
  • 72
  • 109
  • Thanks for the response EitanT! – hkf Sep 12 '12 at 17:53
  • I am having trouble with colfilt. I get the error: "Attempt to reference field of non-structure array." The same function works well with blockproc. Not sure what I am doing wrong. – hkf Sep 12 '12 at 21:49
  • Does your function use the `block_struct.data`? Unlike `blockproc`, `colfilt` does not need a _block_struct_. Define your function like I did, and I believe it should work. – Eitan T Sep 13 '12 at 07:15
  • The function that I have written requires the input data to be 64 X 10. Colfilt seems to change it to a column during the first process and then in to a matrix. And yes, I am not using the struct format as you suggest. – hkf Sep 13 '12 at 16:58
  • Yes, it seems that `func` is expected to treat each column of its input as an individual block. I need to wrap my head around it a bit and revise `func`. – Eitan T Sep 13 '12 at 17:07
  • Okay, I think it works now, but the processing now takes time because of all the "sliding" blocks that need to be processed. – Eitan T Sep 14 '12 at 02:10
  • 1
    Thanks. The code works but it is 2 orders of magnitude slower than the 'for' loops I had in place! – hkf Sep 14 '12 at 16:53
  • Does it produce the same results? – Eitan T Sep 14 '12 at 18:49
  • Yes. I ran it on a small test image. I am going to have to use blockproc as suggested by gevang. Thanks – hkf Sep 14 '12 at 18:57
  • No problem. Maybe I'll find a workaround, but use what works best for you – Eitan T Sep 14 '12 at 19:36