{
"cells": [
{
"cell_type": "raw",
"id": "1bc5cc4b",
"metadata": {},
"source": [
"---\n",
"title: \"Day 22: Clearing the Runway\"\n",
"disable_content_template: 1\n",
"tags:\n",
" - 'image processing'\n",
" - 'computer vision'\n",
"author: Zaki Mughal\n",
"data:\n",
" bio: zmughal\n",
" description: Hough transform for line extraction\n",
"images:\n",
" banner:\n",
" src: 'Runway_34,_Nagoya_Airfield_(3937428018).jpg'\n",
" alt: 'Runway 34, Nagoya Airfield'\n",
" data:\n",
" attribution: |-\n",
" Runway 34, Nagoya Airfield, by Kentaro Iemoto from Tokyo, Japan, CC BY-SA 2.0, via Wikimedia Commons\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "9d225ec4",
"metadata": {},
"source": [
"The North Pole Workshop's logistics department has been busy this year with\n",
"upgrading the sleigh's dash-cam with a new more powerful computer. They've\n",
"asked the research and development department if there is anything that can be\n",
"done with the extra cycles to make the journey safer. After bouncing around a\n",
"couple ideas, they came up with a plan for an autoland system…but since Santa\n",
"and the reindeer do not operate in typical environments (such as roofs), this\n",
"would not be a typical autoland system (which use microwave/radio guidance) so\n",
"they couldn't select anything off-the-shelf (and they do know their shelves).\n",
"\n",
"---\n",
"\n",
"To help iterate over the requirements for this new system, the elves knew that\n",
"PDL could help them process images and try different approaches quickly. This\n",
"iterative process will start with simpler scenarios in order to ensure the\n",
"algorithms are on the right track. We'll start with an airport runway as a test\n",
"image.\n",
"\n",
"# Setup\n",
"\n",
"First we need to setup the environment to do everything we need: read images,\n",
"processing them with some filters, and then plot the results. Since image\n",
"processing requires frequent visualisation, this work is being done in Jupyter\n",
"Notebook via [`Devel::IPerl`](https://p3rl.org/Devel::IPerl), but the setup\n",
"below allows for display of plots both in a notebook, but also as a regular\n",
"Perl script. You can [download the notebook here](20241222-hough-lines.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5665996e-85ed-498e-a21e-ca0056566ead",
"metadata": {},
"outputs": [],
"source": [
"use v5.36;\n",
"use utf8;\n",
"use feature qw(signatures postderef);\n",
"\n",
"use constant IN_IPERL => !! $ENV{PERL_IPERL_RUNNING};\n",
"no if IN_IPERL, warnings => 'redefine'; # fewer messages when re-running cells\n",
"no if ! IN_IPERL, warnings => 'void'; # fewer messages for variable at end of cell outside of IPerl\n",
"\n",
"use PDL 2.095;\n",
"use PDL::Constants qw(PI);\n",
"use PDL::IO::Pic;\n",
"use PDL::Image2D;\n",
"use PDL::ImageRGB;\n",
"\n",
"use SVG;\n",
"use MooX::Struct ();\n",
"use Path::Tiny qw(path);\n",
"use List::Util ();\n",
"use Encode qw(encode_utf8);\n",
"\n",
"use Data::Printer;\n",
"use Data::Printer::Filter::PDL ();\n",
"\n",
"use PDL::Graphics::Gnuplot qw(gpwin);\n",
"\n",
"# Make PDL::Graphics::Gnuplot compatible with IPerl\n",
"use if IN_IPERL, 'Devel::IPerl::Plugin::PDLGraphicsGnuplot';\n",
"if( IN_IPERL ) {\n",
"\tIPerl->load_plugin('PDLGraphicsGnuplot');\n",
"}\n",
"\n",
"sub PDL::Graphics::Gnuplot::fancy_display {\n",
"\tif(IN_IPERL){\n",
"\t\tIPerl->display($_[0]);\n",
"\t} else {\n",
"\t\t#sleep 1; $_[0]->close;\n",
"\t\t$_[0]->pause_until_close;\n",
"\t}\n",
"}\n",
"\n",
"my $gp = gpwin();"
]
},
{
"cell_type": "markdown",
"id": "aad7adf7",
"metadata": {},
"source": [
"# Image processing for feature extraction\n",
"\n",
"The elves know that runways are characterized by strong straight lines: the\n",
"edges of the runway, centerline markings, and threshold stripes. If they can\n",
"reliably detect these lines in images from the sleigh's camera, they'll have a\n",
"key piece of their autoland system.\n",
"\n",
"But how do you find straight lines in an image? Let's break this down into steps:\n",
"\n",
"1. First, we need to find places in the image where there are sharp changes in\n",
" brightness — these are edges that might be part of the lines we're looking\n",
" for:\n",
"\n",
" ```mermaid\n",
" graph LR\n",
" A[Color Image] -->|Convert to greyscale| B[Greyscale]\n",
" B -->|Smooth noise| C[Blurred]\n",
" C -->|Find edges| D[Edge points]\n",
" ```\n",
"\n",
"2. But just having edge points isn't enough - we need to figure out which ones\n",
" form straight lines. This is tricky because:\n",
" - edges might be broken (gaps in runway lines);\n",
" - there's usually noise (not all edges are runway lines); and\n",
" - lines might be partial (only part of the runway visible).\n",
"\n",
"3. This is where the Hough transform comes in. Instead of trying to connect\n",
" edge points directly, it:\n",
" 1. Takes each edge point.\n",
" 2. Finds all possible lines through that point.\n",
" 3. Lets the points \"vote\" on which lines they might be part of.\n",
" 4. Finds the lines that got the most votes.\n",
"\n",
"The rest of this notebook shows how we implement this pipeline using PDL, step\n",
"by step.\n",
"\n",
"# Prepare image\n",
"\n",
"After selecting a test image\n",
"(Runway 34, Nagoya Airfield,\n",
"by Kentaro Iemoto from Tokyo, Japan,\n",
"CC BY-SA 2.0, via Wikimedia Commons),\n",
"we need to read it in using `rpic()`\n",
"and take a look at the dimensions."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e9522290",
"metadata": {
"tags": [
"hide-output-execute_result"
]
},
"outputs": [],
"source": [
"# https://commons.wikimedia.org/wiki/File:Runway_34,_Nagoya_Airfield_(3937428018).jpg\n",
"my $runway = rpic('Runway_34,_Nagoya_Airfield_(3937428018).jpg');\n",
"\n",
"say np $runway;"
]
},
{
"cell_type": "markdown",
"id": "ac139252",
"metadata": {},
"source": [
"and then view the image in Gnuplot:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3abe81c7",
"metadata": {},
"outputs": [],
"source": [
"$gp->image( $runway );\n",
"$gp->fancy_display;"
]
},