Computationally efficient algorithms for statistical image processing. Implementation in R

In the series of our earlier papers on the subject, we proposed a novel statistical hypothesis testing method for detection of objects in noisy images. The method uses results from percolation theory and random graph theory. We developed algorithms that allowed to detect objects of unknown shapes in the presence of nonparametric noise of unknown level and of unknown distribution. No boundary shape constraints were imposed on the objects, only a weak bulk condition for the object’s interior was required. Our algorithms have linear complexity and exponential accuracy. In the present paper, we describe an implementation of our nonparametric hypothesis testing method. We provide a program that can be used for statistical experiments in image processing. This program is written in the statistical programming language R.

💡 Research Summary

The paper presents a concrete implementation in the statistical programming language R of a non‑parametric hypothesis‑testing framework for detecting objects in noisy digital images. The underlying methodology, introduced in a series of earlier works by the authors, leverages percolation theory and random graph theory to treat an image as a lattice graph whose vertices correspond to pixels. Each pixel is binarized (0 = background, 1 = potential object) and the set of active vertices forms clusters according to a chosen neighborhood structure (4‑, 6‑ or 8‑neighbor connectivity).

The central statistical problem is cast as a hypothesis test: under the null hypothesis the image contains only noise, so the distribution of the size of the largest percolation cluster can be derived from the theory of site percolation. If the observed maximal cluster exceeds a critical value obtained from this null distribution, the null is rejected and an object is declared present. Crucially, the method does not require any prior knowledge of the object’s shape, smoothness, or even the noise distribution; only a weak bulk condition on the object interior is assumed. This makes the approach especially suitable for highly irregular, non‑convex, or even disconnected objects that appear in medical imaging, remote sensing, materials science, and urban analysis.

From an algorithmic perspective the paper details three main components:

1. Image conversion and neighbor‑list construction – Using the rimage package, JPEG or other raster formats are read and converted into a vector of gray‑scale values (restricted to 0/1 for percolation). Helper functions (img2mat, mat2img, vec2mat, mat2vec) enable seamless switching between matrix and vector representations. Neighbor lists for rectangular grids are pre‑computed once via the functions adj4mat, adj6mat, and adj8mat. Each list stores, for every pixel index, the indices of its neighboring pixels according to the chosen connectivity. This preprocessing runs in O(N) time where N is the total number of pixels and requires only O(N) memory.

2. Depth‑First Search (DFS) for cluster extraction – The authors implement Tarjan’s DFS algorithm (1972) in the R function spanning_tree. Starting from a seed pixel, the algorithm explores all reachable active pixels, assigning a sequential label (rank) that also encodes the distance from the seed in the spanning tree. A “mother” vector records the parent of each visited vertex, enabling back‑tracking when no further unlabelled neighbors exist. The DFS proceeds until the entire connected component is labelled, then repeats with a new seed to discover all clusters. Because each pixel is visited at most once, the algorithm achieves linear time complexity O(N) and linear space usage.

3. Newman‑Ziff algorithm for null‑distribution simulation – To obtain the distribution of the maximal cluster size M under the null hypothesis, the paper adopts the Newman‑Ziff algorithm (2000). The algorithm treats the number of active sites N as a binomial random variable with parameters (S, p), where S is the total number of sites and p the activation probability. Conditional on a fixed N = n, the algorithm simulates the random permutation of sites, activates the first n sites, and records the size Xₙ of the largest cluster formed. Repeating this for all n = 0,…,S yields the conditional probabilities pₘₙ = P(M = m | N = n). The unconditional distribution follows from the law of total probability: P(M = m) = ∑ₙ pₘₙ · Binom(S, n; p). The authors provide both a pure R implementation and a faster C++ version (called from R) to handle larger lattices.

The paper emphasizes that the entire pipeline—image loading, neighbor list generation, DFS clustering, and Newman‑Ziff simulation—has linear computational complexity with respect to the number of pixels, while the statistical power grows exponentially with the object’s pixel count. This “exponential accuracy” claim is supported by theoretical results proved in the earlier works and by empirical experiments reported in the current manuscript.

A notable practical feature is the built‑in data‑driven stopping rule: the DFS and Newman‑Ziff procedures automatically terminate when no further progress is possible, eliminating the need for user‑specified iteration limits. Consequently, the method can be regarded as an unsupervised learning technique within the broader machine‑learning taxonomy.

The authors also discuss performance considerations. Pure R code, while convenient for prototyping and statistical experimentation, suffers from interpreter overhead and becomes slow for images larger than a few hundred thousand pixels. The C++ implementation dramatically reduces runtime, suggesting that production‑level applications should rely on the compiled version while retaining R for orchestration, visualization, and statistical analysis.

In terms of applicability, the method shines when objects have unknown, highly irregular shapes and when the noise may be heavy‑tailed or even lack a density function. Because the hypothesis test only depends on the percolation threshold and not on parametric assumptions about the noise, it is robust to a wide variety of real‑world imaging conditions. Potential domains include cancer detection in histopathology slides, crack detection in civil‑engineering materials, urban feature extraction from satellite imagery, and detection of faint astronomical sources.

The paper concludes by outlining future directions: (i) extending the implementation to GPU‑accelerated or parallel C++ code for massive images, (ii) incorporating spatially correlated noise models, (iii) adapting the framework to multi‑scale or multi‑channel data (e.g., color or hyperspectral images), and (iv) integrating the method into larger machine‑learning pipelines for downstream tasks such as segmentation or classification.

Overall, the manuscript delivers a thorough translation of sophisticated percolation‑based statistical theory into a usable R package, demonstrates its algorithmic efficiency, and provides concrete code that can be directly employed for experimental studies in statistical image analysis.