msu_iqa_features.py 16 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
12
from . import galbally_iqm_features as iqm
from . import tan_specular_highlights as tsh
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
13

Amir Mohammadi's avatar
Amir Mohammadi committed
14
15
''' Utility functions '''

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
16
17
18

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

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

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
24
25
    hsv = np.zeros_like(rgbImage)
    bob.ip.color.rgb_to_hsv(rgbImage, hsv)
Amir Mohammadi's avatar
Amir Mohammadi committed
26
27
28
29
    h = hsv[0, :, :]
    s = hsv[1, :, :]
    v = hsv[2, :, :]
#
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
30
31
32
33
34
35
    return (h, s, v)


def imshow(image):
    import matplotlib
    from matplotlib import pyplot as plt
Amir Mohammadi's avatar
Amir Mohammadi committed
36
37
38
    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
39
        outImg = image.tolist()
40
        print(len(outImg))
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
41
42
        result = np.dstack((outImg[0], outImg[1]))
        outImg = np.dstack((result, outImg[2]))
Amir Mohammadi's avatar
Amir Mohammadi committed
43
44
45
        # [:,:,1], cmap=mpl.cm.gray)
        plt.imshow((outImg * 255.0).astype(np.uint8))

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

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
51
52
53
    plt.show()


Amir Mohammadi's avatar
Amir Mohammadi committed
54
55
'''Auxilliary functions'''

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
56
57

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

    # 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
70
    if orientation is 'horizontal':
Amir Mohammadi's avatar
Amir Mohammadi committed
71
        refEdge = iqm.edge_thinning(refSobelX[:, :], refSobelX[:, :], thinning)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
72
73
    else:
        if orientation is 'vertical':
Amir Mohammadi's avatar
Amir Mohammadi committed
74
75
76
77
            refEdge = iqm.edge_thinning(
                refSobelY[
                    :, :], refSobelY[
                    :, :], thinning)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
78
        else:
Amir Mohammadi's avatar
Amir Mohammadi committed
79
80
81
82
            refEdge = iqm.edge_thinning(
                refSobelX[
                    :, :], refSobelY[
                    :, :], thinning)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
83
84
85

    return refEdge

86

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

    Parameters:

Amir Mohammadi's avatar
Amir Mohammadi committed
93
94
95
    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).

96
97
    Returns:

Amir Mohammadi's avatar
Amir Mohammadi committed
98
99
100
101
102
103
    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.
104
    """
Amir Mohammadi's avatar
Amir Mohammadi committed
105
106
107
108
109
110
111
    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
112
113
    grayImage = np.zeros_like(h, dtype='uint8')
    bob.ip.color.rgb_to_gray(rgbImage, grayImage)
Amir Mohammadi's avatar
Amir Mohammadi committed
114

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

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

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

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
125
126
127
    # 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
128

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

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

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

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

    return fv


151
def compute_iqa_specularity_features(rgbImage, startEps=0.05):
Amir Mohammadi's avatar
Amir Mohammadi committed
152
153
154
    """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
155
156
    """

Amir Mohammadi's avatar
Amir Mohammadi committed
157
158
159
160
    # separate the specular and diffuse components of input color image.
    speckleFreeImg, diffuseImg, speckleImg = tsh.remove_highlights(
        rgbImage.astype(float), startEps, verboseFlag=False)
    # speckleImg contains the specular-component
161

Amir Mohammadi's avatar
Amir Mohammadi committed
162
163
164
    if len(speckleImg.shape) == 3:
        speckleImg = speckleImg[0]
    speckleImg = speckleImg.clip(min=0)
165

Amir Mohammadi's avatar
Amir Mohammadi committed
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
    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.
    """
197

Amir Mohammadi's avatar
Amir Mohammadi committed
198
199
200
    assert len(
        image.shape) == 2, 'marzilianoBlur():: input image should be ' \
                           'a 2D array (gray level image)'
201

Amir Mohammadi's avatar
Amir Mohammadi committed
202
203
    # compute vertical edge-map of image using sobel
    edgeMap = sobelEdgeMap(image, 'vertical')
204

Amir Mohammadi's avatar
Amir Mohammadi committed
205
206
207
208
209
    # 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
210
211

    blurImg = image
Amir Mohammadi's avatar
Amir Mohammadi committed
212
213
214
215
    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
216
217
218
    blurMetric = 0
    if len(row) > 0:

Amir Mohammadi's avatar
Amir Mohammadi committed
219
220
221
222
223
224
225
226
227
228
229
        # 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
230
        edgeWidths = np.zeros_like(row, dtype=int)
Amir Mohammadi's avatar
Amir Mohammadi committed
231

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
232
233
234
        for i in range(len(row)):
            rEdge = row[i]
            cEdge = col[i]
Amir Mohammadi's avatar
Amir Mohammadi committed
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
            # 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
254
                    cStart = j
Amir Mohammadi's avatar
Amir Mohammadi committed
255
256
257
258
259
                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
260
                    cStart = j
Amir Mohammadi's avatar
Amir Mohammadi committed
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277

            # 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
