Best Fit Sphere of a revolution conic

Question

Good morning,

I'm trying to write a program that finds the best fit sphere R and distance to the apex A for a corneal profile modeled as a conic section for apical rafius r and asphericity Q.

graphic representation of the problem

This paper ( Corneal Elevation Topography: Best Fit Sphere, Elevation Distance, Asphericity, Toricity, and Clinical Implications ) explains how to do it.

Given the following formulas for the conic section where usually $r$ $\in$ [6,8.5] and $Q$ $\in$ ]-1,0]:

$$ C(r, Q) = \frac{ r - \sqrt{ r^{2} - (Q + 1)\rho^{2} } }{Q + 1} $$

and the sphere:

$$ S(R, A) = A + R - \sqrt{ R^{2} - \rho^{2} } $$

The problem can be solved by determining the A and B that minimize the following integral:

$$ e(A, R) = \int_{0}^{1} \left ( \frac{ r - \sqrt{ r^{2} - (Q + 1)\rho^{2} } }{Q + 1} - A - R + \sqrt{ R^{2} - \rho^{2}} \right )^{2} \rho d \rho $$

Now, for a math illiterate like me, minimize A and R with this formula is at best a challenge. What I tried up to now is to calculate the integral by dividing the interval [0,1] in many small fractions, and calculate the formula for each R $\in$ [r, r+3] (steps of 0.01) and A $\in$ [0,0.1] (steps of 0.01).

In code the formula I'm using is:

function(A, R, r, Q, dp) {
      return Math.pow( 
              ((r - Math.sqrt( r*r - (Q + 1) * dp*dp )) / (Q + 1) ) - 
              (R + A - Math.sqrt(R*R - dp*dp)),
            2);
}

Which is calculated from 0 to 1 in steps of 0.001.

The problem is that in my solution A is always 0, and R is smaller than the results cited in the paper.

I know I'm missing all the basic theory here, so that's why I wrote this post.

Any help would be really appreciated.

Thanks

Giulio

Answer 1

$$ e(A, R) = \int_{0}^{1} \left ( \frac{ r - \sqrt{ r^{2} - (Q + 1)\rho^{2} } }{Q + 1} - A - R + \sqrt{ R^{2} - \rho^{2}} \right )^{2} \rho d \rho $$

Performing the integration with a suitable symbolic processor we obtain

$$ \left\{ \begin{array}{rcl} s_1 & = & 6 (Q+1)^2 (Q+2)\\ s_2 & = & 16 (Q+1)^2 \left(r^2-1\right)^{3/2} (aQ(r+1))\\ s_3 & = & 16 \left(r^2-Q-1\right)^{3/2} (aQ(r+1))\\ s_4 & = & 12 (Q+1) \left(a^2 (Q+1)^2+2 a Q (Q+1) r+2 \left(Q^2+Q+1\right) r^2\right)\\ s_5 & = & 2 (Q+2) \left(8 a Q (Q+1) \left(r^2\right)^{3/2}+8 Q^2 \sqrt{r^2} r^3+3 (Q+1) r^4\right)\\ s_6 & = & 6 (Q+1) \sqrt{r^2-1} \sqrt{r^2-Q-1} \left((Q+2) \left(r^2-2\right)+2\right)\\ s_7 & = & 3 Q^2 \sqrt{Q+1} r^4 \log \left(\left(2-\frac{Q+2}{\sqrt{Q+1}}\right) r^2\right)\\ s_8 & = & 3 Q^2 \sqrt{Q+1} r^4 \log \left(2 \sqrt{r^2-1} \sqrt{r^2-Q-1}-\frac{(Q+2) \left(r^2-2\right)+2}{\sqrt{Q+1}}\right) \end{array} \right. $$

and then

$$ e(r,a,Q) = \frac{1}{24 (Q+1)^3}\left(-s_1+s_2-s_3+s_4-s_5+s_6-s_7+s_8\right) $$

I hope this helps.

Try this MATHEMATICA script

CC[r_, q_] := (r + Sqrt[r^2 + q rho^2])/q

SS[r_, a_] := a + r - Sqrt[r^2 - rho^2]

int = Integrate[(CC[r, q + 1] - SS[r, a])^2 rho, rho]

int0 = int /. {rho -> 0}

int1 = int /. {rho -> 1}

int01 = int1 - int0

gr1 = Plot3D[(int01 /. {q -> -0.1}), {a, 0, 0.1}, {r, 6, 8.5}];

gr2 = Plot3D[(int01 /. {q -> -0.35}), {a, 0, 0.1}, {r, 6, 8.5}];

gr3 = Plot3D[(int01 /. {q -> -0.7}), {a, 0, 0.1}, {r, 6, 8.5}];

Show[gr1, gr2, gr3, PlotRange -> All]