
{
"cell_type": "markdown",
"id": "a0d0dc13",
"metadata": {},
"source": [
"Note the dimensions are `[3 W H]` where 3 is the number of channels, `W` is the\n",
"width (x-extent), and `H` is the height (y-extent).\n",
"\n",
"We can't process the color channels in the image directly, so we need to\n",
"convert it to a single greyscale (luminance) channel. We can use `rgbtogr`\n",
"which uses a standard formula."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "024b7d7c",
"metadata": {},
"outputs": [],
"source": [
"my $grey = rgbtogr($runway)->double / 255;\n",
"$gp->image($grey,\n",
"\t{ title => 'Greyscale image',\n",
"\t clut => 'grey' });\n",
"$gp->fancy_display;"
]
},
{
"cell_type": "markdown",
"id": "978b7b02",
"metadata": {},
"source": [
"# Gaussian blur\n",
"\n",
"We need to find edges in the image, but all real-world images have some amount\n",
"of noise (e.g., from random variation in brightness from pixel to pixel).\n",
"This noise can interfere with further processing such as edge\n",
"detection such that the noise may be picked up as spurious edges.\n",
"\n",
"To smooth out this noise, we apply a 2D Gaussian filter to the image.\n",
"The filter:\n",
"\n",
"1. Looks at each pixel and its neighbors.\n",
"2. Creates a weighted average where:\n",
" - The center pixel counts most;\n",
" - Nearby pixels count less;\n",
" - Far pixels barely count at all;\n",
" - The weights follow a bell curve (Gaussian) shape.\n",
"\n",
"The parameter σ (sigma) controls how much smoothing we do:\n",
"\n",
"- Small σ: less smoothing, keeps more detail.\n",
"- Large σ: more smoothing, might blur real edges.\n",
"\n",
"The way we can perform this weighted average is by performing convolution\n",
"with the `conv2d()` function.\n",
"\n",
"\n",
"Implementation note:\n",
"The 2D isotropic Gaussian is a separable filter, which means that\n",
"\n",
"$$G(x,y) = e^{-\\frac{x^2 + y^2}{2\\sigma^2}}$$\n",
"\n",
"can be written as the product of two 1D Gaussians:\n",
"\n",
"$$G(x,y) = e^{-\\frac{x^2}{2\\sigma^2}} \\cdot e^{-\\frac{y^2}{2\\sigma^2}}$$\n",
"\n",
"This lets us compute 1D Gaussian values and use this to create the 2D Gaussian\n",
"values for the kernel, which is much more efficient than computing the full 2D\n",
"function directly (though the kernel is typically small).\n",
"It is also possible to use two consecutive `conv1d()` calls here, but we'll\n",
"skip that for now.\n",
"\n",
"Note that the above formulae are written without the normalisation factor.\n",
"Here we approximate the normalisation factor by dividing by the sum.\n",
"\n",
"See more at:\n",
"\n",
"- R. Fisher, S. Perkins, A. Walker, and E. Wolfart,\n",
" [Spatial Filters - Gaussian Smoothing](https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm).\n",
" \"Hypermedia Image Processing Reference (HIPR)\", 2003.\n",
"- [Gaussian function - Wikipedia](https://en.wikipedia.org/wiki/Gaussian_function#Two-dimensional_Gaussian_function).\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d5427ca2",
"metadata": {},
"outputs": [],
"source": [
"sub gaussian_kernel( $sigma = 1.0 ) {\n",
"\tdie 'sigma must be positive' unless $sigma > 0;\n",
"\n",
"\tmy $GAUSSIAN_CUTOFF = 3; # Standard deviations from center (±3σ)\n",
"\n",
"\t# Kernel size: Ensure it's odd for symmetry around the center pixel\n",
"\t# Total width = 2 * 3σ + 1.\n",
"\tmy $size = pdl(2 * ceil($GAUSSIAN_CUTOFF * $sigma) + 1)->sclr;\n",
"\n",
"\t# Create 1D Gaussian values centered at 0\n",
"\t# using symmetric x-coordinates.\n",
"\tmy $x = xvals($size) - ($size-1)/2;\n",
"\tmy $gauss_1d = exp(-($x**2)/(2 * $sigma**2));\n",
"\n",
"\t# Create 2D kernel using outer product\n",
"\tmy $kernel = $gauss_1d->outer($gauss_1d);\n",
"\t$kernel /= $kernel->sum; # normalise\n",
"\n",
"\treturn $kernel;\n",
"}\n",
"\n",
"sub gaussian_blur($img, $sigma = 1) {\n",
"\treturn conv2d($img, gaussian_kernel($sigma));\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "750dae95",
"metadata": {},
"source": [
"## Test Gaussian blur on a point source\n",
"\n",
"We can see how this works by testing on a simple \"point source\" image — a\n",
"single bright pixel surrounded by darkness:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "eefffa53",
"metadata": {
"lines_to_next_cell": 2,
"tags": [
"hide-output-execute_result"
]
},
"outputs": [],
"source": [
"# Create a point source.\n",
"my $point = zeroes(21,21);\n",
"$point->set(10,10, 1);"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1845c98d",
"metadata": {},
"outputs": [],
"source": [
"my $point_blurred = gaussian_blur($point, 1);\n",
"$gp->image($point_blurred,\n",
"\t{ title => encode_utf8('Gaussian Blur of Point Source (σ=1)'),\n",
"\t clut => 'grey', },\n",
");\n",
"$gp->fancy_display;"
]
},
{
"cell_type": "markdown",
"id": "fb2d39e5",
"metadata": {},
"source": [
"We can see this took the single point and \"spread out\" its intensity in a\n",
"radial manner. Furthermore, `$point_blurred` continues to sum up to 1 which\n",
"means that brightness is preserved due to the normalisation factor."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b37f719a",
"metadata": {
"tags": [
"hide-output-execute_result"
]
},
"outputs": [],
"source": [
"say \"Σ p = \", $point_blurred->sum->sclr;"
]
},
{
"cell_type": "markdown",
"id": "3675f1cf",
"metadata": {},
"source": [
"## Use Gaussian blur on image\n",
"\n",
"We can now apply the Gaussian blur filter to the runway image."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "839797b0",
"metadata": {},
"outputs": [],
"source": [
"my $blurred = gaussian_blur($grey, 2);\n",
"$gp->image($blurred,\n",
"\t{ title => encode_utf8('Gaussian Blur (σ=2)'),\n",
"\t clut => 'grey', },\n",
");\n",
"$gp->fancy_display;"
]
},