How Hebog finds sources

This page is for astronomers who want to know what Hebog does to an image and why. It follows one call to hebog.find_sources() from a FITS image to the published catalogue. For field definitions and units, see the output reference. For how the same steps are spread over many machines, see How Hebog distributes work.

Overview

flowchart TD
    image[/"FITS image (Jy/beam)"/]
    subgraph noise["1 · Background and noise"]
        bg["Sigma-clipped statistics on a coarse grid<br/>→ background and RMS maps"]
        fine["Finer RMS grid around bright sources,<br/>with source pixels protected"]
        bg --> fine
    end
    subgraph detect["2 · Detection"]
        snr["Residual = image − background<br/>S/N = residual / RMS"]
        filters["Beam-matched filters and<br/>B3 à trous wavelet scales 1–3"]
        flood["Seed at the detection threshold,<br/>grow to the island threshold"]
        snr --> filters --> flood
    end
    subgraph character["3 · Characterisation"]
        deblend["Deblend: split islands at significant<br/>peaks separated by a saddle"]
        fit["Fit elliptical Gaussians to<br/>background-subtracted pixels"]
        associate["Associate components into sources<br/>using cross-scale evidence"]
        measure["Sum the fitted components; also measure<br/>each source in its own non-overlapping aperture"]
        deblend --> fit --> associate --> measure
    end
    products[/"catalogue.fits · rms.fits<br/>source-mask.fits · diagnostics.json"/]

    image --> noise --> detect --> character --> products

The approach belongs to the same family as PyBDSF and Aegean: estimate local noise, threshold into islands, and fit Gaussians. It differs mainly in how it treats extended emission, in reporting an aperture flux alongside the summed Gaussian flux, and in running every step on independent tiles. See the comparison with other source finders.

Three populations, not one

Hebog keeps three things apart that are easy to conflate:

Term	Meaning	Catalogue table
Island	A connected footprint in the published detection mask	ISLANDS
Gaussian component	One successfully fitted elliptical Gaussian	GAUSSIAN_COMPONENTS
Source	One or more detected components that Hebog associates as a single object in the image	SOURCES

An island can hold several independent sources. A source can span several islands. A detected component has no Gaussian row if its fit failed or was not admissible. None of the three is automatically one astrophysical object.

Step by step

1. Check the image

Hebog needs pixel values in Jy/beam, an ICRS or FK5 J2000 celestial WCS, a restoring beam and a reference frequency. NaN pixels are allowed and are excluded everywhere, as is every pixel of a 3×3 square of one repeated value, so zero padding is excluded to its last pixel (see invalid pixels). Anything else is rejected before analysis, with an error that says what is missing. Capability and status lists the exact requirements, including the current 22,500-pixel size limit.

2. Estimate background and noise

The background (slowly varying offset) and the RMS (local noise) are estimated separately:

• Iteratively sigma-clipped statistics are computed in windows on a coarse grid, then interpolated to every pixel.
• Pixels that belong to candidate sources are protected: they are excluded from the statistics, so bright or extended emission does not inflate the noise estimate around itself.
• Near bright sources, where noise changes quickly because of imaging artefacts, both profiles switch to a finer RMS grid, whose windows the continuum profile also protects from sources. That bright-region grid only raises the noise: where it reads lower than the coarse grid, the coarse RMS stands, much as PyBDSF applies its small box near a bright source only out to where it stops reading well above its large one. Bright sources close together share one region, which in a crowded deep field can cover the whole image; its fine grid is estimated in blocks of bounded size, so the memory one task needs does not grow with the region. Emission cut off at a block's edge is protected in case it belongs to a bright source beyond, so a block can exclude slightly more than one read of the whole region would.
• At image edges, fine RMS values are extended as constants, never extrapolated towards zero. Background and coarse-RMS slopes are taken from grid samples at least as far apart as the distance being extrapolated, so a genuine gradient is preserved without amplifying small errors. A coarse RMS that falls towards the edge is held at the edge cell's value, so it never goes below the edge cell; one that rises, as noise does towards the edge of a primary-beam-corrected image, is extended.

In the reviewed continuum profile the coarse grid uses 150-pixel windows every 50 pixels and the fine RMS grid 35-pixel windows every 7 pixels. Windows that overlap protected sources are dropped and each gap takes the value of its nearest clean window. A window whose clipped spread is no greater than single precision's resolution at its brightest valid value (2⁻²³ of that value), as a window of samples that all hold one value is, measures no noise and is dropped the same way: the published RMS is never zero or rounding.

