Gradient Domain Operations, take 5

[georg.viehoever]

Hi,

here the latest version of the GradientDomain module. The highlight:

- to get rid of the artefacts at stars located between two images, I have implemented a feathering approach.
- The new parameter "Feather Radius" controls the amount of feathering.
- This works very well, as demonstrated below. Only problem: If you have a large feather radius, the process can be slow. Currently, the radius therefore is limited to 20. Maybe Juan can advise how to make this faster.
- The feathering does not attempt to resolve the issues at the seam between image and background. I dont plan to handle these, removing a few pixels there should not be too painful.
- No important changes elsewhere.

The screenshot attached below shows:

- top row: MergeMosaic with Feather Radius 0 (=off). This is the old behaviour.
- bottom row: Feather Readius=10. The problematic star is nicely merged in. On the right you can see how the masks between the images are feathered.

Sources as usual tested on Fedora14-x64, but no significant problems expected for other platforms.

Enjoy,
Georg

[Harry page]

Hi

Excellent Stuff , now we are rocking   8)

Are the stars at the edge a real problem to sort , as it seems you are so near to perfect now

Harry

[Juan Conejero]

Hi Georg,

Really excellent results. I have attached a modified version of GradientMergeMosaic.cpp that uses FFTConvolution and SeparableConvolution to boost the feathering process. Instead of a LinearFilter I use a GaussianFilter, mainly because it is separable. If a Gaussian filter is not acceptable for some reason (which I doubt) we could use another separable 2D function, such as a box average filter. Both FFT-based convolution and separable convolution are O(N) processes, in contrast with regular convolution which is O(N^2), so you should see a very significant performance improvement.

Now the only bottleneck that remains to be fixed is the erosion filter that you apply before convolution. This is also O(N^2) so it is slow for large radii despite the fact that PCL uses a parallel implementation.  There are several approaches. If what you want is just to shrink the mask by a constant amount around its edges, then how about a simple crop? Cropping is O(N) so we could easily have our cake and eat it tonight Let me know.

So please test the new implementation and if it works as it should, I'll compile the whole thing for Linux and Windows. My modifications are enclosed by "BEGIN CODE BY J.C." and "END CODE BY J.C." comments. Let me know also if you want to replace erosion with cropping.

[georg.viehoever]

Wow, I just love this. Helpful responses with just 70 minutes. Incredible...

-Juan, thanks for the code 8). I'll test it this weekend. I am not sure if a gaussian will work. I need some operation that smoothes the edges of a region encoded by pixels with value==1 (this is just what you see in the mask on the lower right of the screenshot), making the pixels of the region go smoothly from 1 to 0 (background). The smoothed region must not be larger than the original one. The smoothed region is then used for feathering.
   I implemented this by moving the region in by n pixels via erosion, and then smoothly extending it again via convololution. The width of the smoothing needs to be reasonably controllable I even thought about wavelets, but I am not sure if this is the right way. Any operation that has the described characteristics would probably be ok.

-Juan, Crop() is not what I need. I need to move the edges of an arbitrarily shaped region by n pixels inwards. Erosion() just does this. What would be the class/process that you have in mind?

- Harry, I just do not have an elegant implementation idea for the border issue yet. I guess I will have to sleep on it a couple of nights 

Georg

[georg.viehoever]

Hi Juan,

