Skip to content

Instantly share code, notes, and snippets.

@LeeKamentsky
Created April 13, 2018 12:41
Show Gist options
  • Save LeeKamentsky/c0298ede9276792d2aef83b6436b2806 to your computer and use it in GitHub Desktop.
Save LeeKamentsky/c0298ede9276792d2aef83b6436b2806 to your computer and use it in GitHub Desktop.
Neuroglancer for alignment inspection
Display the source blob
Display the rendered blob
Raw
{
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment