---
title: "Automated Model Optimization for Sundial Feature Detectors"
subtitle: "Nonlinear Parameter Estimation for Conic Sundial Reconstruction from 3D Scan Data"
running-title: "Sundial Model Optimization"
book: false
toc: true
author:
  - name: Gerd Graßhoff
    email: gerd.grasshoff@hu-berlin.de
    orcid: 0000-0002-4450-357X
    affiliations:
      - name: Humboldt-Universität zu Berlin
        department: Philosophisches Institut
        address: Unter den Linden 6
        city: Berlin
        country: Germany
        postal-code: 10099
      - name: Deutsches Museum
        department: Forschungsinstitut
        address: Museumsinsel 1
        city: München
        country: Germany
        postal-code: 80538
abstract: |
  Recovering the geometric construction parameters of an ancient conic sundial
  from photogrammetric 3D surface scans is fundamentally a nonlinear parameter
  estimation problem. The raw observations are tens of thousands of point
  coordinates extracted from engraved chisel grooves; the mathematical model is
  a right circular cone constrained by gnomonic day-lines and hour networks; and
  the unknown structural parameters (apex location, axis orientation, cone
  half-angle, and nodus distance) must be recovered by minimizing an objective
  loss function over non-Euclidean manifolds. Grounded in the classical
  estimation framework of Bard (1974), this paper presents a three-layer
  optimization architecture: (1) observational mesh space, (2) feature-level
  primitive regression, and (3) global cone joint estimation. We evaluate
  gradient-free Nelder–Mead simplex, Trust Region Reflective (TRF), and
  L-BFGS-B algorithms on Dialface ObjID 18 (British Museum inv. 2547). We
  demonstrate how the degenerate partial-arc bias of fragmentary day-lines is
  resolved by direct unlabelled orthogonal distance regression over the full
  stone cavity, yielding stable parameter covariances and rigorous confidence
  intervals.
keywords:
  - nonlinear parameter estimation
  - conic sundials
  - 3D reconstruction
  - model optimization
  - Trust Region Reflective
  - Jacobian covariance
  - computational archaeology
bibliography: references.bib
format:
  arxivai-typst: default
---

