msu_iqa_features.py 16.4 KB
Newer Older
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
1 2 3 4 5 6 7 8 9 10
'''
Created on 9 Feb 2016

@author: sbhatta
'''

import numpy as np
import scipy.signal as ssg
import bob.ip.base
import bob.ip.color
11
from . import galbally_iqm_features as iqm
12 13 14

#from . import tan_specular_highlights as tsh
from ._library import remove_highlights
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
15

Amir Mohammadi's avatar
Amir Mohammadi committed
16 17
''' Utility functions '''

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
18 19 20

def matlab_rgb2hsv(rgbImage):
    # first normalize the range of values to 0-1
Amir Mohammadi's avatar
Amir Mohammadi committed
21

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
22
    isUint8 = True
Amir Mohammadi's avatar
Amir Mohammadi committed
23 24 25
    if isUint8:
        rgbImage = rgbImage.astype(np.float64) / 255.0

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
26 27
    hsv = np.zeros_like(rgbImage)
    bob.ip.color.rgb_to_hsv(rgbImage, hsv)
Amir Mohammadi's avatar
Amir Mohammadi committed
28 29 30
    h = hsv[0, :, :]
    s = hsv[1, :, :]
    v = hsv[2, :, :]
David GEISSBUHLER's avatar
David GEISSBUHLER committed
31

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
32 33 34 35 36 37
    return (h, s, v)


def imshow(image):
    import matplotlib
    from matplotlib import pyplot as plt
Amir Mohammadi's avatar
Amir Mohammadi committed
38 39 40
    if len(image.shape) == 3:
        # imshow() expects color image in a slightly different format, so first
        # rearrange the 3d data for imshow...
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
41
        outImg = image.tolist()
42
        print(len(outImg))
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
43 44
        result = np.dstack((outImg[0], outImg[1]))
        outImg = np.dstack((result, outImg[2]))
Amir Mohammadi's avatar
Amir Mohammadi committed
45 46 47
        # [:,:,1], cmap=mpl.cm.gray)
        plt.imshow((outImg * 255.0).astype(np.uint8))

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
48
    else:
Amir Mohammadi's avatar
Amir Mohammadi committed
49 50
        if(len(image.shape) == 2):
            # display gray image.
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
51
            plt.imshow(image.astype(np.uint8), cmap=matplotlib.cm.gray)
Amir Mohammadi's avatar
Amir Mohammadi committed
52

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
53 54 55
    plt.show()


Amir Mohammadi's avatar
Amir Mohammadi committed
56 57
'''Auxilliary functions'''

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
58 59

def sobelEdgeMap(image, orientation='both'):
Amir Mohammadi's avatar
Amir Mohammadi committed
60 61 62 63 64 65 66 67 68 69 70 71

    # bob..sobel returns filter-responses which need to be thresholded to get
    # the edge-map
    thinning = 1
    refImage = image.astype(np.float)

    # compute edge map for reference image
    # returns 3D image. 1st dim is the edge-direction. 1st component is
    # vertical; 2nd component is hor. responses
    refSobel_sep = bob.ip.base.sobel(refImage)
    refSobelX = refSobel_sep[0, :, :]
    refSobelY = refSobel_sep[1, :, :]
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
72
    if orientation is 'horizontal':
Amir Mohammadi's avatar
Amir Mohammadi committed
73
        refEdge = iqm.edge_thinning(refSobelX[:, :], refSobelX[:, :], thinning)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
74 75
    else:
        if orientation is 'vertical':
Amir Mohammadi's avatar
Amir Mohammadi committed
76 77 78 79
            refEdge = iqm.edge_thinning(
                refSobelY[
                    :, :], refSobelY[
                    :, :], thinning)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
80
        else:
Amir Mohammadi's avatar
Amir Mohammadi committed
81 82 83 84
            refEdge = iqm.edge_thinning(
                refSobelX[
                    :, :], refSobelY[
                    :, :], thinning)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
85 86 87

    return refEdge

88

