# -*- coding: utf-8 -*-
# Copyright (C) 2016-2020 by Brendt Wohlberg <brendt@ieee.org>
# All rights reserved. BSD 3-clause License.
# This file is part of the SPORCO package. Details of the copyright
# and user license can be found in the 'LICENSE.txt' file distributed
# with the package.
"""Image quality metrics and related functions
Note that the well-known SSIM metric is not implemented here as it is
available in a number of other Python packages, including:
- `scikit-image <https://github.com/scikit-image/scikit-image>`_
- `PyMetrikz <https://bitbucket.org/kuraiev/pymetrikz>`_
- `Video Quality Metrics <https://github.com/aizvorski/video-quality>`_
- `pyssim <https://github.com/jterrace/pyssim>`_
Some implementations are also available in unpackaged collections of
Python code:
- `src <https://github.com/helderc/src>`_
- `python <https://github.com/mubeta06/python>`_
|
"""
from __future__ import division
import numpy as np
from scipy import ndimage
from scipy import signal
__author__ = """Brendt Wohlberg <brendt@ieee.org>"""
[docs]
def mae(vref, vcmp):
"""
Compute Mean Absolute Error (MAE) between two images.
Parameters
----------
vref : array_like
Reference image
vcmp : array_like
Comparison image
Returns
-------
x : float
MAE between `vref` and `vcmp`
"""
if np.iscomplexobj(vref) or np.iscomplexobj(vcmp):
dtype = np.complex128
else:
dtype = np.float64
r = np.asarray(vref, dtype=dtype).ravel()
c = np.asarray(vcmp, dtype=dtype).ravel()
return np.mean(np.abs(r - c))
[docs]
def mse(vref, vcmp):
"""
Compute Mean Squared Error (MSE) between two images.
Parameters
----------
vref : array_like
Reference image
vcmp : array_like
Comparison image
Returns
-------
x : float
MSE between `vref` and `vcmp`
"""
if np.iscomplexobj(vref) or np.iscomplexobj(vcmp):
dtype = np.complex128
else:
dtype = np.float64
r = np.asarray(vref, dtype=dtype).ravel()
c = np.asarray(vcmp, dtype=dtype).ravel()
return np.mean(np.abs(r - c)**2)
[docs]
def snr(vref, vcmp):
"""
Compute Signal to Noise Ratio (SNR) of two images.
Parameters
----------
vref : array_like
Reference image
vcmp : array_like
Comparison image
Returns
-------
x : float
SNR of `vcmp` with respect to `vref`
"""
dv = np.var(vref)
with np.errstate(divide='ignore'):
rt = dv / mse(vref, vcmp)
return 10.0 * np.log10(rt)
[docs]
def psnr(vref, vcmp, rng=None):
"""
Compute Peak Signal to Noise Ratio (PSNR) of two images. The PSNR
calculation defaults to using the less common definition in terms
of the actual range (i.e. max minus min) of the reference signal
instead of the maximum possible range for the data type
(i.e. :math:`2^b-1` for a :math:`b` bit representation).
Parameters
----------
vref : array_like
Reference image
vcmp : array_like
Comparison image
rng : None or int, optional (default None)
Signal range, either the value to use (e.g. 255 for 8 bit samples) or
None, in which case the actual range of the reference signal is used
Returns
-------
x : float
PSNR of `vcmp` with respect to `vref`
"""
if rng is None:
rng = np.abs(vref.max() - vref.min())
dv = (rng + 0.0)**2
with np.errstate(divide='ignore'):
rt = dv / mse(vref, vcmp)
return 10.0 * np.log10(rt)
[docs]
def isnr(vref, vdeg, vrst):
"""
Compute Improvement Signal to Noise Ratio (ISNR) for reference,
degraded, and restored images.
Parameters
----------
vref : array_like
Reference image
vdeg : array_like
Degraded image
vrst : array_like
Restored image
Returns
-------
x : float
ISNR of `vrst` with respect to `vref` and `vdeg`
"""
msedeg = mse(vref, vdeg)
mserst = mse(vref, vrst)
with np.errstate(divide='ignore'):
rt = msedeg / mserst
return 10.0 * np.log10(rt)
[docs]
def bsnr(vblr, vnsy):
"""
Compute Blurred Signal to Noise Ratio (BSNR) for a blurred and noisy
image.
Parameters
----------
vblr : array_like
Blurred noise free image
vnsy : array_like
Blurred image with additive noise
Returns
-------
x : float
BSNR of `vnsy` with respect to `vblr` and `vdeg`
"""
blrvar = np.var(vblr)
nsevar = np.var(vnsy - vblr)
with np.errstate(divide='ignore'):
rt = blrvar / nsevar
return 10.0 * np.log10(rt)
[docs]
def pamse(vref, vcmp, rescale=True):
"""
Compute Perceptual-fidelity Aware Mean Squared Error (PAMSE) IQA metric
:cite:`xue-2013-perceptual`. This implementation is a translation of the
reference Matlab implementation provided by the authors of
:cite:`xue-2013-perceptual`.
Parameters
----------
vref : array_like
Reference image
vcmp : array_like
Comparison image
rescale : bool, optional (default True)
Rescale inputs so that `vref` has a maximum value of 255, as assumed
by reference implementation
Returns
-------
score : float
PAMSE IQA metric
"""
# Calculate difference, promoting to float if vref and vcmp have integer
# dtype
emap = np.asarray(vref, dtype=np.float64) - \
np.asarray(vcmp, dtype=np.float64)
# Input images in reference code on which this implementation is
# based are assumed to be on range [0,...,255].
if rescale:
emap *= (255.0 / vref.max())
sigma = 0.8
herr = ndimage.filters.gaussian_filter(emap, sigma)
score = np.mean(herr**2)
return score
[docs]
def gmsd(vref, vcmp, rescale=True, returnMap=False):
"""
Compute Gradient Magnitude Similarity Deviation (GMSD) IQA metric
:cite:`xue-2014-gradient`. This implementation is a translation of the
reference Matlab implementation provided by the authors of
:cite:`xue-2014-gradient`.
Parameters
----------
vref : array_like
Reference image
vcmp : array_like
Comparison image
rescale : bool, optional (default True)
Rescale inputs so that `vref` has a maximum value of 255, as assumed
by reference implementation
returnMap : bool, optional (default False)
Flag indicating whether quality map should be returned in addition to
scalar score
Returns
-------
score : float
GMSD IQA metric
quality_map : ndarray
Quality map
"""
# Input images in reference code on which this implementation is
# based are assumed to be on range [0,...,255].
if rescale:
scl = (255.0 / vref.max())
else:
scl = np.float32(1.0)
T = 170.0
dwn = 2
dx = np.array([[1, 0, -1], [1, 0, -1], [1, 0, -1]]) / 3.0
dy = dx.T
ukrn = np.ones((2, 2)) / 4.0
aveY1 = signal.convolve2d(scl * vref, ukrn, mode='same', boundary='symm')
aveY2 = signal.convolve2d(scl * vcmp, ukrn, mode='same', boundary='symm')
Y1 = aveY1[0::dwn, 0::dwn]
Y2 = aveY2[0::dwn, 0::dwn]
IxY1 = signal.convolve2d(Y1, dx, mode='same', boundary='symm')
IyY1 = signal.convolve2d(Y1, dy, mode='same', boundary='symm')
grdMap1 = np.sqrt(IxY1**2 + IyY1**2)
IxY2 = signal.convolve2d(Y2, dx, mode='same', boundary='symm')
IyY2 = signal.convolve2d(Y2, dy, mode='same', boundary='symm')
grdMap2 = np.sqrt(IxY2**2 + IyY2**2)
quality_map = (2*grdMap1*grdMap2 + T) / (grdMap1**2 + grdMap2**2 + T)
score = np.std(quality_map)
if returnMap:
return (score, quality_map)
else:
return score