
Rafal,
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:
x1,D(x1)
x2,D(x2)
x3,D(X3)
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))/(x3x2)
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.
David