Is there a way to parallelize the convolution component of EdgeDetect?

Question

Provided an image like -

test = Import["http://upload.wikimedia.org/wikipedia/commons/d/d5/Sunflowers.jpg"]

We can run EdgeDetect (with the "Canny", "ShenCastan", "Sobel" options) to convolve a kernel of some form with an image to approximate first-order derivatives that indicate borders or significant and sudden variations in an image.

I was wondering if there exists a mechanism for parallelizing EdgeDetect over multiple cores? Certainly it seems like this should be possible, since the "heavy lifting" computationally speaking is more or less just a convolution process? Perhaps we could decompose the input image into smaller chunks and then perform some fast "stitching" and post-processing operation?

Update: Cormullion (in the comments) make the good point that one can decompose an image in as many fragments as you have cores using ImagePartition, run EdgeDetect on each fragment, then stitch the resultant EdgeDetect transformed image back together using ImageAssemble. I suppose one would use ParallelTable for this?

The trouble, as cormullion notes, is dealing with the seams of the reconstituted / assembled image.

I've also noted that ParallelTable doesn't really offer a great speedup, especially for smaller images (I use a larger image on purpose in the below example). Is there a better way of proceeding, or perhaps storing the image in memory?

---

Let's step through cormullion's suggestion for a cropped version of the above image. I have twelve cores, so I'm going to try to crop and chop the image up into twelve fragments, but it should be clear how to adjust this according to the number of cores on your machine.

First, let's grab the image and crop it:

testImage = Import["http://upload.wikimedia.org/wikipedia/commons/d/d5/Sunflowers.jpg"]
testImage = ImageCrop[testImage, {2004, 2004}]

Our variable test is now a {2004, 2004} pixel crop of the original imported image. This lets us break it up into twelve {501, 668} pixel fragments. Remembering that the matrix will be $4 \times 3$, we can flatten the output of ImagePartition:

imageFragmentSet = Flatten[ImagePartition[testImage, {501, 668}]]

Now we can compare:

flattenedImageSet = ParallelTable[EdgeDetect[imageFragmentSet[[i]]], {i, 1, Length[imageFragmentSet]}] // AbsoluteTiming

With:

EdgeDetect[testImage] // AbsoluteTiming

The result is something like a $\approx 2$ to $\approx 3$ fold speedup. Can this be better optimized?

We can now stitch the image together like so:

reconEDImage = ImageAssemble[Partition[flattenedImageSet, 4]]

Note that you need to put a [[2]] in front of flattenedImageSet if you neglect to remove the // AbsoluteTiming command.

The problem now is - how does one fix the seams between the reassembled image fragments in order to match the output of EdgeDetect[testImage]?

Answer 1

It is already the case that the "convolution component of EdgeDetect is multithreaded."

EdgeDetect is composed of a sequence of operations. Some of these operations are indeed multi-threaded, and some are not (because, in short, as the OP says it's tricky).

The CPU usage screenshot is a bit misleading in that single threaded operations tend to last longer in wall clock time than the multi-threaded ones.

Answer 2

This is actually not the answer, but some investigation (I hope helpful).

As in Simon Woods's comments, I thought that EdgeDetect uses multithreading. But actual CPU load is

It is almost always equal to 25%, which corresponds to 1 of my 4 cores. So parallelization may be useful. But there are always parallel overheads. May be you have 6-cores CPU with hyperthreading? For computations hyperthreading gives much less then 2x speedup. So I think that your 3x speedup is great achievement.

Another approach is more rough algorithm:

EdgeDetect[test]; // AbsoluteTiming // First

1.451953

ImageConvolve[test, {{-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1}}]; // AbsoluteTiming // First

0.276083

Or using CUDA:

data = CUDAMemoryLoad@ImageData[test];
edges = CUDAMemoryAllocate[Real, Dimensions@ImageData[test]];

edges = CUDAImageConvolve[data, {{-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1}}]; // 
  AbsoluteTiming // First

result = Image@CUDAMemoryGet[edges];

0.026881

Note, that loading data to CUDA memory requires significant time (~0.2 s).