Interactive isophote fitter and residual visualiser for galaxy images. Fits concentric elliptical isophotes to a galaxy's smooth light distribution using the Jedrzejewski (1987) algorithm via photutils.isophote, synthesises a model image from the fit, and shows the residual — where the smooth model has been subtracted away — alongside the input in a synced dual-view editor. The residual is where the actually-interesting galaxy structure lives: spiral arms, tidal shells, dust lanes, HII regions, stellar streams, star-forming knots, companion galaxies, boxy or discy isophote distortion, bars, rings. GLIMR adds the whole interactive workflow around photutils' fitter: Shift+click centre picking, Ctrl-drag ellipse adjustment, wedge-mask carve-out, 4× quick preview, adaptive ETA with ring-count cap, three-flavour residual export, and FITS save with every seed and fit parameter embedded.
_SyncedView that mirrors both the transform and both scrollbar positions through a viewChanged signal); the Shift+click centre picker + Ctrl-drag ellipse adjustment for setting the seed geometry visually; the live-updating overlay set (cyan seed ellipse, dashed yellow maximum SMA boundary, dotted magenta minimum SMA boundary); the wedge-mask carve-out with on-canvas wedge preview, used to tell the fitter to ignore an angular range so a nearby star, a jet, or a companion galaxy doesn't bias the smooth-light fit; the input normalisation path that stretches linear data to a 0.25 median before fitting so the Jedrzejewski iteration behaves on linear FITS inputs rather than only on pre-stretched views; the 4× quick-preview path (block-mean downsample, geometry scaling c → c/d, a → a/d, Δ → max(Δ/d, 0.5), nearest-neighbour upsample of the model, full-resolution residual recomputation against the upsampled model, isolist entries rescaled back for the boundary overlay); the adaptive ETA display with ring-count cap (seconds-per-ring profile stored per downsample level, exponential moving average learned from each fit, cross-scale extrapolation τ(d_1) ≈ τ(d_2)·(d_2/d_1)^2 when the other level's profile is missing, soft cap at 10 000 rings that auto-raises the step size when the user picks too many); the symmetric residual preview stretch (zero maps to gray 0.5, with S the 99.5th percentile of |R| measured inside the fit boundary and the preview expressed as D = 0.5 + R/(2S), so positive and negative residual deviations are visually equal); the feathered elliptical blend that shows the residual inside the fit boundary and the original outside it so you can see the galaxy in context; the three-flavour residual push (normalised min→0/max→1 with optional pre-shift so values are non-negative; visible variant that pushes the exact right-pane display including the blend; stretched variant with the symmetric zero-at-gray preview applied and the stretch scale S recorded in metadata); and the FITS save path that embeds every seed and fit parameter (ISO_EPS, ISO_PA, ISO_SMA0, ISO_MIN, ISO_MAX, ISO_STEP, ISO_SCLIP, ISO_NCLIP, ISO_FXC, ISO_FPA, ISO_FEPS, ISO_HARM, ISO_WEDGE, ISO_WPA, ISO_WWID, ORIGMIN, ORIGMAX, PEDESTAL, DSFACTOR) plus the pedestal offset applied when residuals are shifted to ≥ 0 on save so a dropped FITS can be read back and inverted without guessing.
The signature move is the symmetric residual preview stretch with zero at mid-gray:
where R = I - M is the raw signed residual and S is the 99.5th percentile of |R| inside the fit boundary. Zero residual maps to gray 0.5, positive residuals (bright structure the model missed — spiral arms, HII regions, foreground stars) brighten toward white, and negative residuals (dark structure the model overshot — dust lanes, holes in the disc, model over-prediction) darken toward black. For a galaxy with real spiral structure sitting on a smooth bulge, the output is close to a two-tone etching with the arms and dust lanes clearly separated from mid-gray — much more legible than clipping either sign to zero.
GLIMR takes the active mono image as its input. For colour frames, convert to luminance first (or run on individual channels and composite). Linear data (direct from a stack) works but benefits from the normalisation option below; already-stretched views work directly. The dialog is non-modal, so you can keep SASpro's main window usable while GLIMR is open.
The dialog is split: input on the left, residual preview on the right. The two views are synced — zoom with the mouse wheel on either side and the other follows exactly, drag-pan on either and the other follows. This is the one trick that makes isophote fitting actually usable: being able to put your eye on exactly the same patch of the galaxy on both sides and compare input to residual at full zoom, without constantly re-navigating.
The right pane shows the input image until you've run the first fit; after that it shows the residual.
Shift+left-click anywhere in the left pane sets the fit centre to that pixel. This is the fastest way to get a good seed on a galaxy whose centre you can see — just click on it. The centre coordinates are shown as text under the sliders.
Alternatively, drag the cyan ellipse itself: hold Ctrl and drag. The cursor turns to the move cursor while Ctrl is down to confirm you're in drag mode. Without Ctrl the ellipse ignores clicks (so you can't move it accidentally while panning the view).
ε = 1 - b/a where a is the semi-major axis and b the semi-minor. 0 = perfect circle, 0.5 = axis ratio 1:2, 0.9 = highly flattened disc.The cyan ellipse on the left pane updates live as you move either slider, so you can see immediately whether the seed matches the galaxy's apparent shape and tilt. If the seed is good, the Jedrzejewski iteration converges quickly and the residual is clean. If the seed is wrong, the fit still usually converges but you'll see a residual signature of the seed error (quadrupole-like patterns inside the galaxy on an ellipticity error, dipole-like patterns on a centre error).
The semi-major axis range sliders set which radii get fitted:
Three ellipses are drawn on the left pane for visual feedback: the cyan solid ellipse at sma0, a yellow dashed ellipse at max sma, and a magenta dotted ellipse at min sma. All three share the ellipticity and PA.
SASpro caps the ring count at 10 000 by auto-raising the step size if your current range and step would exceed it — the fit will still complete, but at a slightly coarser radial sampling than you asked for. The ETA label will tell you.
σ-clip (sclip) (default 3.0) and σ-clip iters (nclip) (default 1) control the clipping of outlier pixels at each isophote. Pixels more than sclip · σ away from the ring median are excluded, iterated up to nclip times. This matters for the fit on real galaxies: a bright foreground star inside the isophote region will drag the fit if not clipped out. Default 3σ / 1 iteration handles most cases; raise iterations to 2 or 3 if you still see a bright star influencing the result.
Three checkboxes constrain which geometric parameters the Jedrzejewski iteration solves for at each radius:
(x0, y0) fixed at the seed. Useful when you know the centre is accurate (bright nucleus) and don't want the fit wandering it on noisy outer rings.ε fixed. Only sensible when you're certain of the ellipticity (unusual).By default none are fixed, so the fit tracks real variations with radius — useful for diagnosing galaxy structure (isophote twists, radial PA variations).
Add a3/b3/a4/b4 in model (default off): include the third- and fourth-order Fourier harmonics in the synthesised model image. The Jedrzejewski fit always measures these as byproducts — they are the departures from pure ellipticity at each radius and are important diagnostics (positive a4 indicates discy isophotes, negative a4 indicates boxy ones) — but by default the model image is built from only the pure elliptical component, so the residual will contain the a3/b3/a4/b4 structure as part of what you see.
Turning this on folds the higher harmonics into the model, so the residual becomes cleaner but also hides any boxy/discy signature. Turn it on when you want the residual to show only non-elliptical, non-harmonic structure (shells, streams, companions); leave it off when you want to see everything the model couldn't capture.
Check Exclude wedge (deg) to tell the fitter to ignore an angular wedge:
width degrees.A translucent dense-pattern overlay appears on the left pane showing exactly what region the fit will ignore. The excluded pixels are set to numpy.ma.masked before they enter the Jedrzejewski iteration, so they neither contribute to the median at any ring nor to the Fourier residuals. Use cases: a bright foreground star on one side of a galaxy, a bright jet, a nearby companion whose light contaminates one side. The wedge mask lets you fit the smooth galaxy light cleanly from the uncontaminated side only.
Check Normalize before fitting (Linear Data) to apply SASpro's shared midtones-transfer stretch to a target median of 0.25 before fitting. Useful when the input is linear (direct from a stack) — the Jedrzejewski iteration converges much more reliably when the galaxy's brightness distribution is in a reasonable dynamic range than on raw linear data where the dynamic range runs from ~0 to maybe 0.02. Leave this off when the input is already stretched.
Normalisation is applied to the fitting input only. The output model and residual are measured against the normalised version, so if you save the residual back to a view and expect it to be subtractable from the original linear data, that won't work directly — you'd need to invert the stretch, which SASpro doesn't do automatically. For most science use cases (visualisation of structure, published figures) the normalised residual is what you want anyway.
Check Quick preview (4× downsample) (default on) to fit at a 4× downsampled resolution and upsample the model back for the residual. The fit runs 16× faster at the cost of a modest loss in finest-scale accuracy — the output model is a sharp 4×4 pixel staircase at the boundaries of each ring, visible under extreme zoom but not at normal viewing. Perfect for interactive iteration while you're dialing in the seed, PA, ε, and SMA range.
Turn off for the final fit when you want the model at full input resolution. The ring estimate label updates to show the ETA at the current downsample level and tags the result “(quick)” or “(full)” so you know which you're looking at.
Clicks run the fit on a background thread. A modal progress dialog shows the current phase (preparing, building mask, fitting isophotes, building model, upsampling/finalizing) and the ETA based on the learned seconds-per-ring profile. When it's done, the right pane updates to show the residual preview.
If photutils raises an error (empty isolist, degenerate geometry, singular matrix on a bad seed), the error is surfaced in a message box and the dialog stays usable so you can adjust and retry.
After the fit, the right pane shows the residual with the symmetric stretch applied (hero equation). The stretch scale S is the 99.5th percentile of |R| measured inside the fit boundary (the last valid isolist entry), so the normalisation is driven by the galaxy's structure rather than by sky noise or by a hot pixel at the frame corner.
Check Show original outside max ellipse (default on) to blend the original image back into the preview outside the fit boundary. The transition is a 4%-feather elliptical alpha mask, so you see the residual inside the galaxy region and the actual image outside it — keeps context so you can navigate.
Each product (model, residual) has an export row with multiple output options:
Writes a 32-bit float FITS with every seed and fit parameter in the primary header. See §15.14 for the full keyword set.
Pushes a new SASpro document with the array normalised min→0, max→1. For residuals, the Shift residuals to ≥ 0 on save checkbox adds a pedestal first so the min is non-negative before normalisation (useful for pipelines that expect non-negative images; the pedestal is recorded in the metadata and in FITS saves).
Pushes exactly what the right pane is showing — symmetric stretch plus optional outside-ellipse blend, feather and all. Takes the preview buffer directly. Useful for publishing “here's what the fit revealed” figures.
Pushes the symmetric-stretch residual without the outside-ellipse blend. The whole image gets the 0.5 + R/(2S) stretch applied, so the sky noise around the galaxy also sits at mid-gray. Useful when the structure you're looking for might extend beyond the fit boundary (shells, tidal streams).
The isophote model, Jedrzejewski fit, residual preview, and quick-preview downsample internals.
Each isophote is a parametric ellipse at semi-major axis a with centre (c_x, c_y), ellipticity ε, and position angle θ. In rotated coordinates:
a point lies on the isophote when:
The galaxy's smooth light distribution is described by a stack of such isophotes, each at its own radius, with its own intensity I_0(a), centre, ellipticity, and PA. The Jedrzejewski algorithm fits these parameters at each radius independently from the image data.
At each radius, the Jedrzejewski (1987) algorithm samples the image around a trial ellipse and decomposes the azimuthal brightness variation into a Fourier series:
Four low-order coefficients drive geometric updates:
Equivalently, the Jacobian relating the four parameters (c_x, c_y, ε, θ) to the four Fourier residuals is inverted to give the Newton step:
and the iteration continues until the Fourier residuals drop below a convergence threshold. The output is the isolist: one IsophoteFit record per ring with the final (a, I_0, c_x, c_y, ε, θ) plus the a3/b3/a4/b4 higher-order harmonics measured against the final geometry (a3/b3 = pure three-fold distortion, a4/b4 = boxy/discy).
GLIMR configures the integration mode as bilinear and sigma-clipping at the user settings; everything else is photutils' responsibility.
Given an isolist, photutils.isophote.build_ellipse_model renders a 2D image M(x, y) that samples I_0(a) along each pixel's elliptical radius according to the isolist:
Internally this interpolates between adjacent rings' I_0 values at each pixel's elliptical radius under the ring's local geometry (centre, ellipticity, PA). When Add a3/b3/a4/b4 in model is on, the azimuthal variation at each ring is enriched:
so the model includes the boxy/discy character of the isophotes. Off (default), the model is strictly pure elliptical, and the a3…b4 structure shows up as part of the residual.
Trivial but critical:
Signed, floating point. Positive wherever the image is brighter than the smooth model (missed structure); negative wherever the image is darker (dust lanes, over-fit). GLIMR carries R as float32 throughout; stretching only happens at display time.
For display, the signed residual is linearly stretched so R = 0 maps to mid-gray and the strongest ±residuals extend to black/white. The stretch scale S is the 99.5th percentile of |R| measured only inside the fit boundary:
where ℰ_max is the last valid isolist ellipse. Measuring S only inside the fit boundary means a hot pixel at the frame corner can't drive the normalisation, and the galaxy's residual structure gets the full dynamic range. The display map is the hero equation:
Why the 99.5th percentile rather than max: a single bright foreground star inside the fit boundary that the mask didn't catch would saturate the stretch and make the galaxy's actual residual structure invisible. The 99.5 cut keeps the top 0.5% of residual pixels clipped (which is where that star lives) while the rest of the structure maps linearly across the dynamic range.
When Show original outside max ellipse is on, GLIMR blends the symmetric-stretched residual back into the normalised input using an elliptical feather. The elliptical radius at a pixel is:
with β = 0.04 (4% feather width in units of normalised elliptical radius). Inside the ellipse m < 1 so α = 1; across the 4% band around the boundary α ramps smoothly to 0; outside α = 0. The final preview is:
No sharp transition, no visible seam between the residual region and the context region.
The wedge mask carves out an angular sector:
where φ = atan2(y - c_y, x - c_x) is the angle from the centre to the pixel and φ_0, w are the wedge's central PA and full width. The arctan2(sin, cos) construction wraps the angular difference to [-π, π] correctly without a discontinuity at 180°. Masked pixels become numpy.ma.masked entries in the fitter's input array, so they're excluded from the Fourier sum at every ring.
4× downsample runs the fit on an image that's 16× smaller in pixel count. Block-mean downsample:
Geometry is scaled consistently:
The max(·, 0.5) floor on step size prevents the ring count from blowing up at the downsample level (a sub-half-pixel step on the small image is wasted; the small image can't resolve rings closer than that). The fit runs on the small image. Then the model is upsampled with nearest-neighbour to the input resolution:
and the residual is computed at full resolution against the upsampled model (not against the downsampled image — so the residual still shows full-resolution structure, just measured against a slightly staircased smooth model). The isolist entries are rescaled back to full-res coordinates for the boundary overlays. Overall speedup: close to 16× on the fit itself.
Ring count estimate:
Soft cap at 10 000 rings: if the user's settings would exceed this, the step size is bumped up:
Per-downsample-level seconds-per-ring profile τ(d) is stored in a dict with EMA update after each fit:
So the ETA improves with each fit at that downsample level. When the current level's profile is missing (first time running at a new downsample), SASpro extrapolates from the other level using the roughly quadratic pixel-count scaling:
so if you've only ever run at 4× downsample and then switch off quick preview, the ETA for the full-res fit is seeded as roughly 16× the quick-preview per-ring time. Rough but correct order-of-magnitude. The ETA displayed on the ring label is τ · N_rings plus a “(quick)” / “(full)” tag.
The fit runs in a QThread-wrapped _FitWorker (QObject). Progress is reported through a progress(int, str) signal at five phases: Preparing, Building mask, Fitting isophotes, Building model, Upsampling / finalizing. The main thread handles only UI updates — the dialog stays responsive, the progress bar advances, you can see the status label changing. On completion, the model, residual, and isolist arrive back on the UI thread via a finished signal; on error, a error(str) signal routes to a message box. Worker and thread are both deleteLater'd after completion so repeated fits don't leak.
When Normalize before fitting (Linear Data) is on, GLIMR runs stretch_mono_image(x, target_median=0.25) before everything else — same midtones-transfer stretch as Autostretch and GHS, same code path. The normalised image is cached as _in01 and used both for display (left pane) and as the input to the fitter. The output model and residual are therefore in the normalised domain. Off: the raw image is clipped to [0, 1] and used directly.
The outside-ellipse context is the normalised input (_in01), not the raw image. This keeps the brightness levels in the context region matched to the brightness scale of the residual stretch, which keeps the transition across the feather smooth. Mixing a raw-domain outside with a stretched-domain inside would print a visible brightness jump at the feather boundary.
The three residual push buttons emit different data arrays with different metadata:
| Variant | Data | Metadata highlights |
|---|---|---|
| Normalized | (R - min) / (max - min), with optional pedestal if Shift residuals to ≥ 0 on save is on | push_variant: "normalized_01", source_range_min, source_range_max |
| Visible | Exact right-pane preview buffer | push_variant: "visible_preview", preview_blend, feather_frac: 0.04 |
| Stretched | Symmetric stretch applied to the whole image (no outside blend) | push_variant: "preview_stretch", stretch_pct: 99.5, stretch_scale_S, zero_maps_to_gray_0p5: true |
All three variants include the full isophote-parameter dict in metadata (isophote_params key), so the pushed document carries a complete record of how it was produced.
FITS saves embed the fit parameters as header keywords. Full set:
| Keyword | Meaning |
|---|---|
CREATOR | “GLIMR” |
PRODUCT | “model” or “resid” |
ORIGMIN, ORIGMAX | min/max of the data before any pedestal is added |
PEDESTAL | pedestal offset added if the residual was shifted to ≥ 0 at save time |
DSFACTOR | downsample factor used for the fit (then upsampled) |
ISO_EPS, ISO_PA, ISO_SMA0 | seed geometry (ellipticity, position angle, initial SMA) |
ISO_MIN, ISO_MAX, ISO_STEP | SMA range and step |
ISO_SCLIP, ISO_NCLIP | sigma-clip settings |
ISO_FXC, ISO_FPA, ISO_FEPS | fix-flags for centre, PA, ellipticity |
ISO_HARM | true if high harmonics were folded into the model |
ISO_WEDGE, ISO_WPA, ISO_WWID | wedge mask (enabled flag, central PA, full width) |
Reading a saved residual back and inverting the pedestal to recover signed residuals: R = data - PEDESTAL. The ORIGMIN / ORIGMAX keywords let you confirm the original range so you can verify the inversion.
0.5 + R/(2S) zero-at-gray display map; the feathered elliptical blend for outside-boundary context; the three-flavour residual push (normalised, visible, stretched) with full metadata capture; the FITS save with every seed/fit parameter embedded and the pedestal-offset round-trip support; the adaptive ETA with per-downsample seconds-per-ring EMA, cross-scale extrapolation, and the 10 000-ring soft cap that auto-raises step size; and the background-threaded fit worker with five-phase progress reporting.
(A_1, B_1, A_2, B_2) are inverted to update the four geometric parameters (c_x, c_y, ε, θ) until convergence. The higher-order harmonics (a_3, b_3, a_4, b_4) are measured as byproducts and are standard galaxy-morphology diagnostics: a_4 > 0 indicates discy isophotes, a_4 < 0 indicates boxy ones, following Bender & Möllenhoff (1987, A&A 177:71). Model synthesis is photutils.isophote.build_ellipse_model, which interpolates the isolist's I_0(a) and (optionally) high-harmonic Fourier coefficients to produce the model image.
stretch_mono_image, same code behind Autostretch and GHS) used by the input normalisation path.
photutils.isophote.Ellipse, EllipseGeometry, build_ellipse_model). NumPy provides array operations, masked arrays (numpy.ma.masked_array for the wedge mask), and linear-algebra helpers for the internal auto-estimate (eigendecomposition of the 2×2 intensity-weighted covariance matrix via linalg.eigh). Astropy provides FITS I/O (astropy.io.fits) for save. PyQt6 provides the UI: QGraphicsView / QGraphicsScene for the dual synced views, the custom DraggableEllipse item for the Ctrl-drag seed, the _SyncedView peer-mirroring, the QThread+QObject worker pattern, and the QProgressDialog for fit progress.