
{
"cell_type": "markdown",
"id": "99a4d67f",
"metadata": {},
"source": [
"# Sobel edge detection\n",
"\n",
"Now that we've smoothed out the noise, we can look for edges. But what exactly\n",
"is an edge in an image? An edge occurs where there's a sharp change in\n",
"brightness — like the transition from the dark runway surface to the bright\n",
"markings.\n",
"\n",
"Mathematically, these changes in brightness are derivatives: how quickly is the\n",
"brightness changing as we move across the image? The Sobel operator helps us\n",
"find these changes by looking at how brightness differs between neighboring\n",
"pixels in both horizontal (x) and vertical (y) directions.\n",
"\n",
"\n",
"Implementation note:\n",
"The Sobel operator uses two 3×3 kernels:\n",
"\n",
"For horizontal changes (x-direction):\n",
"\n",
"$$K_x = \\frac{1}{8}\\begin{bmatrix}\n",
"-1 & 0 & +1 \\\\\n",
"-2 & 0 & +2 \\\\\n",
"-1 & 0 & +1\n",
"\\end{bmatrix}$$\n",
"\n",
"For vertical changes (y-direction):\n",
"\n",
"$$K_y = \\frac{1}{8}\\begin{bmatrix}\n",
"-1 & -2 & -1 \\\\\n",
" 0 & 0 & 0 \\\\\n",
"+1 & +2 & +1\n",
"\\end{bmatrix}$$\n",
"\n",
"We calculate both directions and combine them to get the total edge strength:\n",
"\n",
"\n",
"strength² = (horizontal change)² + (vertical change)²\n",
"\n",
"The factor of ⅛ \n",
"normalizes the kernels so the maximum response to a\n",
"perfect edge is 1.\n",
"\n",
"Note: We store strength² directly rather than taking its square root since:\n",
"\n",
"1. We only care about relative strengths.\n",
"2. Avoiding the square root is more efficient.\n",
"3. The threshold comparison will still work the same.\n",
"\n",
"See more at:\n",
"\n",
"* [Sobel operator - Wikipedia](https://en.wikipedia.org/wiki/Sobel_operator)\n",
"* [Edge detection - Wikipedia](https://en.wikipedia.org/wiki/Edge_detection)\n",
"
\n",
"\n",
"We can threshold this strength to get a **binary edge image**: pixels where the\n",
"strength is above some threshold become 1 (edge), and others become 0\n",
"(non-edge). The threshold is chosen relative to the average edge strength in\n",
"the image — this helps it adapt to different images automatically.\n",
"\n",
"\n",
"Note:\n",
"Often further post-processing of the edge detection result is done such as\n",
"thinning which can help create more clear edges without spurious points, but we\n",
"won't do that here.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0d03eb25",
"metadata": {},
"outputs": [],
"source": [
"sub sobel_operator($img) {\n",
"\t# Define Sobel kernel and normalise it\n",
"\tmy $kx = pdl([\n",
"\t\t[-1, 0, 1],\n",
"\t\t[-2, 0, 2],\n",
"\t\t[-1, 0, 1],\n",
"\t]);\n",
"\t$kx /= $kx->abs->sum;\n",
"\n",
"\t# Stack x and y kernels\n",
"\tmy $k = cat($kx, $kx->transpose);\n",
"\n",
"\t# Add dummy dimension to image for conv2d\n",
"\tmy $img_3d = $img->dummy(-1);\n",
"\n",
"\t# Single convolution and magnitude calculation.\n",
"\t# Returns strength² (no sqrt).\n",
"\treturn conv2d($img_3d, $k)->mv(-1,0)->pow(2)->sumover;\n",
"}\n",
"\n",
"sub default_edge_threshold($strength2) {\n",
"\tmy $mean_strength2 = $strength2->avg;\n",
"\tmy $thresh2 = 4 * $mean_strength2;\n",
"\treturn $strength2 > $thresh2;\n",
"}\n",
"\n",
"sub sobel_edges($img) {\n",
"\treturn default_edge_threshold(sobel_operator($img));\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "b8cb564f",
"metadata": {},
"source": [
"## Test Sobel edge detection on point source\n",
"\n",
"First, let's see how edge detection works on our point source test image. This\n",
"is useful because:\n",
"\n",
"1. We know exactly what the input looks like (single bright pixel)\n",
"2. The edges should appear where brightness changes most rapidly\n",
"3. With a single point, we expect to see edges form a simple pattern\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5ad6b0c8",
"metadata": {},
"outputs": [],
"source": [
"# Test sobel_edges on point source\n",
"my $point_edges = sobel_edges($point);\n",
"$gp->image($point_edges,\n",
"\t{ title => 'Sobel Edges of Point Source',\n",
"\t clut => 'grey', },\n",
");\n",
"$gp->fancy_display;"
]
},
{
"cell_type": "markdown",
"id": "59323348",
"metadata": {},
"source": [
"The result shows edges right next to the original point — this makes sense\n",
"because that's where the brightness changes most sharply.\n",
"\n",
"## Use Sobel edge detection on image\n",
"\n",
"Now let's apply edge detection directly to our blurred runway image:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "574172fd",
"metadata": {},
"outputs": [],
"source": [
"my $binary_edges = sobel_edges($blurred);\n",
"$gp->image($binary_edges,\n",
"\t{ title => 'Sobel Edges Binarised',\n",
"\t clut => 'grey', },\n",
");\n",
"$gp->fancy_display;"
]
},