322 lines
No EOL
11 KiB
C
322 lines
No EOL
11 KiB
C
/*
|
|
* Copyright (c) 2011, Tom Distler (http://tdistler.com)
|
|
* All rights reserved.
|
|
*
|
|
* The BSD License
|
|
*
|
|
* Redistribution and use in source and binary forms, with or without
|
|
* modification, are permitted provided that the following conditions are met:
|
|
*
|
|
* - Redistributions of source code must retain the above copyright notice,
|
|
* this list of conditions and the following disclaimer.
|
|
*
|
|
* - Redistributions in binary form must reproduce the above copyright notice,
|
|
* this list of conditions and the following disclaimer in the documentation
|
|
* and/or other materials provided with the distribution.
|
|
*
|
|
* - Neither the name of the tdistler.com nor the names of its contributors may
|
|
* be used to endorse or promote products derived from this software without
|
|
* specific prior written permission.
|
|
*
|
|
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
|
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
|
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
|
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
|
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
|
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
|
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
|
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|
* POSSIBILITY OF SUCH DAMAGE.
|
|
*/
|
|
|
|
#include "iqa.h"
|
|
#include "convolve.h"
|
|
#include "decimate.h"
|
|
#include "math_utils.h"
|
|
#include "ssim.h"
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
|
|
|
|
/* Forward declarations. */
|
|
IQA_INLINE static double _calc_luminance(float, float, float, float);
|
|
IQA_INLINE static double _calc_contrast(double, float, float, float, float);
|
|
IQA_INLINE static double _calc_structure(float, double, float, float, float, float);
|
|
static int _ssim_map(const struct _ssim_int *, void *);
|
|
static float _ssim_reduce(int, int, void *);
|
|
|
|
/*
|
|
* SSIM(x,y)=(2*ux*uy + C1)*(2sxy + C2) / (ux^2 + uy^2 + C1)*(sx^2 + sy^2 + C2)
|
|
* where,
|
|
* ux = SUM(w*x)
|
|
* sx = (SUM(w*(x-ux)^2)^0.5
|
|
* sxy = SUM(w*(x-ux)*(y-uy))
|
|
*
|
|
* Returns mean SSIM. MSSIM(X,Y) = 1/M * SUM(SSIM(x,y))
|
|
*/
|
|
float iqa_ssim(const unsigned char *ref, const unsigned char *cmp, int w, int h, int stride,
|
|
int gaussian, const struct iqa_ssim_args *args)
|
|
{
|
|
int scale;
|
|
int x,y,src_offset,offset;
|
|
float *ref_f,*cmp_f;
|
|
struct _kernel low_pass;
|
|
struct _kernel window;
|
|
float result;
|
|
double ssim_sum=0.0;
|
|
struct _map_reduce mr;
|
|
|
|
/* Initialize algorithm parameters */
|
|
scale = _max( 1, _round( (float)_min(w,h) / 256.0f ) );
|
|
if (args) {
|
|
if(args->f)
|
|
scale = args->f;
|
|
mr.map = _ssim_map;
|
|
mr.reduce = _ssim_reduce;
|
|
mr.context = (void*)&ssim_sum;
|
|
}
|
|
window.kernel = (float*)g_square_window;
|
|
window.w = window.h = SQUARE_LEN;
|
|
window.normalized = 1;
|
|
window.bnd_opt = KBND_SYMMETRIC;
|
|
if (gaussian) {
|
|
window.kernel = (float*)g_gaussian_window;
|
|
window.w = window.h = GAUSSIAN_LEN;
|
|
}
|
|
|
|
/* Convert image values to floats. Forcing stride = width. */
|
|
ref_f = (float*)malloc(w*h*sizeof(float));
|
|
cmp_f = (float*)malloc(w*h*sizeof(float));
|
|
if (!ref_f || !cmp_f) {
|
|
if (ref_f) free(ref_f);
|
|
if (cmp_f) free(cmp_f);
|
|
return INFINITY;
|
|
}
|
|
for (y=0; y<h; ++y) {
|
|
src_offset = y*stride;
|
|
offset = y*w;
|
|
for (x=0; x<w; ++x, ++offset, ++src_offset) {
|
|
ref_f[offset] = (float)ref[src_offset];
|
|
cmp_f[offset] = (float)cmp[src_offset];
|
|
}
|
|
}
|
|
|
|
/* Scale the images down if required */
|
|
if (scale > 1) {
|
|
/* Generate simple low-pass filter */
|
|
low_pass.kernel = (float*)malloc(scale*scale*sizeof(float));
|
|
if (!low_pass.kernel) {
|
|
free(ref_f);
|
|
free(cmp_f);
|
|
return INFINITY;
|
|
}
|
|
low_pass.w = low_pass.h = scale;
|
|
low_pass.normalized = 0;
|
|
low_pass.bnd_opt = KBND_SYMMETRIC;
|
|
for (offset=0; offset<scale*scale; ++offset)
|
|
low_pass.kernel[offset] = 1.0f/(scale*scale);
|
|
|
|
/* Resample */
|
|
if (_iqa_decimate(ref_f, w, h, scale, &low_pass, 0, 0, 0) ||
|
|
_iqa_decimate(cmp_f, w, h, scale, &low_pass, 0, &w, &h)) { /* Update w/h */
|
|
free(ref_f);
|
|
free(cmp_f);
|
|
free(low_pass.kernel);
|
|
return INFINITY;
|
|
}
|
|
free(low_pass.kernel);
|
|
}
|
|
|
|
result = _iqa_ssim(ref_f, cmp_f, w, h, &window, &mr, args);
|
|
|
|
free(ref_f);
|
|
free(cmp_f);
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
/* _iqa_ssim */
|
|
float _iqa_ssim(float *ref, float *cmp, int w, int h, const struct _kernel *k, const struct _map_reduce *mr, const struct iqa_ssim_args *args)
|
|
{
|
|
float alpha=1.0f, beta=1.0f, gamma=1.0f;
|
|
int L=255;
|
|
float K1=0.01f, K2=0.03f;
|
|
float C1,C2,C3;
|
|
int x,y,offset;
|
|
float *ref_mu,*cmp_mu,*ref_sigma_sqd,*cmp_sigma_sqd,*sigma_both;
|
|
double ssim_sum, numerator, denominator;
|
|
double luminance_comp, contrast_comp, structure_comp, sigma_root;
|
|
struct _ssim_int sint;
|
|
|
|
/* Initialize algorithm parameters */
|
|
if (args) {
|
|
if (!mr)
|
|
return INFINITY;
|
|
alpha = args->alpha;
|
|
beta = args->beta;
|
|
gamma = args->gamma;
|
|
L = args->L;
|
|
K1 = args->K1;
|
|
K2 = args->K2;
|
|
}
|
|
C1 = (K1*L)*(K1*L);
|
|
C2 = (K2*L)*(K2*L);
|
|
C3 = C2 / 2.0f;
|
|
|
|
ref_mu = (float*)malloc(w*h*sizeof(float));
|
|
cmp_mu = (float*)malloc(w*h*sizeof(float));
|
|
ref_sigma_sqd = (float*)malloc(w*h*sizeof(float));
|
|
cmp_sigma_sqd = (float*)malloc(w*h*sizeof(float));
|
|
sigma_both = (float*)malloc(w*h*sizeof(float));
|
|
if (!ref_mu || !cmp_mu || !ref_sigma_sqd || !cmp_sigma_sqd || !sigma_both) {
|
|
if (ref_mu) free(ref_mu);
|
|
if (cmp_mu) free(cmp_mu);
|
|
if (ref_sigma_sqd) free(ref_sigma_sqd);
|
|
if (cmp_sigma_sqd) free(cmp_sigma_sqd);
|
|
if (sigma_both) free(sigma_both);
|
|
return INFINITY;
|
|
}
|
|
|
|
/* Calculate mean */
|
|
_iqa_convolve(ref, w, h, k, ref_mu, 0, 0);
|
|
_iqa_convolve(cmp, w, h, k, cmp_mu, 0, 0);
|
|
|
|
for (y=0; y<h; ++y) {
|
|
offset = y*w;
|
|
for (x=0; x<w; ++x, ++offset) {
|
|
ref_sigma_sqd[offset] = ref[offset] * ref[offset];
|
|
cmp_sigma_sqd[offset] = cmp[offset] * cmp[offset];
|
|
sigma_both[offset] = ref[offset] * cmp[offset];
|
|
}
|
|
}
|
|
|
|
/* Calculate sigma */
|
|
_iqa_convolve(ref_sigma_sqd, w, h, k, 0, 0, 0);
|
|
_iqa_convolve(cmp_sigma_sqd, w, h, k, 0, 0, 0);
|
|
_iqa_convolve(sigma_both, w, h, k, 0, &w, &h); /* Update the width and height */
|
|
|
|
/* The convolution results are smaller by the kernel width and height */
|
|
for (y=0; y<h; ++y) {
|
|
offset = y*w;
|
|
for (x=0; x<w; ++x, ++offset) {
|
|
ref_sigma_sqd[offset] -= ref_mu[offset] * ref_mu[offset];
|
|
cmp_sigma_sqd[offset] -= cmp_mu[offset] * cmp_mu[offset];
|
|
sigma_both[offset] -= ref_mu[offset] * cmp_mu[offset];
|
|
}
|
|
}
|
|
|
|
ssim_sum = 0.0;
|
|
for (y=0; y<h; ++y) {
|
|
offset = y*w;
|
|
for (x=0; x<w; ++x, ++offset) {
|
|
|
|
if (!args) {
|
|
/* The default case */
|
|
numerator = (2.0 * ref_mu[offset] * cmp_mu[offset] + C1) * (2.0 * sigma_both[offset] + C2);
|
|
denominator = (ref_mu[offset]*ref_mu[offset] + cmp_mu[offset]*cmp_mu[offset] + C1) *
|
|
(ref_sigma_sqd[offset] + cmp_sigma_sqd[offset] + C2);
|
|
ssim_sum += numerator / denominator;
|
|
}
|
|
else {
|
|
/* User tweaked alpha, beta, or gamma */
|
|
|
|
/* passing a negative number to sqrt() cause a domain error */
|
|
if (ref_sigma_sqd[offset] < 0.0f)
|
|
ref_sigma_sqd[offset] = 0.0f;
|
|
if (cmp_sigma_sqd[offset] < 0.0f)
|
|
cmp_sigma_sqd[offset] = 0.0f;
|
|
sigma_root = sqrt(ref_sigma_sqd[offset] * cmp_sigma_sqd[offset]);
|
|
|
|
luminance_comp = _calc_luminance(ref_mu[offset], cmp_mu[offset], C1, alpha);
|
|
contrast_comp = _calc_contrast(sigma_root, ref_sigma_sqd[offset], cmp_sigma_sqd[offset], C2, beta);
|
|
structure_comp = _calc_structure(sigma_both[offset], sigma_root, ref_sigma_sqd[offset], cmp_sigma_sqd[offset], C3, gamma);
|
|
|
|
sint.l = luminance_comp;
|
|
sint.c = contrast_comp;
|
|
sint.s = structure_comp;
|
|
|
|
if (mr->map(&sint, mr->context))
|
|
return INFINITY;
|
|
}
|
|
}
|
|
}
|
|
|
|
free(ref_mu);
|
|
free(cmp_mu);
|
|
free(ref_sigma_sqd);
|
|
free(cmp_sigma_sqd);
|
|
free(sigma_both);
|
|
|
|
if (!args)
|
|
return (float)(ssim_sum / (double)(w*h));
|
|
return mr->reduce(w, h, mr->context);
|
|
}
|
|
|
|
|
|
/* _ssim_map */
|
|
int _ssim_map(const struct _ssim_int *si, void *ctx)
|
|
{
|
|
double *ssim_sum = (double*)ctx;
|
|
*ssim_sum += si->l * si->c * si->s;
|
|
return 0;
|
|
}
|
|
|
|
/* _ssim_reduce */
|
|
float _ssim_reduce(int w, int h, void *ctx)
|
|
{
|
|
double *ssim_sum = (double*)ctx;
|
|
return (float)(*ssim_sum / (double)(w*h));
|
|
}
|
|
|
|
|
|
/* _calc_luminance */
|
|
IQA_INLINE static double _calc_luminance(float mu1, float mu2, float C1, float alpha)
|
|
{
|
|
double result;
|
|
float sign;
|
|
/* For MS-SSIM* */
|
|
if (C1 == 0 && mu1*mu1 == 0 && mu2*mu2 == 0)
|
|
return 1.0;
|
|
result = (2.0 * mu1 * mu2 + C1) / (mu1*mu1 + mu2*mu2 + C1);
|
|
if (alpha == 1.0f)
|
|
return result;
|
|
sign = result < 0.0 ? -1.0f : 1.0f;
|
|
return sign * pow(fabs(result),(double)alpha);
|
|
}
|
|
|
|
/* _calc_contrast */
|
|
IQA_INLINE static double _calc_contrast(double sigma_comb_12, float sigma1_sqd, float sigma2_sqd, float C2, float beta)
|
|
{
|
|
double result;
|
|
float sign;
|
|
/* For MS-SSIM* */
|
|
if (C2 == 0 && sigma1_sqd + sigma2_sqd == 0)
|
|
return 1.0;
|
|
result = (2.0 * sigma_comb_12 + C2) / (sigma1_sqd + sigma2_sqd + C2);
|
|
if (beta == 1.0f)
|
|
return result;
|
|
sign = result < 0.0 ? -1.0f : 1.0f;
|
|
return sign * pow(fabs(result),(double)beta);
|
|
}
|
|
|
|
/* _calc_structure */
|
|
IQA_INLINE static double _calc_structure(float sigma_12, double sigma_comb_12, float sigma1, float sigma2, float C3, float gamma)
|
|
{
|
|
double result;
|
|
float sign;
|
|
/* For MS-SSIM* */
|
|
if (C3 == 0 && sigma_comb_12 == 0) {
|
|
if (sigma1 == 0 && sigma2 == 0)
|
|
return 1.0;
|
|
else if (sigma1 == 0 || sigma2 == 0)
|
|
return 0.0;
|
|
}
|
|
result = (sigma_12 + C3) / (sigma_comb_12 + C3);
|
|
if (gamma == 1.0f)
|
|
return result;
|
|
sign = result < 0.0 ? -1.0f : 1.0f;
|
|
return sign * pow(fabs(result),(double)gamma);
|
|
} |