User Guide
==========

.. admonition:: Run this guide
   :class: tip

   .. grid:: 2
      :gutter: 2
      :margin: 2 0 0 0

      .. grid-item::
         .. button-link:: https://colab.research.google.com/github/7jameslondon/MagTrack/blob/master/examples/examples.ipynb
            :color: primary
            :shadow:
            :expand:

            Open in Google Colab

      .. grid-item::
         .. button-link:: https://github.com/7jameslondon/MagTrack/blob/master/examples/examples.ipynb
            :color: secondary
            :shadow:
            :expand:

            Or download and run

Import MagTrack
----------------------

To start, import magtrack. Here we will also import NumPy and a simulation that can generate bead images.

.. code-block:: python

   import magtrack
   from magtrack.simulation import simulate_beads
   import numpy as np

Simulating bead stacks
----------------------

The :mod:`magtrack.simulation` module can synthesize bead image stacks for testing and benchmarking. The snippet below
matches generates a stack of 100 frames where the bead drifts laterally while oscillating
in the vertical direction:

.. code-block:: python

   n_frames = 100
   roi = 64
   nm_per_px = 100.0

   x_true = np.linspace(-200, 200, n_frames)
   y_true = 100.0 * np.sin(np.linspace(1, 6 * np.pi, n_frames))
   z_true = np.zeros_like(x_true)
   xyz_true = np.stack([x_true, y_true, z_true], axis=1)

   stack = simulate_beads(xyz_true, size_px=roi, nm_per_px=nm_per_px)

The result ``stack`` is a three-dimensional NumPy array with shape ``(roi, roi, n_frames)``. Throughout the examples we
convert between pixels and nanometers using ``nm_per_px`` to compare the algorithms against the known truth.

Estimating lateral positions
----------------------------

For a first-pass lateral estimate, use :func:`magtrack.center_of_mass`. The ``background`` keyword controls how the
intensity bias is removed before computing the centroid. By default, ``"median"`` subtraction removes the frame median
for robust centering; ``"none"`` leaves the data unchanged and ``"mean"`` subtracts the frame mean:

.. code-block:: python

   x_com_px, y_com_px = magtrack.center_of_mass(stack)

   x_com_nm = x_com_px * nm_per_px
   y_com_nm = y_com_px * nm_per_px

Because the simulated data provides the ground-truth positions, you can compute the per-frame localization error and
plot it alongside the measured trajectories.

Refining XY positions
---------------------

MagTrack offers multiple refinement strategies after the coarse center-of-mass estimate. Auto-convolution narrows in on
the peak intensity around the bead, and the sub-pixel variant increases accuracy using a localized curve fit:

.. code-block:: python

   x_ac_px, y_ac_px = magtrack.auto_conv(stack, x_com_px, y_com_px)
   x_acsp_px, y_acsp_px = magtrack.auto_conv_sub_pixel(stack, x_com_px, y_com_px)

   error_acsp_nm = np.sqrt((x_acsp_px * nm_per_px - x_true)**2 +
                           (y_acsp_px * nm_per_px - y_true)**2)

For even finer control, Quadrant Interpolation (QI) can be iterated. Starting from the same seed coordinates, refine the
result multiple times to progressively reduce the residual error:

.. code-block:: python

   x_qi1_px, y_qi1_px = magtrack.qi(stack, x_com_px, y_com_px)
   x_qi2_px, y_qi2_px = magtrack.qi(stack, x_qi1_px, y_qi1_px)
   x_qi3_px, y_qi3_px = magtrack.qi(stack, x_qi2_px, y_qi2_px)

Each pass returns a new set of sub-pixel coordinates. Comparing the error traces from the center-of-mass, auto
convolution, and QI routines highlights the accuracy/throughput trade-offs.

Working with radial profiles and axial lookup
---------------------------------------------

To recover axial motion, convert each frame into a radial intensity profile and compare it with a precomputed lookup
library (Z-LUT). Start by simulating a reference stack and scan through known ``z`` offsets to build a table of radial profiles:

.. code-block:: python

   z_reference = np.arange(-10_000, 10_100, 100)
   xyz_reference = np.column_stack([
       np.zeros_like(z_reference),
       np.zeros_like(z_reference),
       z_reference,
   ])
   reference_stack = simulate_beads(xyz_reference, size_px=roi, nm_per_px=nm_per_px)

   reference_profiles = magtrack.radial_profile(
       reference_stack,
       np.full(z_reference.shape, roi / 2, dtype=float),
       np.full(z_reference.shape, roi / 2, dtype=float),
   )

   zlut = np.vstack([z_reference, reference_profiles])

With ``zlut`` prepared, analyze a fresh recording by extracting its radial profiles and calling
:func:`magtrack.lookup_z` to obtain sub-frame axial positions:

.. code-block:: python

   eval_profiles = magtrack.radial_profile(
       stack,
       x_com_px,
       y_com_px,
   )
   z_fit = magtrack.lookup_z(eval_profiles, zlut)

Plotting ``z_fit`` against the simulated truth illustrates how closely the axial estimate tracks the bead.

Running the full pipeline
-------------------------

The :func:`magtrack.stack_to_xyzp` convenience function wraps the lateral and axial steps. Provide a stack of bead images
and the Z-LUT, and it returns the lateral coordinates, axial positions, and fitted radial profiles in one call:

.. code-block:: python

   x_fit, y_fit, z_fit, profiles_fit = magtrack.stack_to_xyzp(stack, zlut=zlut)

This is the quickest path when you already have a reference lookup table and want to process experimental recordings.
Inspect ``profiles_fit`` to debug how well the observed profiles match the lookup data.

Accelerating with GPUs
----------------------

All MagTrack functions can operate on either NumPy arrays (CPU) or CuPy arrays (GPU). To use the GPU, move the stack to
device memory and call the same functions.

.. code-block:: python

   import cupy as cp

   # Move stack to GPU
   stack_gpu = cp.asarray(stack)

   x_com_px_gpu, y_com_px_gpu = magtrack.center_of_mass(stack_gpu, background="median")

   # Move values back to CPU
   x_com_px = cp.asnumpy(x_com_px_gpu)
   y_com_px = cp.asnumpy(y_com_px_gpu)