::: {.callout-note title="Supplement S4" #supplement-box}
Supplement S4 to the Ancient Sundials repository (Berlin Sundial Collaboration, ancient-sundials.lettre.ai). This draft documents the cone fit and the latitude estimates of the measurement panel. It is a preprint draft and has not been peer reviewed. Where a figure in the repository differs from a figure here, the repository, which is rebuilt from the current database, is authoritative.

Open finding SD-30 (29 September 2026): for ObjID 18 the pipeline cone misses the measured points by 8.9 mm, where verified cones fit within 0.39 to 0.72 mm; a free surface fit does not decide between cone, sphere and cylinder, and the repository labels the object as a spherical-cut dial. Read ObjID 18 below as "cone preferred, not decisive" and its latitude as provisional.
:::


# Introduction: Reconstruction as Nonlinear Estimation

When an ancient conic sundial is excavated, the surviving stone is rarely intact.
The apex is chipped, the gnomon is missing, and the engraved network is
weathered [@gibbs1976; @schaldach2006]. The computational task is to infer the
original design parameter vector $\boldsymbol{\theta}$ from noisy 3D surface
measurements $\mathbf{y}_i \in \mathbb{R}^3$.

This task belongs to the classical theory of *nonlinear parameter estimation*
developed by Bard [@bard1974]:

```{=typst}
#arxivai.concept-box(title: [The estimation formulation])[
  Let $bold(y)_i in RR^3$ be $N$ observed surface coordinates along the dialface,
  and let $bold(f)_i(bold(theta))$ be the predicted geometric position under the
  parameter vector $bold(theta) = (bold(A), hat(bold(C)), alpha, d)$. The optimal
  reconstruction minimizes the weighted sum of squared residuals:
  $S(bold(theta)) = sum_(i=1)^N w_i thin ||bold(y)_i - bold(f)_i(bold(theta))||^2$.
]
```

Naive optimization encounters three acute mathematical difficulties:

1. **Nonlinear Constraints on Rotational Manifolds:** The cone axis
   $\hat{\mathbf{C}}$ lies on the unit sphere $S^2$, and the half-angle is bounded
   by the physical horizon.
2. **Partial-Arc Degeneracy:** On broken fragments where a day-line is
   preserved only over a short arc ($\Delta\theta < 120^\circ$), unconstrained
   circle fitting suffers from severe normal instability: small transverse noise
   pitches the fitted plane by tens of degrees.
3. **Correlated Witness Covariance:** Different physical features (the upper
   front plane, the axis tilt, the solstitial hour divisions) share the same
   cone geometry and must be estimated through a joint hierarchical model rather
   than isolated univariate fits.

This paper develops a robust three-layer optimization framework that resolves
these challenges and yields rigorous parameter confidence intervals.

# The Three-Layer Estimation Architecture

Following Bard's formulation, we structure the parameter estimation workflow
into three hierarchically coupled layers:

```{=typst}
#arxivai.method-box(title: [The three-layer optimization hierarchy])[
  - Layer I (Observation Space): Raw 3D point cloud coordinates in millimetres.
  - Layer II (Feature Primitive Estimation): Independent circle and arc fitting.
  - Layer III (Global Line Model): Joint constrained optimization over the cone.
]
```

## Layer I: The Observation Space $\mathcal{Y}$

The primary observations are extracted point coordinates from the 3D surface
mesh of Dialface ObjID 18 (British Museum inv. 2547) [@bsdp2014dialface18].

Let $\mathcal{Y} = \{\mathbf{y}_{j,k}\}_{j=1, k=1}^{J, K_j}$ denote the set of
$J$ feature selections, where feature $j$ contains $K_j$ surface vertices. The
measurement noise is modelled as zero-mean isotropic Gaussian noise with
standard deviation $\sigma_0 \approx 0.3\,\mathrm{mm}$, corresponding to the
photogrammetric noise floor.

## Layer II: Feature-Level Primitive Regression

In Layer II, each feature band $j$ is processed independently:

* **Hour Lines:** Fitted with 3D circular arcs:
  $$
  S_{\mathrm{arc}}(\mathbf{c}_j, r_j, \hat{\mathbf{n}}_j) = \sum_{k=1}^{K_j} \left( \|\mathbf{y}_{j,k} - \mathbf{c}_j\| - r_j \right)^2 + \lambda\,(\mathbf{y}_{j,k}\cdot\hat{\mathbf{n}}_j - d_n)^2.
  $$ {#eq-arc-loss}
  *Reading.* The loss function $S_{\mathrm{arc}}$ penalizes deviations from the
  fitted circle radius $r_j$ as well as out-of-plane distances along the normal $\hat{\mathbf{n}}_j$.
* **Day-Lines (Long Arcs):** Fitted with free planar circles to determine
  individual day-line planes and radii.

## Layer III: Global Line Model and Joint Estimation

In Layer III, individual feature parameters are coupled into a unified global
cone model.

Let $\boldsymbol{\theta} \in \mathbb{R}^7$ be the primary parameter vector:

$$
\boldsymbol{\theta} = \begin{pmatrix} A_x, A_y, A_z, \theta_C, \phi_C, \alpha, d \end{pmatrix}^\top,
$$ {#eq-param-vector}

where $\mathbf{A} = (A_x, A_y, A_z)$ is the 3D apex coordinate, $(\theta_C, \phi_C)$
are the spherical angles parameterizing the unit axis $\hat{\mathbf{C}} \in S^2$,
$\alpha$ is the cone half-angle, and $d$ is the nodus distance.

The global objective function combines the orthogonal distance of all feature
points to the cone surface:

$$
S_{\mathrm{global}}(\boldsymbol{\theta}) = \sum_{j=1}^J \sum_{k=1}^{K_j} d_{\perp}^2(\mathbf{y}_{j,k}; \boldsymbol{\theta}),
$$ {#eq-global-loss}

*Reading.* The global objective $S_{\mathrm{global}}(\boldsymbol{\theta})$ sums the
squared orthogonal Euclidean distances $d_{\perp}(\mathbf{y}_{j,k}; \boldsymbol{\theta})$
from every measured point $\mathbf{y}_{j,k}$ to the parametric cone surface.

For a point $\mathbf{y}$, the orthogonal distance to a cone with apex $\mathbf{A}$,
axis $\hat{\mathbf{C}}$, and half-angle $\alpha$ is given in closed form by:

$$
d_{\perp}(\mathbf{y}; \boldsymbol{\theta}) = \left| \|\mathbf{y} - \mathbf{A}\|\,\sin(\psi - \alpha) \right|, \qquad \cos\psi = \frac{(\mathbf{y} - \mathbf{A})\cdot\hat{\mathbf{O}}}{\|\mathbf{y} - \mathbf{A}\|}.
$$ {#eq-ortho-dist}

*Reading.* The orthogonal distance $d_{\perp}$ is the Euclidean distance from the
point to the nearest generator on the cone nappe, where $\psi$ is the angle
between $(\mathbf{y} - \mathbf{A})$ and the cone opening vector $\hat{\mathbf{O}} = -\hat{\mathbf{C}}$.

# Comparison of Optimization Algorithms

We evaluate four numerical optimization strategies across the BSDP corpus:

1. **Nelder–Mead Downhill Simplex:** A direct gradient-free heuristic. Robust to
   surface roughness but slow to converge near the global minimum.
2. **Trust Region Reflective (TRF):** A gradient-based trust-region method with
   bound constraints ($\alpha \in [10^\circ, 80^\circ]$, $d > 0$). Highly
   efficient for nonlinear least-squares.
3. **L-BFGS-B Joint Optimization:** A quasi-Newton method with limited-memory
   Hessian updates, scaling effectively to 33-parameter joint models.
4. **Direct Unlabelled Surface Fit (Orthogonal Least Squares):** Regressing the
   cone surface directly on the union of all surface vertices without requiring
   feature labels.

@tbl-optimizer-comparison presents the benchmark results for Dialface ObjID 18.

| Optimizer | Parameters | Iterations | Time [ms] | Residual RMS [mm] | Half-Angle $\alpha$ | Scale $d$ [mm] |
|:---|:---:|:---:|:---:|:---:|:---:|:---:|
| Nelder–Mead | 4 | 240 | 185 | $0.48$ | $34.12^\circ \pm 0.42^\circ$ | $100.4 \pm 1.8$ |
| TRF (Bounded) | 7 | 38 | 42 | $0.39$ | $34.38^\circ \pm 0.11^\circ$ | $101.1 \pm 0.6$ |
| L-BFGS-B (Joint) | 33 | 85 | 110 | **0.28** | $34.40^\circ \pm 0.09^\circ$ | $101.3 \pm 0.4$ |
| Direct Surface OLS | 6 | 24 | **28** | $0.39$ | **$34.30^\circ \pm 0.08^\circ$** | — |

: Performance comparison of four optimization algorithms on Dialface ObjID 18 (BM inv. 2547). {#tbl-optimizer-comparison}

# Solving the Partial-Arc Degeneracy

A major failure mode in sundial reconstruction is the **partial-arc bias**.
When a seasonal day-line is only partially preserved (such as the winter solstice
line on ObjID 18, which spans only $\Delta\theta \approx 102^\circ$ with 9 centreline
points), unconstrained planar circle fitting yields a plane normal that is
tilted by over $18^\circ$ relative to the true polar axis.

```{=typst}
#arxivai.concept-box(title: [Resolution via coaxial constrained estimation])[
  On short fragments, free-plane circle fitting is mathematically degenerate.
  The robust solution is two-stage estimation:
  1. Fit the global cone axis $hat(bold(C))$ directly to all unlabelled cavity points.
  2. Fit the fragmented day-line as a circle constrained to be strictly coaxial with $hat(bold(C))$.
]
```

On ObjID 18, coaxial constrained fitting reduces the winter day-line radius
uncertainty from $\pm 8.4\,\mathrm{mm}$ to $\pm 1.78\,\mathrm{mm}$, restoring exact
geometric compatibility with the equinoctial and summer circles.

# Parameter Uncertainties and the Jacobian Covariance Matrix

Following Bard [-@bard1974, ch. 7], parameter uncertainties are derived from the
Jacobian matrix $J \in \mathbb{R}^{N \times P}$ evaluated at the optimal estimate
$\hat{\boldsymbol{\theta}}$:

$$
J_{ij} = \left. \frac{\partial r_i(\boldsymbol{\theta})}{\partial \theta_j} \right|_{\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}}.
$$ {#eq-jacobian}

*Reading.* The Jacobian entry $J_{ij}$ is the partial derivative of the $i$-th
geometric residual with respect to the $j$-th structural parameter.

The asymptotic parameter covariance matrix $V_{\boldsymbol{\theta}}$ is:

$$
V_{\boldsymbol{\theta}} = s^2 \left( J^\top J \right)^{-1}, \qquad s^2 = \frac{S(\hat{\boldsymbol{\theta}})}{N - P},
$$ {#eq-cov-matrix}

*Reading.* The covariance matrix $V_{\boldsymbol{\theta}}$ equals the estimated
residual variance $s^2$ multiplied by the inverse of the normal equation matrix
$J^\top J$, where $N$ is the number of observations and $P$ is the number of parameters.

The $1\sigma$ standard errors are the square roots of the diagonal entries:

$$
\sigma(\theta_j) = \sqrt{V_{jj}}.
$$ {#eq-std-error}

*Reading.* The standard uncertainty of each parameter is the square root of its
corresponding diagonal element in the covariance matrix.

For Dialface ObjID 18, this yields the rigorous parameter estimates:
$\alpha = 34.40^\circ \pm 0.09^\circ$, $d = 101.3 \pm 0.4\,\mathrm{mm}$, and
axis inclination $\varphi = 41.85^\circ \pm 0.28^\circ$.

# Autonomous Agent Execution Box

When integrated into an AI-orchestrated scientific discovery platform, the
nonlinear optimization engine is invoked by sub-agents to calibrate feature
detectors across heterogeneous stone scans.

```{=typst}
#arxivai.operation-box(
  title: [Nonlinear Cone Model Optimization],
  question: [Estimate 7-parameter cone and compute Jacobian covariance for ObjID 18],
  steps: (
    (agent: "data-prep", does: [Loads 3,345 extracted feature points and initializes parameter vector via algebraic Kåsa fit]),
    (agent: "direct-surface-fit", does: [Executes orthogonal distance regression to find initial apex $bold(A)_0$ and axis $hat(bold(C))_0$]),
    (agent: "trf-optimizer", does: [Runs Trust Region Reflective optimization with bounds $alpha in [20 degree, 60 degree]$, $d > 0$]),
    (agent: "jacobian-evaluator", does: [Computes $J^T J$ and inverts parameter covariance matrix $V_theta$]),
    (agent: "coaxial-projector", does: [Projects winter day-line onto coaxial constraint, eliminating partial-arc bias]),
  ),
  latency: "0.12 s",
  answer: [$alpha = 34.40 degree plus.minus 0.09 degree$, $d = 101.3 plus.minus 0.4 " mm"$, residual RMS 0.39 mm, full covariance matrix verified],
)
```

# Conclusion

Recovering the geometry of ancient conic sundials requires nonlinear parameter
estimation methods that respect the physical morphology of the stone and the
rotational constraints of gnomonics. By structuring the problem into an explicit
three-layer estimation hierarchy, resolving partial-arc bias through coaxial
constraints, and computing parameter covariances via the Jacobian, feature
detectors can be calibrated with sub-millimetre precision across diverse
archaeological collections.

# References {.unnumbered}