adapt examples/precentral_gyrus_border: add small cubes instead of radiating triangles

examples/precentral_gyrus_border.ipynb

@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Draw Border Around Precentral Gyrus"
+    "# Add Cubes Around Surface of Precentral Gyrus"
    ]
   },
   {
@@ -58,42 +58,44 @@
   {
    "cell_type": "code",
    "execution_count": 3,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "312074"
-      ]
-     },
-     "execution_count": 3,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
    "source": [
+    "import itertools, numpy\n",
     "from freesurfer_surface import Vertex, Triangle\n",
     "\n",
-    "FACTOR = 1.3\n",
-    "for segment in border_chain.segments():\n",
-    "    surface.add_rectangle((\n",
-    "        segment.vertex_indices[0],\n",
-    "        segment.vertex_indices[1],\n",
-    "        surface.add_vertex(Vertex(\n",
-    "            right=surface.vertices[segment.vertex_indices[1]].right * FACTOR,\n",
-    "            anterior=surface.vertices[segment.vertex_indices[1]].anterior * FACTOR,\n",
-    "            superior=surface.vertices[segment.vertex_indices[1]].superior * FACTOR,\n",
-    "        )),\n",
-    "    ))\n",
-    "len(surface.triangles)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "def slice_offset(iterable, offset):\n",
+    "    return itertools.islice(iterable, offset, len(iterable))\n",
+    "\n",
+    "def unit_vector(vector: numpy.array) -> numpy.array:\n",
+    "    return vector / numpy.linalg.norm(vector)\n",
+    "\n",
+    "BORDER_WIDTH = 0.8 # mm\n",
+    "\n",
+    "for vertex_triplet_indices in zip(*(slice_offset(border_chain.vertex_indices, o) for o in range(3))):\n",
+    "    vertex_triplet_coords = tuple(numpy.array(surface.vertices[idx])\n",
+    "                                  for idx in vertex_triplet_indices)\n",
+    "    backward_vector = vertex_triplet_coords[0] - vertex_triplet_coords[1]\n",
+    "    forward_vector = unit_vector(vertex_triplet_coords[2] - vertex_triplet_coords[1]) * BORDER_WIDTH / 2\n",
+    "    upward_vector = unit_vector(numpy.cross(backward_vector, forward_vector)) * BORDER_WIDTH / 2\n",
+    "    sideward_vector = unit_vector(numpy.cross(upward_vector, forward_vector)) * BORDER_WIDTH / 2\n",
+    "    cross_section_corners = [vertex_triplet_coords[1] + v\n",
+    "                             for v in [-upward_vector -sideward_vector,\n",
+    "                                       -upward_vector +sideward_vector,\n",
+    "                                       +upward_vector +sideward_vector,\n",
+    "                                       +upward_vector -sideward_vector]]\n",
+    "    corner_indices = [[surface.add_vertex(Vertex(*(coords + forward_dir))) for coords in cross_section_corners]\n",
+    "                      for forward_dir in [-forward_vector, +forward_vector]]\n",
+    "    for surface_corner_indices in corner_indices:\n",
+    "        surface.triangles.append(Triangle(surface_corner_indices[:3]))\n",
+    "        surface.triangles.append(Triangle(surface_corner_indices[2:] + surface_corner_indices[:1]))\n",
+    "    for i in range(4):\n",
+    "        j = (i + 1) % 4\n",
+    "        surface.triangles.append(Triangle((corner_indices[0][i], corner_indices[0][j], corner_indices[1][j])))\n",
+    "        surface.triangles.append(Triangle((corner_indices[0][i], corner_indices[1][j], corner_indices[1][i])))\n",
+    "\n",
     "surface.write_triangular('precentral-border.lh.pial')"
    ]
   },

freesurfer_surface/__init__.py

@@ -330,6 +330,8 @@ class Surface:
             for vertex in self.vertices:
                 surface_file.write(struct.pack('>fff', *vertex))
             for triangle in self.triangles:
+                assert all(vertex_index < len(self.vertices)
+                           for vertex_index in triangle.vertex_indices)
                 surface_file.write(struct.pack('>III', *triangle.vertex_indices))
             surface_file.write(self._TAG_OLD_USEREALRAS
                                + struct.pack('>I', 1 if self.using_old_real_ras else 0))