
{
"cell_type": "markdown",
"id": "d1e9d179",
"metadata": {},
"source": [
"Looking at the result:\n",
"\n",
"- Strong edges appear along the runway boundaries.\n",
"- The runway markings create clear edges.\n",
"\n",
"These edge points are what we'll use for line detection, but we still have a\n",
"challenge: how do we figure out which edge points belong to the same straight\n",
"line? We could try to trace along contiguous pixels, but this can not deal with\n",
"the less well-defined edges (e.g., broken disconnected edges due to shadows).\n",
"\n",
"That's where the Hough transform comes in.\n",
"\n",
"# Hough transform\n",
"\n",
"As mentioned earlier, the Hough transform works by using points in the image to\n",
"\"vote\" for which line they are in. Well, technically, the Hough transform is\n",
"more general in that it can be used for any model to match against as long as\n",
"you have a valid parameterisation of the model in the image space.\n",
"\n",
"There are widely used parameterisations for lines and circles. We'll be using\n",
"the line parameterisation because we are looking for lines.\n",
"\n",
"## Define line parameterisation\n",
"\n",
"What do we mean by parameterisation? It's just a fancy way of saying\n",
"what variables do you need to define a shape. For example, a 2D line segment\n",
"needs four parameters $(x_0, y_0), (x_1, y_1)$. Similarly, a 2D line can be\n",
"parameterised using a [slope-intercept form](https://en.wikipedia.org/wiki/Linear_equation#Slope%E2%80%93intercept_form_or_Gradient-intercept_form)\n",
"$y = mx + y_0$ with the two parameters $m$ and $y_0$.\n",
"\n",
"We're looking for lines, so that could work, but it has issues when dealing\n",
"with vertical or nearly vertical lines. Our parameter range for $m$ would have\n",
"to be very large in those cases. An alternative parameterisation is\n",
"[Hesse normal form](https://en.wikipedia.org/wiki/Hesse_normal_form) which for\n",
"a 2D line follows the equation\n",
"$$ \\rho = x \\cos \\theta + y \\sin \\theta $$\n",
"where given a target line, $\\rho$ is the shortest distance between that line and the\n",
"origin and $\\theta$ is the angle between the axis and a line of length $\\rho$\n",
"that intersects the target line. Maybe that's easier to see with a picture:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4a4ab46c",
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"sub create_hesse_diagram {\n",
"\tmy $svg = SVG->new(width => 200, height => 200);\n",
"\n",
"\t# Define arrowhead marker\n",
"\tmy $defs = $svg->defs();\n",
"\tmy $marker = $defs->marker(\n",
"\t\tid => 'arrowhead',\n",
"\t\tviewBox => '0 0 10 10',\n",
"\t\trefX => 10,\n",
"\t\trefY => 5,\n",
"\t\tmarkerWidth => 6,\n",
"\t\tmarkerHeight => 6,\n",
"\t\torient => 'auto',\n",
"\t);\n",
"\t$marker->path(\n",
"\t\td => 'M 0 0 L 10 5 L 0 10 z',\n",
"\t\tstyle => { fill => '#999' }\n",
"\t);\n",
"\n",
"\t# Set coordinate system with origin at center\n",
"\tmy $g = $svg->group(\n",
"\t\ttransform => 'translate(25,25) scale(1,1)',\n",
"\t\tstyle => {\n",
"\t\t\t'stroke-width' => '1.5',\n",
"\t\t\t'font-family' => 'sans-serif',\n",
"\t\t\t'text-anchor' => 'middle',\n",
"\t\t}\n",
"\t);\n",
"\n",
"\t# Axes (dashed)\n",
"\t$g->line(\n",
"\t\tx1 => 0, y1 => 0, x2 => 140, y2 => 0,\n",
"\t\tstyle => {\n",
"\t\t\t'stroke' => '#999',\n",
"\t\t\t'stroke-dasharray' => '4,4',\n",
"\t\t\t'marker-end' => 'url(#arrowhead)',\n",
"\t\t}\n",
"\t);\n",
"\t$g->line(\n",
"\t\tx1 => 0, y1 => 0, x2 => 0, y2 => 140,\n",
"\t\tstyle => {\n",
"\t\t\t'stroke' => '#999',\n",
"\t\t\t'stroke-dasharray' => '4,4',\n",
"\t\t\t'marker-end' => 'url(#arrowhead)',\n",
"\t\t}\n",
"\t);\n",
"\n",
"\t# Blue ρ line at 45°\n",
"\t$g->line(\n",
"\t\tx1 => 0, y1 => 0, x2 => 50, y2 => 50,\n",
"\t\tstyle => { stroke => 'blue' }\n",
"\t);\n",
"\n",
"\t# Red perpendicular line at -45°\n",
"\t$g->line(\n",
"\t\tx1 => 0, y1 => 100, x2 => 100, y2 => 0,\n",
"\t\tstyle => { stroke => 'red' }\n",
"\t);\n",
"\n",
"\t# Perpendicular marker (rotated empty square)\n",
"\t$g->rect(\n",
"\t\tx => 50, y => 50, width => 8, height => 8,\n",
"\t\ttransform => 'rotate(225, 50, 50)',\n",
"\t\tstyle => { fill => 'none', stroke => 'black', 'stroke-width' => 1 }\n",
"\t);\n",
"\n",
"\t# θ arc (45°)\n",
"\t$g->path(\n",
"\t\td => 'M 20,0 A 20,20 0 0,1 14.14,14.14',\n",
"\t\tstyle => { stroke => 'green', fill => 'none' }\n",
"\t);\n",
"\n",
"\t# Labels (flipped for correct orientation)\n",
"\tmy $t = $g->group(transform => 'scale(1,1)');\n",
"\t$t->text(x => 25, y => 45, style => { 'font-size' => '16px' })->cdata('ρ');\n",
"\t$t->text(x => 30, y => 20, style => { 'font-size' => '16px' })->cdata('θ');\n",
"\t$t->text(x => -10, y => -5, style => { 'font-size' => '12px' })->cdata('O');\n",
"\t$t->text(x => 145, y => -5, style => { 'font-size' => '12px' })->cdata('x');\n",
"\t$t->text(x => 5, y => 150, style => { 'font-size' => '12px' })->cdata('y');\n",
"\n",
"\t# Origin point\n",
"\t$g->circle(\n",
"\t\tcx => 0, cy => 0, r => 2,\n",
"\t\tstyle => { fill => 'black' }\n",
"\t);\n",
"\n",
"\treturn $svg;\n",
"}\n",
"\n",
"# Save diagram\n",
"if(IN_IPERL) {\n",
"\tIPerl->svg( bytestream => encode_utf8(create_hesse_diagram()->xmlify), width => '50%' );\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "dad99d5b",
"metadata": {},
"source": [
"This diagram uses the top-left as the origin as is common with image processing\n",
"code. For a fixed θ, different values of ρ represent parallel lines (note that\n",
"ρ can be negative). Values of θ can range from -π/2 to π/2 which allows for\n",
"parametrising all possible lines including vertical and horizontal lines.\n",
"\n",
"We can implement an object to hold all this information in a single place.\n",
"By default, the parameter space is discretised to ~1 pixel resolution in ρ\n",
"and 1° resolution in θ."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "31e39f6d",
"metadata": {},
"outputs": [],
"source": [
"use MooX::Struct -retain, HoughParameters => [\n",
"\t# Rho parameters\n",
"\trho_range => [ required => 1 ],\n",
"\tnum_rho => [ required => 1 ],\n",
"\n",
"\t# Theta parameters\n",
"\ttheta_range => [ required => 1 ],\n",
"\tnum_theta => [ required => 1 ],\n",
"\n",
"\t# Lazy builders for actual parameter values\n",
"\trho_res => [ builder => '_build_rho_res' ],\n",
"\ttheta_res => [ builder => '_build_theta_res' ],\n",
"\trhos => [ builder => '_build_rhos' ],\n",
"\tthetas => [ builder => '_build_thetas' ],\n",
"\t_build_rho_res => sub {\n",
"\t\tmy $self = shift;\n",
"\t\tmy ($min, $max) = @{$self->rho_range};\n",
"\t\t($max - $min) / ($self->num_rho - 1);\n",
"\t},\n",
"\t_build_theta_res => sub {\n",
"\t\tmy $self = shift;\n",
"\t\tmy ($min, $max) = @{$self->theta_range};\n",
"\t\t($max - $min) / ($self->num_theta - 1);\n",
"\t},\n",
"\t_build_rhos => sub {\n",
"\t\tmy $self = shift;\n",
"\t\tsequence($self->num_rho)->xlinvals(@{$self->rho_range});\n",
"\t},\n",
"\t_build_thetas => sub {\n",
"\t\tmy $self = shift;\n",
"\t\tsequence($self->num_theta)->xlinvals(@{$self->theta_range});\n",
"\t},\n",
"\n",
"\t# Class method for default parameters from image\n",
"\tfrom_image => sub {\n",
"\t\tmy ($class, $img, %rest) = @_;\n",
"\t\tmy $max_rho = $img->shape->magnover;\n",
"\t\treturn $class->new(\n",
"\t\t\trho_range => [-$max_rho, $max_rho],\n",
"\t\t\tnum_rho => int(2 * $max_rho), # ~1 pixel resolution\n",
"\t\t\ttheta_range => [-PI/2, PI/2],\n",
"\t\t\tnum_theta => 180, # ~1 degree resolution\n",
"\t\t\t%rest,\n",
"\t\t);\n",
"\t},\n",
"];"
]
},