[SOLVED] making sense of sigmoidal-contrast

[NicolasRobidoux]

I'm probably reinventing the wheel, but here is:
Consider the translated sigmoidal function

s(x) = 1/(1+exp(a(d-x))) with a>0 and x in [-inf,inf].

s is an increasing function with range (0,1).
Now, consider the derived function

S(x) = 1/2 + ( s(x)-s(m) ) / ( s(1)-s(0) ) with x in [0,1] and m in (0,1).

S is an increasing function too. (The "1/2 +" was not chosen by accident: see below.)
We would like so choose d so that S(0)=0, S(1)=1 and S(m)=1/2, that is, so that the range of S is [0,1] and the midpoint m is sent to 1/2. It turns out that these three conditions are satisfied if

s(m) = (s(0)+s(1))/2

which is equivalent to

d = (ln D) / a, where
D = ( (A+1)M-2A ) / ( A+1-2M ),
A = exp(a),
M = exp(a*m)

One consequence is that if this is the approach taken, a cannot be too close to zero, since D->-1 as a->0.
P.S. m can't be too close to zero either, since D->-1 as m->0 as well. Since things are symmetric, one expects trouble when m is close to 1, and indeed D->-A as m->1.
P.S. ... and of course, m=1/2 gives d=1/2 irregardless of the value of a, so things check out.
(Checking my algebra/reasoning and figuring what what constraints one must impose to stay away from singularities.)

[NicolasRobidoux]

An alternate formulation:

S(x) = ( s(x) - s(0) ) / ( s(1) - s(0) )

is the only affine function of y=s(x) which satisfies S(0)=0 and S(1)=1.
The fixed point condition S(m)=m is then equivalent to

D = ( A*(M-1)-m*M*(A-1) ) / ( (A-1)*m+1-M )

