Brisque changes

live_python_qa/brisque.py

@@ -1,14 +1,17 @@
 from .utils import compute_image_mscn_transform, extract_subband_feats
-
+from .imresize import imresize
 import numpy as np
 import cv2
 
 
-# TODO: Compare with Matlab
+
+
 def brisque(image):
     y_mscn = compute_image_mscn_transform(image)
-    half_scale = cv2.resize(image, dsize=(0, 0), fx=0.5, fy=0.5, interpolation=cv2.INTER_LANCZOS4)
+    half_scale = imresize(image, scalar_scale = 0.5, method = 'bicubic')
     y_half_mscn = compute_image_mscn_transform(half_scale)
     feats_full = extract_subband_feats(y_mscn)
     feats_half = extract_subband_feats(y_half_mscn)
     return np.concatenate((feats_full, feats_half))
+
+

live_python_qa/imresize.py

@@ -0,0 +1,149 @@
+"""
+Matlab Interpolation in Python taken from https://github.com/fatheral/matlab_imresize
+
+"""
+
+from __future__ import print_function
+import numpy as np
+from math import ceil, floor
+
+def deriveSizeFromScale(img_shape, scale):
+    output_shape = []
+    for k in range(2):
+        output_shape.append(int(ceil(scale[k] * img_shape[k])))
+    return output_shape
+
+def deriveScaleFromSize(img_shape_in, img_shape_out):
+    scale = []
+    for k in range(2):
+        scale.append(1.0 * img_shape_out[k] / img_shape_in[k])
+    return scale
+
+def triangle(x):
+    x = np.array(x).astype(np.float64)
+    lessthanzero = np.logical_and((x>=-1),x<0)
+    greaterthanzero = np.logical_and((x<=1),x>=0)
+    f = np.multiply((x+1),lessthanzero) + np.multiply((1-x),greaterthanzero)
+    return f
+
+def cubic(x):
+    x = np.array(x).astype(np.float64)
+    absx = np.absolute(x)
+    absx2 = np.multiply(absx, absx)
+    absx3 = np.multiply(absx2, absx)
+    f = np.multiply(1.5*absx3 - 2.5*absx2 + 1, absx <= 1) + np.multiply(-0.5*absx3 + 2.5*absx2 - 4*absx + 2, (1 < absx) & (absx <= 2))
+    return f
+
+def contributions(in_length, out_length, scale, kernel, k_width):
+    if scale < 1:
+        h = lambda x: scale * kernel(scale * x)
+        kernel_width = 1.0 * k_width / scale
+    else:
+        h = kernel
+        kernel_width = k_width
+    x = np.arange(1, out_length+1).astype(np.float64)
+    u = x / scale + 0.5 * (1 - 1 / scale)
+    left = np.floor(u - kernel_width / 2)
+    P = int(ceil(kernel_width)) + 2
+    ind = np.expand_dims(left, axis=1) + np.arange(P) - 1 # -1 because indexing from 0
+    indices = ind.astype(np.int32)
+    weights = h(np.expand_dims(u, axis=1) - indices - 1) # -1 because indexing from 0
+    weights = np.divide(weights, np.expand_dims(np.sum(weights, axis=1), axis=1))
+    aux = np.concatenate((np.arange(in_length), np.arange(in_length - 1, -1, step=-1))).astype(np.int32)
+    indices = aux[np.mod(indices, aux.size)]
+    ind2store = np.nonzero(np.any(weights, axis=0))
+    weights = weights[:, ind2store]
+    indices = indices[:, ind2store]
+    return weights, indices
+
+def imresizemex(inimg, weights, indices, dim):
+    in_shape = inimg.shape
+    w_shape = weights.shape
+    out_shape = list(in_shape)
+    out_shape[dim] = w_shape[0]
+    outimg = np.zeros(out_shape)
+    if dim == 0:
+        for i_img in range(in_shape[1]):
+            for i_w in range(w_shape[0]):
+                w = weights[i_w, :]
+                ind = indices[i_w, :]
+                im_slice = inimg[ind, i_img].astype(np.float64)
+                outimg[i_w, i_img] = np.sum(np.multiply(np.squeeze(im_slice, axis=0), w.T), axis=0)
+    elif dim == 1:
+        for i_img in range(in_shape[0]):
+            for i_w in range(w_shape[0]):
+                w = weights[i_w, :]
+                ind = indices[i_w, :]
+                im_slice = inimg[i_img, ind].astype(np.float64)
+                outimg[i_img, i_w] = np.sum(np.multiply(np.squeeze(im_slice, axis=0), w.T), axis=0)        
+    if inimg.dtype == np.uint8:
+        outimg = np.clip(outimg, 0, 255)
+        return np.around(outimg).astype(np.uint8)
+    else:
+        return outimg
+
+def imresizevec(inimg, weights, indices, dim):
+    wshape = weights.shape
+    if dim == 0:
+        weights = weights.reshape((wshape[0], wshape[2], 1, 1))
+        outimg =  np.sum(weights*((inimg[indices].squeeze(axis=1)).astype(np.float64)), axis=1)
+    elif dim == 1:
+        weights = weights.reshape((1, wshape[0], wshape[2], 1))
+        outimg =  np.sum(weights*((inimg[:, indices].squeeze(axis=2)).astype(np.float64)), axis=2)
+    if inimg.dtype == np.uint8:
+        outimg = np.clip(outimg, 0, 255)
+        return np.around(outimg).astype(np.uint8)
+    else:
+        return outimg
+
+def resizeAlongDim(A, dim, weights, indices, mode="vec"):
+    if mode == "org":
+        out = imresizemex(A, weights, indices, dim)
+    else:
+        out = imresizevec(A, weights, indices, dim)
+    return out
+
+def imresize(I, scalar_scale=None, method='bicubic', output_shape=None, mode="vec"):
+    if method is 'bicubic':
+        kernel = cubic
+    elif method is 'bilinear':
+        kernel = triangle
+    else:
+        print ('Error: Unidentified method supplied')
+
+    kernel_width = 4.0
+    # Fill scale and output_size
+    if scalar_scale is not None:
+        scalar_scale = float(scalar_scale)
+        scale = [scalar_scale, scalar_scale]
+        output_size = deriveSizeFromScale(I.shape, scale)
+    elif output_shape is not None:
+        scale = deriveScaleFromSize(I.shape, output_shape)
+        output_size = list(output_shape)
+    else:
+        print ('Error: scalar_scale OR output_shape should be defined!')
+        return
+    scale_np = np.array(scale)
+    order = np.argsort(scale_np)
+    weights = []
+    indices = []
+    for k in range(2):
+        w, ind = contributions(I.shape[k], output_size[k], scale[k], kernel, kernel_width)
+        weights.append(w)
+        indices.append(ind)
+    B = np.copy(I) 
+    flag2D = False
+    if B.ndim == 2:
+        B = np.expand_dims(B, axis=2)
+        flag2D = True
+    for k in range(2):
+        dim = order[k]
+        B = resizeAlongDim(B, dim, weights[dim], indices[dim], mode)
+    if flag2D:
+        B = np.squeeze(B, axis=2)
+    return B
+
+def convertDouble2Byte(I):
+    B = np.clip(I, 0.0, 1.0)
+    B = 255*B
+    return np.around(B).astype(np.uint8)

live_python_qa/utils.py

@@ -41,12 +41,12 @@ def estimateggdparam(vec):
     r_gam = (gamma(1.0/gam)*gamma(3.0/gam))/((gamma(2.0/gam))**2)
     # print(np.mean(vec))
     sigma_sq = np.mean(vec**2)  #-(np.mean(vec))**2
-    sigma = np.sqrt(sigma_sq)
+    #sigma = np.sqrt(sigma_sq)
     E = np.mean(np.abs(vec))
     rho = sigma_sq / (E**2 + 1e-6)
     array_position = (np.abs(rho - r_gam)).argmin()
     alphaparam = gam[array_position]
-    return alphaparam, sigma
+    return alphaparam, sigma_sq
 
 
 def compute_image_mscn_transform(image, C=1, avg_window=None, extend_mode='constant'):