
{
"cell_type": "markdown",
"id": "6ab1e425",
"metadata": {},
"source": [
"## Build the Hough space accumulator\n",
"\n",
"The voting procedure works by varying θ and finding the ρ for a line that a\n",
"given edge point would be part of. These are then accumulated into a 2D ndarray\n",
"using `histogram2d()` which places the votes in the appropriate bin for each\n",
"parameter.\n",
"\n",
"This accumulator is also known as the Hough space."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d5546196",
"metadata": {},
"outputs": [],
"source": [
"sub line_hough($binary_edges, $params) {\n",
"\t# Create accumulator array\n",
"\tmy $accumulator = zeroes($params->num_rho, $params->num_theta);\n",
"\n",
"\t# Find edge points and convert to 2D matrix of [x y] coordinates\n",
"\t# shape [2,N].\n",
"\tmy $edge_points = whichND($binary_edges);\n",
"\n",
"\tsay '[INFO] ', \"Processing \", $edge_points->dim(1), \" edge points\";\n",
"\n",
"\t# Create transformation matrix [180,2]\n",
"\t# for [x,y] coordinate order\n",
"\tmy $transform = cat(\n",
"\t\tcos($params->thetas), # for x\n",
"\t\tsin($params->thetas), # for y\n",
"\t);\n",
"\n",
"\t# Calculate ρ = edge_points × transform giving [180,N]\n",
"\tmy $rho = $edge_points x $transform;\n",
"\n",
"\t# Scale rho values to accumulator indices using params values\n",
"\tmy $rho_min = $params->rho_range->[0];\n",
"\tmy $rho_max = $params->rho_range->[1];\n",
"\tmy $rho_scale = ($params->num_rho - 1) / ($rho_max - $rho_min);\n",
"\tmy $rho_idx = ($rho - $rho_min) * $rho_scale;\n",
"\n",
"\t# Create theta indices [180,N] where each row is 0..179\n",
"\tmy $theta_idx = zeros($params->num_theta, $rho_idx->dim(1));\n",
"\t$theta_idx += sequence($params->num_theta);\n",
"\n",
"\t# Filter valid rho indices\n",
"\tmy $valid = ($rho_idx >= 0) & ($rho_idx < $params->num_rho);\n",
"\t$rho_idx = $rho_idx->where($valid);\n",
"\t$theta_idx = $theta_idx->where($valid);\n",
"\n",
"\t# Use histogram2d to accumulate votes\n",
"\t$accumulator = histogram2d($rho_idx, $theta_idx,\n",
"\t 1, 0, $params->num_rho,\n",
"\t 1, 0, $params->num_theta);\n",
"\n",
"\treturn $accumulator;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "25cb179d",
"metadata": {},
"source": [
"## Find lines (peaks) in the Hough space\n",
"\n",
"A maximum in the Hough space represents a line where many edge points\n",
"participate and those edge points are thus likely collinear. The following\n",
"looks for those peaks iteratively, but suppresses neighbouring values (which\n",
"represent nearby lines) so that they are not detected as peaks in later\n",
"iterations. The `range()` function is the perfect tool for extracting out a\n",
"neighbourhood and modifying it. However, be careful that you use `sever()` so\n",
"that you don't accidentally end up zeroing out the peak value!"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9d806440",
"metadata": {},
"outputs": [],
"source": [
"sub find_lines($accumulator, $params, $threshold_ratio=0.3, $nhood_size=[5,5], $max_peaks=20) {\n",
"\t# Pre-allocate peaks array\n",
"\tmy $peak_indices = zeroes(2, $max_peaks); # [rho_idx, theta_idx] × max_peaks\n",
"\tmy $peak_values = zeroes($max_peaks);\n",
"\tmy $num_peaks = 0;\n",
"\n",
"\t# Calculate deltas for neighbourhood\n",
"\tmy $rho_theta_delta = pdl(@$nhood_size);\n",
"\n",
"\t# Initialize working copy and get flattened view.\n",
"\tmy $H = $accumulator->copy;\n",
"\tmy $H_flat = $H->clump(-1);\n",
"\tmy $threshold = $H_flat->index(maximum_ind($H_flat)) * $threshold_ratio;\n",
"\n",
"\twhile ($num_peaks < $max_peaks) {\n",
"\t\t# Find peak location and value using maximum_ind\n",
"\t\tmy $max_ind = maximum_ind($H_flat);\n",
"\t\tmy $max_val = $H_flat->index($max_ind)->sever;\n",
"\n",
"\t\tlast if $max_val < $threshold;\n",
"\n",
"\t\tmy $peak_idx = pdl(one2nd($H, $max_ind));\n",
"\n",
"\t\t# Record peak\n",
"\t\t$peak_indices->slice(\":,$num_peaks\") .= $peak_idx;\n",
"\t\t$peak_values->index($num_peaks) .= $max_val;\n",
"\t\t$num_peaks++;\n",
"\n",
"\t\t# Suppress neighbourhood using range.\n",
"\t\t# Adjust center point back by half the neighbourhood size.\n",
"\t\tmy $corner = $peak_idx - $rho_theta_delta/2;\n",
"\t\tmy $r = $H->range($corner, $rho_theta_delta, 't');\n",
"\t\t$r .= 0;\n",
"\t}\n",
"\n",
"\t# Truncate results if fewer peaks found\n",
"\tif ($num_peaks < $max_peaks) {\n",
"\t\t$peak_indices = $peak_indices->slice(':,0:' . ($num_peaks-1));\n",
"\t\t$peak_values = $peak_values->slice('0:' . ($num_peaks-1));\n",
"\t}\n",
"\n",
"\t# Create parameter space values directly\n",
"\tmy $peak_params = cat(\n",
"\t\t$params->rhos->index($peak_indices->using(0)),\n",
"\t\t$params->thetas->index($peak_indices->using(1)),\n",
"\t);\n",
"\n",
"\treturn ($peak_params, $peak_values);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "0ec3acba",
"metadata": {},
"source": [
"## Extract lines from peaks\n",
"\n",
"Then using these peaks, we transform back from the Hough space to the binary\n",
"edge image space in order to find the edge points that participate in that line.\n",
"These segments might be disconnected, so there is some tolerance for that in\n",
"the `$fillgap` parameter."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4e9c4aeb",
"metadata": {},
"outputs": [],
"source": [
"sub extract_line_segments($edges, $peak_params, $fillgap=20, $minlength=40) {\n",
"\tmy ($width, $height) = $edges->dims;\n",
"\n",
"\t# Find all edge points in image as [2,N] array of [x,y] coordinates\n",
"\tmy $edge_points = whichND($edges);\n",
"\n",
"\t# Get parameters for each line from [N,2] array\n",
"\tmy $rhos = $peak_params->slice(':,(0)');\n",
"\tmy $thetas = $peak_params->slice(':,(1)');\n",
"\n",
"\tmy @all_segments;\n",
"\tmy $delta_rho = 0.5; # Half-pixel tolerance\n",
"\n",
"\tfor my $i (0..$rhos->nelem-1) {\n",
"\t\tmy $rho = $rhos->at($i);\n",
"\t\tmy $theta = $thetas->at($i);\n",
"\n",
"\t\t# Project points onto normal direction (cos θ, sin θ) for [x,y] coordinates\n",
"\t\tmy $normal_proj = pdl([cos($theta), sin($theta)]);\n",
"\t\tmy $rho_all = inner($edge_points, $normal_proj);\n",
"\n",
"\t\t# Find points belonging to this line\n",
"\t\tmy $points_idx = which(abs($rho_all - $rho) <= $delta_rho);\n",
"\t\tnext unless $points_idx->nelem >= 2;\n",
"\n",
"\t\t# Get coordinates of line points, maintaining [2,N] structure\n",
"\t\tmy $points = $edge_points->slice([],$points_idx);\n",
"\n",
"\t\t# Calculate spans using minmaximum\n",
"\t\tmy $spans = pdl((minmaximum($points->xchg(-1,0)))[0,1])->xchg(0,-1)->diff2->squeeze;\n",
"\n",
"\t\t# Sort based on dominant direction\n",
"\t\tmy $sort_idx = $spans->at(0) > $spans->at(1)\n",
"\t\t\t# More horizontal - sort by x, then y\n",
"\t\t\t? qsorti($points->slice('(0)') * $height + $points->slice('(1)'))\n",
"\t\t\t# More vertical - sort by y, then x\n",
"\t\t\t: qsorti($points->slice('(1)') * $width + $points->slice('(0)'));\n",
"\n",
"\t\t$points = $points->slice([],$sort_idx);\n",
"\n",
"\t\t# Find gaps using point-to-point distances\n",
"\t\tmy $diffs = $points->xchg(0,1)->diff2->xchg(0,1);\n",
"\t\tmy $gap_idx = which($diffs->magnover > $fillgap);\n",
"\n",
"\t\t# Create segment bounds - [2,N] array where dimension 0 is [start,end]\n",
"\t\tmy $seg_bounds = zeroes(2, $gap_idx->nelem + 1);\n",
"\n",
"\t\t# Handle segment bounds based on gaps\n",
"\t\tif ($gap_idx->nelem) {\n",
"\t\t\t$seg_bounds .= cat(\n",
"\t\t\t\tpdl(0)->append($gap_idx + 1),\n",
"\t\t\t\t$gap_idx->append($points->dim(1) - 1)\n",
"\t\t\t)->transpose;\n",
"\t\t} else {\n",
"\t\t\t$seg_bounds .= pdl([0, $points->dim(1) - 1]);\n",
"\t\t}\n",
"\n",
"\t\t# Process each segment\n",
"\t\tfor my $j (0..$seg_bounds->dim(1)-1) {\n",
"\t\t\tmy $endpoints = $points->slice([],$seg_bounds->slice(\",($j)\"));\n",
"\t\t\tmy $length = $endpoints->xchg(0,1)->diff2->xchg(0,1)->magnover->sclr;\n",
"\n",
"\t\t\tnext unless $length >= $minlength;\n",
"\n",
"\t\t\tpush @all_segments, {\n",
"\t\t\t\tpoint1 => $endpoints->slice(',(0)')->unpdl,\n",
"\t\t\t\tpoint2 => $endpoints->slice(',(1)')->unpdl,\n",
"\t\t\t\trho => $rho,\n",
"\t\t\t\ttheta => $theta,\n",
"\t\t\t\tlength => $length,\n",
"\t\t\t\tdominant_axis => $spans->at(0) > $spans->at(1) ? 'x' : 'y',\n",
"\t\t\t};\n",
"\t\t}\n",
"\t}\n",
"\n",
"\treturn \\@all_segments;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "65d03a84",
"metadata": {},
"source": [
"# Pipeline\n",
"\n",
"Now let's put this all together and take a look at what this looks like.\n",
"\n",
"## Test pattern\n",
"\n",
"First to ensure that various parts of this algorithm are working, let's go even\n",
"simpler and create a synthetic test image that only contains various lines\n",
"which are guaranteed to be edges.\n",
"\n",
"An easy way to do this is to create a vector image using SVG and then render it\n",
"to a raster image."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8ada6873",
"metadata": {},
"outputs": [],
"source": [
"sub create_test_pattern($width=400, $height=400) {\n",
"\t# Create SVG object\n",
"\tmy $svg = SVG->new(\n",
"\t\twidth => $width,\n",
"\t\theight => $height,\n",
"\t);\n",
"\n",
"\t# White background\n",
"\t$svg->rectangle(\n",
"\t\tx => 0,\n",
"\t\ty => 0,\n",
"\t\twidth => $width,\n",
"\t\theight => $height,\n",
"\t\tfill => 'white'\n",
"\t);\n",
"\n",
"\t# Create group for all lines with common style\n",
"\tmy $lines = $svg->group(\n",
"\t\tstyle => {\n",
"\t\t\tstroke => 'black',\n",
"\t\t\t'stroke-width' => 2\n",
"\t\t}\n",
"\t);\n",
"\n",
"\t# Complete test lines\n",
"\t$lines->line(x1 => 50, y1 => 50, x2 => 350, y2 => 50); # Horizontal\n",
"\t$lines->line(x1 => 350, y1 => 50, x2 => 350, y2 => 350); # Vertical\n",
"\t$lines->line(x1 => 50, y1 => 50, x2 => 350, y2 => 350); # 45 degrees\n",
"\n",
"\t# Standalone edges and partial segments\n",
"\t$lines->line(x1 => 50, y1 => 150, x2 => 150, y2 => 150); # Broken horizontal 1\n",
"\t$lines->line(x1 => 200, y1 => 150, x2 => 300, y2 => 150); # Broken horizontal 2\n",
"\t$lines->line(x1 => 250, y1 => 200, x2 => 300, y2 => 230); # Short diagonal\n",
"\t$lines->line(x1 => 280, y1 => 280, x2 => 320, y2 => 280); # Isolated short\n",
"\n",
"\t# Merge test cases\n",
"\t# Nearly parallel vertical lines\n",
"\t$lines->line(x1 => 100, y1 => 200, x2 => 100, y2 => 350);\n",
"\t$lines->line(x1 => 102, y1 => 200, x2 => 102, y2 => 350);\n",
"\n",
"\t# Slightly diverging lines\n",
"\t$lines->line(x1 => 150, y1 => 200, x2 => 180, y2 => 350);\n",
"\t$lines->line(x1 => 155, y1 => 200, x2 => 190, y2 => 350);\n",
"\n",
"\t# Closely spaced parallel horizontal segments\n",
"\t$lines->line(x1 => 200, y1 => 200, x2 => 300, y2 => 200);\n",
"\t$lines->line(x1 => 200, y1 => 203, x2 => 300, y2 => 203);\n",
"\n",
"\t# Nearly coincident short segments\n",
"\t$lines->line(x1 => 250, y1 => 100, x2 => 300, y2 => 120);\n",
"\t$lines->line(x1 => 251, y1 => 101, x2 => 301, y2 => 121);\n",
"\n",
"\t# Create temporary files using Path::Tiny\n",
"\tmy $svg_file = Path::Tiny->tempfile(SUFFIX => '.svg');\n",
"\tmy $png_file = Path::Tiny->tempfile(SUFFIX => '.png');\n",
"\n",
"\t# Save SVG\n",
"\t$svg_file->spew_utf8($svg->xmlify);\n",
"\n",
"\t# Convert to PNG using full path strings\n",
"\tsystem('convert', '-density', '150', $svg_file->absolute, $png_file->absolute) == 0\n",
"\t\tor die \"ImageMagick convert failed: $!\";\n",
"\n",
"\t# Read PNG into PDL\n",
"\tmy $img = rpic(\"$png_file\");\n",
"\n",
"\t# Binary (edge = 1, background = 0)\n",
"\tmy $binary = $img->double / 255 < 0.5;\n",
"\n",
"\treturn $binary;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "608a04f4",
"metadata": {},
"source": [
"This looks like:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "70b5a717",
"metadata": {},
"outputs": [],
"source": [
"my $binary_test_image = create_test_pattern();\n",
"$gp->image($binary_test_image,\n",
"\t{ title => 'Test pattern',\n",
"\t clut => 'grey', },\n",
");\n",
"$gp->fancy_display;"
]
},