possible bug in new fx function jinc IM 6.6.0-4

[fmw42]

IM 6.6.0-4 Q16 Mac OSX Tiger

I am testing the new sinc and jinc fx functions.

The sinc seems to be working correctly (mathematically correct). It should have a central peak of amplitude 1 and a first negative peak of amplitude -0.21 (so that it ranges from -0.21 to 1)

convert -size 1x256 gradient: -rotate 90 -fx "sinc(u)" miff:- |\
im_profile - tmp_sinc_profile.gif


convert -size 1x256 gradient: -rotate 90 -fx "sinc(10*u)" miff:- |\
im_profile - tmp_sinc10_profile.gif


#scaled to 0 to 1 so no clipping
convert -size 1x256 gradient: -rotate 90 -fx "(sinc(10*u)+0.21)/1.21" miff:- |\
im_profile - tmp_sinc10a_profile.gif



The jinc does not come out mathematically according to one of two definitions that I have found.

The jinc function ranges from -0.06613 to 0.5. It can be computed as a series approximation to the Bessel function. The formula is obtained from The Handbook Of Mathematical Functions by Abramowitz and Stegun, formula 9.4.1 and 9.4.3, p369-370. I have scripted this successfully in my FFT work. See http://www.fmwconcepts.com/imagemagick/ ... properties

The other definition is basically to scale that by two so that it ranges from -0.1323 to 1 and can be found at http://books.google.com/books?id=d8FMlH ... on&f=false Scroll to the next page for a table of the peak amplitudes.

Testing the implemented jinc, I find that the main peak is at about 80% and so does not fit either of these two models. Perhaps you are using a different model. If so, could you provide me with the main peak amplitude and the first negative peak amplitude (or reference), so that I can do appropriate scaling.
Or (better) perhaps you might reformulate/normalize it so that the main peak is of amplitude of 1 so that it kind of parallels the sinc function in that regard.


convert -size 1x256 gradient: -rotate 90 -fx "jinc(u)" miff:- |\
im_profile - tmp_jinc_profile.gif


convert -size 1x256 gradient: -rotate 90 -fx "jinc(5*u)" miff:- |\
im_profile - tmp_jinc5_profile.gif

[magick]

We added a patch to fix this problem to ImageMagick 6.6.0-5 Beta available by sometime tomorrow. Thanks.

[fmw42]

We added a patch to fix this problem to ImageMagick 6.6.0-5 Beta available by sometime tomorrow. Thanks.

This is still not working and is perhaps worse as it has a spike in it near the origin (left side) and more importantly is not normalized/scaled (to 1 at the origin)

convert -size 1x256 gradient: -rotate 90 -fx "jinc(u)" miff:- |\
im_profile - tmp_jinc_profile2.gif

[magick]

We modified the jinc formula ti match your discussion: J1(πr)/(πr). In C code it looks like this:

• if (LocaleNCompare(expression,"jinc",4) == 0)
  {
  alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,beta,
  exception);
  if (alpha == 0.0)
  return((MagickRealType) (MagickPI/4.0));
  gamma=(MagickRealType) (j1((double) (MagickPI*alpha))/
  (MagickPI*alpha));
  return(gamma);
  }

Can you spot our error?

[fmw42]

We modified the jinc formula ti match your discussion: J1(πr)/(πr). In C code it looks like this:

• if (LocaleNCompare(expression,"jinc",4) == 0)
  {
  alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,beta,
  exception);
  if (alpha == 0.0)
  return((MagickRealType) (MagickPI/4.0));
  gamma=(MagickRealType) (j1((double) (MagickPI*alpha))/
  (MagickPI*alpha));
  return(gamma);
  }

Can you spot our error?

It looks like the function you have now is essentially correct according to the Handbook Of Mathematical Functions by Abramowitz and Stegun for j1 converted to jinc, such that it has a max value of 0.5, but you have put the value at jinc(x) for x=0 to pi/4 rather than 0.5. That is where you get a peak value of pi/4=0.7853975 (which is the approx 0.8 that I observered).


You can leave it at this if you want, by simply correcting the x=0 (alpha=0?) to 0.5 and then it should range from a max value of 0.5 at x=0 to a value of -0.06613 at the first lobe.

Alternately, I would suggest that you double your result, so that it scales the max value to 1 for x near 0 and then set the x=0 value to 1. Namely jinc=2*j1(pi*r)/(pi*/r)

That way it should then range from a max of 1 down to a min of -0.1323 at the first negative lobe.

I can then test that scaling either way.

You might want to consult with Anthony and Rick Mabry to see which form would be preferable. But it really does not matter as one can scale it appropriately if one knows which format is used and the min and max values.

My second suggestion of scaling by 2 was simply so that it resembles the result more closely of the sinc which has a max of 1 at x=0 for sinc(x). That is both sinc and jinc are normalized to 1 at the origin.

