# How to find a closed curve approximation to contour lines that look like ellipses?

I am finding it very difficult to find resources giving a overview of techniques that can be used to approximate a closed curve. For example, the wikipedia article on curve fitting has further links to linear/polynomial regression techniques that I am familiar with but their section for closed curve fitting is very brief and limited to mentions on fitting circles and ellipses. The references I can find here are also similar brief or specific:

Fitting of Closed Curve in the Polar Coordinate.

Fitting a closed curve on the roots of ${x \choose k}-c$

My specific problem is as follows. I have a function $f(x,y)$ for which I have plotted some contour lines:

Each contour looks like part of a closed ellipse to me. To solve for the contour lines I considered solving the equation, in polar coordinates, $f(r \sin(\theta), r \cos(\theta)) = k$. However, this is not solvable as confirmed by Sage math.

Thus, I am looking for general approaches that will let me find an approximate closed curve to these contour lines. The only thing I can think of doing is to write down a regression equation for $r$ in terms of $\theta$, and then solve the least squares problem: i.e. find coefficients of the equation such that $(f(x, y) - k)^2)$ is minimised.

EDIT: As requested in the comment, $$f(x,y) = e^{-x - y} \sum_{m = 0}^\infty \frac{x^my^m}{(m!)^2}$$ and contour lines are obtained by considering rate parameters $x$ and $y$ such that the probability of a Skellam distributed random variable $Z$ taking value zero is held constant. While solutions specific to this equation are appreciated (such as in one of the previously linked posts where Stirling's approximation is used), please note that more generic methods of fitting closed curves would be more instructive for me in the future.

• What is the function $f(x,y)$ you have plotted? Please update the question with its definition. – unseen_rider Apr 17 '17 at 15:36
• Perhaps useful: Fuzzyfied Lissajous Analysis , but not with incomplete closed contours. And perhaps this: Ellipses made Useful on page 7 : Best Fit Ellipse. – Han de Bruijn Apr 20 '17 at 19:52
• How do you know if your curves are closed or not ? – Yves Daoust Apr 22 '17 at 14:14
• What kind of representation of the approximation curves are you looking for ? A simple way is to approximate $f(x,y)$ by a bivariate polynomial and use its iso-lines. – Yves Daoust Apr 22 '17 at 14:16
• @YvesDaoust I do not know whether they are closed or not. However, they look like something in polar coordinates, or approximately an ellipse, would give a good approximation. – Alex Apr 24 '17 at 2:34

Perhaps they are ellipses with a focus at (0,0).
Then you have $$\sqrt{x^2+y^2}+\sqrt{(x-a)^2+(y-a)^2}=C\\ (x-a)^2+(y-a)^2=(C-\sqrt{x^2+y^2})^2$$ Tidy it up, leave just the square-root on one side, square again.
Note that $2x^2+2y^2=(x+y)^2+(x-y)^2$.
This converts to $$x-y=\pm\sqrt{g(x+y)}$$ for some quadratic function of $(x+y)$.
Then you can plot $(x,y)=((x+y)/2+(x-y)/2,(x+y)/2-(x-y)/2)$ and see if they fit.

The question remains, what is the relation between $a$ and $C$?
Suppose an ellipse cuts the line $x=y$ at $(-b,-b)$ and $(d,d)$ where $d>a$. Then the sum of their distances to the foci $(0,0)$ and $(a,a)$ is $$C = \sqrt{2}b+\sqrt{2}(a+b)=\sqrt{2}d+\sqrt{2}(d-a)$$ If you know $b$ and $d$, you can find $a$ and $C$.

• This isn't correct for considered contour lines. – Yuri Negometyanov Apr 24 '17 at 15:47

The most correct way to plot contour lines is analytic one.
Eg, closed form of the OP function is $$f(x,y) = e^{-(x+y)}\sum_{m=0}^\infty\frac{x^my^m}{(m!)^2} = e^{-(x+y)}J_0(2\sqrt{-xy}) = e^{-(x+y)}I_0(2\sqrt{xy}),$$ where
$J_0$ is Bessel Function of the First Kind,
$I_0(x)$ is Modified Bessel Function of the First Kind.

Contour lines of this function has specific form for k=0.01, k=0.1, k=1, k=10, but they are not closed at all.

The equations of ellipses in this case are $$\frac{(x+y-c_k)^2}{a_k^2}+\frac{(x-y-c_k)^2}{b_k^2} = 1,$$ where the requred parameters should be calculated.

To approximate $f(x,y)$ near the bisector of the first quadrant, it can be used the asymptotic series of Abramowitz & Stegun for $I_0(t)$: $$\sqrt{t}e^{-t}I_0(t) \simeq Q_8(t),$$ $$Q_8(t) = 0.39894 228 + \dfrac{0.01328 592}{t} + \dfrac{0.00225 319}{t^2} - \dfrac{0.00157 565}{t^3} + \dfrac{0.00916 281}{t^4} - \dfrac{0.02057 706}{t^5} + \dfrac{0.02635 537}{t^6} - \dfrac{0.01647 633}{t^7} + \dfrac{0.00392 377}{t^8} +\varepsilon,$$ $$-3.75\le x\le3.75,\quad |\varepsilon|<1.9\cdot10^{-7},$$ then $$f(x,y) \simeq e^{\large -(\sqrt x - \sqrt y)^2}\frac{Q_8(2\sqrt{xy})}{\sqrt[4]{4xy}},$$ $$f(x,x) \simeq \frac{Q_8(2x)}{\sqrt{2x}}.$$ Can be considered the end of a segment of the big semiaxis of the ellipse in the intersection point between contour line and bisector, with the equation $$f(ax_k, ax_k) = k.$$

On the other hand, for the contour lines $$d(x,y) = \frac{dy}{dx} = -\frac{f'_x}{f'_{y}} = -\,\frac{e^{-x-y}\left(-I_1(2\sqrt{xy})+I_0(2\sqrt{xy})\large{\frac{\sqrt y}{\sqrt x}}\right)}{e^{-x-y}\left(-I_1(2\sqrt{xy})+I_0(2\sqrt{xy})\large{\frac{\sqrt x}{\sqrt y}}\right)}$$ $$= -\,\frac{yI_0(2\sqrt{xy})-{\sqrt{xy}}I_1(2\sqrt{xy})}{xI_0(2\sqrt{xy}) -\sqrt{xy}I_1(2\sqrt{xy})},$$ and $$\boxed{c_k = \frac{bx_k+by_k}2,\quad a_k = \sqrt2(ax_k-c_k), \quad b_k = \sqrt2|bx_k-c_k|},$$ where $$f(bx_k, by_k) = 0,\quad d(bx_k, by_k) = 1$$ $($there are two ends of the ellipse small axis $(bx_k,by_k), (by_k,bx_k))$.

It can be used the asymptotic series of Abramowitz & Stegun for $I_1(t)$: $$\frac1t e^{-t}I_0(t) \simeq R_8(t),$$ $$R_8(t) = \dfrac12 + 0.87890 594t^2 + 0.51490 869t^4 - 0.15084 934t^6 + 0.02658 833t^8 - 0.00301 532t^{10} + 0.00032 411t^{12} +\varepsilon,$$ $$-3.75\le x\le3.75,\quad |\varepsilon|<8\cdot10^{-9}.$$

• While this is interesting, your remark about the periodic features doesn't seem useful as I am only interested in the top right quadrant for $(x, y)$. To be specific, I am not interested at all in fitting the curve in the other quadrants. – Alex Apr 24 '17 at 2:37
• @Alex The constraints related to the first quadrant have not been marked. At the same time, they determine the real form of the required area. – Yuri Negometyanov Apr 24 '17 at 4:04
• I am afraid I do not understand your comment. In what sense are the constraints not marked? As the x and y variables are parameters of a Skellam distribution they have to be constrained. – Alex Apr 24 '17 at 4:35
• @Alex ready, waiting for notes... – Yuri Negometyanov Apr 24 '17 at 15:45