
{
"cell_type": "markdown",
"id": "74a67a88",
"metadata": {},
"source": [
"Now let's visualise the various parts of this algorithm.\n",
"\n",
"First, lets create a function that let's us see what the Hough space looks like\n",
"and optionally where the peaks are."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b2a13a70",
"metadata": {},
"outputs": [],
"source": [
"sub plot_hough_space($gp, $hough_space, $params, $peak_params=undef, $nhood_size=[15,15], $title='Hough Transform Space [ log(H + 1) ]') {\n",
"\tmy $log_H = log($hough_space + 1);\n",
"\n",
"\t# Base plot remains the same\n",
"\tmy @plot_items = (\n",
"\t\t{\n",
"\t\t\ttitle => $title,\n",
"\t\t\tylabel => encode_utf8('ρ (pixels)'), # rho\n",
"\t\t\txlabel => encode_utf8('θ (radians)'), # theta\n",
"\t\t\tclut => 'sepia',\n",
"\t\t\txtics => {\n",
"\t\t\t\tlabels => [\n",
"\t\t\t\t\tmap { [ encode_utf8($_->[0]), $_->[1] ] }\n",
"\t\t\t\t\t\t['-π/2' , -PI/2],\n",
"\t\t\t\t\t\t['-π/4' , -PI/4],\n",
"\t\t\t\t\t\t[ '0' , 0 ],\n",
"\t\t\t\t\t\t[ 'π/4' , PI/4],\n",
"\t\t\t\t\t\t[ 'π/2' , PI/2],\n",
"\t\t\t\t],\n",
"\t\t\t},\n",
"\t\t},\n",
"\t\twith => 'image',\n",
"\t\t$params->thetas->dummy(0, $params->num_rho),\n",
"\t\t$params->rhos->dummy(1, $params->num_theta),\n",
"\t\t$log_H,\n",
"\t);\n",
"\n",
"\t# Add rectangles around peaks if provided\n",
"\tif (defined $peak_params) {\n",
"\t\tmy $drho = $params->rho_res * $nhood_size->[0]/2;\n",
"\t\tmy $dtheta = $params->theta_res * $nhood_size->[1]/2;\n",
"\n",
"\t\tfor my $i (0..$peak_params->dim(0)-1) {\n",
"\t\t\tmy $theta = $peak_params->slice(\"($i),1\")->sclr;\n",
"\t\t\tmy $rho = $peak_params->slice(\"($i),0\")->sclr;\n",
"\n",
"\t\t\tmy $rect = make_rectangle($theta, $rho, 2*$dtheta, 2*$drho);\n",
"\n",
"\t\t\t$rect->slice('(0)')->inplace->clip($params->theta_range->@*);\n",
"\t\t\t$rect->slice('(1)')->inplace->clip($params->rho_range->@*);\n",
"\n",
"\t\t\tpush @plot_items, (\n",
"\t\t\t\t{ with => 'lines',\n",
"\t\t\t\t lc => 'rgb \"#FFFFFF\"', # white lines\n",
"\t\t\t\t dt => '1',\n",
"\t\t\t\t},\n",
"\t\t\t\t$rect->using(0,1)\n",
"\t\t\t);\n",
"\t\t}\n",
"\t}\n",
"\n",
"\t$gp->plot(@plot_items);\n",
"}\n",
"\n",
"sub make_rectangle($x0, $y0, $width, $height) {\n",
"\t# Create rectangle coordinates (5 points to close the box)\n",
"\treturn pdl([\n",
"\t\t[$x0-$width/2, $y0-$height/2],\n",
"\t\t[$x0+$width/2, $y0-$height/2],\n",
"\t\t[$x0+$width/2, $y0+$height/2],\n",
"\t\t[$x0-$width/2, $y0+$height/2],\n",
"\t\t[$x0-$width/2, $y0-$height/2], # close the rectangle\n",
"\t]);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "b11ea808",
"metadata": {},
"source": [
"And we absolutely need to see the image with the detected lines overlaid."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5013eaf8",
"metadata": {},
"outputs": [],
"source": [
"sub plot_detected_lines($gp, $orig_img, $segments, $title='Line Detection Results') {\n",
"\t$gp->plot(\n",
"\t\t{\n",
"\t\t\ttitle => $title,\n",
"\t\t\tkey => 'off',\n",
"\t\t},\n",
"\t\t{ with => 'image' }, $orig_img,\n",
"\t\t( map {\n",
"\t\t\t (\n",
"\t\t\t { with => 'lines',\n",
"\t\t\t\tlc => 'rgb \"#00FF00\"', # green color\n",
"\t\t\t\tdt => 1, # solid\n",
"\t\t\t\tlw => 2, }, # thicker lines\n",
"\t\t\t\tpdl($_->{point1}[0], $_->{point2}[0]), # x coordinates\n",
"\t\t\t\tpdl($_->{point1}[1], $_->{point2}[1]), # y coordinates\n",
"\t\t\t )\n",
"\t\t } @$segments ),\n",
"\t);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "9344f641",
"metadata": {},
"source": [
"And a helper function to orchestrate the whole process of the Hough transform, peak\n",
"finding, and line segment extraction."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "34a82fa3",
"metadata": {},
"outputs": [],
"source": [
"sub visualise_hough_transform($gp, $orig_img, $binary_edges, $params, %opts) {\n",
"\tmy $max_rho = $params->rho_range->[1];\n",
"\n",
"\tmy %defaults = (\n",
"\t\tthreshold_ratio => 0.3,\n",
"\t\tnhood_size => [15,15],\n",
"\t\tmax_peaks => 20,\n",
"\t\tfillgap => 50,\n",
"\t\tminlength => int($max_rho * 0.2), # 20% of max possible line length\n",
"\t);\n",
"\t%opts = (%defaults, %opts);\n",
"\n",
"\t# Run Hough transform\n",
"\tmy $hough_space = line_hough($binary_edges, $params);\n",
"\n",
"\t# Find peaks\n",
"\tmy ($peak_params, $peak_values) = find_lines($hough_space,\n",
"\t\t\t$params,\n",
"\t\t\t$opts{threshold_ratio},\n",
"\t\t\t$opts{nhood_size},\n",
"\t\t\t$opts{max_peaks});\n",
"\n",
"\tplot_hough_space($gp, $hough_space, $params, $peak_params, $opts{nhood_size});\n",
"\t$gp->fancy_display;\n",
"\n",
"\t# Extract line segments\n",
"\tmy $segments = extract_line_segments($binary_edges, $peak_params,\n",
"\t\t$opts{fillgap}, $opts{minlength});\n",
"\n",
"\t# Visualize detected segments\n",
"\tplot_detected_lines($gp, $orig_img, $segments);\n",
"\t$gp->fancy_display;\n",
"\n",
"\treturn ($hough_space, $peak_params, $segments);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "3b3d281c",
"metadata": {},
"source": [
"## Hough line feature extraction visualised"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b9c3e60f",
"metadata": {
"tags": [
"hide-output-execute_result"
]
},
"outputs": [],
"source": [
"do {\n",
"\t# For test pattern\n",
"\tmy $test_img = create_test_pattern();\n",
"\tmy $test_params = HoughParameters->from_image($test_img);\n",
"\tmy ($test_hough, $test_peaks, $test_segments)\n",
"\t\t= visualise_hough_transform($gp, $test_img, $test_img, $test_params,\n",
"\t\t\tnhood_size => [100,100],\n",
"\t\t\tfillgap => 10,\n",
"\t\t );\n",
"};"
]
},