Coverage for python / lsst / pipe / tasks / prettyPictureMaker / _localContrast.py: 9%
142 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-01 02:17 -0800
1# This file is part of pipe_tasks.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22from __future__ import annotations
24__all__ = ("localContrast",)
26import numpy as np
27from numpy.typing import NDArray
28import cv2
29from numba import njit, prange
30from numba.typed.typedlist import List
31from collections.abc import Sequence
32from itertools import cycle
35@njit(fastmath=True, parallel=True, error_model="numpy", nogil=True)
36def r(
37 img: NDArray, out: NDArray, g: float, sigma: float, shadows: float, highlights: float, clarity: float
38) -> NDArray:
39 """
40 Apply a post-processing effect to an image using the specified parameters.
42 Parameters:
43 img : `NDArray`
44 The input image array of shape (n_images, height, width).
45 out : `NDArray`
46 The output image array where the result will be stored. Should have the same shape as `img`.
47 g : `float`
48 A parameter for gamma correction.
49 sigma : `float`
50 Parameter that defines the scale at which a change should be considered an edge.
51 shadows : `float`
52 Shadow adjustment factor.
53 highlights `float`
54 Highlight adjustment factor. Negative values INCREASE highlights.
55 clarity : `float`
56 Clarity adjustment factor.
58 Returns:
59 result : `NDArray`
60 The processed image array with the same shape as `out`.
61 """
63 h_s = (highlights, shadows)
65 # Iterate over each pixel in the image
66 for i in prange(out.shape[0]):
67 # Get the current image slice
68 imgI = img[i]
69 # Get the corresponding output slice
70 outI = out[i]
72 # Iterate over each pixel in the image
73 for j in prange(out.shape[1]):
74 # Calculate the contrast adjusted by gamma correction
75 c = imgI[j] - g
76 # Determine the sign of the contrast adjustment
77 s = np.sign(c)
79 # Compute the transformation term t based on the signed contrast
80 t = s * c / (2.0 * sigma)
81 # Clamp t to be within [0, 1]
82 t = max(0, min(t, 1))
84 t2 = t * t
85 # Complement of t
86 mt = 1.0 - t
88 # Determine the index based on the sign of c (either 0 or 1)
89 index = np.uint8(np.bool_(1 + s))
91 # Compute the final pixel value using the transformation and
92 # additional terms for shadows/highlights and clarity
93 val = g + s * sigma * 2 * mt * t + t2 * (s * sigma + s * sigma * h_s[index])
94 val = val + clarity * c * np.exp(-(c * c) / (2.0 * sigma * sigma / 3.0))
96 # Assign the computed value to the output image
97 outI[j] = val
99 return out
102def makeGaussianPyramid(
103 img: NDArray, padY: list[int], padX: list[int], out: List[NDArray] | None
104) -> Sequence[NDArray]:
105 """
106 Create a Gaussian Pyramid from an input image.
108 Parameters:
109 img : `NDArray`
110 The input image, which will be processed to create the pyramid.
111 padY : `list` of `int`
112 List containing padding sizes along the Y-axis for each level of the pyramid.
113 padX : `list` of `int`
114 List containing padding sizes along the X-axis for each level of the pyramid.
115 out `numba.typed.typedlist.List` of `NDarray` or `None`
116 Optional list to store the output images of the pyramid levels.
117 If None, a new list is created.
119 Returns:
120 pyramid : `Sequence` of `NDArray`
121 A sequence of images representing the Gaussian Pyramid.
123 Notes:
124 - The function creates a padded version of the input image and then
125 reduces its size using `cv2.pyrDown` to generate each level of the
126 pyramid.
127 - If 'out' is provided, it will be used to store the pyramid levels;
128 otherwise, a new list is dynamically created.
129 - Padding is applied only if specified by non-zero values in `padY` and
130 `padX`.
131 """
132 # Initialize the output pyramid list if not provided
133 if out is None:
134 pyramid = List()
135 else:
136 pyramid = out
138 # Apply padding only if needed, ensuring the type matches the input image
139 if padY[0] or padX[0]:
140 paddedImage = cv2.copyMakeBorder(
141 img, *(0, padY[0]), *(0, padX[0]), cv2.BORDER_REPLICATE, None if out is None else pyramid[0], None
142 ).astype(img.dtype)
143 else:
144 paddedImage = img
146 # Store the first level of the pyramid (padded image)
147 if out is None:
148 pyramid.append(paddedImage)
149 else:
150 # This might not be sound all the time, copy might be needed!
151 # Update the first level in the provided list
152 pyramid[0] = paddedImage
154 # Generate each subsequent level of the Gaussian Pyramid
155 for i in range(1, len(padY)):
156 if padY[i] or padX[i]:
157 paddedImage = cv2.copyMakeBorder(
158 paddedImage, *(0, padY[i]), *(0, padX[i]), cv2.BORDER_REPLICATE, None, None
159 ).astype(img.dtype)
160 # Downsample the image
161 paddedImage = cv2.pyrDown(paddedImage, None if out is None else pyramid[i])
163 # Append to the list if not provided externally
164 if out is None:
165 pyramid.append(paddedImage)
166 return pyramid
169def makeLapPyramid(
170 img: NDArray,
171 padY: list[int],
172 padX: list[int],
173 gaussOut: List[NDArray] | None,
174 lapOut: List[NDArray] | None,
175 upscratch: List[NDArray] | None = None,
176) -> Sequence[NDArray]:
177 """
178 Create a Laplacian pyramid from the input image.
180 This function constructs a Laplacian pyramid from the input image. It first
181 generates a Gaussian pyramid and then, for each level (except the last),
182 subtracts the upsampled version of the next lower level from the current
183 level to obtain the Laplacian levels. If `lapOut` is None, it creates a
184 new list to store the Laplacian pyramid; otherwise, it uses the provided
185 `lapOut`.
187 Parameters
188 ----------
189 img : `NDArray`
190 The input image as a numpy array.
191 padY : `list` of `int`
192 List of padding sizes for rows (vertical padding).
193 padX : `list` of `int`
194 List of padding sizes for columns (horizontal padding).
195 gaussOut : `numba.typed.typedlist.List` of `NDArray` or None
196 Preallocated storage for the output of the Gaussian pyramid function.
197 If `None` new storage is allocated.
198 lapOut : `numba.typed.typedlist.List` of `NDArray` or None
199 Preallocated for the output Laplacian pyramid. If None, a new
200 `numba.typed.typedlist.List` is created.
201 upscratch : `numba.typed.typedlist.List` of `NDarray`, optional
202 List to store intermediate results of pyramids (default is None).
204 Returns
205 -------
206 results : `Sequence` of `NDArray`
207 The Laplacian pyramid as a sequence of numpy arrays.
209 """
210 pyramid = makeGaussianPyramid(img, padY, padX, gaussOut)
211 if lapOut is None:
212 lapPyramid = List()
213 else:
214 lapPyramid = lapOut
215 for i in range(len(pyramid) - 1):
216 upsampled = cv2.pyrUp(pyramid[i + 1], None if upscratch is None else upscratch[i + 1])
217 if padY[i + 1] or padX[i + 1]:
218 upsampled = upsampled[
219 : upsampled.shape[0] - 2 * padY[i + 1], : upsampled.shape[1] - 2 * padX[i + 1]
220 ]
221 if lapOut is None:
222 lapPyramid.append(pyramid[i] - upsampled)
223 else:
224 cv2.subtract(pyramid[i], upsampled, dst=lapPyramid[i])
225 if lapOut is None:
226 lapPyramid.append(pyramid[-1])
227 else:
228 lapPyramid[-1][:, :] = pyramid[-1]
229 return lapPyramid
232@njit(fastmath=True, parallel=True, error_model="numpy", nogil=True)
233def _calculateOutput(
234 out: List[NDArray],
235 pyramid: List[NDArray],
236 gamma: NDArray,
237 pyramidVectorsBottom: List[NDArray],
238 pyramidVectorsTop: List[NDArray],
239):
240 """
241 Computes the output by interpolating between basis vectors at each pixel in
242 a Gaussian pyramid.
244 The function iterates over each pixel in the Gaussian pyramids
245 and interpolates between the corresponding basis vectors from
246 `pyramidVectorsBottom` and `pyramidVectorsTop`. If a pixel value is outside
247 the range defined by gamma, it skips interpolation.
249 Parameters:
250 -----------
251 out : `numba.typed.typedlist.List` of `np.ndarray`
252 A list of numpy arrays representing the output image pyramids.
253 pyramid : `numba.typed.typedlist.List` of `np.ndarray`
254 A list of numpy arrays representing the Gaussian pyramids.
255 gamma : `np.ndarray`
256 A numpy array containing the range for pixel values to be considered in
257 the interpolation.
258 pyramidVectorsBottom : `numba.typed.typedlist.List` of `np.ndarray`
259 A list of numpy arrays representing the basis vectors at the bottom
260 level of each pyramid layer.
261 pyramidVectorsTop : `numba.typed.typedlist.List` of `np.ndarray`
262 A list of numpy arrays representing the basis vectors at the top level
263 of each pyramid layer.
265 """
266 # loop over each pixel in the gaussian pyramid
267 # gammaDiff = gamma[1] - gamma[0]
268 for level in prange(0, len(pyramid) - 1):
269 yshape = pyramid[level].shape[0]
270 xshape = pyramid[level].shape[1]
271 plevel = pyramid[level]
272 outlevel = out[level]
273 basisBottom = pyramidVectorsBottom[level]
274 basisTop = pyramidVectorsTop[level]
275 for y in prange(yshape):
276 plevelY = plevel[y]
277 outLevelY = outlevel[y]
278 basisBottomY = basisBottom[y]
279 basisTopY = basisTop[y]
280 for x in prange(xshape):
281 val = plevelY[x]
282 if not (val >= gamma[0] and val <= gamma[1]):
283 continue
284 a = (plevelY[x] - gamma[0]) / (gamma[1] - gamma[0])
285 outLevelY[x] = (1 - a) * basisBottomY[x] + a * basisTopY[x]
288def levelPadder(numb: int, levels: int) -> list[int]:
289 """Determine if each level of transform will need to be padded by
290 one to make the level divisible by two.
292 Parameters
293 ----------
294 numb : int
295 The size of the input dimension
296 levels : int
297 The number of times the dimensions will be reduced by a factor of two
299 Returns
300 -------
301 padds : list of int
302 A list where the entries are either zero or one depending on if the
303 size will need padded to be a power of two.
305 """
306 pads = []
307 if numb % 2 != 0:
308 pads.append(1)
309 numb += 1
310 else:
311 pads.append(0)
312 for _ in range(levels):
313 numb /= 2
314 if numb % 2 != 0:
315 pads.append(1)
316 numb += 1
317 else:
318 pads.append(0)
319 return pads
322def localContrast(
323 image: NDArray,
324 sigma: float,
325 highlights: float = -0.9,
326 shadows: float = 0.4,
327 clarity: float = 0.15,
328 maxLevel: int | None = None,
329 numGamma: int = 20,
330) -> NDArray:
331 """Enhance the local contrast of an input image.
333 Parameters
334 ----------
335 image : `NDArray`
336 Two dimensional numpy array representing the image to have contrast
337 increased.
338 sigma : `float`
339 The scale over which edges are considered real and not noise.
340 highlights : `float`
341 A parameter that controls how highlights are enhansed or reduced,
342 contrary to intuition, negative values increase highlights.
343 shadows : `float`
344 A parameter that controls how shadows are deepened.
345 clarity : `float`
346 A parameter that relates to the contrast between highlights and
347 shadow.
348 maxLevel : `int` or `None`
349 The maximum number of image pyramid levels to enhanse the contrast over.
350 Each level has a spatial scale of roughly 2^(level) pixles.
351 numGamma : `int`
352 This is an optimization parameter. This algorithm divides up contrast
353 space into a certain numbers over which the expensive computation
354 is done. Contrast values in the image which fall between two of these
355 values are interpolated to get the outcome. The higher the numGamma,
356 the smoother the image is post contrast enhancement, though above
357 some number there is no decerable difference.
359 Returns
360 -------
361 image : `NDArray`
362 Two dimensional numpy array of the input image with increased local
363 contrast.
365 Raises
366 ------
367 ValueError
368 Raised if the max level to enhance to is greater than the image
369 supports.
371 Notes
372 -----
373 This function, and it's supporting functions, spiritually implement the
374 algorithm outlined at
375 https://people.csail.mit.edu/sparis/publi/2011/siggraph/
376 titled "Local Laplacian Filters: Edge-aware Image Processing with Laplacian
377 Pyramid". This is not a 1:1 implementation, it's optimized for the
378 python language and runtime performance. Most notably it transforms only
379 certain levels and linearly interpolates to find other values. This
380 implementation is inspired by the ony done in the darktable image editor:
381 https://www.darktable.org/2017/11/local-laplacian-pyramids/. None of the
382 code is in common, nor is the implementation 1:1, but reading the original
383 paper and the darktable implementation gives more info about this function.
384 Specifically some variable names follow the paper/other implementation,
385 and may be confusing when viewed without that context.
387 """
388 # ensure the supplied values are floats, and not ints
389 highlights = float(highlights)
390 shadows = float(shadows)
391 clarity = float(clarity)
393 # Determine the maximum level over which the image will be inhanced
394 # and the amount of padding that will be needed to be added to the
395 # image.
396 maxImageLevel = int(np.min(np.log2(image.shape)))
397 if maxLevel is None:
398 maxLevel = maxImageLevel
399 if maxImageLevel < maxLevel:
400 raise ValueError(
401 f"The supplied max level {maxLevel} is is greater than the max of the image: {maxImageLevel}"
402 )
403 support = 1 << (maxLevel - 1)
404 padY_amounts = levelPadder(image.shape[0] + support, maxLevel)
405 padX_amounts = levelPadder(image.shape[1] + support, maxLevel)
406 imagePadded = cv2.copyMakeBorder(
407 image, *(0, support), *(0, support), cv2.BORDER_REPLICATE, None, None
408 ).astype(image.dtype)
410 # build a list of intensities
411 gamma = np.linspace(image.min(), image.max(), numGamma)
413 # make gaussian pyramid
414 pyramid = makeGaussianPyramid(imagePadded, padY_amounts, padX_amounts, None)
416 finalPyramid = List()
417 for sample in pyramid[:-1]:
418 finalPyramid.append(np.zeros_like(sample))
419 finalPyramid.append(pyramid[-1])
421 # make a working array for gaussian pyramid in Lap
422 # make two working arrays for laplace as the true value is interpolated
423 # between the endpoints.
424 # This prevents needing re-allocations which can be time consuming.
425 tmpGauss = List()
426 tmpLap1 = List()
427 tmpLap2 = List()
428 upscratch = List()
429 for i, sample in enumerate(pyramid):
430 tmpGauss.append(np.empty_like(sample))
431 tmpLap1.append(np.empty_like(sample))
432 tmpLap2.append(np.empty_like(sample))
433 if i == 0:
434 upscratch.append(np.empty((0, 0), dtype=image.dtype))
435 continue
436 upscratch.append(np.empty((sample.shape[0] * 2, sample.shape[1] * 2), dtype=image.dtype))
437 # cycle between the endpoints, because there is no reason to recalculate both
438 # endpoints as only one changes for each bin.
439 cycler = iter(cycle((tmpLap1, tmpLap2)))
440 # allocate temporary arrays to use for each bin
441 outCycle = iter(cycle((np.copy(imagePadded), np.copy(imagePadded))))
442 prevImg = r(
443 imagePadded, next(outCycle), gamma[0], sigma, shadows=shadows, highlights=highlights, clarity=clarity
444 )
445 prevLapPyr = makeLapPyramid(
446 prevImg, padY_amounts, padX_amounts, tmpGauss, next(cycler), upscratch=upscratch
447 )
449 for value in range(1, len(gamma) - 1):
450 pyramidVectors = List()
451 pyramidVectors.append(prevLapPyr)
452 newImg = r(
453 imagePadded,
454 next(outCycle),
455 gamma[value],
456 sigma,
457 shadows=shadows,
458 highlights=highlights,
459 clarity=clarity,
460 )
461 prevLapPyr = makeLapPyramid(
462 newImg, padY_amounts, padX_amounts, tmpGauss, next(cycler), upscratch=upscratch
463 )
464 pyramidVectors.append(prevLapPyr)
466 _calculateOutput(
467 finalPyramid,
468 pyramid,
469 np.array((gamma[value - 1], gamma[value])),
470 pyramidVectors[0],
471 pyramidVectors[1],
472 )
473 del pyramidVectors
475 # time to reconstruct
476 output = finalPyramid[-1]
477 for i in range(-2, -1 * len(finalPyramid) - 1, -1):
478 upsampled = cv2.pyrUp(output)
479 upsampled = upsampled[
480 : upsampled.shape[0] - 2 * padY_amounts[i + 1], : upsampled.shape[1] - 2 * padX_amounts[i + 1]
481 ]
482 output = finalPyramid[i] + upsampled
483 return output[:-support, :-support]