Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
imfilter speed

Subject: imfilter speed

From: Patrick

Date: 8 Jun, 2011 19:41:04

Message: 1 of 9

I have a 3D volumetric image and wanted to clean it up a bit. It doesn't seem the image processing toolbox can do what I need (consider the points in neighboring slices as well as in the current slice), so I rolled my own averaging algorithm. The problem is its really slow - it takes almost 30 minutes to do even on our really fast 64-bit machine. When using a similar imfilter averaging function on the same image though, it can be done in under a minute. Granted, for 2D filtering the matrix is only 2 dimensions but that doesn't seem like it should slow it downt that much as it still has to do all the slices. I was wondering how the image processing toolbox executes so quickly considering it has to do the same thing I'm doing - loop through every point in the image. Is there some quicker way to do it that I'm missing? Is there another way to do it besides a double for loop -
like some series of matrix multiplications?

Subject: imfilter speed

From: ImageAnalyst

Date: 8 Jun, 2011 21:11:10

Message: 2 of 9

Patrick:
I'd guess that it might use the method of separable kernels (two 1D
filters instead of a 2D filter) in the event that it detects a
separable kernel - that's a huge speed up when it can be applied.

And secondly, it probably also keeps track of what pixels are in the
window and doesn't get them again. After all, if you have a 7 by 7
window and move it over 1 pixel, 42 of those pixels are the same
(they're still inside the window even though it shifted) and only 7
have changed (you've lost the "old" 7 but gained 7 "new" pixels). Why
computer 42 more things than you have to?

Subject: imfilter speed

From: Patrick

Date: 8 Jun, 2011 21:43:04

Message: 3 of 9

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <ec392abe-777a-4e54-bc1d-66a333219047@k16g2000yqm.googlegroups.com>...
> Patrick:
> I'd guess that it might use the method of separable kernels (two 1D
> filters instead of a 2D filter) in the event that it detects a
> separable kernel - that's a huge speed up when it can be applied.
>
> And secondly, it probably also keeps track of what pixels are in the
> window and doesn't get them again. After all, if you have a 7 by 7
> window and move it over 1 pixel, 42 of those pixels are the same
> (they're still inside the window even though it shifted) and only 7
> have changed (you've lost the "old" 7 but gained 7 "new" pixels). Why
> computer 42 more things than you have to?

Hmm, even with a 1x1 matrix as a filter ( essentially an identity operation), its still really slow so it appears the main problem is not the redundant computations ... I think it may be the looping, which I I understand in MATLAB is comparatively slow.

Subject: imfilter speed

From: Patrick

Date: 8 Jun, 2011 21:45:05

Message: 4 of 9

Here is the code.

function avgData = octImgFilterAvg3D(imgData, frameIdx, radius)

sz = size(imgData);
imgWidth = sz(1);
imgHeight = sz(2);
imgFrames = sz(4);
avgData = zeros(imgWidth, imgHeight);

frIdx1 = max(frameIdx-radius, 1);
frIdx2 = min(frameIdx+radius, imgFrames);
for h=1:imgHeight
    hIdx1 = max(h-radius, 1);
    hIdx2 = min(h+radius, imgHeight);
    for w=1:imgWidth
        wIdx1 = max(w-radius, 1);
        wIdx2 = min(w+radius, imgWidth);
        avgData(w, h) = mean(mean(mean(imgData(wIdx1:wIdx2, hIdx1:hIdx2, 1, frIdx1:frIdx2))));
    end
end

Subject: imfilter speed

From: ImageAnalyst

Date: 8 Jun, 2011 22:42:25

Message: 5 of 9

Patrick:
Maybe it's the memory access or the order of the indexes. because I
ran this code

clc;
tic
for sliceZ = 1:300
for x = 1:1000
for y = 1:1000
r=1;
end
end
end
toc
and it ran in 1.3 seconds on a 5 year old laptop computer. A far cry
from 30 minutes. So it's not the looping per se that's the problem,
but maybe the memory access. Maybe you could try flipping the order
of the loops. I think MATLAB goes down rows first and then moves over
to the next column and goes down those rows, etc., so maybe you can
swap the first two indexes in your array and see if that's faster.

In other words, try
for index1 = 1:1000
for index2 = 1:1000
                            whatever = M(index1, index2, sliceZ);
and see if that's faster than
for index2 = 1:1000
for index1 = 1:1000
                            whatever = M(index1, index2, sliceZ);
because maybe you have them in the least efficient order.

By the way, you aren't writing anything to the command window (or
poking into spreadsheets) inside the loop are you? Because I think
that can slow things down.
ImageAnalyst

Subject: imfilter speed

From: Patrick

Date: 8 Jun, 2011 23:38:05

Message: 6 of 9

I'm not doing much inside that loop. I wrote a test script similiar to what I'm doing here and the killer appears to be:

        avgData(w, h) = mean(mean(mean(imgData(wIdx1:wIdx2, hIdx1:hIdx2, 1, frIdx1:frIdx2))));

Using reshape() and then only a single call to mean():

        numPts = (frIdx2-frIdx1+1)*(hIdx2-hIdx1+1)*(wIdx2-wIdx1+1);
        avgData(w, h) = mean(reshape(imgData(wIdx1:wIdx2, hIdx1:hIdx2, 1, frIdx1:frIdx2), numPts, 1));


This improves the speed by a factor of 3. Still not really close imfilter speed but a signifcant improvement. I might try the running average approach to avoid duplicate computations but not sure how much thats going to help as I am typically using small matrices (5x5) anyway, and the additional if/else statements that need to be used might add as much computation to the loop as naively computing the average for all points does.

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <ec9b93ce-c7b2-4b71-b2bf-13d9643bf9e8@g26g2000yqc.googlegroups.com>...
> Patrick:
> Maybe it's the memory access or the order of the indexes. because I
> ran this code
>
> clc;
> tic
> for sliceZ = 1:300
> for x = 1:1000
> for y = 1:1000
> r=1;
> end
> end
> end
> toc
> and it ran in 1.3 seconds on a 5 year old laptop computer. A far cry
> from 30 minutes. So it's not the looping per se that's the problem,
> but maybe the memory access. Maybe you could try flipping the order
> of the loops. I think MATLAB goes down rows first and then moves over
> to the next column and goes down those rows, etc., so maybe you can
> swap the first two indexes in your array and see if that's faster.
>
> In other words, try
> for index1 = 1:1000
> for index2 = 1:1000
> whatever = M(index1, index2, sliceZ);
> and see if that's faster than
> for index2 = 1:1000
> for index1 = 1:1000
> whatever = M(index1, index2, sliceZ);
> because maybe you have them in the least efficient order.
>
> By the way, you aren't writing anything to the command window (or
> poking into spreadsheets) inside the loop are you? Because I think
> that can slow things down.
> ImageAnalyst

Subject: imfilter speed

From: ImageAnalyst

Date: 9 Jun, 2011 02:28:58

Message: 7 of 9

You're simply taking the mean? Well talk about your separable
kernels! Heck, just do 2 scans of 1D convolutions, one in each 2D
direction, then do this for every slice, then do it in the third
dimension and you're done. So if your image is 640x480x300, you do
(640 + 480) * 300 convolutions to blur in the x,y lanes, then another
640*480 convolutions to average along the slice dimension. All told
you'd have about 337,120 1D convolutions. That sure would beat a 3D
convolution of 92 million voxels.

http://www.mathworks.com/matlabcentral/fileexchange/27957-separable-n-dimensional-convolution
http://blogs.mathworks.com/steve/2006/11/28/separable-convolution-part-2/
http://www.dspguide.com/ch24/3.htm
http://songho.ca/dsp/convolution/convolution2d_separable.html

Subject: imfilter speed

From: Rune Allnor

Date: 9 Jun, 2011 04:26:36

Message: 8 of 9

On Jun 8, 9:41 pm, "Patrick " <praph...@gmail.com> wrote:
> I have a 3D volumetric image and wanted to clean it up a bit.  It doesn't seem the image processing toolbox can do what I need (consider the points in neighboring slices as well as in the current slice), so I rolled my own averaging algorithm.   The problem is its really slow - it takes almost 30 minutes to do even on our really fast 64-bit machine.   When using a similar imfilter averaging function on the same image though, it can be done in under a minute.   Granted, for 2D filtering the matrix is only 2 dimensions but that doesn't seem like it should slow it downt that much as it still has to do all the slices.  I was wondering  how the image processing toolbox executes so quickly considering it has to do the same thing I'm doing - loop through every point in the image.   Is there some quicker way to do it that I'm missing?  Is there another way to do it besides a double for loop -
> like some series of matrix multiplications?

Large-volume data processing is *hard* to get to run fast.
For high-speed work one needs to understand the details of
how the computer hardware works, as well as how the software
interacts with that hardware. Then one needs a programming
language that is efficient to usea along with a compiler
that is able to optimize what you've done so far.

If you, by rolling your own matlab code, got withing 50x
run-time of what might be a compiled function, then your
efforts were very acceptable. Don't expect to get much
closer to that if you stick with matlab.

Rune

Subject: imfilter speed

From: Steve Eddins

Date: 9 Jun, 2011 12:02:18

Message: 9 of 9

On 6/8/2011 3:41 PM, Patrick wrote:
> I have a 3D volumetric image and wanted to clean it up a bit. It doesn't
> seem the image processing toolbox can do what I need (consider the
> points in neighboring slices as well as in the current slice), so I
> rolled my own averaging algorithm.
 > [snip]

imfilter can do three-dimensional filtering.

--
Steve Eddins
http://blogs.mathworks.com/steve/

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us