89
def compute_msu_iqa_features(rgbImage):
90
    """Computes image-quality features for the given input color (RGB) image.
Amir Mohammadi's avatar
Amir Mohammadi committed
91
    This is the main function to call.
92 93 94

    Parameters:

Amir Mohammadi's avatar
Amir Mohammadi committed
95 96 97
    rgbImage (:py:class:`numpy.ndarray`): A ``uint8`` array with 3 dimensions,
        representing the RGB input image of shape [3,M,N] (M rows x N cols).

98 99
    Returns:

Amir Mohammadi's avatar
Amir Mohammadi committed
100 101 102 103 104 105
    featSet (:py:class:`numpy.ndarray`): a 1D numpy array of 121 float32
        scalars. This function returns the image-quality features (for face
        anti- spoofing) that have been described by Wen et al. in their paper:
        "Face  spoof  detection  with  image distortion analysis", IEEE Trans.
        on Information Forensics and Security, vol. 10(4), pp. 746-761, April
        2015.
106
    """
Amir Mohammadi's avatar
Amir Mohammadi committed
107 108 109 110 111 112 113
    rgbImage = rgbImage.copy()
    assert len(
        rgbImage.shape) == 3, 'compute_msu_iqa_features():: image should be ' \
        'a 3D array (containing a rgb image)'
    # defined above. Calls Bob's rgb_to_hsv() after rescaling the input image.
    h, s, v = matlab_rgb2hsv(rgbImage)

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
114 115
    grayImage = np.zeros_like(h, dtype='uint8')
    bob.ip.color.rgb_to_gray(rgbImage, grayImage)
Amir Mohammadi's avatar
Amir Mohammadi committed
116

117
    # compute blur-features
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
118
    blurFeat = blurriness(grayImage)
Amir Mohammadi's avatar
Amir Mohammadi committed
119

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
120 121
    pinaBlur = marzilianoBlur(grayImage)
    pinaBlur /= 30.0
Amir Mohammadi's avatar
Amir Mohammadi committed
122

123
    # compute color-diversity features
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
124
    colorHist, totNumColors = calColorHist(rgbImage)
Amir Mohammadi's avatar
Amir Mohammadi committed
125 126
    totNumColors /= 2000.0  # as done in Matlab code provided by MSU.

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
127 128 129
    # calculate mean, deviation and skewness of each channel
    # use histogram shifting for the hue channel
    momentFeatsH = calmoment_shift(h)
Amir Mohammadi's avatar
Amir Mohammadi committed
130

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
131 132 133 134 135
    momentFeats = momentFeatsH.copy()
    momentFeatsS = calmoment(s)
    momentFeats = np.hstack((momentFeats, momentFeatsS))
    momentFeatsV = calmoment(v)
    momentFeats = np.hstack((momentFeats, momentFeatsV))
136

137
    # compute the image-specularity features
138
    speckleFeats = compute_iqa_specularity_features(rgbImage, startEps=0.06)
Amir Mohammadi's avatar
Amir Mohammadi committed
139 140 141

    # stack the various feature-values in the same order as in MSU's matlab
    # code.
142
    fv = speckleFeats.copy()
Amir Mohammadi's avatar
Amir Mohammadi committed
143

144
    fv = np.hstack((fv, momentFeats))
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
145 146 147 148 149 150 151 152
    fv = np.hstack((fv, colorHist))
    fv = np.hstack((fv, totNumColors))
    fv = np.hstack((fv, blurFeat))
    fv = np.hstack((fv, pinaBlur))

    return fv


153
def compute_iqa_specularity_features(rgbImage, startEps=0.05):
Amir Mohammadi's avatar
Amir Mohammadi committed
154 155 156
    """Returns three features characterizing the specularity present in input
    color image. First the specular and diffuse components of the input image
    are separated using the
157 158
    """

Amir Mohammadi's avatar
Amir Mohammadi committed
159
    # separate the specular and diffuse components of input color image.
160 161 162 163 164 165 166
    # original python version
    #speckleFreeImg, diffuseImg, speckleImg = tsh.remove_highlights(
    #    rgbImage.astype(float), startEps, verboseFlag=False)

    speckleFreeImg, diffuseImg, speckleImg = remove_highlights(
            rgbImage.astype(np.float32), startEps)

