PixInsight Forum (historical)
Software Development => New Scripts and Modules => Topic started by: georg.viehoever on 2011 July 28 12:10:14
-
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
-
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
-
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.
-
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
-
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
-
is there an implementation of a Box Average Filter in PCL?
Sure:
/*
* 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:
/*
* 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 :) )
-
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.
-
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
-
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.
-
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 :)
-
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.
-
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:
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:
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 :)
-
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
-
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
-
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.
/// 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
-
Try this one:
/*
* Shrink the mask with a fast approximation to erosion of a binary image
* with a square structuring element of size shrinkCount_p.
*/
void GradientsMergeMosaic::erodeMaskFast( weightImageType_t& rMask_p, int shrinkCount_p )
{
int w = rMask_p.Width();
int h = rMask_p.Height();
if ( w <= 2*shrinkCount_p || h <= 2*shrinkCount_p )
{
rMask_p.Zero();
return;
}
weightImageType_t tmp( rMask_p );
for ( int i = 0; i < shrinkCount_p; ++i )
{
weightImageType_t::sample* r0 = rMask_p.ScanLine( i );
for ( int j = 0; j < w; ++j )
*r0++ = weightImageType_t::pixel_traits::MinSampleValue();
weightImageType_t::sample* r1 = rMask_p.ScanLine( h-i-1 );
for ( int j = 0; j < w; ++j )
*r1++ = weightImageType_t::pixel_traits::MinSampleValue();
}
for ( int i = shrinkCount_p, i0 = 0, i1 = 2*shrinkCount_p; i1 < h; ++i, ++i0, ++i1 )
{
const weightImageType_t::sample* m0 = tmp.ScanLine( i0 );
const weightImageType_t::sample* m = tmp.ScanLine( i );
const weightImageType_t::sample* m1 = tmp.ScanLine( i1 );
weightImageType_t::sample* r = rMask_p.ScanLine( i );
for ( int j = 0; j < shrinkCount_p; ++j )
*r++ = weightImageType_t::pixel_traits::MinSampleValue();
for ( int j = shrinkCount_p, j0 = 0, j1 = 2*shrinkCount_p; j1 < w; ++j, ++j0, ++j1, ++m0, ++m1 )
*r++ = (*m0 && *m1 && m[j0] && m[j1]) ?
weightImageType_t::pixel_traits::MaxSampleValue() :
weightImageType_t::pixel_traits::MinSampleValue();
for ( int j = 0; j < shrinkCount_p; ++j )
*r++ = weightImageType_t::pixel_traits::MinSampleValue();
}
}
I haven't compared it, but although your convolution+binarization solution takes advantage of PCL's parallel implementations, I think this routine should be faster even without parallelization. Parallelizing it would be very easy, in fact. Worth trying.
Let me know when you have your code ready for release.
-
Hi Juan,
thanks for the code, will try it later.
Thinking a little bit about the code in http://pixinsight.com/forum/index.php?topic=3289.msg22654#msg22654, I feel that certain morphological transforms probably are separable just as convolutions. I think for a square erosion that is certainly the case. I guess someone more versed in group theory will probably be able to explain that we can see max/min and and/or operations as equivalents to +* in this context, so the use of an outer product as in separable convolutions is possible here as well.
Maybe an opportunity to implement SeparableMorphologicalTransforms in PCL?
I hope to have some code later this weekend.
Georg
-
Hi Georg,
AFAIK there is no correspondence between separability in linear operations, where an outer product can be used to transform a 2D matrix operation into two 1D vector operations, and separability in nonlinear operations such as erosion and dilation (and, by generalization, any rank-order operator).
Right now PCL implements 'brute force' morphological transformations, whose complexity is roughly O((n*N)^2), n being the size of the structuring element and N representing the size of the image. My implementations have the advantage that one can use arbitrary structuring elements; for example, circular and star-like elements, which are very useful to process astronomical images. This confers great flexibility but the problem is that these operations are extremely slow for moderately sized filters.
There exist separable morphological transformations, but AFAIK they are always restricted to square and similar kernels. The most well-known one is the van Herk/Gil-Werman algorithm (Marcel van Herk, A fast algorithm for local minimum and maximum filters on rectangular and octagonal kernels, Pattern Recognition Letters, 13, pp. 517-521, 1992). This algorithm has complexity independent on the size of the structuring element and as the paper title tells it works for rectangular and octogonal SEs. Octogonal is not too far from circular, so it sounds as a quite useful alternative for DS images.
Since your work has shown that we need better/faster morphological transformations, and I need them also for several projects (one of them is a morphological wavelet transform tool that I want to publish this Fall-Winter), I'll try to implement some variants of the van Herk/Gil-Werman algorithm in PCL ASAP.
Looking forward to your code to publish it as a development update.
-
Something that I've seen in your code and I'd like to comment a bit. You use to write code to generate convolution filters such as:
SeparableFilter filter( sideLength,1.0/nElements );
Note that in PCL, proper weighting of filter kernels is automatically carried out by the Convolution, SeparableConvolution and FFTConvolution classes. So the code above is equivalent to (the much easier to understand version):
SeparableFilter filter( sideLength,1.0 );
As long as all kernel elements are equal, a filter constructed this way is equivalent to a box average filter.
-
Older messages on this topic in http://pixinsight.com/forum/index.php?topic=3283.0
Newer messages on this topic in http://pixinsight.com/forum/index.php?topic=3295.0