diff --git a/demo/segmentation/level_set_segmentation.ipynb b/demo/segmentation/level_set_segmentation.ipynb new file mode 100644 index 00000000..3e76147d --- /dev/null +++ b/demo/segmentation/level_set_segmentation.ipynb @@ -0,0 +1,354 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Morphological snakes algorithm for cell segmentation\n", + "This notebook demonstrates how to use the GPU-implementation of the `morphological_chan_vese` segmentation algorithm. It is useful for sementing objects that are surrounded by intensity gradients, such as cells surrounded by marked membranes in fluorescence microscopy images.\n", + "\n", + "See also\n", + "* [scikit-image documentation of morphological_chan_vese](https://scikit-image.org/docs/stable/api/skimage.segmentation.html#skimage.segmentation.morphological_chan_vese)\n", + "* [A Morphological Approach to Curvature-Based Evolution of Curves and Surfaces, Márquez-Neila et al](https://ieeexplore.ieee.org/document/6529072)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import numpy as np\n", + "import pyclesperanto_prototype as cle\n", + "from skimage.data import cells3d\n", + "from skimage.segmentation import morphological_chan_vese\n", + "import timeit\n", + "import stackview\n", + "from skimage.io import imread\n", + "cle.select_device(\"TX\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "image = cells3d()[30,0,...]\n", + "#cle.asarray(image)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is recommended to apply the algorithm to an image with equalized background intensity." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "\n", + "cle._ image
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
shape(256, 256)
dtypefloat32
size256.0 kB
min0.14973111
max13.251159
\n", + "\n", + "
" + ], + "text/plain": [ + "cl.OCLArray([[0.9178023 , 1.090274 , 1.4726648 , ..., 0.7241896 , 0.6729306 ,\n", + " 0.6711849 ],\n", + " [0.917227 , 1.0674586 , 1.6082848 , ..., 0.70604837, 0.71757126,\n", + " 0.69825596],\n", + " [0.93129945, 1.2519282 , 1.7861464 , ..., 0.83825123, 0.76400876,\n", + " 0.79552275],\n", + " ...,\n", + " [0.68131644, 0.71315765, 0.6904286 , ..., 1.7797264 , 1.7004015 ,\n", + " 1.8670223 ],\n", + " [0.7588168 , 0.8026982 , 0.8340448 , ..., 1.5616127 , 1.7600938 ,\n", + " 1.6124313 ],\n", + " [0.71979123, 0.7125934 , 0.85827726, ..., 1.2010278 , 1.3635693 ,\n", + " 1.3926189 ]], dtype=float32)" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "background_subtracted = cle.divide_by_gaussian_background(image, sigma_x=15, sigma_y=15)\n", + "background_subtracted" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "
\n", + "\n", + "\n", + "cle._ image
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
shape(256, 256)
dtypefloat32
size256.0 kB
min0.0
max1.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=float32)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result_gpu = cle.morphological_chan_vese(background_subtracted, num_iter=3, smoothing=5)\n", + "result_gpu" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When executing this notebook locally, you can tune parameters manually here:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "97eee12823da4277bddbeaa87ac08366", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(interactive(children=(IntSlider(value=100, description='num_iter'), IntSlider(value=1, descript…" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stackview.interact(cle.morphological_chan_vese, background_subtracted)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Benchmarking\n", + "The GPU-implementation is supposed to be faster than the CPU-version in scikit-image. This depends on the input data and the applied parameters. You can use the code section below to determine speedup for your specific use-case." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "image = cells3d()[30,0,...]" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "def level_set_gpu():\n", + " result_gpu = cle.morphological_chan_vese(image, num_iter=20, smoothing=10)\n", + "\n", + "def level_set_cpu():\n", + " result_cpu = morphological_chan_vese(image, num_iter=20, smoothing=10)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "GPU time = 1.839248200000001\n" + ] + } + ], + "source": [ + "time_in_s = timeit.timeit(level_set_gpu, number=3)\n", + "print(\"GPU time = \", time_in_s)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU time = 1.938188799999999\n" + ] + } + ], + "source": [ + "time_in_s = timeit.timeit(level_set_cpu, number=3)\n", + "print(\"CPU time = \", time_in_s)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "result_gpu = cle.morphological_chan_vese(image, num_iter=20, smoothing=10)\n", + "#result_gpu" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "result_cpu = morphological_chan_vese(image, num_iter=20, smoothing=10)\n", + "#stackview.insight(result_cpu)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "3b476c74083f46958693634b4632bf2e", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(HBox(children=(VBox(children=(ImageWidget(height=256, width=256),)),)), IntSlider(value=128, de…" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stackview.curtain(result_cpu, np.asarray(result_gpu))" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cle.array_equal(result_cpu, result_gpu)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "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.15" + }, + "vscode": { + "interpreter": { + "hash": "f7e7e84a546d85c5beef5f7108b3ea7afa16a9e8e3877bd0c0d54e4df90a2b54" + } + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/demo/smoothing.ipynb b/demo/smoothing.ipynb new file mode 100644 index 00000000..fb2f0dc7 --- /dev/null +++ b/demo/smoothing.ipynb @@ -0,0 +1,408 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "eb73de20-5c0d-4cae-836c-ee3e240cad68", + "metadata": {}, + "outputs": [], + "source": [ + "import pyclesperanto_prototype as cle\n", + "import numpy as np\n", + "from skimage.segmentation.morphsnakes import sup_inf, inf_sup" + ] + }, + { + "cell_type": "markdown", + "id": "ff3849b7", + "metadata": {}, + "source": [ + "### 2D test" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "abc572bf", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAGgCAYAAABonJYyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAXQklEQVR4nO3df2xV9f3H8delVw4ta+8mpPdyR8tK0qRKNbJWyUpDSZQuky0xLkb5IRj/gQHKHYkUphuMjF5gWWO0AkIM2YYEsoxkuGwZnbpG0mw0ddWuLOBiB41607mRe+uQNtDP9w++nO1aanuh7X2vfT6S80c/53NvP/1w8/Tk3mMbcM45AQBMmpLtBQAAhkakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwLAxi/TevXtVUlKiadOmqaKiQm+99dZYfSsAmLCCY/Gkx44dUywW0969e7Vw4UK9/PLL+sY3vqEzZ86ouLj4cx87MDCgDz/8UPn5+QoEAmOxPADIOuecent7FY1GNWXK51wvuzFw3333ubVr16aNlZWVuS1btgz72O7ubieJg4ODY1Ic3d3dn9vEUX+7o7+/X21tbaqtrU0br62tVUtLy6D5fX19SqVS/uH4pXwAJpH8/PzPPT/qkf7444919epVhcPhtPFwOKxEIjFofjweVygU8o/h3g4BgIlkuLd1x+yDw89+Y+fcDRezdetWJZNJ/+ju7h6rJQHA/5xR/+Bw5syZysnJGXTV3NPTM+jqWpI8z5PneaO9DACYEEb9Snrq1KmqqKhQU1NT2nhTU5OqqqpG+9sBwIQ2Jrfgbdq0SY8//rgqKyv1ta99TQcOHNCFCxe0du3asfh2ADBhjUmkH330Uf3zn//Ujh079NFHH6m8vFy/+c1vNGfOnLH4dgAwYQWcsXveUqmUQqFQtpcBAOMimUyqoKBgyPP87g4AMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOC2V7AZOScy/YSAGRZKpVSKBQadh5X0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYFhGkY7H47r33nuVn5+vwsJCPfTQQzp79mzaHOectm/frmg0qtzcXC1evFidnZ2jumgAmCwyinRzc7PWr1+vP/7xj2pqatKVK1dUW1urf//73/6cPXv2qKGhQY2NjWptbVUkEtGSJUvU29s76osHgIku4G7hD+794x//UGFhoZqbm7Vo0SI55xSNRhWLxVRXVydJ6uvrUzgc1u7du7VmzZpBz9HX16e+vj7/61QqpaKioptd0v8E/sYhgOt/4zCZTKqgoGDIebf0nnQymZQk3X777ZKkrq4uJRIJ1dbW+nM8z1NNTY1aWlpu+BzxeFyhUMg/JnqgASATNx1p55w2bdqk6upqlZeXS5ISiYQkKRwOp80Nh8P+uc/aunWrksmkf3R3d9/skgBgwgne7AM3bNigd999V6dOnRp0LhAIpH3tnBs0dp3nefI872aXAQAT2k1dST/11FM6ceKE3nzzTc2ePdsfj0QikjToqrmnp2fQ1TUAYHgZRdo5pw0bNuj48eN64403VFJSkna+pKREkUhETU1N/lh/f7+am5tVVVU1OisGgEkko7c71q9fryNHjuhXv/qV8vPz/SvmUCik3NxcBQIBxWIx1dfXq7S0VKWlpaqvr1deXp6WL18+Jj8AAExkGd2CN9T7yocOHdITTzwh6drV9g9/+EO9/PLLunjxohYsWKCXXnrJ/3BxONdvS5nIuAUPwEhvwbul+6THApEGMBmMy33SAICxRaQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYcFsLwAYSiAQyPYScJOcc9lewoTBlTQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABh2S5GOx+MKBAKKxWL+mHNO27dvVzQaVW5urhYvXqzOzs5bXScATEo3HenW1lYdOHBAd999d9r4nj171NDQoMbGRrW2tioSiWjJkiXq7e295cUCwGRzU5H+5JNPtGLFCh08eFBf+tKX/HHnnJ5//nk9++yzevjhh1VeXq6f/vSnunTpko4cOXLD5+rr61MqlUo7AADX3FSk169fr6VLl+qBBx5IG+/q6lIikVBtba0/5nmeampq1NLScsPnisfjCoVC/lFUVHQzSwKACSnjSB89elRvv/224vH4oHOJREKSFA6H08bD4bB/7rO2bt2qZDLpH93d3ZkuCQAmrGAmk7u7u7Vx40adPHlS06ZNG3JeIBBI+9o5N2jsOs/z5HleJssAgEkjoyvptrY29fT0qKKiQsFgUMFgUM3NzXrhhRcUDAb9K+jPXjX39PQMuroGAAwvo0jff//96ujoUHt7u39UVlZqxYoVam9v19y5cxWJRNTU1OQ/pr+/X83Nzaqqqhr1xQPARJfR2x35+fkqLy9PG5s+fbpmzJjhj8diMdXX16u0tFSlpaWqr69XXl6eli9fPnqrBoBJIqNIj8TmzZv16aefat26dbp48aIWLFigkydPKj8/f7S/FQBMeAHnnMv2Iv5bKpVSKBTK9jLGlLEtN2uoD5thH6/x4V1vXTKZVEFBwZDz+N0dAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYRlH+oMPPtDKlSs1Y8YM5eXl6Z577lFbW5t/3jmn7du3KxqNKjc3V4sXL1ZnZ+eoLhoAJouMIn3x4kUtXLhQt912m37729/qzJkz+slPfqIvfvGL/pw9e/aooaFBjY2Nam1tVSQS0ZIlS9Tb2zvaaweAic9loK6uzlVXVw95fmBgwEUiEbdr1y5/7PLlyy4UCrn9+/ff8DGXL192yWTSP7q7u52kCX1gZLL978TBa3wsJZNJJ8klk8nPnZfRlfSJEydUWVmpRx55RIWFhZo/f74OHjzon+/q6lIikVBtba0/5nmeampq1NLScsPnjMfjCoVC/lFUVJTJkgBgQsso0u+//7727dun0tJS/e53v9PatWv19NNP62c/+5kkKZFISJLC4XDa48LhsH/us7Zu3apkMukf3d3dN/NzAMCEFMxk8sDAgCorK1VfXy9Jmj9/vjo7O7Vv3z6tWrXKnxcIBNIe55wbNHad53nyPC/TdQPApJDRlfSsWbN05513po3dcccdunDhgiQpEolI0qCr5p6enkFX1wCA4WUU6YULF+rs2bNpY+fOndOcOXMkSSUlJYpEImpqavLP9/f3q7m5WVVVVaOwXACYZDL5NPL06dMuGAy6nTt3uvfee8+9+uqrLi8vzx0+fNifs2vXLhcKhdzx48ddR0eHW7ZsmZs1a5ZLpVIZfeI5kQ+MTLb/nTh4jY+lkd7dkfFuvvbaa668vNx5nufKysrcgQMH0s4PDAy4bdu2uUgk4jzPc4sWLXIdHR0ZL3wiHxiZbP87cfAaH0sjjXTAOedkSCqVUigUyvYyxpSxLTdrqA+bYR+v8eFdb10ymVRBQcGQ8/jdHQBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADMvot+AB44n/IQLgShoATCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIZlFOkrV67oueeeU0lJiXJzczV37lzt2LFDAwMD/hznnLZv365oNKrc3FwtXrxYnZ2do75wAJgMMor07t27tX//fjU2Nuqvf/2r9uzZox//+Md68cUX/Tl79uxRQ0ODGhsb1draqkgkoiVLlqi3t3fUFw8AE13AOedGOvmb3/ymwuGwXnnlFX/s29/+tvLy8vTzn/9czjlFo1HFYjHV1dVJkvr6+hQOh7V7926tWbNm0HP29fWpr6/P/zqVSqmoqOhWfibzMthyABNUKpVSKBRSMplUQUHBkPMyupKurq7W66+/rnPnzkmS3nnnHZ06dUoPPvigJKmrq0uJREK1tbX+YzzPU01NjVpaWm74nPF4XKFQyD8meqABIBPBTCbX1dUpmUyqrKxMOTk5unr1qnbu3Klly5ZJkhKJhCQpHA6nPS4cDuv8+fM3fM6tW7dq06ZN/teT4UoaAEYqo0gfO3ZMhw8f1pEjRzRv3jy1t7crFospGo1q9erV/rxAIJD2OOfcoLHrPM+T53k3sXQAmPgyivQzzzyjLVu26LHHHpMk3XXXXTp//rzi8bhWr16tSCQi6doV9axZs/zH9fT0DLq6BgAML6P3pC9duqQpU9IfkpOT49+CV1JSokgkoqamJv98f3+/mpubVVVVNQrLBYDJJaMr6W9961vauXOniouLNW/ePP35z39WQ0ODnnzySUnX3uaIxWKqr69XaWmpSktLVV9fr7y8PC1fvnxMfgAAmMgyivSLL76o73//+1q3bp16enoUjUa1Zs0a/eAHP/DnbN68WZ9++qnWrVunixcvasGCBTp58qTy8/NHffEAMNFldJ/0eLh+7+BEZmzLAWTBmNwnDQAYX0QaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMCyj/+MQo2Oo3wgIAJ/FlTQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwzF2nnXLaXAADjZrjmmYt0b29vtpcAAONmuOYFnLFL14GBAX344Ydyzqm4uFjd3d0qKCjI9rLMSqVSKioqYp+GwT6NDPs0MqOxT8459fb2KhqNasqUoa+Xgze7yLEyZcoUzZ49W6lUSpJUUFDAi2UE2KeRYZ9Ghn0amVvdp1AoNOwcc293AAD+g0gDgGFmI+15nrZt2ybP87K9FNPYp5Fhn0aGfRqZ8dwncx8cAgD+w+yVNACASAOAaUQaAAwj0gBgGJEGAMPMRnrv3r0qKSnRtGnTVFFRobfeeivbS8qaeDyue++9V/n5+SosLNRDDz2ks2fPps1xzmn79u2KRqPKzc3V4sWL1dnZmaUV2xCPxxUIBBSLxfwx9umaDz74QCtXrtSMGTOUl5ene+65R21tbf559km6cuWKnnvuOZWUlCg3N1dz587Vjh07NDAw4M8Zl31yBh09etTddttt7uDBg+7MmTNu48aNbvr06e78+fPZXlpWfP3rX3eHDh1yf/nLX1x7e7tbunSpKy4udp988ok/Z9euXS4/P9/98pe/dB0dHe7RRx91s2bNcqlUKosrz57Tp0+7r3zlK+7uu+92Gzdu9MfZJ+f+9a9/uTlz5rgnnnjC/elPf3JdXV3u97//vfvb3/7mz2GfnPvRj37kZsyY4X7961+7rq4u94tf/MJ94QtfcM8//7w/Zzz2yWSk77vvPrd27dq0sbKyMrdly5YsrciWnp4eJ8k1Nzc755wbGBhwkUjE7dq1y59z+fJlFwqF3P79+7O1zKzp7e11paWlrqmpydXU1PiRZp+uqaurc9XV1UOeZ5+uWbp0qXvyySfTxh5++GG3cuVK59z47ZO5tzv6+/vV1tam2tratPHa2lq1tLRkaVW2JJNJSdLtt98uSerq6lIikUjbM8/zVFNTMyn3bP369Vq6dKkeeOCBtHH26ZoTJ06osrJSjzzyiAoLCzV//nwdPHjQP88+XVNdXa3XX39d586dkyS98847OnXqlB588EFJ47dP5n4L3scff6yrV68qHA6njYfDYSUSiSytyg7nnDZt2qTq6mqVl5dLkr8vN9qz8+fPj/sas+no0aN6++231draOugc+3TN+++/r3379mnTpk363ve+p9OnT+vpp5+W53latWoV+/T/6urqlEwmVVZWppycHF29elU7d+7UsmXLJI3f68lcpK8LBAJpXzvnBo1NRhs2bNC7776rU6dODTo32fesu7tbGzdu1MmTJzVt2rQh5032fRoYGFBlZaXq6+slSfPnz1dnZ6f27dunVatW+fMm+z4dO3ZMhw8f1pEjRzRv3jy1t7crFospGo1q9erV/ryx3idzb3fMnDlTOTk5g66ae3p6Bv0Xa7J56qmndOLECb355puaPXu2Px6JRCRp0u9ZW1ubenp6VFFRoWAwqGAwqObmZr3wwgsKBoP+Xkz2fZo1a5buvPPOtLE77rhDFy5ckMTr6bpnnnlGW7Zs0WOPPaa77rpLjz/+uL773e8qHo9LGr99MhfpqVOnqqKiQk1NTWnjTU1NqqqqytKqsss5pw0bNuj48eN64403VFJSkna+pKREkUgkbc/6+/vV3Nw8qfbs/vvvV0dHh9rb2/2jsrJSK1asUHt7u+bOncs+SVq4cOGgWzjPnTunOXPmSOL1dN2lS5cG/cWUnJwc/xa8cdunUfsIchRdvwXvlVdecWfOnHGxWMxNnz7d/f3vf8/20rLiO9/5jguFQu4Pf/iD++ijj/zj0qVL/pxdu3a5UCjkjh8/7jo6OtyyZcsm3S1TN/Lfd3c4xz45d+32xGAw6Hbu3Onee+899+qrr7q8vDx3+PBhfw775Nzq1avdl7/8Zf8WvOPHj7uZM2e6zZs3+3PGY59MRto551566SU3Z84cN3XqVPfVr37Vv91sMpJ0w+PQoUP+nIGBAbdt2zYXiUSc53lu0aJFrqOjI3uLNuKzkWafrnnttddceXm58zzPlZWVuQMHDqSdZ5+cS6VSbuPGja64uNhNmzbNzZ071z377LOur6/PnzMe+8TvkwYAw8y9Jw0A+A8iDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAw7P8AGRg9QTU+PcYAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "array = np.zeros((100, 85), dtype=np.uint8)\n", + "array[15:85, 15:85] = 1\n", + "array[35:65, 35:65] = 0\n", + "cle.imshow(array)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "560a1217", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pyclesperanto\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAGgCAYAAABonJYyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAXdUlEQVR4nO3df2xV9f3H8delV05b1t5NCPdyx4+VpEmVamStkpWGkihdJltiXIzyQzD+AwOUOxIpTDcYWXsLy4hRBIQYsg0JZBnJcNkyOnWNpNlo6qpdWcDFDhr1pnMj99YhbaCf7x98OfNSKr3Q9r7tfT6S8wfnfO7tp59enx7OPdwGnHNOAACTJmR7AgCAoRFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMGzUIr17926VlJQoPz9fFRUVeuutt0brSwHAuBUcjSc9cuSIYrGYdu/erfnz5+vll1/Wt771LZ06dUozZ8783McODAzoww8/VFFRkQKBwGhMDwCyzjmn3t5eRaNRTZjwOefLbhTcd999bvXq1Wn7ysrK3KZNm2742O7ubieJjY2NLSe27u7uz23iiF/u6O/vV1tbm2pra9P219bWqqWlZdD4vr4+pVIpf3N8KB+AHFJUVPS5x0c80h9//LEuX76scDictj8cDiuRSAwaH4/HFQqF/O1Gl0MAYDy50WXdUXvj8Nov7Jy77mQ2b96sZDLpb93d3aM1JQD4whnxNw6nTJmivLy8QWfNPT09g86uJcnzPHmeN9LTAIBxYcTPpCdOnKiKigo1NTWl7W9qalJVVdVIfzkAGNdG5Ra8DRs26PHHH1dlZaW+8Y1vaN++fTp37pxWr149Gl8OAMatUYn0o48+qn//+9/atm2bPvroI5WXl+t3v/udZs2aNRpfDgDGrYAzds9bKpVSKBTK9jQAYEwkk0kVFxcPeZzP7gAAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMCyY7QnkOudctqcAYIwEAoGMH8OZNAAYRqQBwDAiDQCGEekscM75G4Dc8dn/9pPJ5LAeQ6QBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMyyjS8Xhc9957r4qKijR16lQ99NBDOn36dNoY55y2bt2qaDSqgoICLVy4UJ2dnSM6aQDIFRlFurm5WWvXrtWf//xnNTU16dKlS6qtrdV///tff8yOHTu0c+dO7dq1S62trYpEIlq0aJF6e3tHfPIAMN4F3C18qPG//vUvTZ06Vc3NzVqwYIGcc4pGo4rFYqqrq5Mk9fX1KRwOa/v27Vq1atWg5+jr61NfX5//51QqpRkzZtzslL4Q+BxpAKlUSqFQSMlkUsXFxUOOu6Vr0lc/tPr222+XJHV1dSmRSKi2ttYf43meampq1NLSct3niMfjCoVC/jbeAw0AmbjpSDvntGHDBlVXV6u8vFySlEgkJEnhcDhtbDgc9o9da/PmzUomk/7W3d19s1MCgHEneLMPXLdund59912dOHFi0LFrf225c27IX2XueZ48z7vZaQDAuHZTZ9JPPfWUjh07pjfffFPTp0/390ciEUkadNbc09Mz6OwaAHBjGUXaOad169bp6NGjeuONN1RSUpJ2vKSkRJFIRE1NTf6+/v5+NTc3q6qqamRmDAA5JKPLHWvXrtWhQ4f0m9/8RkVFRf4ZcygUUkFBgQKBgGKxmBoaGlRaWqrS0lI1NDSosLBQS5cuHZVvAADGs4xuwRvquvKBAwf0xBNPSLpytv3jH/9YL7/8ss6fP6958+bppZde8t9cvJGrt6WMZ9yCB2C4t+Dd0n3So4FIA8gFY3KfNABgdBFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYFgw2xMAhiMQCGR7CvgczrlsT2Hc4kwaAAwj0gBgGJEGAMO4Jg2zuA79xXHtz4pr1COHM2kAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDbinS8XhcgUBAsVjM3+ec09atWxWNRlVQUKCFCxeqs7PzVucJADnppiPd2tqqffv26e67707bv2PHDu3cuVO7du1Sa2urIpGIFi1apN7e3lueLADkmpuK9CeffKJly5Zp//79+spXvuLvd87p+eef17PPPquHH35Y5eXl+vnPf64LFy7o0KFD132uvr4+pVKptA0AcMVNRXrt2rVavHixHnjggbT9XV1dSiQSqq2t9fd5nqeamhq1tLRc97ni8bhCoZC/zZgx42amBADjUsaRPnz4sN5++23F4/FBxxKJhCQpHA6n7Q+Hw/6xa23evFnJZNLfuru7M50SAIxbGf36rO7ubq1fv17Hjx9Xfn7+kOOu96t0hvpVSJ7nyfO8TKYBADkjozPptrY29fT0qKKiQsFgUMFgUM3NzXrhhRcUDAb9M+hrz5p7enoGnV0DAG4so0jff//96ujoUHt7u79VVlZq2bJlam9v1+zZsxWJRNTU1OQ/pr+/X83NzaqqqhrxyQPAeJfR5Y6ioiKVl5en7Zs0aZImT57s74/FYmpoaFBpaalKS0vV0NCgwsJCLV26dORmDQA5IqNID8fGjRv16aefas2aNTp//rzmzZun48ePq6ioaKS/FACMewHnnMv2JD4rlUopFAplexqjytiSmzXUm82wj9f4jV1tXTKZVHFx8ZDj+OwOADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAsIwj/cEHH2j58uWaPHmyCgsLdc8996itrc0/7pzT1q1bFY1GVVBQoIULF6qzs3NEJw0AuSKjSJ8/f17z58/Xbbfdpt///vc6deqUfvazn+nLX/6yP2bHjh3auXOndu3apdbWVkUiES1atEi9vb0jPXcAGP9cBurq6lx1dfWQxwcGBlwkEnGNjY3+vosXL7pQKOT27t173cdcvHjRJZNJf+vu7naSxvWG4cn2z4mN1/hoSiaTTpJLJpOfOy6jM+ljx46psrJSjzzyiKZOnaq5c+dq//79/vGuri4lEgnV1tb6+zzPU01NjVpaWq77nPF4XKFQyN9mzJiRyZQAYFzLKNLvv/++9uzZo9LSUv3hD3/Q6tWr9fTTT+sXv/iFJCmRSEiSwuFw2uPC4bB/7FqbN29WMpn0t+7u7pv5PgBgXApmMnhgYECVlZVqaGiQJM2dO1ednZ3as2ePVqxY4Y8LBAJpj3PODdp3led58jwv03kDQE7I6Ex62rRpuvPOO9P23XHHHTp37pwkKRKJSNKgs+aenp5BZ9cAgBvLKNLz58/X6dOn0/adOXNGs2bNkiSVlJQoEomoqanJP97f36/m5mZVVVWNwHQBIMdk8m7kyZMnXTAYdPX19e69995zr776qissLHQHDx70xzQ2NrpQKOSOHj3qOjo63JIlS9y0adNcKpXK6B3P8bxheLL9c2LjNT6ahnt3R8ar+dprr7ny8nLneZ4rKytz+/btSzs+MDDgtmzZ4iKRiPM8zy1YsMB1dHRkPPHxvGF4sv1zYuM1PpqGG+mAc87JkFQqpVAolO1pjCpjS27WUG82wz5e4zd2tXXJZFLFxcVDjuOzOwDAMCINAIYRaQAwjEgDgGFEGgAMy+ifhQNj6bN3CHCnh23czTF6OJMGAMOINAAYxuUOfCHw12nkKs6kAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgWEaRvnTpkp577jmVlJSooKBAs2fP1rZt2zQwMOCPcc5p69atikajKigo0MKFC9XZ2TniEweAXJBRpLdv3669e/dq165d+vvf/64dO3bopz/9qV588UV/zI4dO7Rz507t2rVLra2tikQiWrRokXp7e0d88gAw3gWcc264g7/97W8rHA7rlVde8fd997vfVWFhoX75y1/KOadoNKpYLKa6ujpJUl9fn8LhsLZv365Vq1YNes6+vj719fX5f06lUpoxY8atfE/mZbDkAMapVCqlUCikZDKp4uLiIcdldCZdXV2t119/XWfOnJEkvfPOOzpx4oQefPBBSVJXV5cSiYRqa2v9x3iep5qaGrW0tFz3OePxuEKhkL+N90ADQCaCmQyuq6tTMplUWVmZ8vLydPnyZdXX12vJkiWSpEQiIUkKh8NpjwuHwzp79ux1n3Pz5s3asGGD/+dcOJMGgOHKKNJHjhzRwYMHdejQIc2ZM0ft7e2KxWKKRqNauXKlPy4QCKQ9zjk3aN9VnufJ87ybmDoAjH8ZRfqZZ57Rpk2b9Nhjj0mS7rrrLp09e1bxeFwrV65UJBKRdOWMetq0af7jenp6Bp1dAwBuLKNr0hcuXNCECekPycvL82/BKykpUSQSUVNTk3+8v79fzc3NqqqqGoHpAkBuyehM+jvf+Y7q6+s1c+ZMzZkzR3/961+1c+dOPfnkk5KuXOaIxWJqaGhQaWmpSktL1dDQoMLCQi1dunRUvgEAGM8yivSLL76oH/7wh1qzZo16enoUjUa1atUq/ehHP/LHbNy4UZ9++qnWrFmj8+fPa968eTp+/LiKiopGfPIAMN5ldJ/0WLh67+B4ZmzJAWTBqNwnDQAYW0QaAAwj0gBgGJEGAMOINAAYltEteBgZn/0n8tzpAeSOoT4e4/NwJg0AhhFpADCMyx1ZdjN//QGQOziTBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYeYi7ZzL9hQAYMzcqHnmIt3b25vtKQDAmLlR8wLO2KnrwMCAPvzwQznnNHPmTHV3d6u4uDjb0zIrlUppxowZrNMNsE7DwzoNz0isk3NOvb29ikajmjBh6PPl4M1OcrRMmDBB06dPVyqVkiQVFxfzYhkG1ml4WKfhYZ2G51bXKRQK3XCMucsdAID/IdIAYJjZSHuepy1btsjzvGxPxTTWaXhYp+FhnYZnLNfJ3BuHAID/MXsmDQAg0gBgGpEGAMOINAAYRqQBwDCzkd69e7dKSkqUn5+viooKvfXWW9meUtbE43Hde++9Kioq0tSpU/XQQw/p9OnTaWOcc9q6daui0agKCgq0cOFCdXZ2ZmnGNsTjcQUCAcViMX8f63TFBx98oOXLl2vy5MkqLCzUPffco7a2Nv846yRdunRJzz33nEpKSlRQUKDZs2dr27ZtGhgY8MeMyTo5gw4fPuxuu+02t3//fnfq1Cm3fv16N2nSJHf27NlsTy0rvvnNb7oDBw64v/3tb669vd0tXrzYzZw5033yySf+mMbGRldUVOR+/etfu46ODvfoo4+6adOmuVQqlcWZZ8/Jkyfd1772NXf33Xe79evX+/tZJ+f+85//uFmzZrknnnjC/eUvf3FdXV3uj3/8o/vHP/7hj2GdnPvJT37iJk+e7H7729+6rq4u96tf/cp96Utfcs8//7w/ZizWyWSk77vvPrd69eq0fWVlZW7Tpk1ZmpEtPT09TpJrbm52zjk3MDDgIpGIa2xs9MdcvHjRhUIht3fv3mxNM2t6e3tdaWmpa2pqcjU1NX6kWacr6urqXHV19ZDHWacrFi9e7J588sm0fQ8//LBbvny5c27s1snc5Y7+/n61tbWptrY2bX9tba1aWlqyNCtbksmkJOn222+XJHV1dSmRSKStmed5qqmpyck1W7t2rRYvXqwHHnggbT/rdMWxY8dUWVmpRx55RFOnTtXcuXO1f/9+/zjrdEV1dbVef/11nTlzRpL0zjvv6MSJE3rwwQcljd06mfsUvI8//liXL19WOBxO2x8Oh5VIJLI0Kzucc9qwYYOqq6tVXl4uSf66XG/Nzp49O+ZzzKbDhw/r7bffVmtr66BjrNMV77//vvbs2aMNGzboBz/4gU6ePKmnn35anudpxYoVrNP/q6urUzKZVFlZmfLy8nT58mXV19dryZIlksbu9WQu0lcFAoG0PzvnBu3LRevWrdO7776rEydODDqW62vW3d2t9evX6/jx48rPzx9yXK6v08DAgCorK9XQ0CBJmjt3rjo7O7Vnzx6tWLHCH5fr63TkyBEdPHhQhw4d0pw5c9Te3q5YLKZoNKqVK1f640Z7ncxd7pgyZYry8vIGnTX39PQM+j9Wrnnqqad07Ngxvfnmm5o+fbq/PxKJSFLOr1lbW5t6enpUUVGhYDCoYDCo5uZmvfDCCwoGg/5a5Po6TZs2TXfeeWfavjvuuEPnzp2TxOvpqmeeeUabNm3SY489prvuukuPP/64vv/97ysej0sau3UyF+mJEyeqoqJCTU1NafubmppUVVWVpVlll3NO69at09GjR/XGG2+opKQk7XhJSYkikUjamvX396u5uTmn1uz+++9XR0eH2tvb/a2yslLLli1Te3u7Zs+ezTpJmj9//qBbOM+cOaNZs2ZJ4vV01YULFwb9xpS8vDz/FrwxW6cRewtyBF29Be+VV15xp06dcrFYzE2aNMn985//zPbUsuJ73/ueC4VC7k9/+pP76KOP/O3ChQv+mMbGRhcKhdzRo0ddR0eHW7JkSc7dMnU9n727wznWybkrtycGg0FXX1/v3nvvPffqq6+6wsJCd/DgQX8M6+TcypUr3Ve/+lX/FryjR4+6KVOmuI0bN/pjxmKdTEbaOedeeuklN2vWLDdx4kT39a9/3b/dLBdJuu524MABf8zAwIDbsmWLi0QizvM8t2DBAtfR0ZG9SRtxbaRZpytee+01V15e7jzPc2VlZW7fvn1px1kn51KplFu/fr2bOXOmy8/Pd7Nnz3bPPvus6+vr88eMxTrxedIAYJi5a9IAgP8h0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw/4PsfpbC3ly8HwAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "print(\"Pyclesperanto\")\n", + "array_smooth = cle.superior_inferior(cle.inferior_superior(array))\n", + "cle.imshow(array_smooth)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "480d0014", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Scikit-image\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAGgCAYAAABonJYyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAXdUlEQVR4nO3df2xV9f3H8delV05b1t5NCPdyx4+VpEmVamStkpWGkihdJltiXIzyQzD+AwOUOxIpTDcYWXsLy4hRBIQYsg0JZBnJcNkyOnWNpNlo6qpdWcDFDhr1pnMj99YhbaCf7x98OfNSKr3Q9r7tfT6S8wfnfO7tp59enx7OPdwGnHNOAACTJmR7AgCAoRFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMGzUIr17926VlJQoPz9fFRUVeuutt0brSwHAuBUcjSc9cuSIYrGYdu/erfnz5+vll1/Wt771LZ06dUozZ8783McODAzoww8/VFFRkQKBwGhMDwCyzjmn3t5eRaNRTZjwOefLbhTcd999bvXq1Wn7ysrK3KZNm2742O7ubieJjY2NLSe27u7uz23iiF/u6O/vV1tbm2pra9P219bWqqWlZdD4vr4+pVIpf3N8KB+AHFJUVPS5x0c80h9//LEuX76scDictj8cDiuRSAwaH4/HFQqF/O1Gl0MAYDy50WXdUXvj8Nov7Jy77mQ2b96sZDLpb93d3aM1JQD4whnxNw6nTJmivLy8QWfNPT09g86uJcnzPHmeN9LTAIBxYcTPpCdOnKiKigo1NTWl7W9qalJVVdVIfzkAGNdG5Ra8DRs26PHHH1dlZaW+8Y1vaN++fTp37pxWr149Gl8OAMatUYn0o48+qn//+9/atm2bPvroI5WXl+t3v/udZs2aNRpfDgDGrYAzds9bKpVSKBTK9jQAYEwkk0kVFxcPeZzP7gAAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMCyY7QnkOudctqcAYIwEAoGMH8OZNAAYRqQBwDAiDQCGEekscM75G4Dc8dn/9pPJ5LAeQ6QBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMyyjS8Xhc9957r4qKijR16lQ99NBDOn36dNoY55y2bt2qaDSqgoICLVy4UJ2dnSM6aQDIFRlFurm5WWvXrtWf//xnNTU16dKlS6qtrdV///tff8yOHTu0c+dO7dq1S62trYpEIlq0aJF6e3tHfPIAMN4F3C18qPG//vUvTZ06Vc3NzVqwYIGcc4pGo4rFYqqrq5Mk9fX1KRwOa/v27Vq1atWg5+jr61NfX5//51QqpRkzZtzslL4Q+BxpAKlUSqFQSMlkUsXFxUOOu6Vr0lc/tPr222+XJHV1dSmRSKi2ttYf43meampq1NLSct3niMfjCoVC/jbeAw0AmbjpSDvntGHDBlVXV6u8vFySlEgkJEnhcDhtbDgc9o9da/PmzUomk/7W3d19s1MCgHEneLMPXLdund59912dOHFi0LFrf225c27IX2XueZ48z7vZaQDAuHZTZ9JPPfWUjh07pjfffFPTp0/390ciEUkadNbc09Mz6OwaAHBjGUXaOad169bp6NGjeuONN1RSUpJ2vKSkRJFIRE1NTf6+/v5+NTc3q6qqamRmDAA5JKPLHWvXrtWhQ4f0m9/8RkVFRf4ZcygUUkFBgQKBgGKxmBoaGlRaWqrS0lI1NDSosLBQS5cuHZVvAADGs4xuwRvquvKBAwf0xBNPSLpytv3jH/9YL7/8ss6fP6958+bppZde8t9cvJGrt6WMZ9yCB2C4t+Dd0n3So4FIA8gFY3KfNABgdBFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYFgw2xMAhiMQCGR7CvgczrlsT2Hc4kwaAAwj0gBgGJEGAMO4Jg2zuA79xXHtz4pr1COHM2kAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDbinS8XhcgUBAsVjM3+ec09atWxWNRlVQUKCFCxeqs7PzVucJADnppiPd2tqqffv26e67707bv2PHDu3cuVO7du1Sa2urIpGIFi1apN7e3lueLADkmpuK9CeffKJly5Zp//79+spXvuLvd87p+eef17PPPquHH35Y5eXl+vnPf64LFy7o0KFD132uvr4+pVKptA0AcMVNRXrt2rVavHixHnjggbT9XV1dSiQSqq2t9fd5nqeamhq1tLRc97ni8bhCoZC/zZgx42amBADjUsaRPnz4sN5++23F4/FBxxKJhCQpHA6n7Q+Hw/6xa23evFnJZNLfuru7M50SAIxbGf36rO7ubq1fv17Hjx9Xfn7+kOOu96t0hvpVSJ7nyfO8TKYBADkjozPptrY29fT0qKKiQsFgUMFgUM3NzXrhhRcUDAb9M+hrz5p7enoGnV0DAG4so0jff//96ujoUHt7u79VVlZq2bJlam9v1+zZsxWJRNTU1OQ/pr+/X83NzaqqqhrxyQPAeJfR5Y6ioiKVl5en7Zs0aZImT57s74/FYmpoaFBpaalKS0vV0NCgwsJCLV26dORmDQA5IqNID8fGjRv16aefas2aNTp//rzmzZun48ePq6ioaKS/FACMewHnnMv2JD4rlUopFAplexqjytiSmzXUm82wj9f4jV1tXTKZVHFx8ZDj+OwOADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAsIwj/cEHH2j58uWaPHmyCgsLdc8996itrc0/7pzT1q1bFY1GVVBQoIULF6qzs3NEJw0AuSKjSJ8/f17z58/Xbbfdpt///vc6deqUfvazn+nLX/6yP2bHjh3auXOndu3apdbWVkUiES1atEi9vb0jPXcAGP9cBurq6lx1dfWQxwcGBlwkEnGNjY3+vosXL7pQKOT27t173cdcvHjRJZNJf+vu7naSxvWG4cn2z4mN1/hoSiaTTpJLJpOfOy6jM+ljx46psrJSjzzyiKZOnaq5c+dq//79/vGuri4lEgnV1tb6+zzPU01NjVpaWq77nPF4XKFQyN9mzJiRyZQAYFzLKNLvv/++9uzZo9LSUv3hD3/Q6tWr9fTTT+sXv/iFJCmRSEiSwuFw2uPC4bB/7FqbN29WMpn0t+7u7pv5PgBgXApmMnhgYECVlZVqaGiQJM2dO1ednZ3as2ePVqxY4Y8LBAJpj3PODdp3led58jwv03kDQE7I6Ex62rRpuvPOO9P23XHHHTp37pwkKRKJSNKgs+aenp5BZ9cAgBvLKNLz58/X6dOn0/adOXNGs2bNkiSVlJQoEomoqanJP97f36/m5mZVVVWNwHQBIMdk8m7kyZMnXTAYdPX19e69995zr776qissLHQHDx70xzQ2NrpQKOSOHj3qOjo63JIlS9y0adNcKpXK6B3P8bxheLL9c2LjNT6ahnt3R8ar+dprr7ny8nLneZ4rKytz+/btSzs+MDDgtmzZ4iKRiPM8zy1YsMB1dHRkPPHxvGF4sv1zYuM1PpqGG+mAc87JkFQqpVAolO1pjCpjS27WUG82wz5e4zd2tXXJZFLFxcVDjuOzOwDAMCINAIYRaQAwjEgDgGFEGgAMy+ifhQNj6bN3CHCnh23czTF6OJMGAMOINAAYxuUOfCHw12nkKs6kAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgWEaRvnTpkp577jmVlJSooKBAs2fP1rZt2zQwMOCPcc5p69atikajKigo0MKFC9XZ2TniEweAXJBRpLdv3669e/dq165d+vvf/64dO3bopz/9qV588UV/zI4dO7Rz507t2rVLra2tikQiWrRokXp7e0d88gAw3gWcc264g7/97W8rHA7rlVde8fd997vfVWFhoX75y1/KOadoNKpYLKa6ujpJUl9fn8LhsLZv365Vq1YNes6+vj719fX5f06lUpoxY8atfE/mZbDkAMapVCqlUCikZDKp4uLiIcdldCZdXV2t119/XWfOnJEkvfPOOzpx4oQefPBBSVJXV5cSiYRqa2v9x3iep5qaGrW0tFz3OePxuEKhkL+N90ADQCaCmQyuq6tTMplUWVmZ8vLydPnyZdXX12vJkiWSpEQiIUkKh8NpjwuHwzp79ux1n3Pz5s3asGGD/+dcOJMGgOHKKNJHjhzRwYMHdejQIc2ZM0ft7e2KxWKKRqNauXKlPy4QCKQ9zjk3aN9VnufJ87ybmDoAjH8ZRfqZZ57Rpk2b9Nhjj0mS7rrrLp09e1bxeFwrV65UJBKRdOWMetq0af7jenp6Bp1dAwBuLKNr0hcuXNCECekPycvL82/BKykpUSQSUVNTk3+8v79fzc3NqqqqGoHpAkBuyehM+jvf+Y7q6+s1c+ZMzZkzR3/961+1c+dOPfnkk5KuXOaIxWJqaGhQaWmpSktL1dDQoMLCQi1dunRUvgEAGM8yivSLL76oH/7wh1qzZo16enoUjUa1atUq/ehHP/LHbNy4UZ9++qnWrFmj8+fPa968eTp+/LiKiopGfPIAMN5ldJ/0WLh67+B4ZmzJAWTBqNwnDQAYW0QaAAwj0gBgGJEGAMOINAAYltEteBgZn/0n8tzpAeSOoT4e4/NwJg0AhhFpADCMyx1ZdjN//QGQOziTBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYUQaAAwj0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw4g0ABhGpAHAMCINAIYRaQAwjEgDgGFEGgAMI9IAYBiRBgDDiDQAGEakAcAwIg0AhhFpADCMSAOAYeYi7ZzL9hQAYMzcqHnmIt3b25vtKQDAmLlR8wLO2KnrwMCAPvzwQznnNHPmTHV3d6u4uDjb0zIrlUppxowZrNMNsE7DwzoNz0isk3NOvb29ikajmjBh6PPl4M1OcrRMmDBB06dPVyqVkiQVFxfzYhkG1ml4WKfhYZ2G51bXKRQK3XCMucsdAID/IdIAYJjZSHuepy1btsjzvGxPxTTWaXhYp+FhnYZnLNfJ3BuHAID/MXsmDQAg0gBgGpEGAMOINAAYRqQBwDCzkd69e7dKSkqUn5+viooKvfXWW9meUtbE43Hde++9Kioq0tSpU/XQQw/p9OnTaWOcc9q6daui0agKCgq0cOFCdXZ2ZmnGNsTjcQUCAcViMX8f63TFBx98oOXLl2vy5MkqLCzUPffco7a2Nv846yRdunRJzz33nEpKSlRQUKDZs2dr27ZtGhgY8MeMyTo5gw4fPuxuu+02t3//fnfq1Cm3fv16N2nSJHf27NlsTy0rvvnNb7oDBw64v/3tb669vd0tXrzYzZw5033yySf+mMbGRldUVOR+/etfu46ODvfoo4+6adOmuVQqlcWZZ8/Jkyfd1772NXf33Xe79evX+/tZJ+f+85//uFmzZrknnnjC/eUvf3FdXV3uj3/8o/vHP/7hj2GdnPvJT37iJk+e7H7729+6rq4u96tf/cp96Utfcs8//7w/ZizWyWSk77vvPrd69eq0fWVlZW7Tpk1ZmpEtPT09TpJrbm52zjk3MDDgIpGIa2xs9MdcvHjRhUIht3fv3mxNM2t6e3tdaWmpa2pqcjU1NX6kWacr6urqXHV19ZDHWacrFi9e7J588sm0fQ8//LBbvny5c27s1snc5Y7+/n61tbWptrY2bX9tba1aWlqyNCtbksmkJOn222+XJHV1dSmRSKStmed5qqmpyck1W7t2rRYvXqwHHnggbT/rdMWxY8dUWVmpRx55RFOnTtXcuXO1f/9+/zjrdEV1dbVef/11nTlzRpL0zjvv6MSJE3rwwQcljd06mfsUvI8//liXL19WOBxO2x8Oh5VIJLI0Kzucc9qwYYOqq6tVXl4uSf66XG/Nzp49O+ZzzKbDhw/r7bffVmtr66BjrNMV77//vvbs2aMNGzboBz/4gU6ePKmnn35anudpxYoVrNP/q6urUzKZVFlZmfLy8nT58mXV19dryZIlksbu9WQu0lcFAoG0PzvnBu3LRevWrdO7776rEydODDqW62vW3d2t9evX6/jx48rPzx9yXK6v08DAgCorK9XQ0CBJmjt3rjo7O7Vnzx6tWLHCH5fr63TkyBEdPHhQhw4d0pw5c9Te3q5YLKZoNKqVK1f640Z7ncxd7pgyZYry8vIGnTX39PQM+j9Wrnnqqad07Ngxvfnmm5o+fbq/PxKJSFLOr1lbW5t6enpUUVGhYDCoYDCo5uZmvfDCCwoGg/5a5Po6TZs2TXfeeWfavjvuuEPnzp2TxOvpqmeeeUabNm3SY489prvuukuPP/64vv/97ysej0sau3UyF+mJEyeqoqJCTU1NafubmppUVVWVpVlll3NO69at09GjR/XGG2+opKQk7XhJSYkikUjamvX396u5uTmn1uz+++9XR0eH2tvb/a2yslLLli1Te3u7Zs+ezTpJmj9//qBbOM+cOaNZs2ZJ4vV01YULFwb9xpS8vDz/FrwxW6cRewtyBF29Be+VV15xp06dcrFYzE2aNMn985//zPbUsuJ73/ueC4VC7k9/+pP76KOP/O3ChQv+mMbGRhcKhdzRo0ddR0eHW7JkSc7dMnU9n727wznWybkrtycGg0FXX1/v3nvvPffqq6+6wsJCd/DgQX8M6+TcypUr3Ve/+lX/FryjR4+6KVOmuI0bN/pjxmKdTEbaOedeeuklN2vWLDdx4kT39a9/3b/dLBdJuu524MABf8zAwIDbsmWLi0QizvM8t2DBAtfR0ZG9SRtxbaRZpytee+01V15e7jzPc2VlZW7fvn1px1kn51KplFu/fr2bOXOmy8/Pd7Nnz3bPPvus6+vr88eMxTrxedIAYJi5a9IAgP8h0gBgGJEGAMOINAAYRqQBwDAiDQCGEWkAMIxIA4BhRBoADCPSAGAYkQYAw/4PsfpbC3ly8HwAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "print(\"Scikit-image\")\n", + "array_smooth_ref = sup_inf(inf_sup(array))\n", + "cle.imshow(array_smooth_ref)" + ] + }, + { + "cell_type": "markdown", + "id": "48e384da", + "metadata": {}, + "source": [ + "### 3D test" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "d418a958", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 0. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0.]]]\n" + ] + } + ], + "source": [ + "array3d = np.zeros((5, 5, 5))\n", + "array3d[1:4, 1:4, 1:4] = 1\n", + "array3d[2, 2, 2] = 0\n", + "print(array3d)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "3d2b5f2f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pyclesperanto\n", + "\n", + "[[[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 1 1 1 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]]\n" + ] + } + ], + "source": [ + "print(\"Pyclesperanto\\n\")\n", + "print(cle.superior_inferior(cle.inferior_superior(array3d)))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "9c25d0db", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Scikit-image\n", + "\n", + "[[[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 1 1 1 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]]\n" + ] + } + ], + "source": [ + "print(\"Scikit-image\\n\")\n", + "print(sup_inf(inf_sup(array3d)))" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "9d4e616a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[[0. 0. 0. 0. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 0. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. 0. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 1. 1. 1. 0.]\n", + " [0. 0. 0. 0. 0.]]]\n" + ] + } + ], + "source": [ + "array3d = np.zeros((3, 5, 5))\n", + "array3d[:, 1:4, 1:4] = 1\n", + "array3d[1, 2, 2] = 0\n", + "print(array3d)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "90c85d91", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pyclesperanto\n", + "\n", + "[[[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 1 1 1 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]]\n" + ] + } + ], + "source": [ + "print(\"Pyclesperanto\\n\")\n", + "print(cle.superior_inferior(cle.inferior_superior(array3d)))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "226f18fe", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Scikit-image\n", + "\n", + "[[[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 1 1 1 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]]\n", + "\n", + " [[0 0 0 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 1 0 0]\n", + " [0 0 0 0 0]\n", + " [0 0 0 0 0]]]\n" + ] + } + ], + "source": [ + "print(\"Scikit-image\\n\")\n", + "print(sup_inf(inf_sup(array3d)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f0b63e2d", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "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.15" + }, + "vscode": { + "interpreter": { + "hash": "a77a85fb3f9351cc5c809eebe289f69647bd4ccab6f00793c7ad5e93a0cdb95e" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyclesperanto_prototype/_tier1/__init__.py b/pyclesperanto_prototype/_tier1/__init__.py index 7323b6e8..fea0556a 100644 --- a/pyclesperanto_prototype/_tier1/__init__.py +++ b/pyclesperanto_prototype/_tier1/__init__.py @@ -68,6 +68,7 @@ from ._greater_or_equal import greater_or_equal from ._greater_or_equal_constant import greater_or_equal_constant from ._hessian_eigenvalues import hessian_eigenvalues +from ._inferior_superior import inferior_superior from ._laplace_box import laplace_box from ._laplace_diamond import laplace_diamond from ._logarithm import logarithm @@ -161,6 +162,7 @@ from ._sum_x_projection import sum_x_projection from ._sum_y_projection import sum_y_projection from ._sum_z_projection import sum_z_projection +from ._superior_inferior import superior_inferior from ._transpose_xy import transpose_xy from ._transpose_xz import transpose_xz from ._transpose_yz import transpose_yz diff --git a/pyclesperanto_prototype/_tier1/_inferior_superior.py b/pyclesperanto_prototype/_tier1/_inferior_superior.py new file mode 100644 index 00000000..bb9c9418 --- /dev/null +++ b/pyclesperanto_prototype/_tier1/_inferior_superior.py @@ -0,0 +1,40 @@ +from .._tier0 import execute, create_binary_like +from .._tier0 import plugin_function +from .._tier0 import Image + +@plugin_function(categories=['binary processing'], output_creator=create_binary_like) +def inferior_superior(source : Image, destination : Image = None) -> Image: + """Dilates the image respectively with a number of kernels and takes the minimum + value across all results for each pixel. + + Parameters + ---------- + source : Image + destination : Image, optional + + Returns + ------- + destination + + Examples + -------- + >>> import pyclesperanto_prototype as cle + >>> cle.inferior_superior(source, destination) + + References + ---------- + Implemented in inf_sup function in scikit morphological snakes: + .. [1] https://github.com/scikit-image/scikit-image/blob/00177e14097237ef20ed3141ed454bc81b308f82/skimage/segmentation/morphsnakes.py + """ + import numpy as np + if source.dtype != np.uint8: + # the read function in the kernel below is custom and only supports images of type uint8 + source = source.astype(np.uint8) + + parameters = { + "src":source, + "dst":destination + } + + execute(__file__, './inferior_superior_' + str(len(destination.shape)) + 'd_x.cl', 'inferior_superior_' + str(len(destination.shape)) + 'd', destination.shape, parameters) + return destination diff --git a/pyclesperanto_prototype/_tier1/_superior_inferior.py b/pyclesperanto_prototype/_tier1/_superior_inferior.py new file mode 100644 index 00000000..7adba43c --- /dev/null +++ b/pyclesperanto_prototype/_tier1/_superior_inferior.py @@ -0,0 +1,40 @@ +from .._tier0 import execute, create_binary_like +from .._tier0 import plugin_function +from .._tier0 import Image + +@plugin_function(categories=['binary processing'], output_creator=create_binary_like) +def superior_inferior(source : Image, destination : Image = None) -> Image: + """Erodes the image respectively with a number of kernels and takes the maximum + value across all results for each pixel. + + Parameters + ---------- + source : Image + destination : Image, optional + + Returns + ------- + destination + + Examples + -------- + >>> import pyclesperanto_prototype as cle + >>> cle.inferior_superior(source, destination) + + References + ---------- + Implemented in sup_inf function in scikit morphological snakes: + .. [1] https://github.com/scikit-image/scikit-image/blob/00177e14097237ef20ed3141ed454bc81b308f82/skimage/segmentation/morphsnakes.py + """ + import numpy as np + if source.dtype != np.uint8: + # the read function in the kernel below is custom and only supports images of type uint8 + source = source.astype(np.uint8) + + parameters = { + "src":source, + "dst":destination + } + + execute(__file__, './superior_inferior_' + str(len(destination.shape)) + 'd_x.cl', 'superior_inferior_' + str(len(destination.shape)) + 'd', destination.shape, parameters) + return destination diff --git a/pyclesperanto_prototype/_tier1/inferior_superior_2d_x.cl b/pyclesperanto_prototype/_tier1/inferior_superior_2d_x.cl new file mode 100644 index 00000000..8db97f94 --- /dev/null +++ b/pyclesperanto_prototype/_tier1/inferior_superior_2d_x.cl @@ -0,0 +1,83 @@ +__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST; + +#define READ_IMAGE_ZERO_OUTSIDE(a,b,c) read_buffer2duc_zero_outside(GET_IMAGE_WIDTH(a),GET_IMAGE_HEIGHT(a),GET_IMAGE_DEPTH(a),a,b,c) + +inline uchar2 read_buffer2duc_zero_outside(int read_buffer_width, int read_buffer_height, int read_buffer_depth, __global uchar * buffer_var, sampler_t sampler, int2 position ) +{ + int2 pos = (int2){position.x, position.y}; + int pos_in_buffer = pos.x + pos.y * read_buffer_width; + if (pos.x < 0 || pos.x >= read_buffer_width || pos.y < 0 || pos.y >= read_buffer_height) { + return (uchar2){0, 0}; + } + return (uchar2){buffer_var[pos_in_buffer],0}; +} + +__kernel void inferior_superior_2d ( + IMAGE_src_TYPE src, + IMAGE_dst_TYPE dst +) +{ + const int x = get_global_id(0); + const int y = get_global_id(1); + + const int2 pos = (int2){x,y}; + + float value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, pos).x; + + // if value is already 1, dilate will return 1 + if (value == 1) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + /* Dilate with kernel [[1, 0, 0], + [0, 1, 0], + [0, 0, 1]] */ + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, 1})).x; + if (value == 0) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, -1})).x; + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + } + + /* Dilate with kernel [[0, 1, 0], + [0, 1, 0], + [0, 1, 0]] */ + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){0, 1})).x; + if (value == 0) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){0, -1})).x; + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + } + + /* Dilate with kernel [[0, 0, 1], + [0, 1, 0], + [1, 0, 0]] */ + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, 1})).x; + if (value == 0) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, -1})).x; + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + } + + /* Dilate with kernel [[0, 0, 0], + [1, 1, 1], + [0, 0, 0]] */ + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, 0})).x; + if (value == 0) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, 0})).x; + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + } + + // If all dilates are 1 then return 1 + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); +} \ No newline at end of file diff --git a/pyclesperanto_prototype/_tier1/inferior_superior_3d_x.cl b/pyclesperanto_prototype/_tier1/inferior_superior_3d_x.cl new file mode 100644 index 00000000..367df83d --- /dev/null +++ b/pyclesperanto_prototype/_tier1/inferior_superior_3d_x.cl @@ -0,0 +1,194 @@ +__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST; + + +#define READ_IMAGE_ZERO_OUTSIDE(a,b,c) read_buffer3duc_zero_outside(GET_IMAGE_WIDTH(a),GET_IMAGE_HEIGHT(a),GET_IMAGE_DEPTH(a),a,b,c) + +inline uchar2 read_buffer3duc_zero_outside(int read_buffer_width, int read_buffer_height, int read_buffer_depth, __global uchar * buffer_var, sampler_t sampler, int4 position ) +{ + int4 pos = (int4){position.x, position.y, position.z, 0}; + int pos_in_buffer = pos.x + pos.y * read_buffer_width + pos.z * read_buffer_width * read_buffer_height; + if (pos.x < 0 || pos.x >= read_buffer_width || pos.y < 0 || pos.y >= read_buffer_height || pos.z < 0 || pos.z >= read_buffer_depth) { + return (uchar2){0, 0}; + } + return (uchar2){buffer_var[pos_in_buffer],0}; +} + +__kernel void inferior_superior_3d ( + IMAGE_src_TYPE src, + IMAGE_dst_TYPE dst +) +{ + const int x = get_global_id(0); + const int y = get_global_id(1); + const int z = get_global_id(2); + + const int4 pos = (int4){x, y, z, 0}; + + float value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, pos).x; + + // if value is already 0, erode will return 0 + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // printf("pixel at coord x%i y%i z%i\n", x, y, z); + + + // P0 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, 0, 0})).x; + // printf("value %i\n", value); + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // P1 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, 0, j, 0})).x; + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // P2 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){0, i, j, 0})).x; + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // P3 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, j, 0})).x; + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // P4 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){j, i, -i, 0})).x; + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // P5 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, i, 0})).x; + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // P6 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, -i, 0})).x; + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // P7 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, i, j, 0})).x; + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // P8 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, -i, j, 0})).x; + if (value != 0) { + break; + } + } + if (value != 0) { + break; + } + } + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // If all erodes are 0 then return 0 + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); +} \ No newline at end of file diff --git a/pyclesperanto_prototype/_tier1/superior_inferior_2d_x.cl b/pyclesperanto_prototype/_tier1/superior_inferior_2d_x.cl new file mode 100644 index 00000000..aae7c711 --- /dev/null +++ b/pyclesperanto_prototype/_tier1/superior_inferior_2d_x.cl @@ -0,0 +1,85 @@ +__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST; + + +#define READ_IMAGE_ZERO_OUTSIDE(a,b,c) read_buffer2duc_zero_outside(GET_IMAGE_WIDTH(a),GET_IMAGE_HEIGHT(a),GET_IMAGE_DEPTH(a),a,b,c) + +inline uchar2 read_buffer2duc_zero_outside(int read_buffer_width, int read_buffer_height, int read_buffer_depth, __global uchar * buffer_var, sampler_t sampler, int2 position ) +{ + int2 pos = (int2){position.x, position.y}; + int pos_in_buffer = pos.x + pos.y * read_buffer_width; + if (pos.x < 0 || pos.x >= read_buffer_width || pos.y < 0 || pos.y >= read_buffer_height) { + return (uchar2){0, 0}; + } + return (uchar2){buffer_var[pos_in_buffer],0}; +} + + +__kernel void superior_inferior_2d ( + IMAGE_src_TYPE src, + IMAGE_dst_TYPE dst +) +{ + const int x = get_global_id(0); + const int y = get_global_id(1); + + const int2 pos = (int2){x,y}; + + float value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, pos).x; + + // if value is already 0, erode will return 0 + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + /* Erode with kernel [[1, 0, 0], + [0, 1, 0], + [0, 0, 1]] */ + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, 1})).x; + if (value != 0) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, -1})).x; + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + } + + /* Erode with kernel [[0, 1, 0], + [0, 1, 0], + [0, 1, 0]] */ + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){0, 1})).x; + if (value != 0) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){0, -1})).x; + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + } + + /* Erode with kernel [[0, 0, 1], + [0, 1, 0], + [1, 0, 0]] */ + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, 1})).x; + if (value != 0) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, -1})).x; + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + } + + /* Erode with kernel [[0, 0, 0], + [1, 1, 1], + [0, 0, 0]] */ + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){1, 0})).x; + if (value != 0) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int2){-1, 0})).x; + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + } + + // If all erodes are 0 then return 0 + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); +} \ No newline at end of file diff --git a/pyclesperanto_prototype/_tier1/superior_inferior_3d_x.cl b/pyclesperanto_prototype/_tier1/superior_inferior_3d_x.cl new file mode 100644 index 00000000..d997a227 --- /dev/null +++ b/pyclesperanto_prototype/_tier1/superior_inferior_3d_x.cl @@ -0,0 +1,193 @@ +__constant sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST; + +#define READ_IMAGE_ZERO_OUTSIDE(a,b,c) read_buffer3duc_zero_outside(GET_IMAGE_WIDTH(a),GET_IMAGE_HEIGHT(a),GET_IMAGE_DEPTH(a),a,b,c) + +inline uchar2 read_buffer3duc_zero_outside(int read_buffer_width, int read_buffer_height, int read_buffer_depth, __global uchar * buffer_var, sampler_t sampler, int4 position ) +{ + int4 pos = (int4){position.x, position.y, position.z, 0}; + int pos_in_buffer = pos.x + pos.y * read_buffer_width + pos.z * read_buffer_width * read_buffer_height; + if (pos.x < 0 || pos.x >= read_buffer_width || pos.y < 0 || pos.y >= read_buffer_height || pos.z < 0 || pos.z >= read_buffer_depth) { + return (uchar2){0, 0}; + } + return (uchar2){buffer_var[pos_in_buffer],0}; +} + +__kernel void superior_inferior_3d ( + IMAGE_src_TYPE src, + IMAGE_dst_TYPE dst +) +{ + const int x = get_global_id(0); + const int y = get_global_id(1); + const int z = get_global_id(2); + + const int4 pos = (int4){x, y, z, 0}; + + float value = READ_src_IMAGE(src, sampler, pos).x; + + // if value is already 0, erode will return 0 + if (value == 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); + return; + } + + // printf("pixel at coord x%i y%i z%i\n", x, y, z); + + + // P0 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, 0, 0})).x; + // printf("value %i\n", value); + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // P1 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, 0, j, 0})).x; + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // P2 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){0, i, j, 0})).x; + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // P3 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, j, 0})).x; + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // P4 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){j, i, -i, 0})).x; + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // P5 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, i, 0})).x; + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // P6 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, j, -i, 0})).x; + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // P7 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, i, j, 0})).x; + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // P8 + for (int i = -1; i <= 1; i++) { + for (int j = -1; j <= 1; j++) { + value = READ_IMAGE_ZERO_OUTSIDE(src, sampler, (pos + (int4){i, -i, j, 0})).x; + if (value == 0) { + break; + } + } + if (value == 0) { + break; + } + } + if (value != 0) { + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(1)); + return; + } + + // If all erodes are 0 then return 0 + WRITE_dst_IMAGE(dst, pos, CONVERT_dst_PIXEL_TYPE(0)); +} \ No newline at end of file diff --git a/pyclesperanto_prototype/_tier3/__init__.py b/pyclesperanto_prototype/_tier3/__init__.py index 5036b154..2d806aef 100644 --- a/pyclesperanto_prototype/_tier3/__init__.py +++ b/pyclesperanto_prototype/_tier3/__init__.py @@ -56,6 +56,7 @@ from ._standard_deviation_of_touching_neighbors_map import standard_deviation_of_touching_neighbors_map from ._standard_deviation_of_proximal_neighbors_map import standard_deviation_of_proximal_neighbors_map from ._standard_deviation_of_proximal_neighbors_map import standard_deviation_of_proximal_neighbors_map as standard_deviation_of_distal_neighbors_map +from ._morphological_chan_vese import morphological_chan_vese from ._subtract_gaussian_background import subtract_gaussian_background from ._z_position_range_projection import z_position_range_projection \ No newline at end of file diff --git a/pyclesperanto_prototype/_tier3/_morphological_chan_vese.py b/pyclesperanto_prototype/_tier3/_morphological_chan_vese.py new file mode 100644 index 00000000..cfdf5b11 --- /dev/null +++ b/pyclesperanto_prototype/_tier3/_morphological_chan_vese.py @@ -0,0 +1,170 @@ +import numpy as np +from .._tier0 import Image, plugin_function + +@plugin_function +def morphological_chan_vese(input_image: Image, + contour_image : Image = np.zeros((0, 0)), + output_image : Image = None, + num_iter : int = 100, + smoothing : int = 1, + lambda1 : float = 1, + lambda2 : float = 1) -> Image: + """Morphological Active Contours without Edges (MorphACWE) + + Active contours without edges implemented with morphological operators. It + can be used to segment objects in images and volumes without well defined + borders. It is required that the inside of the object looks different on + average than the outside (i.e., the inner area of the object should be + darker or lighter than the outer area on average). + + This is a re-implementation of the algorithm in scikit-image [2]. + + Parameters + ---------- + input_image: Image + contour_image: Image, optional + output_image: Image, optional + num_iter : uint + Number of iterations to run + smoothing: int, optional + lambda1 : float, optional + Weight parameter for the outer region. If `lambda1` is larger than + `lambda2`, the outer region will contain a larger range of values than + the inner region. + lambda2 : float, optional + Weight parameter for the inner region. If `lambda2` is larger than + `lambda1`, the inner region will contain a larger range of values than + the outer region. + + Returns + ------- + Final segmentation + + Notes + ----- + This is a version of the Chan-Vese algorithm that uses morphological + operators instead of solving a partial differential equation (PDE) for the + evolution of the contour. The set of morphological operators used in this + algorithm are proved to be infinitesimally equivalent to the Chan-Vese PDE + (see [1]_). However, morphological operators are do not suffer from the + numerical stability issues typically found in PDEs (it is not necessary to + find the right time step for the evolution), and are computationally + faster. + + The algorithm and its theoretical derivation are described in [1]_. + + References + ---------- + .. [1] A Morphological Approach to Curvature-based Evolution of Curves and + Surfaces, Pablo Márquez-Neila, Luis Baumela, Luis Álvarez. In IEEE + Transactions on Pattern Analysis and Machine Intelligence (PAMI), + 2014, :DOI:`10.1109/TPAMI.2013.106` + .. [2] https://github.com/scikit-image/scikit-image/blob/5e74a4a3a5149a8a14566b81a32bb15499aa3857/skimage/segmentation/morphsnakes.py#L212-L312 + """ + from .._tier0 import create_like + from .._tier1 import superior_inferior, inferior_superior, set + + from .._tier1 import absolute, binary_or, binary_not, gradient_x, gradient_y + from .._tier1 import gradient_z, copy, power + from .._tier1 import greater_constant, smaller_constant + from .._tier1 import mask, add_image_and_scalar, add_images_weighted + + if contour_image.size == 0: + contour_image = checkerboard_level_set(input_image.shape) + + greater_constant(contour_image, destination=output_image, constant=0) + + temp_1 = create_like(output_image) + temp_2 = create_like(output_image) + temp_3 = create_like(output_image) + temp_4 = create_like(output_image) + temp_5 = create_like(output_image) + + for _ in range(num_iter): + set(temp_3, 0) + + # define invert image + binary_not(output_image, destination=temp_1) + + # compute outside contour score + sum_image_value = mask(input_image, mask=temp_1).sum() + sum_contour_value = temp_1.sum() + 1e-8 + c0 = - (sum_image_value / sum_contour_value) + + # compute inside contour score + sum_image_value = mask(input_image, mask=output_image).sum() + sum_contour_value = output_image.sum() + 1e-8 + c1 = - (sum_image_value / sum_contour_value) + + # Image attachment + # compute gradient on contour in all direction + for d in range(input_image.ndim): + if d == 0: + gradient_x(output_image, destination=temp_1) + absolute(temp_1, destination=temp_2) + temp_3 += temp_2 + if d == 1: + gradient_y(output_image, destination=temp_1) + absolute(temp_1, destination=temp_2) + temp_3 += temp_2 + if d == 2: + gradient_z(output_image, destination=temp_1) + absolute(temp_1, destination=temp_2) + temp_3 += temp_2 + + # compute contour evolution according to gradient and score on contour + add_image_and_scalar(input_image, destination=temp_1, scalar=c1) + add_image_and_scalar(input_image, destination=temp_2, scalar=c0) + power(temp_1, destination=temp_4, exponent=2) + power(temp_2, destination=temp_5, exponent=2) + temp_2 = temp_3 * (lambda1 * temp_4 - lambda2 * temp_5) + + # apply contour update on contour image + greater_constant(temp_2, destination=temp_1, constant=0) + smaller_constant(temp_2, destination=temp_3, constant=0) + + binary_or(temp_1, temp_3, destination=temp_4) + binary_not(temp_4, destination=temp_1) + mask(output_image, mask=temp_1, destination=temp_2) + add_images_weighted(temp_2, temp_3, destination=temp_1, factor1=1, factor2=1) + + # smooth contour + for s in range(smoothing): + if s % 2 == 0: + inferior_superior(temp_1, temp_2) + superior_inferior(temp_2, temp_1) + else: + superior_inferior(temp_1, temp_2) + inferior_superior(temp_2, temp_1) + + copy(temp_1, destination=output_image) + + return output_image + + +def checkerboard_level_set(image_shape, square_size=5): + """Create a checkerboard level set with binary values. + Parameters + ---------- + image_shape : tuple of positive integers + Shape of the image. + square_size : int, optional + Size of the squares of the checkerboard. It defaults to 5. + Returns + ------- + out : array with shape `image_shape` + Binary level set of the checkerboard. + See Also + -------- + disk_level_set + """ + + grid = np.mgrid[[slice(i) for i in image_shape]] + grid = (grid // square_size) + + # Alternate 0/1 for even/odd numbers. + grid = grid & 1 + + checkerboard = np.bitwise_xor.reduce(grid, axis=0) + res = np.int8(checkerboard) + return res \ No newline at end of file diff --git a/tests/test_inferior_superior.py b/tests/test_inferior_superior.py new file mode 100644 index 00000000..5d89c3f0 --- /dev/null +++ b/tests/test_inferior_superior.py @@ -0,0 +1,90 @@ +import pyclesperanto_prototype as cle +import numpy as np + + +def test_inferior_superior_2d(): + test = cle.push(np.asarray([ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + ])) + + reference = cle.push(np.asarray([ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 0, 0, 1, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 1, 0, 0, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + ])) + + result = cle.create(test) + cle.inferior_superior(test, result) + + print(result) + + a = cle.pull(result) + b = cle.pull(reference) + assert (np.array_equal(a, b)) + + +def test_inferior_superior_3d(): + test = cle.push(np.asarray([[ + [0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], + + [[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 0, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], + + [[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0] + ]])) + + reference = cle.push(np.asarray([[ + [0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], + + [[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], + + [[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0] + ]])) + + result = cle.create(test) + cle.inferior_superior(test, result) + + print(result) + + a = cle.pull(result) + b = cle.pull(reference) + assert (np.array_equal(a, b)) \ No newline at end of file diff --git a/tests/test_morphological_snakes.py b/tests/test_morphological_snakes.py new file mode 100644 index 00000000..0930b992 --- /dev/null +++ b/tests/test_morphological_snakes.py @@ -0,0 +1,52 @@ +import pyclesperanto_prototype as cle +from skimage import morphology +from skimage import segmentation +import numpy as np + +def test_morphological_chan_vese(): + + image = generate_disk((50, 50), 5) + + reference = cle.push(segmentation.morphological_chan_vese(image, num_iter=1, smoothing=1, lambda1=1, lambda2=1)) + + result = cle.morphological_chan_vese(image, num_iter=1, smoothing=1, lambda1=1, lambda2=1) + + a = cle.pull(result) + b = cle.pull(reference) + + print(a) + print(b) + + assert (np.array_equal(a, b)) + + +def test_morphological_chan_vese_on_membranes_2d(): + from skimage.data import cells3d + from skimage.segmentation import morphological_chan_vese + import pyclesperanto_prototype as cle + + image2d = cells3d()[30, 0, ...] + + result_gpu = cle.morphological_chan_vese(image2d, num_iter=20, smoothing=10) + result_cpu = morphological_chan_vese(image2d, num_iter=20, smoothing=10) + + assert cle.array_equal(result_cpu, result_gpu) + +def test_morphological_chan_vese_on_membranes_3d(): + + from skimage.data import cells3d + from skimage.segmentation import morphological_chan_vese + import pyclesperanto_prototype as cle + + image3d = cells3d()[10:50:, 0, :40, :40] + + result_gpu = cle.morphological_chan_vese(image3d, num_iter=20, smoothing=10) + result_cpu = morphological_chan_vese(image3d, num_iter=20, smoothing=10) + + cle.array_equal(result_cpu, result_gpu) + +def generate_disk(shape, radius): + image = np.zeros(shape) + image[image.shape[0] // 2 - radius:image.shape[0] // 2 + radius + 1, + image.shape[1] // 2 - radius:image.shape[1] // 2 + radius + 1] = morphology.disk(radius) + return image diff --git a/tests/test_superior_inferior.py b/tests/test_superior_inferior.py new file mode 100644 index 00000000..28e69721 --- /dev/null +++ b/tests/test_superior_inferior.py @@ -0,0 +1,158 @@ +import pyclesperanto_prototype as cle +import numpy as np +from skimage.segmentation.morphsnakes import sup_inf, inf_sup + +def test_superior_inferior_2d(): + test = cle.push(np.asarray([ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + ])) + + reference = cle.push(np.asarray([ + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + ])) + + result = cle.create(test) + cle.superior_inferior(test, result) + + print(result) + + a = cle.pull(result) + b = cle.pull(reference) + assert (np.array_equal(a, b)) + + reference2 = sup_inf(test) + assert cle.array_equal(result, reference2) + +def test_superior_inferior_2d_compare_with_skimage_x(): + + array = np.zeros((100, 85), dtype=np.uint8) + array[15:85, 15:85] = 1 + array[35:65, 35:65] = 0 + + result = cle.superior_inferior(cle.inferior_superior(array)) + reference = sup_inf(inf_sup(array)) + + print("result", result) + print("reference", reference) + + assert cle.array_equal(result, reference) + + +def test_superior_inferior_2d_compare_with_skimage_y(): + + array = np.zeros((85, 100), dtype=np.uint8) + array[15:85, 15:85] = 1 + array[35:65, 35:65] = 0 + + result = cle.superior_inferior(cle.inferior_superior(array)) + reference = sup_inf(inf_sup(array)) + + print("result", result) + print("reference", reference) + + assert cle.array_equal(result, reference) + + +def test_superior_inferior_3d(): + test = cle.push(np.asarray([ + [[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]], + + [[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], + + [[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 0, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], + + [[0, 0, 0, 0, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 1, 1, 1, 0], + [0, 0, 0, 0, 0]], + + [[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]], + ])) + + reference = inf_sup(test) + + result = cle.inferior_superior(test) + + print("result", result) + print("reference", reference) + + a = cle.pull(result) + b = cle.pull(reference) + assert (np.array_equal(a, b)) + +def test_superior_inferior_3d_compare_with_skimage_x(): + + array = np.zeros((5, 5, 4)) + array[1:4, 1:4, 1:4] = 1 + array[2, 2, 2] = 0 + + result = cle.superior_inferior(cle.inferior_superior(array)) + reference = sup_inf(inf_sup(array)) + + print("result", result) + print("reference", reference) + + assert cle.array_equal(result, reference) + +def test_superior_inferior_3d_compare_with_skimage_y(): + + array = np.zeros((5, 4, 5)) + array[1:4, 1:4, 1:4] = 1 + array[2, 2, 2] = 0 + + result = cle.superior_inferior(cle.inferior_superior(array)) + reference = sup_inf(inf_sup(array)) + + print("result", result) + print("reference", reference) + + assert cle.array_equal(result, reference) + +def test_superior_inferior_3d_compare_with_skimage_z(): + + array = np.zeros((4, 5, 5)) + array[1:4, 1:4, 1:4] = 1 + array[2, 2, 2] = 0 + + result = cle.superior_inferior(cle.inferior_superior(array)) + reference = sup_inf(inf_sup(array)) + + print("result", result) + print("reference", reference) + + assert cle.array_equal(result, reference)