On images at least 150 pixels on a side, the continuum profile also measures the fine grid across the whole image (local noise), and wherever one of its windows came clean, its RMS is the one published. Three rules decide what a dropped window publishes there:

• A gap with a clean window within half a coarse window never reads below 0.8 of the coarse RMS at its centre, or of the largest such clean window if that is lower. Beside a sharp step in the noise the nearest clean window can lie on the quieter side, and without this floor it would lower the noise beside the step.
• A gap beyond every clean window's reach, as inside a source wider than a coarse window, keeps its nearest clean window's RMS, because the unprotected coarse window there holds the source's emission: on images whose shorter side is 600 pixels or more the coarse grid is not protected, so over extended emission its RMS keeps emission that no clean window nearby reads.
• A window that was measured keeps its own value, so a quieter region still reads its own noise. Within about half a fine window of a step, windows that straddle it read a mixture of both levels, as any windowed estimate does.

Around bright sources a gap takes a clean window's RMS only within one fine window (35 pixels) of it; farther away the coarse RMS stands, so in a crowded field, where few fine windows are clean, one window's noise is not carried across the image or across a step in the noise. Beside such a step in a crowded field the published noise is then the coarse estimate, which mixes the two sides over about half a coarse window. In a field so crowded that clean fine windows lie within about half a coarse window of fewer than half its grid cells, filling the rest from them would copy a few distant windows across the field. There a clean window's RMS stands only within that reach; beyond it the coarse RMS stands, as PyBDSF's does, unless it reads more than three times the nearest clean window's noise, as over a source wider than a coarse window, where that noise stands instead. A field so crowded that no fine window anywhere is clean is too crowded to protect: it keeps the unprotected sigma-clipped coarse statistics, which, like PyBDSF's, include the sources' wings. The bright-region grid still refines the background around bright sources there, but not the RMS: its windows protect only the bright source's own emission, so in such a field they would read a neighbouring extended source as noise and cut into its side nearest the bright source.

Background refinement around bright sources is triggered at 75σ, or at your detection threshold if your island threshold is higher than that. This never changes your detection or island thresholds. Noise structure finer than the grid is not measured.

If no pixel has a finite positive RMS, a sigma threshold has no meaning. Hebog then returns an empty catalogue, a zero mask and an all-NaN RMS image marked unavailable. This is not evidence of an empty sky. An image with too few valid pixels for any coarse window to measure is a cause, as is a noiseless simulation whose windows all fall under the noise floor, as a window that holds a source and its tails does. A noiseless simulation is not reliably so: a window of a source's far tails alone has a spread close to their own values, which is measured as a positive but tiny noise, so add noise to a simulation.

3. Detect compact and extended emission

Hebog forms the residual image − background and the signal-to-noise ratio residual / RMS. Detection then uses two thresholds, as PyBDSF and Aegean do:

• a detection threshold (for example 5σ): an island must contain evidence at this level; and
• a lower island threshold (for example 3σ): an accepted island grows over eight-connected pixels down to this level.

To find emission that is faint per pixel but significant over a larger area, Hebog also evaluates the residual through a bank of beam-matched filters and a B3-spline à trous wavelet transform. These filtered images only help decide whether a region is significant. Islands still grow on the original residual pixels, and all fluxes are measured on the original image.

Two area rules apply. A region promoted only by a filter must cover a minimum area in beams. A region with a pixel directly above the detection threshold can be smaller. Your minimum_island_pixels and optional maximum_island_pixels are applied afterwards.

Diffuse multiscale emission can enlarge the region used to measure a nearby detection, but it cannot merge two detections merely because their faint wings touch. Each pixel of it goes to the detection nearest along the emission itself, within half a beam, so a detection's region is connected through its own pixels; a pixel no detection reaches that way is left out. Multiscale support is kept only when it persists at an adjacent wavelet scale.

4. Deblend and fit Gaussians

Within each island Hebog looks for significant peaks and splits the island into one region per peak. Two peaks stay apart when the fainter lies at least 1.5σ above the pass between them, the highest level at which one connected part of the island still holds both, so a faint source beside a much brighter one keeps its own region. Each region holds the pixels that rise to its own peak.

