{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deformations\n",
    "\n",
    "In PyCA, we have functionality for v-field (vector field) and h-field (deformation field) transformations. We set up a few dummy Image3Ds and Field3Ds to show the syntax."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import PyCA.Core as ca\n",
    "grid = ca.GridInfo(ca.Vec3Di(10, 10, 10), ca.Vec3Df(1.0, 1.0, 1.0), ca.Vec3Df(1.0, 1.0, 1.0))\n",
    "mType = ca.MEM_HOST\n",
    "\n",
    "# image/field to be deformed\n",
    "im = ca.Image3D(grid, mType)\n",
    "ca.SetMem(im, 1.0)\n",
    "f = ca.Field3D(grid, mType)\n",
    "ca.SetMem(f, ca.Vec3Df(1.0, 2.0, 3.0))\n",
    "\n",
    "# Deformation field (v-field)\n",
    "v = ca.Field3D(grid, mType)\n",
    "ca.SetToZero(v)\n",
    "# Deformation field (h-field)\n",
    "h = ca.Field3D(grid, mType)\n",
    "ca.SetToIdentity(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we just defined an image `im` and a field `f` that we will be deformed.\n",
    "The v-field `v` is set to 0, and the h-field `h` is set to x (identity).\n",
    "Each h-field has a corresponding v-field, and they are related by $h(x) = x + v(x)$.\n",
    "We can convert v-fields to h-fields and vice versa:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ca.VtoH_I(v) # v(x) = v(x) + x\n",
    "ca.HtoV_I(v) # v(x) = v(x) - x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying a deformation to an image/field is done via a trilinear interpolation (or bilinear, for 2D images/fields).  Note that as in typical image registration, the deformation being applied is the inverse deformation, i.e. function composition: $I_{def}(x) = I(h(x)) = I(x+v(x))$.\n",
    "The PyCA functions are straightforward:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alpha = 1.0 #scalar multiple for applying a v-field\n",
    "im_def = ca.Image3D(grid, mType)\n",
    "f_def = ca.Field3D(grid, mType)\n",
    "bg = ca.BACKGROUND_STRATEGY_PARTIAL_ZERO\n",
    "\n",
    "ca.ApplyV(im_def, im, v, alpha, bg) # im_def = im(x+alpha*v(x))\n",
    "ca.ApplyVInv(im_def, im, v, alpha, bg) # im_def = im(x-alpha*v(x))\n",
    "ca.ApplyH(im_def, im, h, bg) #im_def = im(h(x))\n",
    "ca.ApplyV(f_def, f, v, bg) #f_def = f(x+v(x))\n",
    "ca.ApplyVInv(f_def, f, v, bg) #f_def = f(x-v(x))\n",
    "ca.ApplyH(f_def, f, h, bg) #f_def = f(h(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, these deformations don't actually do anything (since we have `v` and `h` set to zero and identity); this is just to demonstrate the syntax.\n",
    "\n",
    "You will notice a few things here. First, ApplyV has a sometimes optional `alpha` as a constant multiple for the vector field. Second, we can choose a background strategy (what value to use for interpolations that lie outside the image boundary).  The options are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bg = ca.BACKGROUND_STRATEGY_CLAMP # set background to value at the image boundary\n",
    "bg = ca.BACKGROUND_STRATEGY_PARTIAL_ZERO # background = 0\n",
    "bg = ca.BACKGROUND_STRATEGY_PARTIAL_ID # background is identity (for fields only)\n",
    "bg = ca.BACKGROUND_STRATEGY_WRAP # wrap image (interpolation on the torus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deformation Composition\n",
    "You could do a composition of deformations using the previous functions, but since there is an intrinsic meaning to a v-field and h-field, we have composition functions that work specifically on these objects. This way, you don't need to set the appropriate background strategy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#initialize some variables\n",
    "h_out = ca.Field3D(grid, mType)\n",
    "h = ca.Field3D(grid, mType)\n",
    "h0 = ca.Field3D(grid, mType)\n",
    "h1 = ca.Field3D(grid, mType)\n",
    "v = ca.Field3D(grid, mType)\n",
    "ca.SetToIdentity(h)\n",
    "ca.SetToIdentity(h0)\n",
    "ca.SetToIdentity(h1)\n",
    "ca.SetToZero(v)\n",
    "\n",
    "#Apply Compositions\n",
    "ca.ComposeHH(h_out, h0, h1) # h_out = h0(h1(x))\n",
    "ca.ComposeHV(h_out, h, v, 1.0) # h_out = h(x+v(x))\n",
    "ca.ComposeHVInv(h_out, h, v, 1.0) # h_out = h(x-v(x))\n",
    "ca.ComposeVH(h_out, v, h, 1.0) # h_out = h(x)+v(h(x)) = (x+v(x))(h(x))\n",
    "ca.ComposeVInvH(h_out, v, h, 1.0) # h_out = h(x)-v(h(x)) = (x-v(x))(h(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have a list of deformations that we need composed, we can do this by going forward through the list or backward through the list; both methods are essentially equivalent. I will show this using `ComposeHH`, but doing it with a list of vector fields (or negative vector fields) is similar. \n",
    "\n",
    "For a list of 10 fields, our goal is to solve for $h(x) = h_0 \\circ h_1 ... \\circ h_9(x)$.\n",
    "\n",
    "Going forward though the fields, we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# set up list of h fields\n",
    "hlist = [ca.Field3D(grid, mType) for _ in xrange(10)]\n",
    "scratchF = ca.Field3D(grid, mType)\n",
    "for ht in hlist:\n",
    "    ca.SetToIdentity(ht)\n",
    "    \n",
    "# Compose list        \n",
    "ca.SetToIdentity(h)\n",
    "for ht in hlist:\n",
    "    ca.ComposeHH(scratchF, h, ht)   \n",
    "    h.swap(scratchF)           # h = scratchF, scratchF = h  \n",
    "                                                             \n",
    "# h = x                         : initializiation        \n",
    "# h = h(h0(x)) = h0(x)          : first iteration        \n",
    "# h = h(h1(x)) = h0(h1(x))                               \n",
    "# h = h(h2(x)) = h0(h1(h2(x)))                           \n",
    "# ...                                                    \n",
    "# h = h(h9(x)) = h0(h1(...(h8(h9))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice here that we need a scratch variable `scratchF` because `ComposeHH` must take three distinct Field3Ds. We then just call the swap method.\n",
    "\n",
    "To go backward through the list, we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ca.SetToIdentity(h)\n",
    "for ht in reversed(hlist):\n",
    "    ca.ComposeHH(scratchF, ht, h)\n",
    "    h.swap(scratchF)\n",
    "\n",
    "# h = x                    : initialization\n",
    "# h = h9(h(x)) = h9(x)     : first iteration\n",
    "# h = h8(h(x)) = h8(h9(x)\n",
    "# h = h7(h(x)) = h7(h8(h9(x)))\n",
    "# ...\n",
    "# h = h0(h(x)) = h0(h1(...(h8(h9))))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}