Neuroglancer for alignment inspection

neuroglancer-alignment.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "032dfaec014146bd9acdbe8e4e38ec95",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>VBox</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "VBox(children=(Text(value='/home/leek/data/2018_04_12_point_clouds/fixed_pts.npy', layout=Layout(width='80%')), Text(value='/home/leek/data/2018_04_12_point_clouds/moving_pts.npy', layout=Layout(width='80%')), IntText(value=100)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import neuroglancer\n",
    "from ipywidgets import interact, FloatText, Text, VBox, IntText, Layout\n",
    "#\n",
    "# Path to N x 3 arrays of keypoints for fixed and moving images\n",
    "#\n",
    "layout = Layout(width=\"80%\")\n",
    "fixed_pts_path = Text(\"/home/leek/data/2018_04_12_point_clouds/fixed_pts.npy\",\n",
    "                     layout=layout)\n",
    "moving_pts_path = Text(\"/home/leek/data/2018_04_12_point_clouds/moving_pts.npy\",\n",
    "                      layout=layout)\n",
    "#\n",
    "# Dimensions of the cube enclosing them\n",
    "#\n",
    "hdim = IntText(100)\n",
    "widget = VBox([fixed_pts_path, moving_pts_path, hdim])\n",
    "display(widget)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Read in the points\n",
    "#\n",
    "fixed_pts = np.load(fixed_pts_path.value)\n",
    "moving_pts = np.load(moving_pts_path.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Estimate an X, Y and Z offset from density\n",
    "#\n",
    "fcenter = np.mean(fixed_pts, 0)\n",
    "mcenter = np.mean(moving_pts, 0)\n",
    "off = fcenter - mcenter\n",
    "\n",
    "cmoving_pts = moving_pts + off[np.newaxis, :]\n",
    "#\n",
    "# Find the bounds for the histogram\n",
    "#\n",
    "((zmin, zmax), (ymin, ymax), (xmin, xmax)) =\\\n",
    " [(min(np.min(fixed_pts[:, _]), np.min(cmoving_pts[:, _])), \n",
    "  max(np.max(fixed_pts[:, _]), np.max(cmoving_pts[:, _]))) for _ in range(3)]\n",
    "hrange = ((zmin, zmax), (ymin, ymax), (xmin, xmax))\n",
    "hbins = tuple([hdim.value] * 3)\n",
    "hfixed = np.histogramdd(fixed_pts, bins=hbins, range=hrange)[0]\n",
    "hfixed = hfixed / np.max(hfixed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Fire up the Neuroglancer viewer. Pick a random port.\n",
    "# To choose a port:\n",
    "# neuroglancer.set_server_bind_address(bind_port=###)\n",
    "#\n",
    "viewer = neuroglancer.Viewer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7d6c30ec73e74ca89b4b0078008ad0d2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatText(value=40.518392548531665, description='dx'), FloatText(value=-345.6909436715282, description='dy'), FloatText(value=-355.085833881951, description='dz'), FloatText(value=0.0, description='rot'), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<a href=\"http://127.0.0.1:35621/v/e274d4f9b9079fe61bc6c70bf364192b22db6d8a/\" target=\"_blank\">Viewer</a>"
      ],
      "text/plain": [
       "http://127.0.0.1:35621/v/e274d4f9b9079fe61bc6c70bf364192b22db6d8a/"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "# Here is all the magic. Translate, then rotate the histogram of the moving\n",
    "# point cloud around the center\n",
    "#\n",
    "def disp_rmoving(dx, dy, dz, rot):\n",
    "    #\n",
    "    # np.newaxis is used to choose the broadcast axis for the\n",
    "    # offset array\n",
    "    #\n",
    "    t = moving_pts + np.array((dz, dy, dx))[np.newaxis, :]\n",
    "    #\n",
    "    # Rotate around the center of mass\n",
    "    #\n",
    "    centerx = np.mean(t[:, 2])\n",
    "    centery = np.mean(t[:, 1])\n",
    "    tc = t - np.array((0, centery, centerx))[np.newaxis, :]\n",
    "    rrot = rot * np.pi / 180\n",
    "    #\n",
    "    # Do the rotation, then add the center of mass back in\n",
    "    #\n",
    "    trx = -np.sin(rrot) * tc[:, 1] + np.cos(rrot) * tc[:, 2] + centerx\n",
    "    tryy = np.cos(rrot) * tc[:, 1] + np.sin(rrot) * tc[:, 2] + centery\n",
    "    tr = np.column_stack((t[:, 0], tryy, trx))\n",
    "    hmoving = np.histogramdd(tr, bins=hbins, range=hrange)[0]\n",
    "    hmoving = hmoving / np.max(hmoving)\n",
    "    with viewer.txn() as s:\n",
    "        s.layers[\"fixed\"] = neuroglancer.ImageLayer(\n",
    "            source=neuroglancer.LocalVolume(hfixed),\n",
    "            shader=\"\"\"\n",
    "void main() {\n",
    "  emitRGB(vec3(toNormalized(getDataValue()),0, 0));\n",
    "}\n",
    "            \"\"\")\n",
    "        s.layers[\"rmoving\"] = neuroglancer.ImageLayer(\n",
    "            source=neuroglancer.LocalVolume(hmoving),\n",
    "            shader = \"\"\"\n",
    "void main() {\n",
    "  emitRGB(vec3(0, toNormalized(getDataValue()), 0));\n",
    "}\n",
    "            \"\"\")\n",
    "interact(disp_rmoving, dx=FloatText(off[2]), dy=FloatText(off[1]), dz=FloatText(off[0]),\n",
    "        rot = FloatText(0.0))\n",
    "viewer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Environment (conda_nuggt)",
   "language": "python",
   "name": "conda_nuggt"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}