I'm trying to find a way to numerically extract or estimate the axis of rotation for an arbitrarily shaped helix, which may be curved and have varying pitch and radius. (See here for an example.)
For instance, let's say I have a line (described by a position vector $\textbf{x}$) and an ortho-normal basis at every point (the line tangent $\textbf{t}$, normal $\textbf{n}$, and binormal $\textbf{b}$). If we parameterize the line by $s$, a helix $\textbf{h}$ can be generated using the following equation:
$$ \mathbf{h}(s) = \mathbf{x}(s) + R(s)[\cos(2\pi k s)\mathbf{n}(s) + \sin(2\pi k s)\mathbf{b}(s)] \tag{1} $$
where $R(s)$ can be considered the radius of curvature of the helix at $s$ projected onto the plane spanned by $\mathbf{n}$ and $\mathbf{b}$, and $k$ is a frequency related to the pitch.
Then is it possible, given only $\mathbf{h}(s)$, to get a "good" estimate of $\mathbf{x}(s)$ back?
Initial approach
My initial approach is to numerically calculate the helix's own tangent, normal, and binormal vectors using the following equations (and a finite difference approximation):
$$ \mathbf{t_{h}} = \frac{\mathbf{h}'}{\|{\mathbf{h}'}\|} \quad,\quad \mathbf{n_{h}} = \frac{\mathbf{t_{h}}'}{\|\mathbf{t_{h}}'\|} \quad,\quad \mathbf{b_{h}} = \mathbf{t_{h}}\times\mathbf{n_{h}} $$
$\|\mathbf{t_{h}}'\|$ is also the curvature $\kappa$, which is the inverse of the helix's radius of curvature $R_{h}$. Together, I can roughly estimate the helical axis as:
$$ \mathbf{x_{est}(s)} = \mathbf{h}(s) + R_{h}(s)\mathbf{n_h}(s) \tag{2} $$
To test $(2)$, I generated a helix using Equation $(1)$ and the following definitions: $$ \mathbf{x}(s) = (0, s^2,s)$$ $$ R(s) = 1.0 + 0.5 \sin(2s)$$ $$ k =1 $$
The results are provided by this link, where the purple line is $\mathbf{x_{est}}$ and the white line is $\mathbf{x}$ (the helix is colored by curvature). I don't consider this a "good" approximation of $\mathbf{x}$ and I'd like to do better.
Also, Equation $(2)$ will fail completely if the pitch of the helix is large relative to the radius. I believe this is because the term $R[\cos(\ldots)\mathbf{n}+\sin(\ldots)\mathbf{b}]$ in equation $(1)$ is essentially the projection of $R_{h} \mathbf{n_h}$ onto the plane $\mathbf{t}$, so as the pitch increases, $\mathbf{n}_h \cdot \mathbf{n}$ tends to $0$ and $R_{h}$ more strongly overestimates $R$.
So I'm looking for an equation, algorithm, or technique to reduce the error between $\mathbf{x_{est}}(s)$ and $\mathbf{x}(s)$. I think the answer may be in estimating $R(s)$ and $\mathbf{n}(s)$ better, but it could be something completely different. Any suggestions are welcome.