167 168 169 170
    speckleImg[np.where(np.isinf(speckleImg))] = 0
    speckleImg[np.where(np.isnan(speckleImg))] = 0
    speckleImg[np.where(speckleImg < 0)] = 0
    speckleImg[np.where(speckleImg > 255)] = 255
171

Amir Mohammadi's avatar
Amir Mohammadi committed
172
    # speckleImg contains the specular-component
173

Amir Mohammadi's avatar
Amir Mohammadi committed
174 175 176
    if len(speckleImg.shape) == 3:
        speckleImg = speckleImg[0]
    speckleImg = speckleImg.clip(min=0)
177

Amir Mohammadi's avatar
Amir Mohammadi committed
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
    speckleMean = np.mean(speckleImg)
    # factors 1.5 and 4.0 are proposed by Wen et al. in their paper and matlab
    # code.
    lowSpeckleThresh = speckleMean * 1.5
    hiSpeckleThresh = speckleMean * 4.0
    specklePixels = speckleImg[
        np.where(
            np.logical_and(
                speckleImg >= lowSpeckleThresh,
                speckleImg < hiSpeckleThresh))]

    # percentage of specular pixels in image
    r = float(specklePixels.flatten().shape[
              0]) / (speckleImg.shape[0] * speckleImg.shape[1])
    m = np.mean(specklePixels)  # mean-specularity (of specular-pixels)
    s = np.std(specklePixels)  # std. of specularity (of specular-pixels)

    # scaling by factor of 150 is as done by Wen et al. in their matlab code.
    return np.asarray((r, m / 150.0, s / 150.0), dtype=np.float32)


def marzilianoBlur(image):
    """Method proposed by Marziliano et al. for determining the average width
    of vertical edges, as a measure of blurredness in an image. (Reimplemented
    from the Matlab code provided by MSU.)

    :param image: 2D gray-level (face) image
    :param regionMask: (optional) 2D matrix (binary image), where 1s mark the
        pixels belonging to a region of interest, and 0s indicate pixels
        outside ROI.
    """
209

Amir Mohammadi's avatar
Amir Mohammadi committed
210 211 212
    assert len(
        image.shape) == 2, 'marzilianoBlur():: input image should be ' \
                           'a 2D array (gray level image)'
213

Amir Mohammadi's avatar
Amir Mohammadi committed
214 215
    # compute vertical edge-map of image using sobel
    edgeMap = sobelEdgeMap(image, 'vertical')
216

Amir Mohammadi's avatar
Amir Mohammadi committed
217 218 219 220 221
    # There will be some difference between the result of this function and the
    # Matlab version, because the edgeMap produced by sobelEdgeMap() is not
    # exactly the same as that produced by Matlab's edge() function. Test edge-
    # map generated in Matlab produces the same result as the matlab version of
    # MarzilianoBlur().
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
222 223

    blurImg = image
Amir Mohammadi's avatar
Amir Mohammadi committed
224 225 226 227
    C = blurImg.shape[1]  # number of cols in image
    # row, col contain the indices of the pixels that comprise edge-map.
    (row, col) = edgeMap.nonzero()

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
228 229 230
    blurMetric = 0
    if len(row) > 0:

Amir Mohammadi's avatar
Amir Mohammadi committed
231 232 233 234 235 236 237 238 239 240 241
        # to make the following code work in a similar fashion to the original
        # matlab code, sort the cols in ascending order, and sort the rows
        # according to the cols.
        #     ind = np.lexsort((row,col))
        #     row = row[ind]
        #     col = col[ind]
        # print('lexsort_col: %d' % (1+col))
        # print('lexsort_row: %d' % (1+row))
        # This was only used for debugging (to compare with Matlab code). In
        # fact it is not necessary, so it is commented out.

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
242
        edgeWidths = np.zeros_like(row, dtype=int)
Amir Mohammadi's avatar
Amir Mohammadi committed
243

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
244 245 246
        for i in range(len(row)):
            rEdge = row[i]
            cEdge = col[i]