Let me know when changed and I will test again. But as it is late for you now, it can certainly wait until tomorrow or whenever as it is not critical for me right now. Eventually, I will make use of it in my scripts.

Fred

P.S.

In my FFT use, I needed the pi factors in both numerator and denominator jinc=j1(pi*r)/(pi*r) AND also in the sinc=sin(pi*x)/(pi*x), so I am wondering if you scaled the sinc by pi or not. I guess it does not matter as sinc(pi*x) would take care of it if you had sinc=sin(x)/(x).

Likewise, if you had jinc=j1(r)/(r), I could use jinc(pi*r) to get j1(pi*r)/(pi*r).

It would be nice if both the sinc and jinc had the same arguments of using pi or not. Because that affects the horizontal scale or compression of the results. That is where the zeros occur and where the lobes have min and max values along the horizontal axis.

But you might want to check with Rick (and Anthony) and see if they have any comments or suggestions regarding these issues of using pi or not and of the factor of 2 for the jinc.

I guess my point is to have both the sinc and jinc be scaled to have max value of 1 at the origin (alpha=0) and use the same horizontal scaling factors whether that is pi or 1 times the horizontal (argument) coordinate. See horizontal normalized sinc=sin(pi*x)/(pi*x) vs unnormalize sinc=sin(x)/(x) at http://en.wikipedia.org/wiki/Sinc_function Also see horizontal unnormalized jinc=j1(x)/(x) at http://mathworld.wolfram.com/JincFunction.html. Normalized would be either jinc=j1(pi*x)/(pi*x) or 2*j1(pi*x)/(pi*x) if you also want the amplitude to be normalized to 1.

But I would like to know what you did or will do for the sinc as well, if you don't mind, so that I can scale it appropriately as needed for my FFT use.

P.S. 2

In looking at the code, it appears that you use sinc=sin(pi*x)/(pi*x). So that is fine with me. Thus using jinc=j1(pi*r)/(pi*r) or normalizing it to jinc=2*j1(pi*r)/(pi*r) would be consistent.

[magick]

In looking at the code, it appears that you use sinc=sin(pi*x)/(pi*x). So that is fine with me. Thus using jinc=j1(pi*r)/(pi*r) or normalizing it to jinc=2*j1(pi*r)/(pi*r) would be consistent.

We choose 2*j1(pi*r)/(pi*r). Look for the patch in ImageMagick 6.6.0-5 by end of day. Thanks,

[fmw42]

In looking at the code, it appears that you use sinc=sin(pi*x)/(pi*x). So that is fine with me. Thus using jinc=j1(pi*r)/(pi*r) or normalizing it to jinc=2*j1(pi*r)/(pi*r) would be consistent.

We choose 2*j1(pi*r)/(pi*r). Look for the patch in ImageMagick 6.6.0-5 by end of day. Thanks,

That is good with me. I will test once I see the next beta. Thanks very much.

[fmw42]

Jinc is still not quite right (if I have not tested prematurely). Small change needed. From the code:

if (LocaleNCompare(expression,"jinc",4) == 0)
{
alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,beta,
exception);
gamma=(MagickRealType) (2.0*j1((double) (MagickPI*alpha))/
(MagickPI*alpha));
return(gamma);
}

If I read this correctly, you removed the test for alpha=0. I think you need to put back that test, such that at alpha=0, jinc=1. Sorry if my explanation above was not clear on this point.

convert -size 1x256 gradient: -rotate 90 -fx "jinc(u)" miff:- |\
> im_profile - tmp_jinc_profile3.gif



This diagram goes from x near 0 with value of 1 to x=0 with value probably of 0. It should be smooth all the way to x=0 where the value should be 1. No spike downward or upward near the x origin.

[magick]

Done. Look for the update to ImageMagick 6.6.0-5 Beta within an hour.

[fmw42]

Done. Look for the update to ImageMagick 6.6.0-5 Beta within an hour.

It is working well now.

convert -size 1x256 gradient: -rotate 90 -fx "jinc(u)" miff:- |\
im_profile - tmp_jinc_profile4.gif


Scaled Half Jinc (so range 0 to 1 and no clipping):
convert -size 1x256 gradient: -rotate 90 -fx "(jinc(10*u)+0.1323)/1.1323" miff:- |\
im_profile - tmp_jinc10a_profile4.gif


Scaled Full Jinc (so range 0 to 1 and no clipping):
convert -size 1x256 gradient: -rotate 90 -fx "(jinc(20*(u-0.5))+0.1323)/1.1323" miff:- |\
im_profile - tmp_jinc10b_profile4.gif


Thanks.

Fred

[fmw42]

I have edited the fx documentation to show these two formulae

sinc(x)=sin(pi*x)/(pi*x)

jinc(x)=2*j1(pi*x)/(pi*x)

[magick]

We just released ImageMagick 6.6.0-5. It should be available in just a few minutes.