
{
"cell_type": "markdown",
"id": "230abd3f",
"metadata": {},
"source": [
"Here we see that some of the lines in the test image that are nearby each other\n",
"are not considered a separate edge. This is a consequence of multiple parts of\n",
"the algorithm resolution, delta rho tolerance, and neighbourhood suppression.\n",
"Shorter lines are not picked up because the minimum length parameter.\n",
"\n",
"Vertical and horizontal line segments are found, but do not extend the full\n",
"length of the image elements.\n",
"\n",
"---\n",
"\n",
"And finally for the runway image:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4c52288f",
"metadata": {
"tags": [
"hide-output-execute_result"
]
},
"outputs": [],
"source": [
"do {\n",
"\t# For runway image\n",
"\tmy $params = HoughParameters->from_image($binary_edges);\n",
"\tmy ($hough, $peaks, $segments)\n",
"\t\t= visualise_hough_transform($gp, $runway, $binary_edges, $params,\n",
"\t\t\tthreshold_ratio => 0.1,\n",
"\t\t\tnhood_size => [30,30],\n",
"\t\t\tmax_peaks => 50,\n",
"\t\t\tfillgap => 20,\n",
"\t\t\tminlength => 150,\n",
"\t\t );\n",
"};"
]
},
{
"cell_type": "markdown",
"id": "2675b22c",
"metadata": {},
"source": [
"Line segments in both a more vertical direction and a more horizontal direction\n",
"are found. There are some curious line segments that extend the full length of\n",
"the top and bottom of the image. These could be filtered out using a region of\n",
"interest (ROI) filter (i.e., essentially a mask) as it is likely that edges\n",
"near the border of the image do not represent objects.\n",
"\n",
"# What now?\n",
"\n",
"With this, the elves were satisfied with their progress on the autoland\n",
"system…for now. But there are so many more things to explore. PDL let them\n",
"quickly process and visualise an algorithm, making it easy to iterate on their\n",
"design. For the next iteration they could swap out one edge detection algorithm\n",
"for another, come up with custom ways of processing and post-processing the\n",
"image, implement more robust line segment extraction, or use a probabilistic\n",
"approach to quickly extract more complex objects\n",
"([especially those to avoid](https://www.youtube.com/watch?v=MgIwLeASnkw))."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "IPerl",
"language": "perl",
"name": "iperl"
},
"language_info": {
"file_extension": ".pl",
"mimetype": "text/x-perl",
"name": "perl",
"version": "5.38.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
![Runway 34, Nagoya Airfield](\"https://commons.wikimedia.org/wiki/File:Runway_34,_Nagoya_Airfield_(3937428018).jpg\"){kind=link}