Amir Mohammadi's avatar
Amir Mohammadi committed
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
            # instead of setting them to 'inf' as in MSU's Matlab version
            cStart = 0
            cEnd = 0

            # we want to estimate the edge-width, which will be cEnd - cStart.

            # search for start of edge in horizontal direction
            if cEdge > 0:  # i.e., edge is not on the left-border
                # 2.1: search left of current pixel (backwards)
                if blurImg[
                        rEdge,
                        cEdge] > blurImg[
                        rEdge,
                        cEdge -
                        1]:  # edge corresponds to a local peak; estimate left-
                        # extent of peak
                    j = cEdge - 1
                    while j > 0 and blurImg[rEdge, j] >= blurImg[rEdge, j - 1]:
                        j -= 1
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
266
                    cStart = j
Amir Mohammadi's avatar
Amir Mohammadi committed
267 268 269 270 271
                else:  # edge corresponds to a local valley; determine left-
                    # extent of valley
                    j = cEdge - 1
                    while j > 0 and blurImg[rEdge, j] <= blurImg[rEdge, j - 1]:
                        j -= 1
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
272
                    cStart = j
Amir Mohammadi's avatar
Amir Mohammadi committed
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289

            # search for end of edge in horizontal direction
            cEnd = C - 1  # initialize to right-border of image -- the max.
            # possible position for cEnd
            if cEdge < C - 1:
                if blurImg[
                        rEdge,
                        cEdge] > blurImg[
                        rEdge,
                        cEdge +
                        1]:
                    # edge corresponds to a local peak; estimate right-extent
                    # of peak
                    j = cEdge + 1
                    while j < C - 1 and blurImg[rEdge,
                                                j] >= blurImg[rEdge, j + 1]:
                        j += 1
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
290
                    cEnd = j
Amir Mohammadi's avatar
Amir Mohammadi committed
291 292 293 294 295 296 297
                else:
                    # edge corresponds to a local valley; determine right-
                    # extent of valley
                    j = cEdge + 1
                    while j < C - 1 and blurImg[rEdge,
                                                j] <= blurImg[rEdge, j + 1]:
                        j += 1
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
298
                    cEnd = j
Amir Mohammadi's avatar
Amir Mohammadi committed
299

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
300
            edgeWidths[i] = cEnd - cStart
Amir Mohammadi's avatar
Amir Mohammadi committed
301 302 303 304 305 306

            # sanity-check (edgeWidths should not have negative values)
            negSum = np.sum(edgeWidths[np.where(edgeWidths < 0)])
            assert negSum == 0, 'marzilianoBlur():: edgeWidths[] contains ' \
                                'negative values. YOU CANNOT BE SERIOUS!'

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
307 308
        # Final metric computation
        blurMetric = np.mean(edgeWidths)
Amir Mohammadi's avatar
Amir Mohammadi committed
309 310 311 312 313

    # compute histogram of edgeWidths ...(later)
        # binnum = 100;
        # t = ((1:binnum) - .5) .* C ./ binnum;
        # whist = hist(width_array, t) ./ length(width_array);
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
314 315 316

    return blurMetric

Amir Mohammadi's avatar
Amir Mohammadi committed
317 318 319 320
def calmoment(channel, regionMask=None):
    """ returns the first 3 statistical moments (mean, standard-dev., skewness)
    and 2 other first-order statistical measures of input image
    :param channel: 2D array containing gray-image-like data
321 322
    """

Amir Mohammadi's avatar
Amir Mohammadi committed
323 324 325
    assert len(
        channel.shape) == 2, 'calmoment():: channel should be ' \
                             'a 2D array (a single color-channel)'
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
326 327 328

    t = np.arange(0.05, 1.05, 0.05) + 0.025                 # t = 0.05:0.05:1;

Amir Mohammadi's avatar
Amir Mohammadi committed
329 330
    # pixnum = length(channel(:));
    nPix = np.prod(channel.shape)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
331
    m = np.mean(channel)                            # m = mean(channel(:));