is there an implementation of a Box Average Filter in PCL? The way I am using LinearFilter (v0=V1, all coefficients equal 1/n) in effect is a box filter. It is separable (see for example http://blogs.mathworks.com/steve/2006/11/28/separable-convolution-part-2/), but apparently  the implementation of IsSeparable() does not notice....

If BoxFilter is not available, how would I implement it? Looking at KernelFilter, my guess is that it should not be more than a dozen lines or so.

I also think that I can substitute the Erosion by a BoxFilter followed by Binarize(). Will experiment with it this weekend.

Thanks for the support.
Georg

[Juan Conejero]

is there an implementation of a Box Average Filter in PCL?

Sure:

Code: [Select]

/*
 * Create a box average filter of size 15x15
 */
KernelFilter B( 15, 1.0F ):
bool b = B.IsSeparable(); // b = true
However, you can define the equivalent separable filter directly:

Code: [Select]

/*
 * Create a separable filter equivalent to a box average filter of size 15x15
 */
SeparableFilter B( 15, 1.0F );

The way I am using LinearFilter (v0=V1, all coefficients equal 1/n) in effect is a box filter.

Nope. LinearFilter used that way defines a circular filter with constant value equal to one. This is not a box average filter and isn't separable.

I'm thinking in an O(N)  way to shrink the mask without erosion. Later I'll comment on this (when I sort my ideas )

[Carlos Milovic]

Since it is a binary image, a logical operation may do the trick: !(pix1 && pix2 && ... && pixN). Furthermore, if it is a black pixel, we already know the result, so there is no need to calculate again; only perform that for white pixels. Also, it could be coded with an loop exit, if it ever finds a black pixel...
Anyway, it is not an O(N^2) operation (N the number of pixels). It is O(N*M), M the size of the kernel.

Georg, I'm 99% sure that the problems with the stars at the boundaries will banish with a 3px zero padding, before the gradient calculation/merging. It will slow down the process a little bit, but is not a great deal. Furthemore, since we are using a solver based on FFT calculations (and I would say that this is one of the narrowest bottlenecks), this padding will be smaller in many cases than the padding internally performed, to achieve optimized sizes.

[georg.viehoever]

Hi Carlos,

first of all thanks for your advice regarding feathering. Your recommendation encouraged me to implement it, and apparently this is a successful procedure for handling the seams between images.

...Georg, I'm 99% sure that the problems with the stars at the boundaries will banish with a 3px zero padding, before the gradient calculation/merging. It will slow down the process a little bit, but is not a great deal. Furthemore, since we are using a solver based on FFT calculations (and I would say that this is one of the narrowest bottlenecks), this padding will be smaller in many cases than the padding internally performed, to achieve optimized sizes.

Regarding the problem with stars at the border to the background of the merged image, I am not sure if I understand you directly. The problem happens there because the gradients of the image regions (some very strong, especially at stars) meet with the gradients of the background (all 0). So the effect we see there is just the kind of balancing between two gradient fields that we would expect. Adding a couple more pixels with gradient 0 would not really help, it would just make the background region larger (that already has gradients of 0). So what exactly do you want to pad there? The Laplace field, the dx,dy gradient field, and with which values?

Thanks for your advice,

Georg

[Carlos Milovic]

I'm saying that you should zero pad in the image domain, and then calculate the gradient. The effect is that you'll have an extra line with non zero coefficients. The next zero values  will guarantee a higher order continuum.
The effect is that strange gradients will be much softer than before.
Mixing mirror padding and zero padding I was able to get rid of almost any boundary artefact in the HDRC implementation. Only those nasty large scale gradients remain (I need to write something more complex to deal with them).
Give it a try.

[Juan Conejero]

Hi Georg,

I have devised a simple algorithm to implement a fast erosion for binary images with a block structuring element:

- Let n be the size of the structuring element in pixels.
- Create 4 duplicates of the original mask. Call them L, R, U, D.
- Translate L by n pixels leftward.
- Translate R by n pixels rightward.
- Translate U by n pixels upward.
- Translate D by n pixels downward.
- Compute the bitwise logical AND of L, R, U and D. This is the result of the erosion operation.

The algorithm is approx. O(5*N) and its complexity is independent on n. It is very easy to implement with PCL classes. Later I'll put an example with the Crop and PixelMath tools that demonstrates the way it works. There are no bottlenecks now

[Juan Conejero]

Here is an example for n=100:

http://pixinsight.org/images/forum/20110729/fast-erosion.jpg

With a little work the algorithm can be implemented optimally by computing the translated pixels and the result in a single loop. This and a fast convolution (via FFT-based convolution or separable convolution, depending on the feathering length) can make up a very fast feathering procedure.

[Juan Conejero]

A note of caution (to prevent misunderstandings): this algorithm is not strictly equivalent to a true erosion with a square structure, but is a sufficiently good approximation for this application (shrinkage of a binary mask that consists approximately of a rotated rectangular structure).

Here is a PixelMath expression that implements the same algorithm in a single operation, without requiring the Crop tool:

Code: [Select]

a = pixel( $T, x()-n, y() );
b = pixel( $T, x()+n, y() );
c = pixel( $T, x(), y()-n );
d = pixel( $T, x(), y()+n );
a & b & c & d
with symbols:

Code: [Select]

n=100, a, b, c, d
If instead of AND we use OR (a | b | c | d) the result is an approximation to a dilation filter, but it has serious problems with corners.

Finally, if we use the XOR operator (a &| b &| c &| d) and n=1, the result is a fast edge detector that works quite well

[georg.viehoever]

Hi Georg,

I have devised a simple algorithm to implement a fast erosion for binary images with a block structuring element:

- Let n be the size of the structuring element in pixels.
- Create 4 duplicates of the original mask. Call them L, R, U, D.
- Translate L by n pixels leftward.
- Translate R by n pixels rightward.
- Translate U by n pixels upward.
- Translate D by n pixels downward.
- Compute the bitwise logical AND of L, R, U and D. This is the result of the erosion operation.

The algorithm is approx. O(5*N) and its complexity is independent on n. It is very easy to implement with PCL classes. Later I'll put an example with the Crop and PixelMath tools that demonstrates the way it works. There are no bottlenecks now

Hmm, I am not so sure. Suppose you have a mask with vertical stripes of width n/2. Then this erosion will do nothing, because the translations will map stripes to stripes. The kind of patterns that can be handled this way is certainly limited. It may work properly if the mask is always a single rotated rectangle - which would probably be sufficient for the typical merge of aligned images. Have to think about it.

Thanks,
Georg

[georg.viehoever]

A note of caution (to prevent misunderstandings): this algorithm is not strictly equivalent to a true erosion with a square structure, but is a sufficiently good approximation for this application (shrinkage of a binary mask that consists approximately of a rotated rectangular structure).
...

Just what  I thought. Thanks!
Georg

[georg.viehoever]

Hi,
just tried separable filters. They are just what I needed. With a FeatherRadius of 20 (meaning a sidelength for erosion and convolution of 41), a convolve of Harry's NGC7000 took 50 seconds. With separable filters, it is just 2.3 seconds 8). I also changed erosion, and it now improved from roughly 45 seconds with the PC Lstandard erosion process to 2.2 seconds using the code below.

Code: [Select]

  /// do erosion via convolution. Works only for binary masks with values of 1.0 and 0.0.
 //static
  void GradientsMergeMosaic::erodeMaskConvolve(weightImageType_t & rMask_p,int32 shrinkCount_p)
  {
    int sideLength=1+2*shrinkCount_p;
    int nElements=sideLength*sideLength;

    SeparableFilter filter( sideLength,1.0/nElements);
    SeparableConvolution convolve(filter);
    convolve >> rMask_p;
    rMask_p.ResetSelections();
    rMask_p.Binarize(0.999);
  }
Thanks for your help!
Georg