-fx airy, j0, and j1 functions

[jlettvin]

I failed in my attempt to make changes local to my machine.
I am running ubuntu and libraries (specifically lperl) and other resources are not where they ought to be.
May I ask for the code proposed below to be inserted and made available for ubuntu packages?

Here are the changes needed in "magick/fx.c" function MagickRealType FxEvaluateSubexpression.

/********************************************/
case 'A':
case 'a':
/* abs and acos code */
if (LocaleNCompare(expression,"airy",4) == 0 )
{
alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,beta,
exception);
MagickRealType half = (MagickRealType) 0.5;
MagickRealType zero = (MagickRealType) 0.0;
if( alpha == zero ) return half;
else return((MagickRealType) ( j0((double) alpha)) / (double) alpha);
}
/* asin, alt, and other following code */

/********************************************/

case 'J':
case 'j':
/* j code */
if (LocaleNCompare(expression,"j0",2) == 0 )
{
alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,beta,
exception);
return((MagickRealType) j0((double) alpha));
}
if (LocaleNCompare(expression,"j1",2) == 0 )
{
alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+3,beta,
exception);
return((MagickRealType) j1((double) alpha));
}
break;

[magick]

We'll add support for airy(), j0(), and j1() to the -fx option in the next point release. Thanks.

[jlettvin]

Thank you too.
I may have future suggestions, but those changes cover all my current image needs.
Over and out.

[fmw42]

Don't you need to specify the order of the j0 and j1 functions? Is that the argument to these functions?

As an aside, I believe Anthony has some intention to add sinc, and jinc to the list of -function (due to their use in FFT filters), but I am not sure when he intends to get to those. He has already implemented some of the bessel functions in his -resize filters.

Another aside, IM has B e s s e l O r d e r O n e code already which is j1 order zero.

I wrote a script to implement some of these:

# Bessel function of first kind order 0 from
# formula 9.4.1 and 9.4.3 from Abramowitz and Stegun, p369-370

# Bessel function of first kind order 1 from
# formula 9.4.4 and 9.4.6 from Abramowitz and Stegun, p369-370

Don't know if these will be useful.

[magick]

Type

man j0
man j1

to get a description of these functions. We added sinc() and jinc() to the -fx option.

[fmw42]

thanks

please consider adding these (esp rotated 2D sinc and jinc) to -function at some later time when details of scaling factors can be worked out between Magick, Anthony, Rick Mabry, myself and any others. In my fft (HDRI) filters, the x arguments of jinc (and rotated 2D sinc) needs to be scaled by pi, so some kind of scaling factor(s) would be needed. Also these go positive and negative and so would need to be different for HDRI and non-HDRI as the non_HDRI would need a bias added to shift the negatives to the positive range and scale the max to 1.

[magick]

The best path forward would be to work out these functions with -fx and post them here. We can them code them for the -function option. An alternative would be to wait for Anthony to implement these functions.

[fmw42]

The best path forward would be to work out these functions with -fx and post them here. We can them code them for the -function option. An alternative would be to wait for Anthony to implement these functions.

There is no rush. I suggest we wait for Anthony's input. I can formulate something for discussion and we can all get our input before either you or Anthony create something. Having them in -fx will help in formulating what is needed. So do not hesitate to proceed on that front.

Thanks

Fred

[jlettvin]

My apologies.
My earlier note used the wrong Bessel function, j0, in Airy
where it should have used j1.

Please allow me to document the correct calculation.

http://scienceworld.wolfram.com/physics/AiryDisk.html

http://en.wikipedia.org/wiki/Airy_disk
(See "Mathematical Details about halfway down)

The formula below is the correct mathematical description of
the PSF (Point Spread Function) of a monochromatic point source
on the optic axis of an in focus optical system with a circular aperture.

Airy: I( r ) = I0 * ( j1( r ) * j1( r ) ) / ( r * r )
where
r: ( 2 * pi * a * q ) / ( lambda * R )
a: radius of the aperture
q: radial distance from the optic axis
j1: Bessel function of the first kind of order 1
R: observation distance
lambda: wavelength

Note: The peak value of Airy is 0.25 since
the peak value of ( j1( r ) / r ) is 0.5
and this is correctly handled in the code you released.
It is probably appropriate to leave this intact
rather than conveniently normalizing Airy to 1.0.

I write C++ code to do this and write out .ppm files all the time.
It is actually fairly trivial though tricky code and the results are viewable here:
http://www.aneuran.com/aneuran/wiki
The images on that page illustrate the PSF of the human eye as the pupil changes diameter.

I have attempted to make the fix in your source, recompile, and test.
But the results seem to be incorrect.
In particular, this is the correction I attempted to make in magick/fx.c:

if (LocaleNCompare(expression,"airy",4) == 0)
{
alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,beta, exception);
if (alpha == 0.0)
return(0.25);
double local = j1( (double) alpha ) / alpha;
return( (MagickRealType) ( local * local ) );
}

This is the test command that should produce a monochromatic Airy pattern:
convert WhitePointOnBlackBackground.gif -fx "airy( 0.1 * u )" MonochromaticPointSpreadFunction.gif

To be more consistent with how the human eye PSF converts a TV RGB pixel at one pupil size, use this command:
convert WhitePointOnBlackBackground.gif -channel red -fx "4*airy( 0.56 * u )" -channel green -fx "4*airy( 0.53 * u )" -channel blue -fx "4*airy( 0.43 * u )" HumanPSF.gif

This would account for the longitudinal chromatic aberration of the human optical apparatus.
Human vision also has transverse chromatic aberration where eccentricity from the optic center
displaces the centers of the Airy patterns so they are no longer concentric.
In fact, their relative displacements are in proportion to the difference in lambda
along a radial line from the optic center.
I do not yet seek to address this additional characteristic,
but hope one day to contribute workable fixes to the team.