Amir Mohammadi's avatar
Amir Mohammadi committed
332 333 334 335 336
    # d = sqrt(sum((channel(:) - m) .^ 2) / pixnum);
    d = np.std(channel)
    # s = sum(((channel(:) - m) ./ d) .^ 3) / pixnum;
    s = np.sum(np.power(((channel - m) / d), 3)) / nPix
    # print(t)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
337
    myHH = np.histogram(channel, t)[0]
Amir Mohammadi's avatar
Amir Mohammadi committed
338 339 340 341 342 343 344
    # print(myHH)
    # hh = hist(channel(:),t) / pixnum;
    hh = myHH.astype(float) / nPix

    # H = np.array([m,d,s, np.sum(hh[0:1]), np.sum(hh[-2:-1])])  # H = [m d s
    # sum(hh(1:2)) sum(hh(end-1:end))];
    H = np.array([m, d, s])
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
345
    s0 = np.sum(hh[0:2])
Amir Mohammadi's avatar
Amir Mohammadi committed
346 347
    # print(s0)
    H = np.hstack((H, s0))
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
348
    s1 = np.sum(hh[-2:])
Amir Mohammadi's avatar
Amir Mohammadi committed
349 350 351
    # print(s1)
    H = np.hstack((H, s1))

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
352 353 354
    return H


Amir Mohammadi's avatar
Amir Mohammadi committed
355 356 357 358 359 360 361 362
def calmoment_shift(channel):
    assert len(
        channel.shape) == 2, 'calmoment_shift():: channel should be a ' \
                             '2D array (a single color-channel)'

    channel = channel + 0.5
    channel[np.where(channel > 1.0)] -= 1.0

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
363 364
    H = calmoment(channel)

Amir Mohammadi's avatar
Amir Mohammadi committed
365
    return H
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
366 367 368


def calColorHist(image, m=100):
Amir Mohammadi's avatar
Amir Mohammadi committed
369 370 371 372 373 374 375 376 377 378
    """
    function returns the top 'm' most popular colors in the input image
        :param image: RGB color-image represented in a 3D array
        :param m: integer specifying how many 'top' colors to be counted (e.g.
            for m=10 the function will return the pixel-counts for the top 10
            most popular colors in image)
        :return cHist: counts of the top 100 most popular colors in image
        :return numClrs: total number of distinct colors in image
    """
    # 1. compute color histogram of image (normalized, if specified)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
379
    numBins = 32
Amir Mohammadi's avatar
Amir Mohammadi committed
380
    maxval = 255
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
381 382
    cHist = rgbhist(image, maxval, numBins, 1)

Amir Mohammadi's avatar
Amir Mohammadi committed
383 384 385
    # 2. determine top 100 colors of image from histogram
    y = sorted(cHist, reverse=True)  # [Y, I] = sort(H,'descend');
    cHist = y[0:m]                    # H = Y(1:m)';
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
386 387

    c = np.cumsum(y)                # C = cumsum(Y);
Amir Mohammadi's avatar
Amir Mohammadi committed
388 389
    numClrs = np.where(c > 0.99)[0][0]    # clrnum = find(C>.99,1,'first') - 1;

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
390 391 392 393 394 395 396
    cHist = np.array(cHist)
    return cHist, numClrs


'''
computes 3d color histogram of image
'''
Amir Mohammadi's avatar
Amir Mohammadi committed
397 398


Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
399
def rgbhist(image, maxval, nBins, normType=0):
Amir Mohammadi's avatar
Amir Mohammadi committed
400 401 402 403 404 405 406 407 408 409 410 411 412
    assert len(image.shape) == 3, 'image should be a 3D (rgb) array of shape' \
        ' (3, m,n) where m is no. of rows, and n is no. if cols in image.c$'
    assert normType > -1 and normType < 2, 'rgbhist():: normType should ' \
        ' be only 0 or 1'
    # zeros([nBins nBins nBins]);
    H = np.zeros((nBins, nBins, nBins), dtype=np.uint32)

    decimator = (maxval + 1) / nBins
    numPix = image.shape[1] * image.shape[2]
    # im = reshape(I,[size(I,1)*size(I,2) 3]);
    im = image.reshape(3, numPix).copy()
    im = im.T

