I am the fool who started the first of the APUG threads you listed, in which I suggested a simple mathematical function for fitting densitometric data. Suffice it to say that my methods were not instantly adopted by thousands of photographers, leading to the fame and fortune that I had imagined for myself. But, I have used that function to calculate the Kodak contrast index, and I think that I can suggest a general way to formulate the problem.

The key thing that is needed is a relationship between the exposure and density, which is presumably based on your measured data. I used the function I described, with fit values for the 3 parameters that define that function. Having an analytical function makes things a bit easier, but one could also use a table of closely spaced points.

Let's call the exposure values "x" and the density at a given value of x "D(x)". To use the Kodak definition, the exposure values should be given in log base 2 units (i.e. doublings or "stops"). The CI definition specifies 3 points that lie on a common line:

1. A point on the baseline, where density - (base + fog) = 0.
2. A point to the right of point 1 that lies on the curve and is 0.2 units away from point 1.
3. A point to the right of point 2 that lies on the curve and is 2 units away from point 2.

Call these points:

From Pythagoras, the relationships among the points can be written as:

(x2 - x1)^2 + (D(x2)-D(x1))^2 = 0.2^2
(x3 - x2)^2 + (D(x3)-D(x2))^2 = 2.0^2
(x3 - x1)^2 + (D(x3)-D(x1))^2 = 2.2^2

This gives three equations in three unknowns. Once the values for the three x values are known, then the contrast index can be calculated from the slope of the line, using any pair of points:

CI = (D(x3)-D(x2))/(x3-x2)

Solving the equations has to be done numerically, i.e. iteratively, at least for my analytical function or a table of data points. This requires a starting guess for the three x values, but that isn't too difficult.

I have done this in Python, using the optimize.fsolve function in the scipy package, but I'm sure that similar functionality can found in R.

I hope that this is a bit of help. I'm happy to send you the Python code if that would be useful.