{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Subpixel interpolation errors\n\nImage correlation can easily be more accurate than pixel in displacements.\nHowever for measuring very small displacements, systematic interpolation errors\ncan appear that have a small bias towards integer or 0.5 pixel displacements.\n\nThis example illustrates and studies this systematic error.\n\nThis error is typically reduced using higher order interpolation (order 3 instead of order 1).\n\nWe will also show an elegant mitigation from [Wantz2025]_ called \"Shift-DVC\" that\nwe have implemented in 2D and 3D.\n\n\nWe will apply synthetic displacements to an existing pattern, and it is essential\nto use a different interpolator than what will be used for the correlation, for this \nwe will Fourier shifts.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Import modules\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy\nimport matplotlib.pyplot as plt\nimport spam.deformation\nimport spam.DIC\nimport spam.datasets"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Define function to shift images in Fourier space\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def fourier_shift(im, shift_z=0, shift_y=0, shift_x=0):\n    nz, ny, nx = im.shape\n\n    # FFT freq\n    kz = numpy.fft.fftfreq(nz)\n    ky = numpy.fft.fftfreq(ny)\n    kx = numpy.fft.fftfreq(nx)\n\n    # put them in the freq grid\n    kz, ky, kx = numpy.meshgrid(kz, ky, kx, indexing=\"ij\")\n\n    # phase shift\n    phase = numpy.exp(-2j * numpy.pi * (shift_z * kz + shift_y * ky + shift_x * kx))\n\n    # returb shifted\n    return numpy.real(numpy.fft.ifftn(numpy.fft.fftn(im) * phase))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we will load the data and synthetically\napply subpixel displacements from [0, 1] px\nin N steps\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Load data\nim = spam.datasets.loadSnow()[0:50, 0:50, 0:50]\n\nN = 20\nsteps = numpy.linspace(0, 1, N + 1)\n\n# N x im series of displaced images\nims = numpy.zeros((len(steps), *im.shape))\n\nfor n, step in enumerate(steps):\n    ims[n] = fourier_shift(im, shift_x=step)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we will do image correlations, varying\nthe shift option and the interpolation order\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for shift in [True, False]:\n    for order in [1, 3]:\n        displacements = numpy.zeros((N + 1, 3))\n        for i in range(N + 1):\n            reg = spam.DIC.register(\n                ims[0],\n                ims[i],\n                # margin=4,\n                interpolationOrder=order,\n                # deltaPhiMin=0.001,\n                shift=shift,\n                # PhiInit = spam.deformation.computePhi({'t': [0.,0.,-1.]})\n                # verbose=True\n            )\n            displacements[i] = reg[\"Phi\"][0:3, -1]\n        print(displacements[:, 2])\n        plt.plot(steps, displacements[:, 2] - steps, \".-\", label=f\"{order = }, {shift = }\")\nplt.xlabel(\"Applied displacement (px)\")\nplt.ylabel(\"Error (px)\\n Measured displacement - Applied Displacement\")\nplt.legend()\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}