| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275 |
- import numpy as np
- from scipy.stats import entropy
- from ..util._backends import dispatchable
- from ..util.dtype import dtype_range
- from .._shared.utils import _supported_float_type, check_shape_equality, warn
- __all__ = [
- 'mean_squared_error',
- 'normalized_root_mse',
- 'peak_signal_noise_ratio',
- 'normalized_mutual_information',
- ]
- def _as_floats(image0, image1):
- """
- Promote im1, im2 to nearest appropriate floating point precision.
- """
- float_type = _supported_float_type((image0.dtype, image1.dtype))
- image0 = np.asarray(image0, dtype=float_type)
- image1 = np.asarray(image1, dtype=float_type)
- return image0, image1
- @dispatchable
- def mean_squared_error(image0, image1):
- """
- Compute the mean-squared error between two images.
- Parameters
- ----------
- image0, image1 : ndarray
- Images. Any dimensionality, must have same shape.
- Returns
- -------
- mse : float
- The mean-squared error (MSE) metric.
- Notes
- -----
- .. versionchanged:: 0.16
- This function was renamed from ``skimage.measure.compare_mse`` to
- ``skimage.metrics.mean_squared_error``.
- """
- check_shape_equality(image0, image1)
- image0, image1 = _as_floats(image0, image1)
- return np.mean((image0 - image1) ** 2, dtype=np.float64)
- @dispatchable
- def normalized_root_mse(image_true, image_test, *, normalization='euclidean'):
- """
- Compute the normalized root mean-squared error (NRMSE) between two
- images.
- Parameters
- ----------
- image_true : ndarray
- Ground-truth image, same shape as im_test.
- image_test : ndarray
- Test image.
- normalization : {'euclidean', 'min-max', 'mean'}, optional
- Controls the normalization method to use in the denominator of the
- NRMSE. There is no standard method of normalization across the
- literature [1]_. The methods available here are as follows:
- - 'euclidean' : normalize by the averaged Euclidean norm of
- ``im_true``::
- NRMSE = RMSE * sqrt(N) / || im_true ||
- where || . || denotes the Frobenius norm and ``N = im_true.size``.
- This result is equivalent to::
- NRMSE = || im_true - im_test || / || im_true ||.
- - 'min-max' : normalize by the intensity range of ``im_true``.
- - 'mean' : normalize by the mean of ``im_true``
- Returns
- -------
- nrmse : float
- The NRMSE metric.
- Notes
- -----
- .. versionchanged:: 0.16
- This function was renamed from ``skimage.measure.compare_nrmse`` to
- ``skimage.metrics.normalized_root_mse``.
- References
- ----------
- .. [1] https://en.wikipedia.org/wiki/Root-mean-square_deviation
- """
- check_shape_equality(image_true, image_test)
- image_true, image_test = _as_floats(image_true, image_test)
- # Ensure that both 'Euclidean' and 'euclidean' match
- normalization = normalization.lower()
- if normalization == 'euclidean':
- denom = np.sqrt(np.mean((image_true * image_true), dtype=np.float64))
- elif normalization == 'min-max':
- denom = image_true.max() - image_true.min()
- elif normalization == 'mean':
- denom = image_true.mean()
- else:
- raise ValueError("Unsupported norm_type")
- return np.sqrt(mean_squared_error(image_true, image_test)) / denom
- def peak_signal_noise_ratio(image_true, image_test, *, data_range=None):
- """
- Compute the peak signal to noise ratio (PSNR) for an image.
- Parameters
- ----------
- image_true : ndarray
- Ground-truth image, same shape as im_test.
- image_test : ndarray
- Test image.
- data_range : int, optional
- The data range of the input image (distance between minimum and
- maximum possible values). By default, this is estimated from the image
- data-type.
- Returns
- -------
- psnr : float
- The PSNR metric.
- Notes
- -----
- .. versionchanged:: 0.16
- This function was renamed from ``skimage.measure.compare_psnr`` to
- ``skimage.metrics.peak_signal_noise_ratio``.
- References
- ----------
- .. [1] https://en.wikipedia.org/wiki/Peak_signal-to-noise_ratio
- """
- check_shape_equality(image_true, image_test)
- if data_range is None:
- if image_true.dtype != image_test.dtype:
- warn(
- "Inputs have mismatched dtype. Setting data_range based on "
- "image_true."
- )
- dmin, dmax = dtype_range[image_true.dtype.type]
- true_min, true_max = np.min(image_true), np.max(image_true)
- if true_max > dmax or true_min < dmin:
- raise ValueError(
- "image_true has intensity values outside the range expected "
- "for its data type. Please manually specify the data_range."
- )
- if true_min >= 0:
- # most common case (255 for uint8, 1 for float)
- data_range = dmax
- else:
- data_range = dmax - dmin
- image_true, image_test = _as_floats(image_true, image_test)
- err = mean_squared_error(image_true, image_test)
- data_range = float(data_range) # prevent overflow for small integer types
- return 10 * np.log10((data_range**2) / err)
- def _pad_to(arr, shape):
- """Pad an array with trailing zeros to a given target shape.
- Parameters
- ----------
- arr : ndarray
- The input array.
- shape : tuple
- The target shape.
- Returns
- -------
- padded : ndarray
- The padded array.
- Examples
- --------
- >>> _pad_to(np.ones((1, 1), dtype=int), (1, 3))
- array([[1, 0, 0]])
- """
- if not all(s >= i for s, i in zip(shape, arr.shape)):
- raise ValueError(
- f'Target shape {shape} cannot be smaller than input'
- f'shape {arr.shape} along any axis.'
- )
- padding = [(0, s - i) for s, i in zip(shape, arr.shape)]
- return np.pad(arr, pad_width=padding, mode='constant', constant_values=0)
- def normalized_mutual_information(image0, image1, *, bins=100):
- r"""Compute the normalized mutual information (NMI).
- The normalized mutual information of :math:`A` and :math:`B` is given by:
- .. math::
- Y(A, B) = \frac{H(A) + H(B)}{H(A, B)}
- where :math:`H(X) := - \sum_{x \in X}{p(x) \log p(x)}` is the entropy,
- :math:`X` is the set of image values, and :math:`p(x)` is the probability
- of occurrence of value :math:`x \in X`.
- It was proposed to be useful in registering images by Colin Studholme and
- colleagues [1]_. It ranges from 1 (perfectly uncorrelated image values)
- to 2 (perfectly correlated image values, whether positively or negatively).
- Parameters
- ----------
- image0, image1 : ndarray
- Images to be compared. The two input images must have the same number
- of dimensions.
- bins : int or sequence of int, optional
- The number of bins along each axis of the joint histogram.
- Returns
- -------
- nmi : float
- The normalized mutual information between the two arrays, computed at
- the granularity given by ``bins``. Higher NMI implies more similar
- input images.
- Raises
- ------
- ValueError
- If the images don't have the same number of dimensions.
- Notes
- -----
- If the two input images are not the same shape, the smaller image is padded
- with zeros.
- References
- ----------
- .. [1] C. Studholme, D.L.G. Hill, & D.J. Hawkes (1999). An overlap
- invariant entropy measure of 3D medical image alignment.
- Pattern Recognition 32(1):71-86
- :DOI:`10.1016/S0031-3203(98)00091-0`
- """
- if image0.ndim != image1.ndim:
- raise ValueError(
- f'NMI requires images of same number of dimensions. '
- f'Got {image0.ndim}D for `image0` and '
- f'{image1.ndim}D for `image1`.'
- )
- if image0.shape != image1.shape:
- max_shape = np.maximum(image0.shape, image1.shape)
- padded0 = _pad_to(image0, max_shape)
- padded1 = _pad_to(image1, max_shape)
- else:
- padded0, padded1 = image0, image1
- hist, bin_edges = np.histogramdd(
- [np.reshape(padded0, -1), np.reshape(padded1, -1)],
- bins=bins,
- density=True,
- )
- H0 = entropy(np.sum(hist, axis=0))
- H1 = entropy(np.sum(hist, axis=1))
- H01 = entropy(np.reshape(hist, -1))
- return (H0 + H1) / H01
|