278
                    cEnd = j
Amir Mohammadi's avatar
Amir Mohammadi committed
279
280
281
282
283
284
285
                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
286
                    cEnd = j
Amir Mohammadi's avatar
Amir Mohammadi committed
287

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
288
            edgeWidths[i] = cEnd - cStart
Amir Mohammadi's avatar
Amir Mohammadi committed
289
290
291
292
293
294

            # 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
295
296
        # Final metric computation
        blurMetric = np.mean(edgeWidths)
Amir Mohammadi's avatar
Amir Mohammadi committed
297
298
299
300
301

    # 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
302
303
304
305

    return blurMetric


Amir Mohammadi's avatar
Amir Mohammadi committed
306
307
308
309
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
310
311
    """

Amir Mohammadi's avatar
Amir Mohammadi committed
312
313
314
    assert len(
        channel.shape) == 2, 'calmoment():: channel should be ' \
                             'a 2D array (a single color-channel)'
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
315
316
317

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

Amir Mohammadi's avatar
Amir Mohammadi committed
318
319
    # pixnum = length(channel(:));
    nPix = np.prod(channel.shape)
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
320
    m = np.mean(channel)                            # m = mean(channel(:));
Amir Mohammadi's avatar
Amir Mohammadi committed
321
322
323
324
325
    # 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
326
    myHH = np.histogram(channel, t)[0]
Amir Mohammadi's avatar
Amir Mohammadi committed
327
328
329
330
331
332
333
    # 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
334
    s0 = np.sum(hh[0:2])
Amir Mohammadi's avatar
Amir Mohammadi committed
335
336
    # print(s0)
    H = np.hstack((H, s0))
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
337
    s1 = np.sum(hh[-2:])
Amir Mohammadi's avatar
Amir Mohammadi committed
338
339
340
    # print(s1)
    H = np.hstack((H, s1))

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
341
342
343
    return H


Amir Mohammadi's avatar
Amir Mohammadi committed
344
345
346
347
348
349
350
351
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
352
353
    H = calmoment(channel)

Amir Mohammadi's avatar
Amir Mohammadi committed
354
    return H
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
355
356
357


def calColorHist(image, m=100):
Amir Mohammadi's avatar
Amir Mohammadi committed
358
359
360
361
362
363
364
365
366
367
    """
    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
368
    numBins = 32
Amir Mohammadi's avatar
Amir Mohammadi committed
369
    maxval = 255
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
370
371
    cHist = rgbhist(image, maxval, numBins, 1)

Amir Mohammadi's avatar
Amir Mohammadi committed
372
373
374
    # 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
375
376

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

Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
379
380
381
382
383
384
385
    cHist = np.array(cHist)
    return cHist, numClrs


'''
computes 3d color histogram of image
'''
Amir Mohammadi's avatar
Amir Mohammadi committed
386
387


Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
388
def rgbhist(image, maxval, nBins, normType=0):
Amir Mohammadi's avatar
Amir Mohammadi committed
389
390
391
392
393
394
395
396
397
398
399
400
401
    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

402
403
404
405
406
407
408
409
410
411
412
    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
413
    H = H.ravel()  # H = H(:);
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
414
    # Un-Normalized histogram
Amir Mohammadi's avatar
Amir Mohammadi committed
415
416
417

    if normType == 1:
        H = H.astype(np.float32) / np.sum(H)    # l1 normalization
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
418
419
420
421
    return H


def blurriness(image):
Amir Mohammadi's avatar
Amir Mohammadi committed
422
423
424
425
426
427
428
429
430
431
432
433
434
    """
    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
435
    khor = kver.T
Amir Mohammadi's avatar
Amir Mohammadi committed
436
437
438
439
440
441
442
443
444
445
446
447
448
449

    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
450
    DFver = np.diff(image.astype('int16'), axis=0)
Amir Mohammadi's avatar
Amir Mohammadi committed
451
    DFver[np.where(DFver < 0)] = 0
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
452
453

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

Amir Mohammadi's avatar
Amir Mohammadi committed
456
457
    DBver = np.abs(np.diff(Bver, axis=0))
    DBhor = np.abs(np.diff(Bhor, axis=1))
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
458
459
460

    Vver = DFver.astype(float) - DBver.astype(float)
    Vhor = DFhor.astype(float) - DBhor.astype(float)
Amir Mohammadi's avatar
Amir Mohammadi committed
461
462
    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
463
464

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

Amir Mohammadi's avatar
Amir Mohammadi committed
467
468
    SVver = np.sum(Vver)  # sum(Vver(:));
    SVhor = np.sum(Vhor)  # sum(Vhor(:));
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
469

Amir Mohammadi's avatar
Amir Mohammadi committed
470
471
    BFver = (SFver - SVver) / SFver
    BFhor = (SFhor - SVhor) / SFhor
Sushil BHATTACHARJEE's avatar
Sushil BHATTACHARJEE committed
472

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

Amir Mohammadi's avatar
Amir Mohammadi committed
475
    return blurF