{ "cells": [ { "cell_type": "code", "execution_count": 0, "metadata": { "colab": { "autoexec": { "startup": false, "wait_interval": 0 } }, "colab_type": "code", "id": "mmKBfDbTa_4M" }, "outputs": [], "source": [ "import numpy as np\n", "import scipy.ndimage " ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": { "autoexec": { "startup": false, "wait_interval": 0 } }, "colab_type": "code", "id": "-r22-U0ZnrRz" }, "outputs": [], "source": [ "# neighbour_code_to_normals is a lookup table.\n", "# For every binary neighbour code \n", "# (2x2x2 neighbourhood = 8 neighbours = 8 bits = 256 codes) \n", "# it contains the surface normals of the triangles (called \"surfel\" for \n", "# \"surface element\" in the following). The length of the normal \n", "# vector encodes the surfel area.\n", "#\n", "# created by compute_surface_area_lookup_table.ipynb using the \n", "# marching_cube algorithm, see e.g. https://en.wikipedia.org/wiki/Marching_cubes\n", "#\n", "neighbour_code_to_normals = [\n", " [[0,0,0]],\n", " [[0.125,0.125,0.125]],\n", " [[-0.125,-0.125,0.125]],\n", " [[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n", " [[0.125,-0.125,0.125]],\n", " [[-0.25,-0.0,-0.25],[0.25,0.0,0.25]],\n", " [[0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[0.5,0.0,-0.0],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n", " [[-0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,0.125,0.125]],\n", " [[-0.25,0.0,0.25],[-0.25,0.0,0.25]],\n", " [[0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n", " [[0.5,0.0,0.0],[0.25,-0.25,0.25],[-0.125,0.125,-0.125]],\n", " [[-0.5,0.0,0.0],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n", " [[0.5,0.0,0.0],[0.5,0.0,0.0]],\n", " [[0.125,-0.125,-0.125]],\n", " [[0.0,-0.25,-0.25],[0.0,0.25,0.25]],\n", " [[-0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.0,-0.5,0.0],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n", " [[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.0,0.0,-0.5],[0.25,0.25,0.25],[-0.125,-0.125,-0.125]],\n", " [[-0.125,-0.125,0.125],[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n", " [[-0.125,0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[-0.125,0.125,0.125]],\n", " [[-0.25,0.0,0.25],[-0.25,0.0,0.25],[0.125,-0.125,-0.125]],\n", " [[0.125,0.125,0.125],[0.375,0.375,0.375],[0.0,-0.25,0.25],[-0.25,0.0,0.25]],\n", " [[0.125,-0.125,-0.125],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n", " [[0.375,0.375,0.375],[0.0,0.25,-0.25],[-0.125,-0.125,-0.125],[-0.25,0.25,0.0]],\n", " [[-0.5,0.0,0.0],[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25],[0.125,0.125,0.125]],\n", " [[-0.5,0.0,0.0],[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25]],\n", " [[0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n", " [[0.0,-0.5,0.0],[0.125,0.125,-0.125],[0.25,0.25,-0.25]],\n", " [[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.125,-0.125,0.125],[-0.25,-0.0,-0.25],[0.25,0.0,0.25]],\n", " [[0.0,-0.25,0.25],[0.0,0.25,-0.25],[0.125,-0.125,0.125]],\n", " [[-0.375,-0.375,0.375],[-0.0,0.25,0.25],[0.125,0.125,-0.125],[-0.25,-0.0,-0.25]],\n", " [[-0.125,0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n", " [[-0.0,0.0,0.5],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[0.25,0.25,-0.25],[0.25,0.25,-0.25],[0.125,0.125,-0.125],[-0.125,-0.125,0.125]],\n", " [[0.125,-0.125,0.125],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n", " [[0.5,0.0,0.0],[0.25,-0.25,0.25],[-0.125,0.125,-0.125],[0.125,-0.125,0.125]],\n", " [[0.0,0.25,-0.25],[0.375,-0.375,-0.375],[-0.125,0.125,0.125],[0.25,0.25,0.0]],\n", " [[-0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[0.25,-0.25,0.0],[-0.25,0.25,0.0]],\n", " [[0.0,0.5,0.0],[-0.25,0.25,0.25],[0.125,-0.125,-0.125]],\n", " [[0.0,0.5,0.0],[0.125,-0.125,0.125],[-0.25,0.25,-0.25]],\n", " [[0.0,0.5,0.0],[0.0,-0.5,0.0]],\n", " [[0.25,-0.25,0.0],[-0.25,0.25,0.0],[0.125,-0.125,0.125]],\n", " [[-0.375,-0.375,-0.375],[-0.25,0.0,0.25],[-0.125,-0.125,-0.125],[-0.25,0.25,0.0]],\n", " [[0.125,0.125,0.125],[0.0,-0.5,0.0],[-0.25,-0.25,-0.25],[-0.125,-0.125,-0.125]],\n", " [[0.0,-0.5,0.0],[-0.25,-0.25,-0.25],[-0.125,-0.125,-0.125]],\n", " [[-0.125,0.125,0.125],[0.25,-0.25,0.0],[-0.25,0.25,0.0]],\n", " [[0.0,0.5,0.0],[0.25,0.25,-0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.375,0.375,-0.375],[-0.25,-0.25,0.0],[-0.125,0.125,-0.125],[-0.25,0.0,0.25]],\n", " [[0.0,0.5,0.0],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n", " [[0.25,-0.25,0.0],[-0.25,0.25,0.0],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n", " [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0],[-0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.25,-0.25,0.0],[-0.25,-0.25,0.0]],\n", " [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0]],\n", " [[-0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n", " [[0.0,-0.25,0.25],[0.0,-0.25,0.25]],\n", " [[0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125]],\n", " [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[0.375,-0.375,0.375],[0.0,-0.25,-0.25],[-0.125,0.125,-0.125],[0.25,0.25,0.0]],\n", " [[-0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[-0.25,0.0,0.25],[-0.25,0.0,0.25]],\n", " [[0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.0,0.5,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125]],\n", " [[-0.25,0.25,-0.25],[-0.25,0.25,-0.25],[-0.125,0.125,-0.125],[-0.125,0.125,-0.125]],\n", " [[-0.25,0.0,-0.25],[0.375,-0.375,-0.375],[0.0,0.25,-0.25],[-0.125,0.125,0.125]],\n", " [[0.5,0.0,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125]],\n", " [[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n", " [[-0.0,0.0,0.5],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n", " [[-0.25,-0.0,-0.25],[-0.375,0.375,0.375],[-0.25,-0.25,0.0],[-0.125,0.125,0.125]],\n", " [[0.0,0.0,-0.5],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n", " [[-0.0,0.0,0.5],[0.0,0.0,0.5]],\n", " [[0.125,0.125,0.125],[0.125,0.125,0.125],[0.25,0.25,0.25],[0.0,0.0,0.5]],\n", " [[0.125,0.125,0.125],[0.25,0.25,0.25],[0.0,0.0,0.5]],\n", " [[-0.25,0.0,0.25],[0.25,0.0,-0.25],[-0.125,0.125,0.125]],\n", " [[-0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[-0.25,0.0,0.25],[-0.25,0.0,0.25],[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n", " [[0.125,-0.125,0.125],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n", " [[0.25,0.0,0.25],[-0.375,-0.375,0.375],[-0.25,0.25,0.0],[-0.125,-0.125,0.125]],\n", " [[-0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n", " [[0.25,0.0,0.25],[0.25,0.0,0.25]],\n", " [[-0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n", " [[0.0,-0.5,0.0],[0.125,0.125,-0.125],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n", " [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[0.125,-0.125,0.125]],\n", " [[0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n", " [[0.0,0.25,0.25],[0.0,0.25,0.25],[0.125,-0.125,-0.125]],\n", " [[-0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,-0.125,0.125],[0.125,0.125,0.125]],\n", " [[-0.0,0.0,0.5],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[-0.0,0.5,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.5,0.0,-0.0],[0.25,-0.25,-0.25],[0.125,-0.125,-0.125]],\n", " [[-0.25,0.25,0.25],[-0.125,0.125,0.125],[-0.25,0.25,0.25],[0.125,-0.125,-0.125]],\n", " [[0.375,-0.375,0.375],[0.0,0.25,0.25],[-0.125,0.125,-0.125],[-0.25,0.0,0.25]],\n", " [[0.0,-0.5,0.0],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n", " [[-0.375,-0.375,0.375],[0.25,-0.25,0.0],[0.0,0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[-0.125,0.125,0.125],[-0.25,0.25,0.25],[0.0,0.0,0.5]],\n", " [[0.125,0.125,0.125],[0.0,0.25,0.25],[0.0,0.25,0.25]],\n", " [[0.0,0.25,0.25],[0.0,0.25,0.25]],\n", " [[0.5,0.0,-0.0],[0.25,0.25,0.25],[0.125,0.125,0.125],[0.125,0.125,0.125]],\n", " [[0.125,-0.125,0.125],[-0.125,-0.125,0.125],[0.125,0.125,0.125]],\n", " [[-0.25,-0.0,-0.25],[0.25,0.0,0.25],[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,0.125]],\n", " [[-0.25,-0.25,0.0],[0.25,0.25,-0.0],[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.25,-0.25,0.0],[0.25,0.25,-0.0],[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,0.125]],\n", " [[-0.25,-0.0,-0.25],[0.25,0.0,0.25],[0.125,0.125,0.125]],\n", " [[0.125,-0.125,0.125],[-0.125,-0.125,0.125],[0.125,0.125,0.125]],\n", " [[0.5,0.0,-0.0],[0.25,0.25,0.25],[0.125,0.125,0.125],[0.125,0.125,0.125]],\n", " [[0.0,0.25,0.25],[0.0,0.25,0.25]],\n", " [[0.125,0.125,0.125],[0.0,0.25,0.25],[0.0,0.25,0.25]],\n", " [[-0.125,0.125,0.125],[-0.25,0.25,0.25],[0.0,0.0,0.5]],\n", " [[-0.375,-0.375,0.375],[0.25,-0.25,0.0],[0.0,0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[0.0,-0.5,0.0],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n", " [[0.375,-0.375,0.375],[0.0,0.25,0.25],[-0.125,0.125,-0.125],[-0.25,0.0,0.25]],\n", " [[-0.25,0.25,0.25],[-0.125,0.125,0.125],[-0.25,0.25,0.25],[0.125,-0.125,-0.125]],\n", " [[0.5,0.0,-0.0],[0.25,-0.25,-0.25],[0.125,-0.125,-0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[-0.0,0.5,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[-0.0,0.0,0.5],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,-0.125,0.125],[0.125,0.125,0.125]],\n", " [[-0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[0.0,0.25,0.25],[0.0,0.25,0.25],[0.125,-0.125,-0.125]],\n", " [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[0.0,0.25,0.25],[0.0,0.25,0.25]],\n", " [[0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[0.125,-0.125,0.125]],\n", " [[0.0,-0.5,0.0],[0.125,0.125,-0.125],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.25,0.0,0.25],[0.25,0.0,0.25]],\n", " [[0.125,0.125,0.125],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n", " [[-0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125]],\n", " [[0.25,0.0,0.25],[-0.375,-0.375,0.375],[-0.25,0.25,0.0],[-0.125,-0.125,0.125]],\n", " [[0.125,-0.125,0.125],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n", " [[-0.25,-0.0,-0.25],[0.25,0.0,0.25],[0.25,0.0,0.25],[0.25,0.0,0.25]],\n", " [[-0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[-0.25,0.0,0.25],[0.25,0.0,-0.25],[-0.125,0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.25,0.25,0.25],[0.0,0.0,0.5]],\n", " [[0.125,0.125,0.125],[0.125,0.125,0.125],[0.25,0.25,0.25],[0.0,0.0,0.5]],\n", " [[-0.0,0.0,0.5],[0.0,0.0,0.5]],\n", " [[0.0,0.0,-0.5],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n", " [[-0.25,-0.0,-0.25],[-0.375,0.375,0.375],[-0.25,-0.25,0.0],[-0.125,0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n", " [[-0.0,0.0,0.5],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n", " [[-0.25,0.0,0.25],[0.25,0.0,-0.25]],\n", " [[0.5,0.0,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125]],\n", " [[-0.25,0.0,-0.25],[0.375,-0.375,-0.375],[0.0,0.25,-0.25],[-0.125,0.125,0.125]],\n", " [[-0.25,0.25,-0.25],[-0.25,0.25,-0.25],[-0.125,0.125,-0.125],[-0.125,0.125,-0.125]],\n", " [[-0.0,0.5,0.0],[-0.25,0.25,-0.25],[0.125,-0.125,0.125]],\n", " [[0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[-0.25,0.0,0.25],[-0.25,0.0,0.25]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n", " [[0.375,-0.375,0.375],[0.0,-0.25,-0.25],[-0.125,0.125,-0.125],[0.25,0.25,0.0]],\n", " [[0.0,-0.25,0.25],[0.0,-0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[0.0,0.0,0.5],[0.25,-0.25,0.25],[0.125,-0.125,0.125]],\n", " [[0.0,-0.25,0.25],[0.0,-0.25,0.25]],\n", " [[-0.125,-0.125,0.125],[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n", " [[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.125,-0.125,0.125]],\n", " [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0]],\n", " [[0.125,0.125,0.125],[-0.25,-0.25,0.0],[-0.25,-0.25,0.0]],\n", " [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0],[-0.125,-0.125,0.125]],\n", " [[-0.25,-0.25,0.0],[-0.25,-0.25,0.0],[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n", " [[0.0,0.5,0.0],[0.25,0.25,-0.25],[-0.125,-0.125,0.125]],\n", " [[-0.375,0.375,-0.375],[-0.25,-0.25,0.0],[-0.125,0.125,-0.125],[-0.25,0.0,0.25]],\n", " [[0.0,0.5,0.0],[0.25,0.25,-0.25],[-0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.125,0.125,0.125],[0.25,-0.25,0.0],[-0.25,0.25,0.0]],\n", " [[0.0,-0.5,0.0],[-0.25,-0.25,-0.25],[-0.125,-0.125,-0.125]],\n", " [[0.125,0.125,0.125],[0.0,-0.5,0.0],[-0.25,-0.25,-0.25],[-0.125,-0.125,-0.125]],\n", " [[-0.375,-0.375,-0.375],[-0.25,0.0,0.25],[-0.125,-0.125,-0.125],[-0.25,0.25,0.0]],\n", " [[0.25,-0.25,0.0],[-0.25,0.25,0.0],[0.125,-0.125,0.125]],\n", " [[0.0,0.5,0.0],[0.0,-0.5,0.0]],\n", " [[0.0,0.5,0.0],[0.125,-0.125,0.125],[-0.25,0.25,-0.25]],\n", " [[0.0,0.5,0.0],[-0.25,0.25,0.25],[0.125,-0.125,-0.125]],\n", " [[0.25,-0.25,0.0],[-0.25,0.25,0.0]],\n", " [[-0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[0.0,0.25,-0.25],[0.375,-0.375,-0.375],[-0.125,0.125,0.125],[0.25,0.25,0.0]],\n", " [[0.5,0.0,0.0],[0.25,-0.25,0.25],[-0.125,0.125,-0.125],[0.125,-0.125,0.125]],\n", " [[0.125,-0.125,0.125],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n", " [[0.25,0.25,-0.25],[0.25,0.25,-0.25],[0.125,0.125,-0.125],[-0.125,-0.125,0.125]],\n", " [[-0.0,0.0,0.5],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,0.125],[-0.125,0.125,0.125]],\n", " [[-0.125,0.125,0.125],[0.125,-0.125,0.125]],\n", " [[-0.375,-0.375,0.375],[-0.0,0.25,0.25],[0.125,0.125,-0.125],[-0.25,-0.0,-0.25]],\n", " [[0.0,-0.25,0.25],[0.0,0.25,-0.25],[0.125,-0.125,0.125]],\n", " [[0.125,-0.125,0.125],[-0.25,-0.0,-0.25],[0.25,0.0,0.25]],\n", " [[0.125,-0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.0,-0.5,0.0],[0.125,0.125,-0.125],[0.25,0.25,-0.25]],\n", " [[0.0,-0.25,0.25],[0.0,0.25,-0.25]],\n", " [[0.125,0.125,0.125],[0.125,-0.125,0.125]],\n", " [[0.125,-0.125,0.125]],\n", " [[-0.5,0.0,0.0],[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25]],\n", " [[-0.5,0.0,0.0],[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25],[0.125,0.125,0.125]],\n", " [[0.375,0.375,0.375],[0.0,0.25,-0.25],[-0.125,-0.125,-0.125],[-0.25,0.25,0.0]],\n", " [[0.125,-0.125,-0.125],[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n", " [[0.125,0.125,0.125],[0.375,0.375,0.375],[0.0,-0.25,0.25],[-0.25,0.0,0.25]],\n", " [[-0.25,0.0,0.25],[-0.25,0.0,0.25],[0.125,-0.125,-0.125]],\n", " [[0.0,-0.25,-0.25],[0.0,0.25,0.25],[-0.125,0.125,0.125]],\n", " [[-0.125,0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[-0.125,-0.125,-0.125],[-0.25,-0.25,-0.25],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.0,0.0,-0.5],[0.25,0.25,0.25],[-0.125,-0.125,-0.125]],\n", " [[0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.0,-0.5,0.0],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n", " [[-0.125,-0.125,0.125],[0.125,-0.125,-0.125]],\n", " [[0.0,-0.25,-0.25],[0.0,0.25,0.25]],\n", " [[0.125,-0.125,-0.125]],\n", " [[0.5,0.0,0.0],[0.5,0.0,0.0]],\n", " [[-0.5,0.0,0.0],[-0.25,0.25,0.25],[-0.125,0.125,0.125]],\n", " [[0.5,0.0,0.0],[0.25,-0.25,0.25],[-0.125,0.125,-0.125]],\n", " [[0.25,-0.25,0.0],[0.25,-0.25,0.0]],\n", " [[0.5,0.0,0.0],[-0.25,-0.25,0.25],[-0.125,-0.125,0.125]],\n", " [[-0.25,0.0,0.25],[-0.25,0.0,0.25]],\n", " [[0.125,0.125,0.125],[-0.125,0.125,0.125]],\n", " [[-0.125,0.125,0.125]],\n", " [[0.5,0.0,-0.0],[0.25,0.25,0.25],[0.125,0.125,0.125]],\n", " [[0.125,-0.125,0.125],[-0.125,-0.125,0.125]],\n", " [[-0.25,-0.0,-0.25],[0.25,0.0,0.25]],\n", " [[0.125,-0.125,0.125]],\n", " [[-0.25,-0.25,0.0],[0.25,0.25,-0.0]],\n", " [[-0.125,-0.125,0.125]],\n", " [[0.125,0.125,0.125]],\n", " [[0,0,0]]]" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": { "autoexec": { "startup": false, "wait_interval": 0 } }, "colab_type": "code", "id": "-kruspp5Y1Ip" }, "outputs": [], "source": [ "def compute_surface_distances(mask_gt, mask_pred, spacing_mm):\n", " \"\"\"Compute closest distances from all surface points to the other surface.\n", "\n", " Finds all surface elements \"surfels\" in the ground truth mask `mask_gt` and\n", " the predicted mask `mask_pred`, computes their area in mm^2 and the distance\n", " to the closest point on the other surface. It returns two sorted lists of\n", " distances together with the corresponding surfel areas. If one of the masks\n", " is empty, the corresponding lists are empty and all distances in the other\n", " list are `inf` \n", " \n", " Args:\n", " mask_gt: 3-dim Numpy array of type bool. The ground truth mask.\n", " mask_pred: 3-dim Numpy array of type bool. The predicted mask.\n", " spacing_mm: 3-element list-like structure. Voxel spacing in x0, x1 and x2\n", " direction \n", "\n", " Returns:\n", " A dict with \n", " \"distances_gt_to_pred\": 1-dim numpy array of type float. The distances in mm\n", " from all ground truth surface elements to the predicted surface, \n", " sorted from smallest to largest\n", " \"distances_pred_to_gt\": 1-dim numpy array of type float. The distances in mm\n", " from all predicted surface elements to the ground truth surface, \n", " sorted from smallest to largest \n", " \"surfel_areas_gt\": 1-dim numpy array of type float. The area in mm^2 of \n", " the ground truth surface elements in the same order as \n", " distances_gt_to_pred\n", " \"surfel_areas_pred\": 1-dim numpy array of type float. The area in mm^2 of \n", " the predicted surface elements in the same order as \n", " distances_pred_to_gt\n", " \n", " \"\"\"\n", " \n", " # compute the area for all 256 possible surface elements \n", " # (given a 2x2x2 neighbourhood) according to the spacing_mm\n", " neighbour_code_to_surface_area = np.zeros([256])\n", " for code in range(256):\n", " normals = np.array(neighbour_code_to_normals[code])\n", " sum_area = 0\n", " for normal_idx in range(normals.shape[0]):\n", " # normal vector\n", " n = np.zeros([3])\n", " n[0] = normals[normal_idx,0] * spacing_mm[1] * spacing_mm[2]\n", " n[1] = normals[normal_idx,1] * spacing_mm[0] * spacing_mm[2]\n", " n[2] = normals[normal_idx,2] * spacing_mm[0] * spacing_mm[1]\n", " area = np.linalg.norm(n)\n", " sum_area += area\n", " neighbour_code_to_surface_area[code] = sum_area\n", "\n", " # compute the bounding box of the masks to trim\n", " # the volume to the smallest possible processing subvolume\n", " mask_all = mask_gt | mask_pred\n", " bbox_min = np.zeros(3, np.int64)\n", " bbox_max = np.zeros(3, np.int64)\n", "\n", " # max projection to the x0-axis\n", " proj_0 = np.max(np.max(mask_all, axis=2), axis=1)\n", " idx_nonzero_0 = np.nonzero(proj_0)[0]\n", " if len(idx_nonzero_0) == 0:\n", " return {\"distances_gt_to_pred\": np.array([]), \n", " \"distances_pred_to_gt\": np.array([]), \n", " \"surfel_areas_gt\": np.array([]), \n", " \"surfel_areas_pred\": np.array([])}\n", " \n", " bbox_min[0] = np.min(idx_nonzero_0)\n", " bbox_max[0] = np.max(idx_nonzero_0)\n", "\n", " # max projection to the x1-axis\n", " proj_1 = np.max(np.max(mask_all, axis=2), axis=0)\n", " idx_nonzero_1 = np.nonzero(proj_1)[0]\n", " bbox_min[1] = np.min(idx_nonzero_1)\n", " bbox_max[1] = np.max(idx_nonzero_1)\n", "\n", " # max projection to the x2-axis\n", " proj_2 = np.max(np.max(mask_all, axis=1), axis=0)\n", " idx_nonzero_2 = np.nonzero(proj_2)[0]\n", " bbox_min[2] = np.min(idx_nonzero_2)\n", " bbox_max[2] = np.max(idx_nonzero_2)\n", "\n", " print(\"bounding box min = {}\".format(bbox_min))\n", " print(\"bounding box max = {}\".format(bbox_max))\n", "\n", " # crop the processing subvolume.\n", " # we need to zeropad the cropped region with 1 voxel at the lower, \n", " # the right and the back side. This is required to obtain the \"full\" \n", " # convolution result with the 2x2x2 kernel\n", " cropmask_gt = np.zeros((bbox_max - bbox_min)+2, np.uint8)\n", " cropmask_pred = np.zeros((bbox_max - bbox_min)+2, np.uint8)\n", "\n", " cropmask_gt[0:-1, 0:-1, 0:-1] = mask_gt[bbox_min[0]:bbox_max[0]+1,\n", " bbox_min[1]:bbox_max[1]+1,\n", " bbox_min[2]:bbox_max[2]+1]\n", "\n", " cropmask_pred[0:-1, 0:-1, 0:-1] = mask_pred[bbox_min[0]:bbox_max[0]+1,\n", " bbox_min[1]:bbox_max[1]+1,\n", " bbox_min[2]:bbox_max[2]+1]\n", "\n", " # compute the neighbour code (local binary pattern) for each voxel\n", " # the resultsing arrays are spacially shifted by minus half a voxel in each axis.\n", " # i.e. the points are located at the corners of the original voxels\n", " kernel = np.array([[[128,64],\n", " [32,16]],\n", " [[8,4],\n", " [2,1]]])\n", " neighbour_code_map_gt = scipy.ndimage.filters.correlate(cropmask_gt.astype(np.uint8), kernel, mode=\"constant\", cval=0) \n", " neighbour_code_map_pred = scipy.ndimage.filters.correlate(cropmask_pred.astype(np.uint8), kernel, mode=\"constant\", cval=0) \n", "\n", " # create masks with the surface voxels\n", " borders_gt = ((neighbour_code_map_gt != 0) & (neighbour_code_map_gt != 255))\n", " borders_pred = ((neighbour_code_map_pred != 0) & (neighbour_code_map_pred != 255))\n", "\n", " # compute the distance transform (closest distance of each voxel to the surface voxels)\n", " if borders_gt.any():\n", " distmap_gt = scipy.ndimage.morphology.distance_transform_edt(~borders_gt, sampling=spacing_mm)\n", " else:\n", " distmap_gt = np.Inf * np.ones(borders_gt.shape)\n", "\n", " if borders_pred.any(): \n", " distmap_pred = scipy.ndimage.morphology.distance_transform_edt(~borders_pred, sampling=spacing_mm)\n", " else:\n", " distmap_pred = np.Inf * np.ones(borders_pred.shape)\n", "\n", " # compute the area of each surface element\n", " surface_area_map_gt = neighbour_code_to_surface_area[neighbour_code_map_gt]\n", " surface_area_map_pred = neighbour_code_to_surface_area[neighbour_code_map_pred]\n", "\n", " # create a list of all surface elements with distance and area\n", " distances_gt_to_pred = distmap_pred[borders_gt]\n", " distances_pred_to_gt = distmap_gt[borders_pred]\n", " surfel_areas_gt = surface_area_map_gt[borders_gt]\n", " surfel_areas_pred = surface_area_map_pred[borders_pred]\n", "\n", " # sort them by distance\n", " if distances_gt_to_pred.shape != (0,):\n", " sorted_surfels_gt = np.array(sorted(zip(distances_gt_to_pred, surfel_areas_gt)))\n", " distances_gt_to_pred = sorted_surfels_gt[:,0]\n", " surfel_areas_gt = sorted_surfels_gt[:,1]\n", "\n", " if distances_pred_to_gt.shape != (0,):\n", " sorted_surfels_pred = np.array(sorted(zip(distances_pred_to_gt, surfel_areas_pred)))\n", " distances_pred_to_gt = sorted_surfels_pred[:,0]\n", " surfel_areas_pred = sorted_surfels_pred[:,1]\n", "\n", "\n", " return {\"distances_gt_to_pred\": distances_gt_to_pred, \n", " \"distances_pred_to_gt\": distances_pred_to_gt, \n", " \"surfel_areas_gt\": surfel_areas_gt, \n", " \"surfel_areas_pred\": surfel_areas_pred}\n" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "colab": { "autoexec": { "startup": false, "wait_interval": 0 } }, "colab_type": "code", "id": "y8shBcxlrUtJ" }, "outputs": [], "source": [ "def compute_average_surface_distance(surface_distances):\n", " distances_gt_to_pred = surface_distances[\"distances_gt_to_pred\"]\n", " distances_pred_to_gt = surface_distances[\"distances_pred_to_gt\"]\n", " surfel_areas_gt = surface_distances[\"surfel_areas_gt\"]\n", " surfel_areas_pred = surface_distances[\"surfel_areas_pred\"]\n", " average_distance_gt_to_pred = np.sum( distances_gt_to_pred * surfel_areas_gt) / np.sum(surfel_areas_gt)\n", " average_distance_pred_to_gt = np.sum( distances_pred_to_gt * surfel_areas_pred) / np.sum(surfel_areas_pred)\n", " return (average_distance_gt_to_pred, average_distance_pred_to_gt)\n", "\n", "def compute_robust_hausdorff(surface_distances, percent):\n", " distances_gt_to_pred = surface_distances[\"distances_gt_to_pred\"]\n", " distances_pred_to_gt = surface_distances[\"distances_pred_to_gt\"]\n", " surfel_areas_gt = surface_distances[\"surfel_areas_gt\"]\n", " surfel_areas_pred = surface_distances[\"surfel_areas_pred\"]\n", " if len(distances_gt_to_pred) > 0:\n", " surfel_areas_cum_gt = np.cumsum(surfel_areas_gt) / np.sum(surfel_areas_gt)\n", " idx = np.searchsorted(surfel_areas_cum_gt, percent/100.0)\n", " perc_distance_gt_to_pred = distances_gt_to_pred[min(idx, len(distances_gt_to_pred)-1)]\n", " else:\n", " perc_distance_gt_to_pred = np.Inf\n", " \n", " if len(distances_pred_to_gt) > 0:\n", " surfel_areas_cum_pred = np.cumsum(surfel_areas_pred) / np.sum(surfel_areas_pred)\n", " idx = np.searchsorted(surfel_areas_cum_pred, percent/100.0)\n", " perc_distance_pred_to_gt = distances_pred_to_gt[min(idx, len(distances_pred_to_gt)-1)]\n", " else:\n", " perc_distance_pred_to_gt = np.Inf\n", " \n", " return max( perc_distance_gt_to_pred, perc_distance_pred_to_gt)\n", "\n", "def compute_surface_overlap_at_tolerance(surface_distances, tolerance_mm):\n", " distances_gt_to_pred = surface_distances[\"distances_gt_to_pred\"]\n", " distances_pred_to_gt = surface_distances[\"distances_pred_to_gt\"]\n", " surfel_areas_gt = surface_distances[\"surfel_areas_gt\"]\n", " surfel_areas_pred = surface_distances[\"surfel_areas_pred\"]\n", " rel_overlap_gt = np.sum(surfel_areas_gt[distances_gt_to_pred