proper scaling of the Jinc filter for EWA use

Discuss digital image processing techniques and algorithms. We encourage its application to ImageMagick but you can discuss any software solutions here.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

@imaggie: If the zeros of the denominator are also zeros of the numerator (with a "power" which is no larger in the sense of Taylor expansions) you can fix things so that you get what you "want" by getting under the hood. If this "matching zeros" condition does not hold, your computation as a whole is unstable, which means that unless your input image satisfies some pretty strong conditions a small change in the input may lead to a large change in the output.

There is a quick and dirty way however: when you detect that the denominator is zero, just move the input by "epsilon" (1.e-14, for example; the safe way is to multiply the input by 1+epsilon P.S. unless the input is zero, in which case you replace by epsilon). This is the quick and dirty way of taking limits, which is what I believe you want.

Do I make sense? (#1: This is brutally terse. #2: I don't really know what you're doing.)

P.S. For maximum accuracy, you ideally should also multiply the numerator's input by 1+epsilon when you do that with the denominator's input, although my guess is that no-one will see the difference if you don't.
Last edited by NicolasRobidoux on 2012-01-07T14:27:56-07:00, edited 4 times in total.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

@imaggie: I've not had time to update the ones that are built into resize.c, but I have state-of-the-art ways of approximating extremely accurately fixed combinations of complex functions (trig, Bessel, etc) with polynomials or rational functions which have few terms. The approach is based on a hot rodded version of the Boost C++ Minimax program. If your ratios are fixed (in terms of which functions show up) and indeed you only evaluate where the ratio is defined in the limit (that is, the 0/0 value must converge if evaluated as a limit), my novel approach works.

If you are curious, see chapters 20, 21 and 22 of the masters thesis of my student Chantal Racette which is found here: http://web.cs.laurentian.ca/nrobidoux/m ... is2011.pdf.
What's in her thesis is already outdated (even though she defended December 2), but it gives the idea.

Note: This kind of tweak you do LAST: It's a speed of computation thing, only worth it if you don't use LUTs. If you use LUTs, the "division by zero" problems I alluded to in my previous post can be fixed by replacing 0s by epsilon in the LUT.
Last edited by NicolasRobidoux on 2012-01-07T14:27:45-07:00, edited 1 time in total.
imaggie
Posts: 88
Joined: 2011-12-19T04:15:36-07:00
Authentication code: 8675308

Re: proper scaling of the Jinc filter for EWA use

Post by imaggie »

thanks I'll take a look.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

@imaggie: There is another way of avoiding divisions by zero. Not super accurate if most of your sampling locations are at zeros, but will keep things safe if there are more than one but few compared to the total number of locations you evaluate the function at.

Modify the resize.c code adding a branch at the end of the functions you use, checking if the square of the result is less than a small epsilon^2, and replacing the computed value by the small epsilon if it is.

This is way less difficult than it looks.

Of course, these changes can't be reverted back to the main release of ImageMagick, but it allows you to explore various possibilities without worrying about divisions by zero.

Warning: If a large fraction of the locations sampled happen to be at zeros of the denominator, you'll need to carefully modify the numerators as well.
Last edited by NicolasRobidoux on 2012-01-07T14:27:33-07:00, edited 1 time in total.
User avatar
fmw42
Posts: 25562
Joined: 2007-07-02T17:14:51-07:00
Authentication code: 1152
Location: Sunnyvale, California, USA

Re: proper scaling of the Jinc filter for EWA use

Post by fmw42 »

Nicolas,

I believe he is trying to deconvolve a blurred image (using FFT techniques), whose blurring function in the frequency domain is a sinc (motion blur) and a jinc (lens defocus) in the presence of image noise. All the methods of avoiding divide by zero you mentioned are standard techniques called pseudo-inverse filter. There is also a least square solution called Wiener filtering. If there were no noise, one could get a very good if not perfect reconstruction simply by dividing by the appropriate sinc or jinc, since in the frequency domain, the blurring is simply a multiplication of the sinc or jinc with the transform of the image. This corresponds to convolution with a straight line or circular rect function in the spatial domain. However in the presence of noise, the regularization process of these techniques produces ringing around high contrast areas (which you don't get if there is no noise in the image). By increasing the epsilon, you can get a smoother result, but get more ringing. If you decrease the epsilon, you get a sharper reconstruction, but have amplified the noise.

He is trying to use sinc or jinc windowing applied to the sinc or jinc reconstruction filter to reduce the ringing.

Hope that helps clarify the background for all this.

Fred
imaggie
Posts: 88
Joined: 2011-12-19T04:15:36-07:00
Authentication code: 8675308

Re: proper scaling of the Jinc filter for EWA use

Post by imaggie »

Are the internal jinc zero-crossing available anywhere ?

Fred sent me a link with the first 5 but I'd like at least 6 or even 7.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

@imaggie: There is an app which computes any number of them at higher precision (you need to divide by pi): http://cose.math.bas.bg/webMathematica/ ... lZeros.jsp

The ImageMagick internal values are found here:http://trac.imagemagick.org/browser/Ima ... e/resize.c. Find "jinc_zeros[16]".
User avatar
fmw42
Posts: 25562
Joined: 2007-07-02T17:14:51-07:00
Authentication code: 1152
Location: Sunnyvale, California, USA

Re: proper scaling of the Jinc filter for EWA use

Post by fmw42 »

NicolasRobidoux wrote:@imaggie: There is an app which computes any number of them at higher precision (you need to divide by pi): http://cose.math.bas.bg/webMathematica/ ... lZeros.jsp

The ImageMagick internal values are found here:http://trac.imagemagick.org/browser/Ima ... e/resize.c. Find "jinc_zeros[16]".
@imaggie

I sent you all the above and the 16 numbers from resize.c already in past email.

Fred
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

fmw42 wrote:...
Hope that helps clarify the background for all this.

Fred
It does. Thank you (again).
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

I stated at some point that I'd take the time to compute the optimal vertical/horizontal line preserving blur for Lanczos2Sharp with the new approach. Here it is: .9580278036312191.
This is the Axiom code that computes it:

Code: Select all

)cl a

digits 100

R1 := 3.831705970207512315614435886308160766564545274287801928762298989918839309519011470214112874757423127

r1 := R1 / %pi

R2 := 7.0155866698156187535370499814765247432763115029113138960553778269854960155020186630727149301794664578066071840984767118820251703769720326548673375642560309232023970190386543841800591275430452193373925

r2 := R2 / %pi

jinc x == besselJ(1,%pi*x)/x

wind x == besselJ(1,x*(R1/r2))/x

l x == if (x<r2) then ( wind(x) * jinc(x) ) else 0.

f r == [ r, sqrt(1.+1)*r, 2*r, sqrt(1.+4)*r, sqrt(4.+4)*r, 3*r, sqrt(9.+1.)*r, sqrt(9.+4)*r, 4*r, sqrt(16.+1)*r, sqrt(9.+9)*r, sqrt(16.+4)*r, 5*r ]

z := l 1.e-128

totheright r == _
( _
  l((f(r)).1) _
+ l((f(r)).3) _
+ l((f(r)).6) _
+ l((f(r)).9) _
+ l((f(r)).13) _
+ 2 * _
  ( _
    l((f(r)).2) _
  + l((f(r)).5) _
  + l((f(r)).11) _
  + 2 * _
    (   l((f(r)).4) _
      + l((f(r)).7) _
      + l((f(r)).8) _
      + l((f(r)).10) _
      + l((f(r)).12) ) ) ) / ( z + _
4 * _
( l((f(r)).1) _
+ l((f(r)).2) _
+ l((f(r)).3) _
+ l((f(r)).5) _
+ l((f(r)).6) _
+ l((f(r)).9) _
+ l((f(r)).11) _
+ l((f(r)).13) _
+ 2 * _
( l((f(r)).4) _
+ l((f(r)).7) _
+ l((f(r)).8) _
+ l((f(r)).10) _
+ l((f(r)).12) ) ) )

-- To find the one with sum closest to 0 (this is a fully converged result):

[ totheright(1.04381104202789666520062836+.000000000000000000000000001*i) for i in 0..10 ]

-- The optimal blur
-- (to minimize the worst change in a vertical or horizontal line under no-op):

1/1.04381104202789666520062836
With 4 lobes, it's .9870395083298263, obtained with similar code.
Last edited by NicolasRobidoux on 2012-04-11T06:42:29-07:00, edited 1 time in total.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

@imaggie: RE: motion blur and EWA: K.Mueller, T.Möller, J.E.Swan II, R.Crawfis, N.Shareef, R. Yagel, "Splatting Errors and Antialiasing", IEEE Transactions on Visualization and Computer Graphics, ITVCG 4(2): 178-191, June 1998.
http://www.cs.sfu.ca/~torsten/Publicati ... iasing.pdf
This is about ADDING motion blur. But maybe you'll find something useful in there.

P.S. Interesting: it quotes an article that apparently optimized a piecewise cubic for the accurate reconstruction of vector fields: http://www.computer.org/portal/web/csdl ... /38.310726. I'll have to get myself a copy to find out if it relates at all to the Robidoux and RobidouxSharp Keys splines.
henrywho
Posts: 188
Joined: 2011-08-17T06:46:40-07:00
Authentication code: 8675308

Re: proper scaling of the Jinc filter for EWA use

Post by henrywho »

Does it mean......

JincJinc2: -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.9580278036312191 -define filter:lobes=2 -distort resize ....
JincJinc3: -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.9891028367558475 -define filter:lobes=3 -distort resize ....
JincJinc4: -define filter:filter=Jinc -define filter:window=Jinc -define filter:blur=0.9870395083298263 -define filter:lobes=4 -distort resize ....
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

@henrywho: Yes. Here are (hopefully correct) "short" versions of what you call JincJinc2, 3 and 4:

Code: Select all

convert INPUT.IMG \
-filter Lanczos2 \
-define filter:blur=.9580278036312191 \
-distort Resize WIDTHxHEIGHT OUTPUT.IMG
convert INPUT.IMG \
-filter Lanczos \
-define filter:blur=.9891028367558475 \
-distort Resize WIDTHxHEIGHT OUTPUT.IMG
convert INPUT.IMG \
-filter Lanczos \
-define filter:lobes=4 \
-define filter:blur=.9870395083298263 \
-distort Resize WIDTHxHEIGHT OUTPUT.IMG
One way to describe them is "radial" Jinc-windowed Jinc filtering with built-in "awareness" of the fact that they are used on a grid.
Last edited by NicolasRobidoux on 2012-04-19T10:00:43-07:00, edited 1 time in total.
NicolasRobidoux
Posts: 1944
Joined: 2010-08-28T11:16:00-07:00
Authentication code: 8675308
Location: Montreal, Canada

Re: proper scaling of the Jinc filter for EWA use

Post by NicolasRobidoux »

Suggestion: If you want more sharpness (at the expense of more aliasing/jaggies/deep halos), use the RobidouxSharp Keys cubic:

Code: Select all

convert INPUT.IMG \
-filter Cubic -define filter:C=.3689927438004929 \
-distort Resize WIDTHxHEIGHT OUTPUT.IMG
As far as I can figure, it's the Catmull-Rom of the EWA world.
User avatar
anthony
Posts: 8883
Joined: 2004-05-31T19:27:03-07:00
Authentication code: 8675308
Location: Brisbane, Australia

Re: proper scaling of the Jinc filter for EWA use

Post by anthony »

NicolasRobidoux wrote:Suggestion: If you want more sharpness (at the expense of more aliasing/jaggies/deep halos), use the RobidouxSharp Keys cubic:

Code: Select all

convert INPUT.IMG \
-filter Cubic -define filter:C=.3689927438004929 \
-distort Resize WIDTHxHEIGHT OUTPUT.IMG
As far as I can figure, it's the Catmull-Rom of the EWA world.
But does this make a good general purpose default for 'distorts'!
EG: magnification, magnification, and especially no-op and near no-op distortions.

Thing in terms of a animation that is zooming though the 'no-op' distortion. You do not want a sudden discontinuity of effects on the resulting images. But you what the image to be as close to the original at the no-op point.

As previously seen at the beginning of this discussion (first page) I am not good at critically reviewing results of a distortion filter. One day I can think one filter is better, the next another filter is better. It depend too much on the image and the situation.

But the no-op case is very important as many users (like me) are not dealing with real-world images (photos) but more commonly with cartoon like images (computer graphics, text, line drawing), so a good no-op default is vital.
Anthony Thyssen -- Webmaster for ImageMagick Example Pages
https://imagemagick.org/Usage/
Post Reply