Each component is then fitted with an elliptical Gaussian, initialised from image moments and fitted to the original background-subtracted pixels. Neighbouring components whose fitting regions touch are fitted jointly. A fit is admitted only if it converged, stayed within physical bounds, is well conditioned, leaves acceptable residuals and has usable uncertainties. Depending on the data, the admitted model is a free ellipse or a beam-shaped Gaussian. When one or more components of a joint fit cannot be fitted with a free shape, for instance because their own pixels lie along a line, those components take the beam shape and the rest are refitted free; any of these then not significantly extended take the beam too.

Fitting work is bounded. In a crowded field the touching fitting regions can chain across the whole image. A chain with more components or pixels than one joint fit admits (at most sixteen components), or wider than the bounded read window, is fitted island by island instead, as PyBDSF fits every island; a component's fit then does not model a neighbouring island's faint wings. An island too large or too complex for the bounded deblend or joint fit on its own is kept as a detection and recorded as deferred in the diagnostics. It is never silently dropped.

5. Associate and measure sources

In the continuum profile, Hebog decides which components form one source using a hierarchy of wavelet-scale features, the compact Gaussian models and the emission left after subtracting them. Every merge is recorded with its evidence in diagnostics.json. Components without positive evidence stay independent.

A fit's Gaussian models describe its components when the emission left after subtracting them holds no residual feature of its own. As in PyBDSF's search of its residual image, such a feature must be positive, grow from a detection-threshold seed on the original pixels or one of the wavelet or matched-filter scales, and touch the detected support of the fitted components. It counts against only the groups of overlapping Gaussians whose own support it touches, not every component fitted jointly with them.

Residual emission the components share can join them into one source. That join is stricter for a component whose model describes it: the emission must lie beneath its own fitted core, the half-maximum ellipse of its Gaussian. Overlapping Gaussians and resolved arcs can join it too. Resolved arcs form a loop only where they lie on the rim of the hole they enclose.

Each source is measured twice. Its catalogue INTEGRATED_FLUX is the sum of its fitted Gaussian components, which is how PyBDSF defines a source's total flux. Each source also owns a non-overlapping aperture, and the signed sum of background-subtracted pixels in it is published as ASSOCIATION_APERTURE_FLUX; a source with no admitted fit uses this aperture flux as its INTEGRATED_FLUX and is flagged. The two agree for isolated compact sources but part for extended or edge-clipped emission, because a fit integrates sky the image does not cover and misses diffuse emission no component describes. If the aperture sum is not positive, as it can be where neighbouring sources' wings raise the background, a source with an admitted fit keeps its summed flux and leaves ASSOCIATION_APERTURE_FLUX empty, while a source without one gets no catalogue row: Hebog does not substitute a positive-only estimate. A source whose aperture or position cannot be measured at all gets no row either. A source of one fitted component is published as that Gaussian: its position, position errors, peak and shapes are the fit's, as PyBDSF publishes a single-Gaussian source. Any other source's position comes from the detection footprint, not from faint measurement-only wings, so its centroid can lie between two peaks or inside a ring.

In the compact profile, association is skipped: every fitted component is its own source and carries its Gaussian measurement. Diagnostics then declare extended-emission-incomplete.

6. Publish products

The mask is the retained detection footprint, and its connected regions are the catalogue's islands. Publication smooths each detection's noisy edge with a 3×3 opening and keeps sparse edge pixels only at 6σ, but it never removes an accepted detection: one the smoothing would remove entirely, such as a 5–6σ point source in a beam only a few pixels wide, keeps its whole thresholded footprint. Multiscale emission is published with a detection only where it joins that detection's footprint, so no detection is published in parts. A component whose own pixels fall outside the mask, such as a faint companion on a brighter source's rim, belongs to the islands of the detection it was split from. A Gaussian row is published only with its parent source row. All four files are written to a private directory, validated, and moved into place together, so a partial result is never visible.

Hebog does not publish its background map, residual or model images, wavelet planes, or a filtered sky model.

Detailed decision flow

The diagram below shows every accept or reject decision in the current finder. Input and configuration failures stop before any scientific work. A scientific rejection, such as an unseeded region or a rejected fit, does not stop the run: the item is left out of the catalogue and, where it had an identity, recorded in the diagnostics.

