I'm generating set of contour curves that I would like to convert to a parametric equation. I have no idea how to do it other than polynomial fit or elliptical equation, neither of which seems unsuitable. I could even use guidance on how to research this problem. When I lookup "squarish circle" on the Interwebs, all I mostly find discussion of rectangles with rounded corners, which is NOT what I'm looking for.
The below is an idealized example of the type of contour: 
I want to parameterize the curve for 1) correcting of quantization error (not shown in above image) and 2) data compression.
In the above example, like most of the of contour sets I'm working with, the inner most contours lines could easly fit to a circle/ellipse, but the same cannot be said of the outer curves, as they become more squarish.
For my application, the contours can be expected to be smooth, well-behaved, convex curves, and will have at least one axis of symmetry. But the data I'm working have finite sampling, producing quantization error in the calculated contours (particularly along angles away from $0, \pm\pi/2, \pi$). I want to correct for this quantization error by fitting the data to an a priori parametric expression that can capture the squarish nature of the curve.
I've already tried simply fitting calculated $r$ vs. $\theta$ to polynomial curves but even high-order (5+ order) polynomial produce very bad results. Equally important, standard polynomial fitting will not take advantage of the fact that, in my case (due to selective non-uniform sampling), $f(\theta)$, $f'(\theta)$ and $f''(\theta)$ can be made to have very little quantization error at $0, \pm\pi/2, \pi$.
Best Answer
I have developed a generalization of the superellipse, called superconics. The general form is given by
$$f(X) = b(1-c^2|X/a|^q)^{1/p}$$
or its canonical form
$$f(x) = (1-c^2|x|^q)^{1/p}$$
Here, $a$ and $b$ scale the $x$ and $y$ axes, resp. and $c^2=\pm1$. $c^2=1$ corresponds to elliptic and parabolic types and $c^2=-1$ corresponds to hyperbolic types. More to the point, $c^2$ can vary smoothly between.
In addition to the above intrinsic equation for the superconics, we can derive the following $parametric$ equation for superconics in the complex plane:
$$z=|\cos^{2/q}(u)|\text{sgn}(\cos(u))+i|(1-c^2\cos(u))^p|\text{sgn}(\sin(u))$$
where $u = [0,2\pi]$ for a closed curve.
The plot and animation below show the range of the superconics for closed curves (you can also do hyperbolic types). In the plot the number above each image are the values of $p$ and $q$, resp. It seems to me that the parameter range you are interested in is $p\le.5$ and $q\ge2$. Probably the easiest thing to do is to set $q=1/p$; this corresponds to the more familiar superellipse.