This, of course, is different: It's sigmoidal-contrast keeping m in place instead of sending it to 1/2.
(Again, I'll need to triple check: Not a good algebra day.)

[NicolasRobidoux]

The top post points out that things need to be clamped when D<epsilon.

[fmw42]

If you do not already have these or seen them, see

http://www.cs.dartmouth.edu/farid/downl ... ls/fip.pdf p44, eq 6.2
http://osdir.com/ml/video.image-magick. ... 00006.html

[NicolasRobidoux]

Thank you Fred.
----
First step: Rewrite as s(x) = 1/(1+C*exp(-a*x)) (or C/(C+exp(-a*x)) or 1/(1+exp(d-a*x)) and see if this cleans things up some and leads to "clean" clamping.
Second step: Always use S(x) = ( s(x) - s(0) ) / ( s(1) - s(0) ) because having a range equal to [0,1] is not negotiable. Having m sent to (exactly) 1/2 is negotiable.

[fmw42]

Thank you Fred.
----
First step: Rewrite as s(x) = 1/(1+C*exp(-a*x)) (or 1/(C+exp(-a*x) or 1/(1+exp(d-a*x)) and see if this cleans things up some and leads to "clean" clamping.
Second step: Always use S(x) = ( s(x) - s(0) ) / ( s(1) - s(0) ) because having a range equal to [0,1] is not negotiable. Having m sent to (exactly) 1/2 is negotiable.

I think you need a constant term in the exponential so that the curve can be shifted left or right as in the current implementation of -sigmoidal-contrast, i.e, exp(-a*x+b). That is, unless you are just interested in a new formula that is symmetric and centered at 0.5,0.5.

Your second step can be written as S(x) = ( s(x) - s(0) ) / ( s(1) from the first paper I reference above and still have it normalize to the range of [0,1]. But please explore your other normalization technique.

[NicolasRobidoux]

@Fred: With your b,

C=exp(b)

in the first rewrite, and

C=exp(-b)

in the second (which had a typo, now corrected).
All these versions are equivalent (in exact arithmetic) unless one clamps something, which is what I'm going to have to do. But yes, C will be tightly related to the midpoint value.
Indeed, what's tricky is allowing a wide range of "midpoints" to be shifted to 1/2 while allowing low and high contrast. And getting this to be accurately reversible.

[fmw42]

@Fred: With your b,

C=exp(b)

in the first rewrite, and

C=exp(-b)

in the second (which had a typo, now corrected).

DOH! I completely missed the C coefficient that you added to your expression. Of course they are equivalent with the right choices of sign. Sorry for my oversight!

[NicolasRobidoux]

It's starting to look like I have to consider two cases: 0<m<=1/2, and 1/2<=m<1, which of course should give the same thing when m=1/2 (the standard case).

[NicolasRobidoux]

Now I see that fudging things by using ( s(0) + s(1) ) / 2 computed from the actual m was a pretty nice and expedient solution.

[NicolasRobidoux]

Maybe what "you" want is that the inflexion point of ( s(x) - s(0) ) / ( s(1) - s(0) ) be at x=m? This is an alternate meaning of "m is the midpoint" as one understands things in the context of S-shaped curves. This may give a well regularized way of getting the gist of what is desired when m is not the usual 1/2. If

s(x) = 1/( 1 + C*exp(-a*x) ),
S(x) = ( s(x) - s(0) ) / ( s(1) - s(0) ),
A = exp(-a),
M = exp(-a*m),

then

C=( sqrt( (M-A)^2 + 4*M*A ) - (M-A) ) / ( 2*M*A )

puts the inflexion point of S exactly at x=m. Note that C is always > 0, even for a=0, and that, of course, one gets the usual scaled sigmoid when m=1/2.
I need to figure out what the limit of S is when a->0^+.
P.S. The limit of S as a->0 from the right is S(x)=x. That is: contrast = 0 gives the identity map S(x) = x, not matter what m is, and the singularity of S is removable.
P.S. The above corresponds to the formulation given in Hany Farid's book Fundamentals of Image Processing linked above by Fred) with beta = ln(C).

[NicolasRobidoux]

Here is the Axiom code (FLOSS computer algebra system with elegant language) that verifies that the above choice of C puts m at an inflection point:

Code: Select all

)cl a
A := exp(-a)
M := exp(-a*m)
C := ( sqrt( (A-M)^2 + 4*A*M ) + (A-M) ) / (2*A*M)
s x == 1/(1+C*exp(-a*x))
S x == ( s x - s 0 ) / ( s 1 - s 0 )
eval( D(S x,x,2), x=m )

The output of the last line is 0.
Here is the code that verifies that, for every value of m in (0,1), S(x)->x as a->0:

Code: Select all

limit(S x,a=0)

The output of this is x.
-----
This verifies things from an algebraic viewpoint, but I'm still not sure this is really the way to go.

[NicolasRobidoux]

You see, the other (current in IM) formulation also gives that there is an inflection point at x=m:

Code: Select all

)cl a
s x == 1/(1+exp(a*(m-x)))
S x == ( s x - s 0 ) / ( s 1 - s 0 )
eval( D(S x,x,2), x=m )

This gives 0.
It also gives S(x)->x in the limit as a->0:

Code: Select all

limit(S x,a=0)

This gives x.
Yet, there is no question that the above formulation is simpler to compute and understand: If you translate and scale a function (which is the relationship between S and s), the inflection point(s) move in predictable ways.

[NicolasRobidoux]

Unless I hear protests, this is what I will do:
Since the current version (the "old" version, with code slightly cleaned up) is totally fine except that the removable singularity that occurs when a->0 is not handled correctly, I'll put a branch that replaces S(x) by its limit when a is very small (MagickEpsilon to the rescue!), this limit being the identity S(x)=x (which makes complete sense as what one should get when the contrast a is set to 0).
I'll also try to make the code a bit more transparent, and will add comments.
I'll remove the "+0.5" I can't figure what for.
I'll make sure the inverse is exactly that.
And I'll remove all traces of the "new" version which I believe was put in in February or so.
Note: The notation roughly matches the source code, which is why it may seem somewhat opaque. a is the contrast in [0,infinity), m is the midpoint (inflection point location) in (0,1), s is the standard sigmoidal, S is the scaled sigmoidal (fixed so that S(0)=0 and S(1)=1).
P.S. For my eventual sigmoidal halo reduction enlargement method (with locally determined min/max) purposes, the S with m=1/2 which has the largest second derivative at 0 and 1 occurs for a approximately equal to 4.79871456103093534, a value which works quite well with text, with orthogonal Lanczos 3 for example. If a local min/max pair is used instead of the limits of the gamut, less contrast will go a long way, so 5ish, which is less than what I generally recommend, will actually work well (hopefully).
P.S. No: contrast 4.8 would appear not to be quite enough for best looking "sigmoidized enlargement", even if replacing global min/max by local ones. It's better than enlarging through linear light. But more contrast is better.

[NicolasRobidoux]

Cristy and Anthony and Fred (and Glenn?):
I just svn checked in a major reorganization of the sigmoidal-contrast code. It basically gives the same values as the "old" formulas.
(I unfortunately forgot to fully use svn --changelist again, so 10 files were changed besides enhance.c. Apologies.)
I only changed it in IM7 to make comparison easier, pending triple check and comments.
Two major changes:
1.The limit as contrast->0 of S(x) (namely the identity map S(x)=x) is hardwired as soon as contrast < 4*MagickEpsilon.
2.There was an OpenMP "parallelized for" in a loop which appeared to me be producing a LUT, not processing an image. I removed it. I know how and where to put it back (not quite as in the old code, because I moved the for loop outside of the branches having to do with values of contrast; I understand that when using OpenMP, it may or may not be preferable to do otherwise). If you tell me you want an OpenMP call when producing a 1D LUT, let me know, and I'll put the "parallelized for" back. It would appear to me to be unnecessary. But I've not timed it. (You can ask me to do benchmarks.)
-----
I'll add something to ChangeLog when I push this into IM6.