Eventually, I hope to contribute code that completely emulates
characteristics of human optical transforms when presented with RGB data such as from a TV or computer.
Continuous spectrum characteristics may be beyond what is possible with ImageMagick.

Just to complete the information about this specialized Airy diffraction pattern:
the PSFs of multiple point sources are NOT ADDITIVE!
Airy ought not be used as an image filter.
The correct use of j1 in multiple source systems (pretty much all images) is
Intensity = square( w1 * j1( u1 ) + w2 * j1( u2 ) + ... wN * j1( uN ) );
where w1, w2, ... wN are wave heights for each of the contributing points.
In other words, the wave functions must be added (superposition) before squaring.
I would enjoy contributing to ImageMagick in ways like this
once I understand how to make j1 and Airy work at all.

As I said above, the sample fix I attempted failed to produce the desired results.
The image produced by "convert" does not resemble the images I produce from my own programs
and I know the results from those programs are correct from many years of testing.
Can you advise me on my error?

[jlettvin]

Slight revision on note from last night.

For multiple point sources, the correct formula is:
Intensity = square( w1 * j1( u1 ) / u1 + w2 * j1( u2 ) / u2 + ... + wN * j1( uN ) / uN );

In the formula submitted last night I forgot the divisions.

The morning is wiser than the evening.

[fmw42]

The airy function is [j1(r)/(r)]^2 (which by the way is the square of the jinc function, although the jinc function like the sinc function has been normalized in IM as 2*j1(pi*r)/(pi*r) so that its max at r=0 is 1)

In the code currently we have (IM 6.6.0-10):

if (LocaleNCompare(expression,"airy",4) == 0)
{
alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,beta,
exception);
if (alpha == 0.0)
return(0.5);
return((MagickRealType) (j0((double) alpha)/alpha));
}


Mathematically, this should be changed simply by replacing j0 with j1, squaring and also replacing 0.5 with 0.25

That is change these two lines


return(0.5);
return((MagickRealType) (j0((double) alpha)/alpha));


with these two lines


return(0.25);
return((MagickRealType) (j1((double) alpha)*j1((double) alpha) / (alpha*alpha)));


(hope I have the correct balance on parens)

One could normalize by making it 4*[j1(pi*r)/(pi*r)] (so that it has max value at r=0 of 1) to be consistent with jinc and sinc or leave as it is. I am slightly partial to the normalized form for consistency and maximum dynamic range, but don't really care that much either way. Either way, the user can insert the appropriate normalization himself.

If user jlettvin or any one else, disagrees with my suggestions or analysis, please comment.

[jlettvin]

Thank you for examining my code and pointing out jinc to me.

I tried your code without success:
if (alpha == 0.0)
return(0.25);
return( (MagickRealType)(
(j1((double)alpha)*j1((double)alpha)) /
(alpha*alpha)));

I also tried convert WhitePixelonBlack.gif -fx "jinc(u)*jinc(u)" PSF.gif
also without success even when I use multipliers on u
between 0.001 and 1000.0 by powers of 10.

The images are just simply NOT the diffraction patterns they are supposed to be.
Remember to use my wiki images for reference to see the correct character for these images:
http://www.aneuran.com/aneuran/wiki/

[fmw42]

images created by fx must have a range of 0 to 1 unless using HDRI. I would need to review your scaling values above and try jinc^2. However, note that jinc=2*j1(pi*r)/(pi*r) as defined by IM. The way to test this first is to use a gradient for input so that it ranges from 0 to quantumrange (65535 on Q16). When this goes into fx it will be converted to range 0 to 1. The output should then also be in the range 0 to 1 due to the scaling in the jinc by 2, and also the pi factor will affect the position of the peaks and troughs. IM will then convert it back to quantum range when creating the resulting image. The use Anthony's im_profile script to graph it.

Alternately, feed your fx expression to Anthony's im_graph script and it will avoid the need for the gradient.

I will try to do this, but may not get to it til tomorrow as I have to leave in about an hour.

But Magick has pointed out that the correct expression is coded in the latest Beta or SVN, but not sure which. He sent me the following expression:

#if defined(MAGICKCORE_HAVE_J1)
if (LocaleNCompare(expression,"airy",4) == 0)
{
alpha=FxEvaluateSubexpression(fx_info,channel,x,y,expression+4,beta,
exception);
if (alpha == 0.0)
return(0.25);
gamma=j1((double) alpha )/alpha;
return(gamma*gamma);
}
#endif


which appears to me to be correct. So you can try that by downloading the beta.

But in the mean time, let me see what I can produce. I will get back as soon as I have an example.

[fmw42]

Here are some examples using jinc^2 where jinc(x)=2*j1(pi*x)/(pi*x)

im_graph '-fx jinc(u)*jinc(u)' airy1.gif


The above is equivalent to [2*j1(x)/x]^2 where x ranges from 0 to pi. This looks about right from the graph on http://scienceworld.wolfram.com/physics/AiryDisk.html

im_graph '-fx jinc(2*u)*jinc(2*u)' airy2.gif


im_graph '-fx jinc(4*(u-0.5))*jinc(4*(u-0.5))' airy4s.gif



convert -size 128x128 radial-gradient: -negate -fx "jinc(2*u)*jinc(2*u)" airy2_grad.gif



convert airy2_grad.gif -evaluate log 1000 airy2_grad_log1000.gif


I hope this helps.

Fred

[jlettvin]

Thank you.

Using your command-line, I was able to duplicate the correct image character.
I will consider this issue closed.

Do you consider my pursuit of adding code for emulating human vision
a worthwile addition to ImageMagick?