Optimize adaptive histogram equalization

[hakonanes]

The implementation of adaptive histogram equalization (calling scikit-image's equalize_adapthist()) should be optimized.

Potential optimizations include:

• Use HyperSpy's improved map() function (seems like this function handles chunks better than the current use of Dask's map_blocks() in kikuchipy)
• Use the more restricted (and thus faster) private function for rescaling intensities after equalization (https://github.com/pyxem/kikuchipy/blob/develop/kikuchipy/pattern/_pattern.py#L88)
• Use scikit-image's private function _clahe() directly instead of via the public API (dangerous since the function can change whenever)
• Use another implementation? See https://stackoverflow.com/questions/43569566/adaptive-histogram-equalization-in-python, https://gist.github.com/sadimanna/52c320ce5c49e200ce398f800d39a2c1.

The implementation is spread across two files:

kikuchipy/kikuchipy/signals/ebsd.py

Lines 800 to 888 in c5df967

	def adaptive_histogram_equalization(
	self,
	kernel_size: Optional[Union[Tuple[int, int], List[int]]] = None,
	clip_limit: Union[int, float] = 0,
	nbins: int = 128,
	):
	"""Enhance the local contrast in an EBSD scan inplace using
	adaptive histogram equalization.
	This method uses :func:`skimage.exposure.equalize_adapthist`.
	Parameters
	----------
	kernel_size
	Shape of contextual regions for adaptive histogram
	equalization, default is 1/4 of image height and 1/4 of
	image width.
	clip_limit
	Clipping limit, normalized between 0 and 1 (higher values
	give more contrast). Default is 0.
	nbins
	Number of gray bins for histogram ("data range"), default is
	128.
	See also
	--------
	~kikuchipy.signals.EBSD.rescale_intensity
	~kikuchipy.signals.EBSD.normalize_intensity
	Examples
	--------
	To best understand how adaptive histogram equalization works,
	we plot the histogram of the same image before and after
	equalization:
	>>> import numpy as np
	>>> import matplotlib.pyplot as plt
	>>> import kikuchipy as kp
	>>> s = kp.data.nickel_ebsd_small()
	>>> s2 = s.inav[0, 0].deepcopy()
	>>> s2.adaptive_histogram_equalization() # doctest: +SKIP
	>>> hist, _ = np.histogram(
	... s.inav[0, 0].data, bins=255, range=(0, 255)
	... )
	>>> hist2, _ = np.histogram(s2.data, bins=255, range=(0, 255))
	>>> fig, ax = plt.subplots(nrows=2, ncols=2)
	>>> _ = ax[0, 0].imshow(s.inav[0, 0].data)
	>>> _ = ax[1, 0].plot(hist)
	>>> _ = ax[0, 1].imshow(s2.data)
	>>> _ = ax[1, 1].plot(hist2)
	Notes
	-----
	* It is recommended to perform adaptive histogram equalization
	only *after* static and dynamic background corrections,
	otherwise some unwanted darkening towards the edges might
	occur.
	* The default window size might not fit all pattern sizes, so it
	may be necessary to search for the optimal window size.
	"""
	# Determine window size (shape of contextual region)
	sig_shape = self.axes_manager.signal_shape
	if kernel_size is None:
	kernel_size = (sig_shape[0] // 4, sig_shape[1] // 4)
	elif isinstance(kernel_size, numbers.Number):
	kernel_size = (kernel_size,) * self.axes_manager.signal_dimension
	elif len(kernel_size) != self.axes_manager.signal_dimension:
	raise ValueError(f"Incorrect value of `shape`: {kernel_size}")
	kernel_size = [int(k) for k in kernel_size]
	# Create dask array of signal patterns and do processing on this
	dask_array = get_dask_array(signal=self)
	# Local contrast enhancement
	equalized_patterns = dask_array.map_blocks(
	func=chunk.adaptive_histogram_equalization,
	kernel_size=kernel_size,
	clip_limit=clip_limit,
	nbins=nbins,
	dtype=self.data.dtype,
	)
	# Overwrite signal patterns
	if not self._lazy:
	with ProgressBar():
	print("Adaptive histogram equalization:", file=sys.stdout)
	equalized_patterns.store(self.data, compute=True)
	else:
	self.data = equalized_patterns

which calls

kikuchipy/kikuchipy/pattern/chunk.py

Lines 364 to 411 in c5df967

	def adaptive_histogram_equalization(
	patterns: Union[np.ndarray, da.Array],
	kernel_size: Union[Tuple[int, int], List[int]],
	clip_limit: Union[int, float] = 0,
	nbins: int = 128,
	) -> np.ndarray:
	"""Local contrast enhancement of a chunk of EBSD patterns with
	adaptive histogram equalization.
	This method makes use of :func:`skimage.exposure.equalize_adapthist`.
	Parameters
	----------
	patterns
	EBSD patterns.
	kernel_size
	Shape of contextual regions for adaptive histogram equalization.
	clip_limit
	Clipping limit, normalized between 0 and 1 (higher values give
	more contrast). Default is 0.
	nbins
	Number of gray bins for histogram. Default is 128.
	Returns
	-------
	equalized_patterns : numpy.ndarray
	Patterns with enhanced contrast.
	"""
	dtype_in = patterns.dtype.type
	equalized_patterns = np.empty_like(patterns)
	for nav_idx in np.ndindex(patterns.shape[:-2]):
	# Adaptive histogram equalization
	equalized_pattern = equalize_adapthist(
	patterns[nav_idx],
	kernel_size=kernel_size,
	clip_limit=clip_limit,
	nbins=nbins,
	)
	# Rescale intensities
	equalized_patterns[nav_idx] = pattern_processing.rescale_intensity(
	equalized_pattern, dtype_out=dtype_in
	)
	return equalized_patterns

[hakonanes]

This method is also applicable to the EBSD master pattern. Applying it to the master pattern before dictionary generation with EBSDMasterPattern.get_patterns() seems to produce comparable relative intensities within patterns to the case where equalization is done to each individual pattern after dictionary generation (without equalization applied to the master pattern).

[hakonanes]

I don't adaptive histogram equalization myself, and have not encountered someone who finds this method in kikuchipy useful. If you're that person, feel free to say so and I'll re-open this issue!