413 414 415 416 417 418 419 420 421 422 423
    p = np.floor(im.astype(float) / decimator).astype(np.uint32)
    # in future versions of numpy (1.13 and above) you can replace this with:
    # unique_p, count = np.unique(p, return_counts=True, axis=0)
    # the following lines were taken from: https://stackoverflow.com/a/16973510
    p2 = np.ascontiguousarray(p).view(
        np.dtype((np.void, p.dtype.itemsize * p.shape[1])))
    unique_p, count = np.unique(p2, return_counts=True)
    unique_p = unique_p.view(p.dtype).reshape(-1, p.shape[1])
    # till here
    H[unique_p[:, 0], unique_p[:, 1], unique_p[:, 2]] = count

Amir Mohammadi's avatar
Amir Mohammadi committed
424
    H = H.ravel()  # H = H(:);
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
425
    # Un-Normalized histogram
Amir Mohammadi's avatar
Amir Mohammadi committed
426 427 428

    if normType == 1:
        H = H.astype(np.float32) / np.sum(H)    # l1 normalization
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
429 430 431 432
    return H


def blurriness(image):
Amir Mohammadi's avatar
Amir Mohammadi committed
433 434 435 436 437 438 439 440 441 442 443 444 445
    """
    function to estimate blurriness of an image, as computed by Di Wen et al.
    in their IEEE-TIFS-2015 paper.
        :param image: a gray-level image
    """

    assert len(
        image.shape) == 2, 'Input to blurriness() function should ' \
                           'be a 2D (gray) image'

    d = 4
    fsize = 2 * d + 1
    kver = np.ones((1, fsize)) / fsize
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
446
    khor = kver.T
Amir Mohammadi's avatar
Amir Mohammadi committed
447 448 449 450 451 452 453 454 455 456 457 458 459 460

    Bver = ssg.convolve2d(
        image.astype(
            np.float32), kver.astype(
            np.float32), mode='same')
    Bhor = ssg.convolve2d(
        image.astype(
            np.float32), khor.astype(
            np.float32), mode='same')

    # implementations of DFver and DFhor below don't look the same as in the
    # Matlab code, but the following implementation produces equivalent
    # results. there might be a bug in Matlab! The 2 commented statements above
    # would correspond to the intent of the Matlab code.
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
461
    DFver = np.diff(image.astype('int16'), axis=0)
Amir Mohammadi's avatar
Amir Mohammadi committed
462
    DFver[np.where(DFver < 0)] = 0
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
463 464

    DFhor = np.diff(image.astype('int16'), axis=1)
Amir Mohammadi's avatar
Amir Mohammadi committed
465
    DFhor[np.where(DFhor < 0)] = 0
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
466

Amir Mohammadi's avatar
Amir Mohammadi committed
467 468
    DBver = np.abs(np.diff(Bver, axis=0))
    DBhor = np.abs(np.diff(Bhor, axis=1))
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
469 470 471

    Vver = DFver.astype(float) - DBver.astype(float)
    Vhor = DFhor.astype(float) - DBhor.astype(float)
Amir Mohammadi's avatar
Amir Mohammadi committed
472 473
    Vver[Vver < 0] = 0  # Vver(find(Vver<0)) = 0;
    Vhor[Vhor < 0] = 0  # Vhor(find(Vhor<0)) = 0;
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
474 475

    SFver = np.sum(DFver)
Amir Mohammadi's avatar
Amir Mohammadi committed
476
    SFhor = np.sum(DFhor)  # sum(DFhor(:));
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
477

Amir Mohammadi's avatar
Amir Mohammadi committed
478 479
    SVver = np.sum(Vver)  # sum(Vver(:));
    SVhor = np.sum(Vhor)  # sum(Vhor(:));
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
480

Amir Mohammadi's avatar
Amir Mohammadi committed
481 482
    BFver = (SFver - SVver) / SFver
    BFhor = (SFhor - SVhor) / SFhor
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
483

Amir Mohammadi's avatar
Amir Mohammadi committed
484
    blurF = max(BFver, BFhor)  # max([BFver BFhor]);
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
485

Amir Mohammadi's avatar
Amir Mohammadi committed
486
    return blurF