{ "cells": [ { "cell_type": "markdown", "id": "e9f9be63-1543-4808-872b-7a8e3ad6ea45", "metadata": {}, "source": [ "# Measure distance along a center line\n", "A common question is how to determine distances of points along the centerline of an object. For this we can skeletonize the object, determine distances along the center line and then find the closes center line pixel to given point coordinates. \n", "\n", "See also:\n", "* https://scikit-image.org/docs/stable/auto_examples/edges/plot_skeleton.html" ] }, { "cell_type": "code", "execution_count": 1, "id": "faa6ce47-e33f-4e9f-8c55-29a73b7db43e", "metadata": {}, "outputs": [], "source": [ "import bia_bob" ] }, { "cell_type": "code", "execution_count": 2, "id": "d1d23dfa", "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "from skimage.io import imread\n", "import napari_segment_blobs_and_things_with_membranes as nsbatwm\n", "import napari_simpleitk_image_processing as nsitk\n", "import numpy as np\n", "import skan\n", "import stackview\n", "import pyclesperanto_prototype as cle\n", "from scipy.spatial.distance import euclidean" ] }, { "cell_type": "markdown", "id": "c7ee349b", "metadata": {}, "source": [ "## Starting point: a binary image\n", "We start using a binary image that looks like an arm." ] }, { "cell_type": "code", "execution_count": 3, "id": "435b72bd", "metadata": { "lines_to_next_cell": 2 }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
shape(100, 100)
dtypeuint16
size19.5 kB
min0
max1
\n", "\n", "
" ], "text/plain": [ "StackViewNDArray([[0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " ...,\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "binary_arm = imread(\"../../data/binary_arm.tif\")\n", "stackview.insight(binary_arm)" ] }, { "cell_type": "markdown", "id": "fbaed36b-f7c5-40a7-bef2-6ddad8747269", "metadata": {}, "source": [ "Furthermore, we continue with a list of coordinates in X/Y format:" ] }, { "cell_type": "code", "execution_count": 4, "id": "81cbd0ce-d30e-4a73-b6a6-b2d41f521ffc", "metadata": {}, "outputs": [], "source": [ "coordinates_xy = np.asarray([\n", " [80, 70],\n", " [70, 70],\n", " [60, 70],\n", " [40, 40]]).T" ] }, { "cell_type": "markdown", "id": "3a8b1fd9-c9b7-4533-814b-bef2de93c8a6", "metadata": {}, "source": [ "We next produce a label image where the given coordinates are labeled. The first coordinate (index=0 in the list) will be labeled with 1, the second with 2, and so on. Background pixels are 0. \n", "We use this label image for visualization and further down, we will also use this image to do the measurement." ] }, { "cell_type": "code", "execution_count": 5, "id": "95a916e7-21ba-402b-9335-d55e167f39ec", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAGgCAYAAAAD9NhnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAZvklEQVR4nO3df0yV5/3/8dcR5gEcnK0azpGKCgnGVtrUgTVDIyytdKufZI1JV39Vu/0xndpCTaoymkiMctQ/DGm+01azqItjmmUsc822yvqD1bDMloaWYYLbypS05UO6mnNstBDl+vzht/c8YpHDD98HzvOR3EnPdV/nnOtcNb687vd1zu1zzjkBAGBgkvUAAADJixACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmBmzEDpw4IDy8vKUlpamoqIivf3222P1VgCAcSp1LF705MmTqqys1IEDB7Ro0SK98sor+t73vqdz585p5syZgz63v79fH3/8sTIzM+Xz+cZieACAMeSc0+XLl5WTk6NJk+6w1nFj4OGHH3YbNmyIaZs7d67bvn37HZ/b1dXlJHFwcHBwjPOjq6vrjn/nj/rluL6+PrW0tKi8vDymvby8XM3NzQP69/b2KhqNeofjR70BYELIzMy8Y59RD6FPP/1U169fVzAYjGkPBoPq7u4e0D8cDisQCHjHnS7XAQDGh6GUVMZsY8Ktb+6cu+2AqqqqFIlEvKOrq2ushgQASDCjvjFh2rRpSklJGbDq6enpGbA6kiS/3y+/3z/awwAAjAOjvhKaPHmyioqK1NjYGNPe2NiokpKS0X47AMA4NiZbtLds2aKnn35axcXF+va3v61Dhw7p4sWL2rBhw1i8HQBgnBqTEHrqqaf0n//8Rzt37tQnn3yiwsJC/eEPf9CsWbPG4u0AAOOUzyXYnuhoNKpAIGA9DADACEUiEWVlZQ3ah9+OAwCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmEm1HgASz49//ONhP/fQoUOjOBIAEx0rIQCAGUIIAGCGEAIAmPE555z1IG4WjUYVCASshzHhjaTuM1zUi4DkEolElJWVNWgfVkIAADOEEADADJfjkoTF5bd4cbkOmFi4HAcASGiEEADATFwhFA6HtWDBAmVmZio7O1tPPPGEOjo6Yvo451RTU6OcnBylp6errKxM7e3tozpoAMDEEFdN6Lvf/a5WrFihBQsW6Nq1a6qurlZbW5vOnTunKVOmSJL27t2r3bt36+jRo5ozZ4527dqlv/zlL+ro6FBmZuYd34Oa0OgYDzWgwVAfAsa/odSE4vrtuD/96U8xj48cOaLs7Gy1tLRoyZIlcs6prq5O1dXVWr58uSTp2LFjCgaDqq+v1/r16we8Zm9vr3p7e73H0Wg0niEBAMaxEdWEIpGIJOmee+6RJHV2dqq7u1vl5eVeH7/fr9LSUjU3N9/2NcLhsAKBgHfk5uaOZEgAgHFk2CHknNOWLVu0ePFiFRYWSpK6u7slScFgMKZvMBj0zt2qqqpKkUjEO7q6uoY7JADAODPsWzls3rxZH3zwgc6cOTPgnM/ni3nsnBvQ9iW/3y+/3z/cYSS18V73Gcytn40aETAxDWsl9Oyzz+rUqVN68803NWPGDK89FApJ0oBVT09Pz4DVEQAAcYWQc06bN29WQ0OD3njjDeXl5cWcz8vLUygUUmNjo9fW19enpqYmlZSUjM6IAQATRlyX4zZt2qT6+nr97ne/U2ZmprfiCQQCSk9Pl8/nU2VlpWpra1VQUKCCggLV1tYqIyNDq1atGpMPMNFM5EtsI8HlOWBiiiuEDh48KEkqKyuLaT9y5IieeeYZSdLWrVt19epVbdy4UZcuXdLChQt1+vTpIX1HCACQXOIKoaF8r9Xn86mmpkY1NTXDHRMAIEnw23EAADPcysEYNaDRQY0ISDzcygEAkNAIIQCAGUIIAGCGmtBdQN3n7qI+BCQGakIAgIRGCAEAzAz7V7Tx1bj8BgBDw0oIAGCGEAIAmCGEAABm2KI9TNR9xg+2bAM22KINAEhohBAAwAwhBAAww/eEMOFxa3AgcbESAgCYIYQAAGYIIQCAGb4nNER8L2hioj4EjB2+JwQASGiEEADADCEEADBDCAEAzBBCAAAzhBAAwAwhBAAwQwgBAMwQQgAAM4QQAMAMt3JAUuM2D4AtVkIAADOEEADADCEEADBDCAEAzBBCAAAzhBAAwAxbtIGbsGUbuLtYCQEAzBBCAAAzhBAAwAw1IWAQN9eIqA8Bo4+VEADADCEEADBDCAEAzBBCAAAzhBAAwAwhBAAwwxZtYIj4SR9g9LESAgCYIYQAAGYIIQCAGWpCQzTY9f871QpuPQ8AuIGVEADADCEEADBDCAEAzPicc856EDeLRqMKBALWwxiRwWpE1IcmLr43BMSKRCLKysoatA8rIQCAGUIIAGCGEAIAmCGEAABmCCEAgJkRhVA4HJbP51NlZaXX5pxTTU2NcnJylJ6errKyMrW3t490nACACWjYP9vzzjvv6NChQ3rwwQdj2vft26f9+/fr6NGjmjNnjnbt2qWlS5eqo6NDmZmZIx7weMS2bFs1wVDs4//tHpP34VYPQPyGtRL6/PPPtXr1ah0+fFjf/OY3vXbnnOrq6lRdXa3ly5ersLBQx44d05UrV1RfX3/b1+rt7VU0Go05AADJYVghtGnTJi1btkyPPvpoTHtnZ6e6u7tVXl7utfn9fpWWlqq5ufm2rxUOhxUIBLwjNzd3OEMCAIxDcYfQiRMn9N577ykcDg8419194zJHMBiMaQ8Gg965W1VVVSkSiXhHV1dXvEMCAIxTcdWEurq6VFFRodOnTystLe0r+/l8vpjHzrkBbV/y+/3y+/3xDCPhxXMrB277MPbGqgZ0Jzf/v6Q+BNxeXCuhlpYW9fT0qKioSKmpqUpNTVVTU5NeeuklpaameiugW1c9PT09A1ZHAADEFUKPPPKI2tra1Nra6h3FxcVavXq1WltblZ+fr1AopMbGRu85fX19ampqUklJyagPHgAwvsV1OS4zM1OFhYUxbVOmTNHUqVO99srKStXW1qqgoEAFBQWqra1VRkaGVq1aNXqjnkC4KyuAZDbqt/feunWrrl69qo0bN+rSpUtauHChTp8+nbTfEQIAfLURh9Bbb70V89jn86mmpkY1NTUjfWkAwATHb8cBAMxwZ9W7jBpPcmKLNpIRd1YFACQ0QggAYIYQAgCYGfUt2hgc3wMCgP9iJQQAMEMIAQDMcDkOuAu46ypwe6yEAABmCCEAgBlCCABghpqQMbZsA0hmrIQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZn3POWQ/iZtFoVIFAwHoYCYOf8Uk+3OYBE0UkElFWVtagfVgJAQDMEEIAADOEEADADDWhcYT6EJA8JkJtkJoQACChEUIAADOEEADADLf3BoAEFE8NeDzXj1gJAQDMEEIAADNs0R7H2LIN4HYS5fIcW7QBAAmNEAIAmCGEAABmqAlNINSIANzKsj5ETQgAkNAIIQCAGUIIAGCGn+2ZQG6+9kt9CMB4wEoIAGCGEAIAmOFy3AR167ZMLs8BSESshAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCG7wklCb43BCARsRICAJghhAAAZgghAIAZakJJits+AEgErIQAAGYIIQCAGS7Hge3bAMywEgIAmCGEAABm4g6hjz76SGvWrNHUqVOVkZGhhx56SC0tLd5555xqamqUk5Oj9PR0lZWVqb29fVQHDQCYGOKqCV26dEmLFi3Sd77zHf3xj39Udna2/vWvf+kb3/iG12ffvn3av3+/jh49qjlz5mjXrl1aunSpOjo6lJmZOdrjxxigRgTgbokrhPbu3avc3FwdOXLEa5s9e7b338451dXVqbq6WsuXL5ckHTt2TMFgUPX19Vq/fv2A1+zt7VVvb6/3OBqNxvsZAADjVFyX406dOqXi4mI9+eSTys7O1vz583X48GHvfGdnp7q7u1VeXu61+f1+lZaWqrm5+bavGQ6HFQgEvCM3N3eYHwUAMN7EFUIffvihDh48qIKCAr322mvasGGDnnvuOf3iF7+QJHV3d0uSgsFgzPOCwaB37lZVVVWKRCLe0dXVNZzPAQAYh+K6HNff36/i4mLV1tZKkubPn6/29nYdPHhQa9eu9fr5fL6Y5znnBrR9ye/3y+/3xztu3EXUiACMlbhWQtOnT9f9998f03bffffp4sWLkqRQKCRJA1Y9PT09A1ZHAADEFUKLFi1SR0dHTNv58+c1a9YsSVJeXp5CoZAaGxu98319fWpqalJJSckoDBcAMJHEdTnu+eefV0lJiWpra/WDH/xAZ8+e1aFDh7zLNT6fT5WVlaqtrVVBQYEKCgpUW1urjIwMrVq1akw+AO4+foEbwGiJK4QWLFig3/72t6qqqtLOnTuVl5enuro6rV692uuzdetWXb16VRs3btSlS5e0cOFCnT59mu8IAQAG8DnnnPUgbhaNRhUIBKyHgSFiJQQktls3Ft1NkUhEWVlZg/bht+MAAGa4lQNG5E7/ymKlBGAwrIQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBm2aGNM8RM/AAbDSggAYIYQAgCYIYQAAGaoCeGu4Q6tAG7FSggAYIYQAgCYIYQAAGaoCcEMNSIArIQAAGYIIQCAGUIIAGCGmhASxp1uFT4Y6knA+MRKCABghhACAJjhchwmBG4ZAYxPrIQAAGYIIQCAGUIIAGCGmhAmnDtt9aZmdHs5z+TEPP746MdGI/mvV/4nLebx+le/MBrJf93z2OaYx5+99v+MRjIxsBICAJghhAAAZgghAIAZn3POWQ/iZtFoVIFAwHoYSBLUhzDRjeTnsEYqEokoKytr0D6shAAAZgghAIAZtmgjqXF3V8AWKyEAgBlCCABghhACAJihJgTcJJ7trNSPgJFjJQQAMEMIAQDMEEIAADOEEADADCEEADBDCAEAzLBFGxgmtnMDI8dKCABghhACAJghhAAAZrizKmCAGhHGkuXdVG/GnVUBAAmNEAIAmCGEAABmqAkBCYZ6EeKVKDWgW1ETAgAkNEIIAGCGEAIAmCGEAABmCCEAgBlCCABghls5AAmGW0QgmbASAgCYIYQAAGbiCqFr167pxRdfVF5entLT05Wfn6+dO3eqv7/f6+OcU01NjXJycpSenq6ysjK1t7eP+sABAONfXDWhvXv36uWXX9axY8c0b948vfvuu/rhD3+oQCCgiooKSdK+ffu0f/9+HT16VHPmzNGuXbu0dOlSdXR0KDMzc0w+BJCsbq0fUSOauBL1p3lGKq4Q+utf/6rvf//7WrZsmSRp9uzZ+tWvfqV3331X0o1VUF1dnaqrq7V8+XJJ0rFjxxQMBlVfX6/169cPeM3e3l719vZ6j6PR6LA/DABgfInrctzixYv1+uuv6/z585Kk999/X2fOnNHjjz8uSers7FR3d7fKy8u95/j9fpWWlqq5ufm2rxkOhxUIBLwjNzd3uJ8FADDOxLUS2rZtmyKRiObOnauUlBRdv35du3fv1sqVKyVJ3d3dkqRgMBjzvGAwqAsXLtz2NauqqrRlyxbvcTQaJYgAIEnEFUInT57U8ePHVV9fr3nz5qm1tVWVlZXKycnRunXrvH4+ny/mec65AW1f8vv98vv9wxg6gFvdqW5AzcjWRK3rjERcIfTCCy9o+/btWrFihSTpgQce0IULFxQOh7Vu3TqFQiFJN1ZE06dP957X09MzYHUEAEBcNaErV65o0qTYp6SkpHhbtPPy8hQKhdTY2Oid7+vrU1NTk0pKSkZhuACACcXFYd26de7ee+91r776quvs7HQNDQ1u2rRpbuvWrV6fPXv2uEAg4BoaGlxbW5tbuXKlmz59uotGo0N6j0gk4iRxcHBwcIzzIxKJ3PHv/LhCKBqNuoqKCjdz5kyXlpbm8vPzXXV1tevt7fX69Pf3ux07drhQKOT8fr9bsmSJa2trG/J7EEIcHBwcE+MYSgj5nHNOCSQajSoQCFgPAwAwQpFIRFlZWYP24bfjAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABghhACAJghhAAAZgghAIAZQggAYIYQAgCYIYQAAGYIIQCAGUIIAGCGEAIAmCGEAABmCCEAgBlCCABgJuFCyDlnPQQAwCgYyt/nCRdCly9fth4CAGAUDOXvc59LsKVHf3+/Pv74YznnNHPmTHV1dSkrK8t6WAkrGo0qNzeXeboD5mlomKehYZ4G55zT5cuXlZOTo0mTBl/rpN6lMQ3ZpEmTNGPGDEWjUUlSVlYW/5OHgHkaGuZpaJinoWGevlogEBhSv4S7HAcASB6EEADATMKGkN/v144dO+T3+62HktCYp6FhnoaGeRoa5mn0JNzGBABA8kjYlRAAYOIjhAAAZgghAIAZQggAYIYQAgCYSdgQOnDggPLy8pSWlqaioiK9/fbb1kMyEw6HtWDBAmVmZio7O1tPPPGEOjo6Yvo451RTU6OcnBylp6errKxM7e3tRiNODOFwWD6fT5WVlV4b83TDRx99pDVr1mjq1KnKyMjQQw89pJaWFu888yRdu3ZNL774ovLy8pSenq78/Hzt3LlT/f39Xh/maRS4BHTixAn3ta99zR0+fNidO3fOVVRUuClTprgLFy5YD83EY4895o4cOeL+/ve/u9bWVrds2TI3c+ZM9/nnn3t99uzZ4zIzM91vfvMb19bW5p566ik3ffp0F41GDUdu5+zZs2727NnuwQcfdBUVFV478+TcZ5995mbNmuWeeeYZ97e//c11dna6P//5z+6f//yn14d5cm7Xrl1u6tSp7tVXX3WdnZ3u17/+tfv617/u6urqvD7M08glZAg9/PDDbsOGDTFtc+fOddu3bzcaUWLp6elxklxTU5Nzzrn+/n4XCoXcnj17vD5ffPGFCwQC7uWXX7YappnLly+7goIC19jY6EpLS70QYp5u2LZtm1u8ePFXnmeebli2bJn70Y9+FNO2fPlyt2bNGucc8zRaEu5yXF9fn1paWlReXh7TXl5erubmZqNRJZZIJCJJuueeeyRJnZ2d6u7ujpkzv9+v0tLSpJyzTZs2admyZXr00Udj2pmnG06dOqXi4mI9+eSTys7O1vz583X48GHvPPN0w+LFi/X666/r/PnzkqT3339fZ86c0eOPPy6JeRotCfcr2p9++qmuX7+uYDAY0x4MBtXd3W00qsThnNOWLVu0ePFiFRYWSpI3L7ebswsXLtz1MVo6ceKE3nvvPb3zzjsDzjFPN3z44Yc6ePCgtmzZop/+9Kc6e/asnnvuOfn9fq1du5Z5+v+2bdumSCSiuXPnKiUlRdevX9fu3bu1cuVKSfx5Gi0JF0Jf8vl8MY+dcwPaktHmzZv1wQcf6MyZMwPOJfucdXV1qaKiQqdPn1ZaWtpX9kv2eerv71dxcbFqa2slSfPnz1d7e7sOHjyotWvXev2SfZ5Onjyp48ePq76+XvPmzVNra6sqKyuVk5OjdevWef2SfZ5GKuEux02bNk0pKSkDVj09PT0D/sWRbJ599lmdOnVKb775pmbMmOG1h0IhSUr6OWtpaVFPT4+KioqUmpqq1NRUNTU16aWXXlJqaqo3F8k+T9OnT9f9998f03bffffp4sWLkvjz9KUXXnhB27dv14oVK/TAAw/o6aef1vPPP69wOCyJeRotCRdCkydPVlFRkRobG2PaGxsbVVJSYjQqW845bd68WQ0NDXrjjTeUl5cXcz4vL0+hUChmzvr6+tTU1JRUc/bII4+ora1Nra2t3lFcXKzVq1ertbVV+fn5zJOkRYsWDdjif/78ec2aNUsSf56+dOXKlQF3BU1JSfG2aDNPo8RwU8RX+nKL9s9//nN37tw5V1lZ6aZMmeL+/e9/Ww/NxE9+8hMXCATcW2+95T755BPvuHLlitdnz549LhAIuIaGBtfW1uZWrlzJVlHnYnbHOcc8OXdj+3pqaqrbvXu3+8c//uF++ctfuoyMDHf8+HGvD/Pk3Lp169y9997rbdFuaGhw06ZNc1u3bvX6ME8jl5Ah5JxzP/vZz9ysWbPc5MmT3be+9S1vO3IyknTb48iRI16f/v5+t2PHDhcKhZzf73dLlixxbW1tdoNOELeGEPN0w+9//3tXWFjo/H6/mzt3rjt06FDMeebJuWg06ioqKtzMmTNdWlqay8/Pd9XV1a63t9frwzyNHPcTAgCYSbiaEAAgeRBCAAAzhBAAwAwhBAAwQwgBAMwQQgAAM4QQAMAMIQQAMEMIAQDMEEIAADOEEADAzP8BY7Z7K4xzPAoAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# draw the coordinates into an image; for visualization purposes\n", "blank_image = cle.create((binary_arm.shape))\n", "labeled_spots = coordinate_visualization = cle.pointlist_to_labelled_spots(coordinates_xy, blank_image)\n", "\n", "# show the labeled pixels on top of the binary image\n", "cle.imshow(binary_arm, continue_drawing=True, max_display_intensity=1)\n", "cle.imshow(labeled_spots, labels=True, alpha=0.6)" ] }, { "cell_type": "markdown", "id": "0ba19ff9-2c9e-4636-8c97-6d23b69c475c", "metadata": {}, "source": [ "## Pre-processing\n", "Before we can skeletonize the image, we need to fill the black holes in the white area." ] }, { "cell_type": "code", "execution_count": 6, "id": "fdd7664a-96b4-456f-9131-fce331a9fa08", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "n-sitk made image
\n", "\n", "\n", "\n", "\n", "\n", "
shape(100, 100)
dtypeuint16
size19.5 kB
min0
max1
\n", "\n", "
" ], "text/plain": [ "StackViewNDArray([[0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " ...,\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0]], dtype=uint16)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "filled_holes = nsitk.binary_fill_holes(binary_arm)\n", "filled_holes" ] }, { "cell_type": "markdown", "id": "16008b80", "metadata": {}, "source": [ "## Skeletonization\n", "The [skeleton of a binary image](https://scikit-image.org/docs/stable/auto_examples/edges/plot_skeleton.html) is a thin line in the center of white areas." ] }, { "cell_type": "code", "execution_count": 7, "id": "7c410305", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "<__array_function__ internals>:200: RuntimeWarning: Converting input from bool to for compatibility.\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "nsbatwm made image
\n", "\n", "\n", "\n", "\n", "\n", "
shape(100, 100)
dtypebool
size9.8 kB
minFalse
maxTrue
\n", "\n", "
" ], "text/plain": [ "StackViewNDArray([[False, False, False, ..., False, False, False],\n", " [False, False, False, ..., False, False, False],\n", " [False, False, False, ..., False, False, False],\n", " ...,\n", " [False, False, False, ..., False, False, False],\n", " [False, False, False, ..., False, False, False],\n", " [False, False, False, ..., False, False, False]])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "skeleton = nsbatwm.skeletonize(filled_holes)\n", "skeleton" ] }, { "cell_type": "markdown", "id": "cf66428b-9a4a-44c6-ac3c-5d60e97c5b02", "metadata": {}, "source": [ "## Distance along the skeleton\n", "We first scan the skeleton and retrieve a list of coordinates." ] }, { "cell_type": "code", "execution_count": 8, "id": "05cc0b63-2963-4bcc-8f66-19c1d9301840", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[33, 12],\n", " [33, 13],\n", " [32, 14],\n", " [32, 15],\n", " [32, 16]], dtype=int64)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def skeleton_to_pointlist(skeleton):\n", " # Convert binary image into a graph representation\n", " graph = skan.Skeleton(skeleton)\n", " \n", " # Extract the coordinates of the polyline from the graph\n", " polyline = graph.path_coordinates(0)\n", "\n", " return polyline\n", "\n", "polyline = skeleton_to_pointlist(skeleton)\n", "polyline[:5]" ] }, { "cell_type": "markdown", "id": "483abf93-e394-4cc5-92db-16adae1107ba", "metadata": {}, "source": [ "We then compute the distance of all pixels within the skeleton to one end of the skeleton. We have little influence to which end the measurement is done." ] }, { "cell_type": "code", "execution_count": 9, "id": "f459ace1-a85e-4fa5-ba26-ad31360791c8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0, 1.0, 2.414213562373095, 3.414213562373095, 4.414213562373095]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def distances_along_pointlist(polyline):\n", " # Initialize a list to store the distances\n", " distances = [0]\n", " \n", " # Compute the distance between each pair of consecutive points in the polyline\n", " # Iterate over pairs of points using zip and calculate the Euclidean distance \n", " # between 'point' and 'follower' for each pair\n", " # Append the distance to the list 'distances'\n", " for point, follower in zip(polyline[:-1], polyline[1:]):\n", " d = euclidean(point, follower)\n", " distances.append(np.sum(distances[-1]) + d)\n", "\n", " return distances\n", "\n", "distances = distances_along_pointlist(polyline)\n", "distances[:5]" ] }, { "cell_type": "markdown", "id": "18904d14-fa00-4668-b427-631ced5c808e", "metadata": {}, "source": [ "We can also visualize the distance along the skeleton as sanity check." ] }, { "cell_type": "code", "execution_count": 10, "id": "7eee496e-07e5-410c-b130-2e130954f1f8", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "cle._ image
\n", "\n", "\n", "\n", "\n", "\n", "
shape(100, 100)
dtypefloat32
size39.1 kB
min0.0
max83.01219
\n", "\n", "
" ], "text/plain": [ "cl.OCLArray([[0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " ...,\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cle.replace_intensities(cle.label_spots(skeleton), [0] + distances)" ] }, { "cell_type": "markdown", "id": "611fc9d6-de62-4ece-a57a-8ad2f3dbc2ac", "metadata": {}, "source": [ "## Distance matrix\n", "For computing a distance matrix from each point coordinate to each pixel on the skeleton, we need to reformat our data." ] }, { "cell_type": "code", "execution_count": 11, "id": "62e69151-1e21-4a83-a542-f1709611bba9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,\n", " 27, 28, 29, 30, 30, 30, 31, 32, 33, 33, 34, 35, 35, 36, 37, 38,\n", " 38, 39, 40, 40, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,\n", " 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,\n", " 68, 69, 70, 71, 72, 73, 74, 75],\n", " [33, 33, 32, 32, 32, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 35,\n", " 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,\n", " 52, 53, 54, 55, 56, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 60,\n", " 60, 60, 61, 61, 61, 62, 62, 62, 62, 62, 62, 63, 63, 64, 64, 64,\n", " 64, 64, 64, 65, 65, 65, 65, 64]], dtype=int64)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# format the [[Y1,X1], [Y2,X2], ...] list into\n", "# [[X1, X2, ...], [Y1, Y2, ...]] clesperanto format\n", "polyline_xy = polyline[:,::-1].T\n", "polyline_xy" ] }, { "cell_type": "markdown", "id": "ad53c610-9580-46fd-aa81-5a950c0560bb", "metadata": {}, "source": [ "A distance matrix gives a 2D represenation of all points in a given list to all points given in another list of points." ] }, { "cell_type": "code", "execution_count": 12, "id": "bfb3d62a-9bab-4ddf-a7dd-0d3222b80c91", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[77.41447 68.7968 60.60528 28.861738]\n", " [76.537575 67.955864 59.816383 27.89265 ]\n", " [76.15773 67.6757 59.665733 27.20294 ]\n", " [75.29276 66.85058 58.898216 26.24881 ]\n", " [74.431175 66.0303 58.137768 25.298222]\n", " [74.094536 65.802734 58.0517 24.698177]\n", " [73.246155 65. 57.31492 23.769728]]\n" ] } ], "source": [ "distance_matrix = cle.generate_distance_matrix(coordinates_xy, polyline_xy)[1:,1:]\n", "print(distance_matrix[:7])" ] }, { "cell_type": "code", "execution_count": 13, "id": "23077c16-88ce-4368-8ef4-42e9457bbf48", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(72, 4)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "distance_matrix.shape" ] }, { "cell_type": "markdown", "id": "0bd2efa2-bf0c-4489-ace3-44d4ea5ff10f", "metadata": {}, "source": [ "The closest point along the skeleton to each of our four points can be determined using argmin projection in the distance matrix." ] }, { "cell_type": "code", "execution_count": 14, "id": "5a05e891-5075-4432-b35a-a1971e548191", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([70., 67., 59., 24.], dtype=float32)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "indices_on_skeleton = np.asarray(cle.arg_minimum_y_projection(distance_matrix))[0]\n", "indices_on_skeleton" ] }, { "cell_type": "markdown", "id": "79e333b6-fb3e-4d8a-9244-c36afcb49eca", "metadata": {}, "source": [ "From these indices, we can read out the distance along the skeleton." ] }, { "cell_type": "code", "execution_count": 15, "id": "b1c48f7b-dac7-438d-9ce0-64a867baa99b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[81.59797974644663, 78.59797974644663, 69.76955262170044, 28.970562748477146]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "distances_along_skeleton = [distances[int(i)] for i in indices_on_skeleton]\n", "distances_along_skeleton" ] }, { "cell_type": "markdown", "id": "f7f9326f-bb05-4d02-bff5-4dfd7e0b4e6e", "metadata": {}, "source": [ "## A helper function to do it all\n", "We now formulate a function that takes the skeleton and the list of coordinates to all the things shown above." ] }, { "cell_type": "code", "execution_count": 16, "id": "93e42e09-a8f4-46f5-817a-59ff127b9412", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[81.59797974644663, 78.59797974644663, 69.76955262170044, 28.970562748477146]" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def distances_of_points_along_skeleton(skeleton, coordinates_xy):\n", " \"\"\"\n", " Takes a binary skeleton image and a list of coordinates in format [[X1, X2, ...], [Y1, Y2, ...]]\n", " It will then determine all distances of the points along the skeleton. \n", " \"\"\"\n", " polyline = skeleton_to_pointlist(skeleton)\n", " distances = distances_along_pointlist(polyline)\n", " \n", " polyline_xy = polyline[:,::-1].T\n", " distance_matrix = cle.generate_distance_matrix(coordinates_xy, polyline_xy)[1:,1:]\n", "\n", " indices_on_skeleton = np.asarray(cle.arg_minimum_y_projection(distance_matrix))[0]\n", " distances_along_skeleton = [distances[int(i)] for i in indices_on_skeleton]\n", "\n", " return distances_along_skeleton\n", "\n", "distances_of_points_along_skeleton(skeleton, coordinates_xy)" ] }, { "cell_type": "markdown", "id": "6abf1047-2c91-4295-92fb-5ddb312fe9e0", "metadata": {}, "source": [ "## A distance map of the binary image\n", "To check that our distance computation makes sense, also visually, we can compute the distance along the center line (skeleton) for all pixels in the binary image.\n", "Therefore we first label all pixels individually." ] }, { "cell_type": "code", "execution_count": 17, "id": "644f59f7-4881-4e32-80a2-e13e1b4073b4", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "cle._ image
\n", "\n", "\n", "\n", "\n", "\n", "
shape(100, 100)
dtypeuint32
size39.1 kB
min0.0
max2725.0
\n", "\n", "
" ], "text/plain": [ "cl.OCLArray([[0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " ...,\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0],\n", " [0, 0, 0, ..., 0, 0, 0]], dtype=uint32)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "all_pixels_labeled = cle.label_spots(binary_arm)\n", "all_pixels_labeled" ] }, { "cell_type": "markdown", "id": "c699816c-19df-4b15-86be-72d45a7beaf0", "metadata": {}, "source": [ "We then compute the distance along the skeleton for all labeled pixel like above and visualize this in an image." ] }, { "cell_type": "code", "execution_count": 18, "id": "0c671db4-76dd-4597-a662-bd2efda832a6", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "cle._ image
\n", "\n", "\n", "\n", "\n", "\n", "
shape(100, 100)
dtypefloat32
size39.1 kB
min0.0
max83.01219
\n", "\n", "
" ], "text/plain": [ "cl.OCLArray([[0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " ...,\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coordinates_all_xy = cle.labelled_spots_to_pointlist(all_pixels_labeled)\n", "distances_all = distances_of_points_along_skeleton(skeleton, coordinates_all_xy)\n", "distance_map = cle.replace_intensities(all_pixels_labeled, [0] + distances_all)\n", "distance_map" ] }, { "cell_type": "markdown", "id": "6b0364aa-aeb8-439a-bbb0-bcb0a7fd493d", "metadata": {}, "source": [ "Using [`stackview.picker`](https://github.com/haesleinhuepf/stackview#pick-intensities) we can hover with the mouse over the image and read out intensities. This only works in a Jupyter-like environment." ] }, { "cell_type": "code", "execution_count": 19, "id": "60603514-1a47-4ca4-b897-70f86e0a141d", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ef07ebf1366046809d0f57c0614aeb33", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(HBox(children=(VBox(children=(ImageWidget(height=400, width=400),)),)), IntSlider(value=50, des…" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "stackview.picker(distance_map, zoom_factor=4)" ] }, { "cell_type": "markdown", "id": "d7a9b8c2-65a5-4572-bdcc-29ef892497c5", "metadata": {}, "source": [ "## Exercise\n", "Compute the distance of the points with respect to the other end of the image. Hint: Flip the input image when computing the path along the skeleton." ] }, { "cell_type": "code", "execution_count": null, "id": "7e7f2575-d9e3-47ee-8919-6a192e8153eb", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "main_language": "python", "notebook_metadata_filter": "-all" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.16" } }, "nbformat": 4, "nbformat_minor": 5 }