A dissolution vug on a static high-resolution borehole image log is a dark blob — a patch of low resistivity where acidic groundwater dissolved the carbonate and left a cavity. The trouble is that the rock around it is not a clean white page. It is a textured, slowly varying field of grey: bedding, cementation contrasts, tool artefacts, the natural intensity gradient of the formation itself. Ask a naive global threshold to "find the dark blobs" and it will either swallow whole intervals of dim matrix as if they were one enormous vug, or, tuned tighter, miss the faint ones entirely. The carbonate background is not noise to be filtered; it is a structured signal that happens to sit on top of the thing you want. In a roughly twenty-month engagement with a mid-sized Middle East carbonate operator we partnered with, the move that made automated vug quantification work was not a bigger neural network. It was a deterministic image-processing step that removes the background by construction: subtract the top intensity modes, then threshold what remains. This piece is about why that works, and why — for this particular problem — a hand-built CV pipeline was the right engineering choice over a learned segmentation mask.
The problem: the matrix is the majority class
Start with the histogram. Take a one-metre patch of a static high-resolution borehole image log, rescale its intensities to the standard 0–255 range, and plot how often each value occurs. The shape is not flat. It is dominated by a tall cluster of values that correspond to the carbonate matrix — the rock that fills most of the borehole wall — with the vugs contributing only a thin tail of darker pixels far to one side. In population terms the matrix is the overwhelming majority class and the vugs are a rare minority, and any segmentation method that does not respect that imbalance will be captured by the majority.
This is precisely why a single global threshold fails. There is no one intensity value that separates "vug" from "not-vug" across a whole well, because the matrix intensity itself drifts with depth, lithology, and tool response. A threshold tuned for one zone clips the next. The conventional fix — local or adaptive thresholding — helps, but on its own it still has to contend with the matrix dominating every local window. You need to get the matrix out of the way first.
Mode subtraction: deleting the background by its own statistics
The insight is that the matrix announces itself in the histogram. The most frequent intensity values — the modes of the distribution — are, almost by definition, the background. So instead of trying to threshold around the matrix, you remove it: identify the top-K most frequent intensity modes in the patch and subtract them out, leaving an image in which the dominant rock fabric has been suppressed and the darker, rarer vug pixels stand in relief.
In the production workflow K was fixed at 5 — the top-five intensity modes (with pixel values on the rescaled 0–255 scale). Each mode is subtracted in turn, producing a small stack of mode-subtracted images (the K, K−1, K−2, K−3, K−4 variants), each of which exposes the vug population against a slightly different residual background. The effect is dramatic in the histogram itself: the most frequent mode in a representative patch carries a frequency of roughly 4,700, while the last selected mode sits near 1,700. Peeling off those high-frequency modes is what flattens the carbonate plateau and lets the dissolution pores survive as resolvable contours.
The second free parameter is the spacing between the modes you pick. Take the five most frequent values literally and you get five nearly identical numbers — consecutive histogram bins differing by fractions of an intensity level — which is useless: each would expose the same pixels and the downstream contour detector would re-trace 90%+ duplicate contours. The fix is a minimum intensity-level separation, δ_m, between selected modes. Setting δ_m = 5 forces the five chosen modes to be genuinely distinct slices of the intensity range, so each mode-subtracted image surfaces a different, complementary slice of the vug population rather than the same blobs five times. Together, K = 5 and δ_m = 5 are the two numbers that define the whole subtraction stage — and, crucially, they were held universal across wells. The same constants applied to a vertical well logged with a high-resolution borehole image log and a horizontal well logged with a compact microresistivity tool kilometres away, which is the property you want from a method you intend to run unattended at scale.
Adaptive thresholding on a flattened field
With the matrix suppressed, thresholding finally has a fair fight. The pipeline applies adaptive (local) thresholding to each mode-subtracted image: rather than one global cutoff, the threshold for each pixel is computed from the mean intensity of a surrounding neighbourhood — a block — with a small offset constant C. On a flattened field this resolves individual vugs cleanly, because the residual matrix is no longer there to bias the local mean.
Block size is the parameter that decides what scale of feature the threshold "sees," and it has a direct physical meaning on an image log. The team swept odd block sizes — 11, 21, 31, 41, 51 pixels — and chose 31 px, because that window covers about 46.58 cm² of borehole wall, comfortably larger than a typical vug but small enough to track the matrix's local variation. The trade-off is concrete and worth internalising: an 11-px block covers only 5.87 cm², too small to establish a stable local background and prone to fragmenting a single vug; a 51-px block covers 126.08 cm², large enough to smear several nearby vugs into one threshold decision. The offset constant C was set to the mean intensity of the current patch by default. These are not magic numbers; they are a window deliberately sized to be bigger than the signal and smaller than the trend.
Turning thresholds into a vug catalog: contours and the circularity gate
A thresholded binary image is not yet a measurement. The next stage extracts contours — the closed outlines of each connected dark region — and treats each contour as a candidate vug. This is where the pipeline stops being generic image processing and starts encoding geology, because not every dark contour is a dissolution pore. Linear features that cut across the image are fractures or bedding; wellbore-parallel streaks are tool artefacts. The discriminating property is shape, and the operative number is circularity.
Vugs are roughly equant — semi-circular to circular — while fractures are elongate. So each contour is scored on a circularity axis and passed through a circularity gate of 0.3–1.0: contours below 0.3 are rejected as fractures or linear artefacts, contours inside the gate enter the vug catalog. Across the validation wells the surviving vug population showed a circularity distribution spanning 0.28–0.85 with a peak in the 0.45–0.7 band — semi-circular dominant, exactly the diagenetic signature you expect — and measured areas of 1–12 cm², concentrated in the 1–3.5 cm² range. The gate is the single geometric decision that separates a vug catalog from an indiscriminate blob count, and the catalog below makes that decision tactile: drag the gate and watch contours fall in and out, including four real vugs that manual picking missed.
Two further filters clean up the catalog before it is counted. Because each of the five mode-subtracted images is thresholded independently, the same physical vug often appears as a contour in more than one image; these duplicates are merged when their centroids fall within a separation threshold (τ = 5) and their bounding boxes overlap past an IoU of 0.2 (20%). The reproducibility of the merge is itself tunable: tightening the duplicate-contour spacing removed roughly 250 redundant contours at the first level, a further ~120 at the second, and stabilised (about 30 removed) by the third — a clean illustration that the cleanup is convergent, not arbitrary. A final contrastive/vicinity step uses a Laplacian-based bounding-box expansion (10%) and a mean-based enlarged radius (2R) to suppress false positives at vug margins. The output, computed every 10 cm down the well, is a full statistical profile per interval: count, total and mean area, an area/porosity spectrum, a circularity spectrum, and an azimuth spectrum — the population a geologist reads by eye, rendered as numbers.
Why deterministic, and why not a learned mask
It is fair to ask why, on an AI engagement, the vug stage is a classical CV pipeline rather than a trained segmentation network. The honest answer is that the problem's constraints favoured determinism, and good applied-ML engineering is knowing when not to reach for the model.
First, supervision is scarce and expensive. Ground truth here is an expert's per-interval vug percentage, not pixel-wise masks. There were not enough hand-segmented vugs to train a segmentation network that would generalise across wells, formations, and two different microresistivity imaging tools. A method whose entire behaviour is governed by a handful of interpretable constants — K, δ_m, block size, circularity gate, IoU — can be specified once and reasoned about, which matters when an interpreter has to defend a number.
Second, the operating point has a physical, conductivity-based gate that is trivial to encode and hard to learn reliably. The pipeline can assert, for instance, that a one-metre zone contains no vugs at all if its minimum pixel value never drops past a conductivity threshold — a deterministic guardrail that prevents the method from hallucinating dissolution porosity in tight rock. Baking that into a learned mask would be an exercise in constraint-injection; in a rule-based pipeline it is one line.
Third, determinism is auditable and fast. Given the same patch and the same constants, the method returns the same contours every time — no stochastic inference, no drift between runs, no GPU. And it is quick: the full pipeline runs at about 15 s/m, against a comparable published morphological approach (Li et al., 2019) that costs roughly 5 minutes per metre — order-of-magnitude faster, which is what turns "interesting demo" into "run it on every well." Against expert picks, the automated profile we shipped reached roughly 85% agreement, with a mean absolute error near 1.21 cm² per interval, and in side-by-side review it actually recovered vugs that the expert's interpretation-software ground truth had missed.
(undefined, undefined) ·The broader lesson for the practitioner is about matching the tool to the structure of the problem. A learned segmentation model is the right answer when you have abundant pixel-level labels and a feature whose appearance varies in ways no fixed rule captures. Vug isolation is the opposite case: the background is statistically self-identifying (it is the histogram's modes), the target's discriminating property is a simple geometric invariant (circularity), and labels are thin. Under those constraints a transparent, deterministic mode-subtraction pipeline is not a fallback — it is the more defensible engineering choice, and the one that survives contact with a petrophysicist who needs to trust the number. It is also a pattern that has travelled well: across the carbonate operators we have worked with in the Middle East and the United States, the right amount of "AI" for a porosity-QC deliverable is sometimes a few well-chosen, legible constants.
Takeaways for the practitioner
Mode subtraction is a small idea with a large payoff: let the background remove itself. The carbonate matrix is the most frequent thing in the histogram, so subtracting the top modes deletes it without ever having to model it, and everything downstream — adaptive thresholding, contour extraction, the circularity gate — becomes a fair fight on a flattened field. The result is a vug detector that is fast, reproducible, auditable, and tuned by numbers a geoscientist can read rather than weights they cannot.
Key takeaways
- The carbonate matrix is the majority class in a high-resolution borehole image log patch's intensity histogram; a single global threshold is captured by it. Subtracting the top-5 intensity modes (K=5) on the rescaled 0-255 range suppresses the background by its own statistics, leaving vugs in relief.
- A minimum mode spacing of delta-m=5 is essential: pick the five literally-most-frequent modes and you get 90%+ duplicate contours; force genuine intensity separation and each mode-subtracted image surfaces a complementary slice of the vug population. First mode frequency ~4,700 vs last selected ~1,700.
- Adaptive thresholding then resolves individual vugs. Block size is sized to the geology: 31 px (46.58 cm²) was chosen over 11 px (5.87 cm², fragments vugs) and 51 px (126.08 cm², smears them), with offset C = patch mean intensity.
- A circularity gate of 0.3-1.0 turns dark contours into a vug catalog, rejecting elongate fractures; the surviving population spans circularity 0.28-0.85 (peak 0.45-0.7) and area 1-12 cm². Cross-mode duplicates merge at IoU 0.2; a conductivity threshold deterministically vetoes vugs in tight zones.
- Determinism was the right call, not a compromise: scarce mask-level labels, an easily-encoded physical guardrail, full auditability, and 15 s/m throughput (vs ~5 min/m for a published baseline) at ~85% agreement and ~1.21 cm² MAE against expert picks — recovering vugs the expert's interpretation-software ground truth missed.
References
[1] Li, X., et al. Automatic characterization of vugs from borehole image logs via path-opening morphology (2019). Used here as the throughput baseline (~5 minutes per metre) against which the mode-subtraction pipeline's 15 s/m is compared.