flowchart TD
    start([SourceFinderRequest, configuration, executor])
    config{Configuration valid?}
    destination{Output directory absent?}
    input{FITS readable and metadata supported?}
    reject_config([Reject invalid thresholds or profile])
    reject_destination([Refuse to overwrite caller-owned output])
    reject_input([Reject invalid unit, WCS, beam, frequency, shape, or size])

    start --> config
    config -- No --> reject_config
    config -- Yes --> destination
    destination -- No --> reject_destination
    destination -- Yes --> input
    input -- No --> reject_input
    input -- Yes --> partition[Plan deterministic background and detection tile cores]

    profile{Continuum or compact profile?}
    continuum[Source-protected background and adaptive local RMS]
    compact[Compact-profile background and RMS policy]
    rms[Reconcile tile estimates into background and RMS planes]
    usable{Any finite, positive RMS?}
    empty[Create empty catalogue, zero mask, and all-NaN unavailable RMS]

    partition --> profile
    profile -- Continuum --> continuum
    profile -- Compact --> compact
    continuum --> rms
    compact --> rms
    rms --> usable
    usable -- No --> empty

    residual[Subtract background; exclude invalid pixels; divide by local RMS]
    filters[Evaluate beam-aware matched filters and residual B3 à trous scales per tile]
    flood[Grow eight-connected candidate islands on original residual pixels at island threshold]
    seed{Candidate has detection-threshold evidence?}
    drop_unseeded[Discard unseeded region]
    area{Meets multiscale area rule or has a direct detection-threshold pixel?}
    drop_area[Discard unsupported sub-area region]
    own[Keep direct labels; assign significant multiscale support within half a beam along it to the nearest seed owner]
    boundary[Refine boundaries with original-pixel S/N and adjacent-scale persistence; preserve owner bridges]
    configured_size{Within caller minimum and optional maximum pixel count?}
    drop_size[Discard component]
    any_component{Any component remains?}

    usable -- Yes --> residual
    residual --> filters
    filters --> flood
    flood --> seed
    seed -- No --> drop_unseeded
    seed -- Yes --> area
    area -- No --> drop_area
    area -- Yes --> own
    own --> boundary
    boundary --> configured_size
    configured_size -- No --> drop_size
    configured_size -- Yes --> any_component
    any_component -- No --> empty

    deblend{Parent is eligible for bounded deblending and has separable peaks?}
    split[Split into deterministic component owners]
    preserve[Keep one owner or record bounded-work deferral]
    fit{Gaussian model is measurable and scientifically admissible?}
    gaussian[Retain admitted Gaussian-component measurement]
    fit_absent[Keep detection identity; record unavailable or deferred fit]
    hierarchy[Build multiscale component hierarchy and association evidence]
    association{Profile allows source association?}
    associate[Choose independent, compact-model, or extended-morphology grouping]
    singleton[Use one source per component; declare extended-emission limitation]
    source_measure[Measure each source once in its signed, non-overlapping owned aperture; sum its fitted components]
    source_valid{Position available, an admitted fit or a positive aperture, and publication support?}
    compact_valid{Admitted component measurement with publication support?}
    source_row[Publish source row]
    source_absent[Keep source disposition without a catalogue row]

    any_component -- Yes --> deblend
    deblend -- Yes --> split
    deblend -- No --> preserve
    split --> fit
    preserve --> fit
    fit -- Yes --> gaussian
    fit -- No --> fit_absent
    gaussian --> hierarchy
    fit_absent --> hierarchy
    hierarchy --> association
    association -- Continuum --> associate
    association -- Compact --> singleton
    associate --> source_measure
    singleton --> compact_valid
    source_measure --> source_valid
    source_valid -- Yes --> source_row
    source_valid -- No --> source_absent
    compact_valid -- Yes --> source_row
    compact_valid -- No --> source_absent

    assemble[Build mask-connected islands and link sources and fitted components]
    bundle[Write and validate catalogue, RMS, mask, and diagnostics in a private bundle]
    publish{All four products valid?}
    success([Atomically rename complete bundle and return SourceFinderResult])
    fail([Leave requested output absent; caller may retry])

    source_row --> assemble
    source_absent --> assemble
    empty --> bundle
    assemble --> bundle
    bundle --> publish
    publish -- Yes --> success
    publish -- No --> fail

Reading the results safely

• Use the catalogue and diagnostics.json together. The diagnostics list every detected component and source, including those without a catalogue row.
• Do not infer a Gaussian fit from a source row unless it carries original-pixel-gaussian-model, as a source of one fitted component does, and do not infer a blank sky from an empty catalogue.
• To compare with a PyBDSF or Aegean component list, use GAUSSIAN_COMPONENTS, not SOURCES.
• A successful run is not scientific qualification; see capability and status.