Legacy ImageMagick Discussions Archive

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.)

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.)

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

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

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.

NicolasRobidoux wrote: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.

@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.

NicolasRobidoux wrote:@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!

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).

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.

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).

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: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: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.

You see, the other (current in IM) formulation also gives that there is an inflection point at x=m:This gives 0.
It also gives S(x)->x in the limit as 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.

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.

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.