BEAT - sbhatta/iqm

This algorithm is a legacy one. The API has changed since its implementation. New versions and forks will need to be updated.

This algorithm is splittable

Endpoint Groups 1

Algorithms have at least one input and one output. All algorithm endpoints are organized in groups. Groups are used by the platform to indicate which inputs and outputs are synchronized together. The first group is automatically synchronized with the channel defined by the block in which the algorithm is deployed.

Group: main

Endpoint Name	Data Format	Nature
video	system/array_4d_uint8/1	Input
iqm_feature_set	system/array_2d_floats/1	Output

import numpy as np
import scipy as sp
import scipy.signal as ssg
import scipy.ndimage.filters as snf

import bob
import bob.io.base
import bob.io.video
import bob.ip.color
import bob.ip.base

'''
functions for computing image-quality features
'''

"""
Main function to be called, to extract a set of image quality-features computed for the input image
@param image: 2d numpy array. Should contain input image of size [M,N] (i.e. M rows x N cols).
@return featSet: a tuple of float-scalars, each representing one image-quality measure.
"""
def computeQualityFeatures(image, gwin):
    
    smoothed = ssg.convolve2d(image, gwin, boundary='symm', mode='same')

return imageQualityMeasures(image, smoothed)

"""
actually computes similarities between the two input images, but also returns some descriptors of the reference-image that are independent of any other image.
Returns a tuple of 18 values, each of which is a float-scalar. 
The quality measures computed in this function correspond to the Image-quality features discussed in Galbally et al., 2014.
"""
def imageQualityMeasures(refImage, testImage):
    assert len(refImage.shape)==2, "refImage should be a 2D array"
    assert len(testImage.shape)==2, "testImage should be a 2D array"
    assert (refImage.shape[0] == testImage.shape[0]), "The two images should have the same width"
    assert (refImage.shape[1] == testImage.shape[1]), "The two images should have the same height"
    
    diffImg = refImage.astype(np.float) - testImage.astype(np.float) 
    diffSq = np.square(diffImg)
    sumDiffSq = np.sum(diffSq)
    absDiffImg = np.absolute(diffImg)
    
    refSq = np.square(refImage.astype(np.float))
    sumRefSq = np.sum(refSq)
    
    numPx = refImage.shape[0]*refImage.shape[1] #number of pixels in each image
    maxPxVal = 255.0;
    
    #1 MSE
    mse = float(sumDiffSq)/float(numPx)
    
    #2 PSNR
    psnr = np.inf
    if mse>0:
        psnr = 10.0*np.log10(maxPxVal*maxPxVal/mse)
    
    #3 AD: Average difference
    ad = float(np.sum(diffImg))/float(numPx)
    
    #4 SC: structural content
    testSq = np.square(testImage.astype(np.float))
    sumTestSq = np.sum(testSq)
    sc=np.inf
    if sumTestSq>0:
        sc = float(sumRefSq)/float(sumTestSq)
    
    #5 NK: normalized cross-correlation
    imgProd = refImage * testImage # element-wise product
    nk = float(np.sum(imgProd))/float(sumRefSq)
    
    #6 MD: Maximum difference
    md = float(np.amax(absDiffImg))
    
    #7 LMSE: Laplacian MSE
    #scipy implementation of laplacian is different from Matlab's version, especially at the image-borders
    # To significant differences between scipy...laplace and Matlab's del2() are:
    #    a. Matlab del2() divides the convolution result by 4, so the ratio (scipy.laplace() result)/(del2()-result) is 4
    #    b. Matlab does a different kind of processing at the boundaries, so the results at the boundaries are different in the 2 calls.
    #In Galbally's Matlab code, there is a factor of 4, which I have dropped (no difference in result), 
    #because this is implicit in scipy.ndimage.filters.laplace()
    op = snf.laplace(refImage, mode='reflect') #mode can be 'wrap', 'reflect', 'nearest', 'mirror', or ['constant' with a specified value]
    opSq = np.square(op)
    sum_opSq = np.sum(opSq)
    tmp1 = (op - (snf.laplace(testImage, mode='reflect')))
    num_op = np.square(tmp1)
    lmse = float(np.sum(num_op))/float(sum_opSq)
    
    #8 NAE: normalized abs. error
    sumRef = np.sum(np.absolute(refImage))
    nae = float(np.sum(absDiffImg))/float(sumRef)
    
    #9 SNRv: SNR in db
    snrv = 10.0*np.log10(float(sumRefSq)/float(sumDiffSq))
    
    #10 RAMDv: R-averaged max diff (r=10)
    #implementation below is same as what Galbally does in Matlab
    r=10
    sorted = np.sort(diffImg.flatten())[::-1] #the [::-1] flips the sorted vector, so that it is in descending order
    topsum = np.sum(sorted[0:r])
    ramdv = np.sqrt(float(topsum)/float(r))
    
    #11,12: MAS: Mean Angle Similarity, MAMS: Mean Angle-Magnitude Similarity
    mas, mams = angleSimilarity(refImage, testImage, diffImg) 
    
    fftRef = np.fft.fft2(refImage)
    fftTest = np.fft.fft2(testImage)
      
    #13, 14: SME: spectral magnitude error; SPE: spectral phase error
    sme, spe = spectralSimilarity(refImage, testImage) 
    
    #15 TED: Total edge difference
    #ted = edgeSimilarity(refImage, testImage)

#16 TCD: Total corner difference
    #tcd = cornerSimilarity(refImage, testImage)
    
    #17, 18: GME: gradient-magnitude error; GPE: gradient phase error  
    gme, gpe = gradientSimilarity(refImage, testImage)
    
    #19 SSIM
    ssim, _ = SSIM(refImage, testImage)
    
    #20 VIF
    vif = VIFp(refImage, testImage)
    
    #21,22,23,24,25: , RRED, BIQI, JQI, NIQE: these parameters are not implemented in this code.
    
    #26 HLFI: high-low frequency index (implemented as done by Galbally in Matlab).
    hlfi=HLFreq_Index(fftRef, refImage.shape[1])
    
    return (mse, psnr, ad, sc, nk, md, lmse, nae, snrv, ramdv, mas, mams, sme, gme, gpe, ssim, vif, hlfi)

"""
SSIM: Structural Similarity index between two gray-level images. The dynamic range is assumed to be 0..255.
Ref:Z. Wang, A.C. Bovik, H.R. Sheikh and E.P. Simoncelli: 
    "Image Quality Assessment: From error measurement to Structural Similarity"
    IEEE Trans. on Image Processing, 13(1), 01/2004
    @param refImage: 2D numpy array (reference image)
    @param testImage: 2D numpy array (test image)
    Both input images should have the same dimensions. This is assumed, and not verified in this function
    @return ssim: float-scalar. The mean structural similarity between the 2 input images.
    @return ssim_map: the SSIM index map of the test image (this map is smaller than the test image).
"""
def SSIM(refImage, testImage):
    M=refImage.shape[0]
    N=refImage.shape[1]
        
    winSz=11 #window size for gaussian filter
    winSgm = 1.5 # sigma for gaussian filter
    
    #input image should be at least 11x11 in size.
    if(M<winSz) or (N<winSz):
        ssim_index = -np.inf
        ssim_map = -np.inf
        
        return ssim_index, ssim_map
    
    #construct the gaussian filter
    gwin = gauss_2d((winSz, winSz), winSgm)
    
    K1 = 0.01 # constants taken from the initial matlab implementation provided by Bovik's lab.
    K2 = 0.03
    L = 255 #dynamic range.
    
    C1 = (K1*L)*(K1*L)
    C2 = (K2*L)*(K2*L)
    
    #ssg is scipy.signal
    mu1 = ssg.convolve2d(refImage, gwin, mode='valid')
    mu2 = ssg.convolve2d(testImage, gwin, mode='valid')
    
    mu1Sq = mu1*mu1
    mu2Sq = mu2*mu2
    mu1_mu2 = mu1*mu2
    
    sigma1_sq = ssg.convolve2d((refImage*refImage), gwin, mode='valid') - mu1Sq
    sigma2_sq = ssg.convolve2d((testImage*testImage), gwin, mode='valid') - mu1Sq
    sigma12 = ssg.convolve2d((refImage*testImage), gwin, mode='valid') - mu1_mu2
    
    assert (C1>0 and C2 > 0), "Conditions for computing ssim with this code are not met. Set the Ks and L to values > 0."
    num1 = (2.0*mu1_mu2 + C1) *(2.0*sigma12 + C2) 
    den1 = (mu1Sq + mu2Sq+C1)*(sigma1_sq + sigma2_sq +C2)
    ssim_map = num1/den1

ssim = np.average(ssim_map)
    
    return ssim, ssim_map

"""
VIF: Visual Information Fidelity measure.
Ref: H.R. Sheikh and A.C. Bovik: "Image Information and Visual Quality", IEEE Trans. Image Processing.
Adapted from Galbally's matlab code, which was provided by Bovik et al's LIVE lab.
    @param refImage: 2D numpy array (reference image)
    @param testImage: 2D numpy array (test image)
    Both input images should have the same dimensions. This is assumed, and not verified in this function
    @return vifp: float-scalar. Measure of visual information fidelity between the 2 input images
"""
def VIFp(refImage, testImage):
    sigma_nsq = 2.0
    num=0
    den=0
    
    #sc is scale, taking values (1,2,3,4)
    for sc in range(1,5):
        N=(2**(4-sc+1))+1
        win = gauss_2d((N,N), (float(N)/5.0))
        
        #ssg is scipy.signal
        
        if sc > 1 :
            refImage = ssg.convolve2d(refImage, win, mode='valid')
            testImage = ssg.convolve2d(testImage, win, mode='valid')
            refImage = refImage[::2, ::2]           #downsample by factor 2 in each direction
            testImage = testImage[::2, ::2]
        
        mu1 = ssg.convolve2d(refImage, win, mode='valid')
        mu2 = ssg.convolve2d(testImage, win, mode='valid')
        mu1Sq = mu1*mu1
        mu2Sq = mu2*mu2
        mu1_mu2 = mu1*mu2
        
        sigma1_sq = ssg.convolve2d((refImage*refImage), win, mode='valid') - mu1Sq
        sigma2_sq = ssg.convolve2d((testImage*testImage), win, mode='valid') - mu2Sq
        sigma12 = ssg.convolve2d((refImage*testImage), win, mode='valid') - mu1_mu2
        
        sigma1_sq[sigma1_sq<0]=0        #set negative filter responses to 0.
        sigma2_sq[sigma2_sq<0]=0
        
        g = sigma12 / (sigma1_sq + 1e-10)
        sv_sq = sigma2_sq - g*sigma12;
        
        g[(sigma1_sq < 1e-10)]=0
        sv_sq[sigma1_sq < 1e-10] = sigma2_sq[sigma1_sq < 1e-10]
        sigma1_sq[sigma1_sq<1e-10]=0
        
        g[(sigma2_sq < 1e-10)]=0
        sv_sq[sigma2_sq < 1e-10] =0
        
        sv_sq[g<0]=sigma2_sq[g<0]
        g[g<0]=0
        sv_sq[sv_sq <= 1e-10] = 1e-10     #sic. As implemented in the original matlab version...
        
        m1 = g*g*sigma1_sq
        m2 = sv_sq+sigma_nsq
        m3 = np.log10(1 + m1/m2)
        
        m4 = np.log10(1 + (sigma1_sq/sigma_nsq))
        
        num += np.sum(m3)
        den += np.sum(m4)
    
    vifp = num/den
    return vifp

"""
HLFI: relative difference between high- and low-frequency energy in image.
Ref: I. Avcibas, N. Memon, B. Sankur: "Steganalysis using image quality metrics", IEEE Trans. Image Processing, 12, 2003.
@param imgFFT: 2D numpy array of complex numbers, representing Fourier transform of test image.
@param ncols: int. Number of columns in image.
@return float-scalar.
"""
def HLFreq_Index(imgFFT, ncols):

N= ncols
    colHalf = int(round(N/2))  #round it up
    freqSel = 0.15
    
    freqCol = round(freqSel*N)
    lowFreqColHalf = int(round(freqCol/2.0)) 
    
    fftRes = imgFFT 
    fftMag = np.abs(fftRes)
    totalEnergy = np.sum(fftMag)

lowIdx = colHalf-lowFreqColHalf
    hiIdx = colHalf + lowFreqColHalf

LowFreqMag = fftMag[:, lowIdx:hiIdx]
    lowFreqMagTotal = np.sum(LowFreqMag)
    
    fftMag[:, lowIdx:hiIdx]=0
    highFreqMagTotal = np.sum(fftMag)

highLowFreqIQ = np.abs(lowFreqMagTotal - highFreqMagTotal)/float(totalEnergy)
    
    return highLowFreqIQ

"""
Image similarity based on gradient. Computes the mean phase and magnitude difference of gradient between input reference and test images.
Ref: I. Avcibas, N. Memon, B. Sankur: "Steganalysis using image quality metrics", IEEE Trans. Image Processing, 12, 2003.
    @param refImage: 2D numpy array (reference image)
    @param testImage: 2D numpy array (test image)
    Both input images should have the same dimensions. This is assumed, and not verified in this function.
    @return difGradMag: float-scalar. Mean difference in gradient-magnitude.
    @return difGradPhase: float-scalar. Mean difference in gradient-phase.
"""
def gradientSimilarity(refImage, testImage):
    
    numPx = refImage.shape[0]*refImage.shape[1] # we assume that testImage is of the same shape as refImage
    
    #compute gradient (a la matlab) for reference image
    refGrad=np.gradient(refImage,5,5) #5: spacing of 5 pixels between 2 sites of grad. evaluation.
    
    refReal = refGrad[0]
    refImag = refGrad[1]
    refGradComplex = refReal + 1j*refImag
    
    refMag=np.abs(refGradComplex)
    refPhase = np.arctan2(refImag, refReal)

#compute gradient for test image
    testGrad=np.gradient(testImage,5) #5: spacing of 5 pixels between 2 sites of grad. evaluation. It applies to both dims.
    testReal = testGrad[0]
    testImag = testGrad[1]
    testGradComplex = testReal + 1j*testImag
    
    testMag=np.abs(testGradComplex)
    testPhase = np.arctan2(testImag, testReal)    
    
    absPhaseDiff = np.abs(refPhase - testPhase) 
    difGradPhase = (np.sum(absPhaseDiff))/float(numPx)
    
    absMagDiff = np.abs(refMag - testMag) 
    difGradMag = float(np.sum(absMagDiff))/float(numPx)
    
    return difGradMag, difGradPhase

"""
find local maxima using 3x3 mask.
Used in cornerSimilarity()
Should produce results very similar to matlabs imregionalmax()
@param img: 2d numpy array. Image-like, containing a 'cornerness'-index for every pixel.
@return regmax: 2d numpy array. Binary image showing corners (which are regions of local maxima in input image).
"""
def regionalmax(img):
    h = img.shape[0]
    w = img.shape[1]
    
    #extend input image borders by repeating border-values
    b = np.zeros((img.shape[0]+2,img.shape[1]+2))
    b[1:-1,1:-1] = img
    b[0,:]=b[1,:]
    b[:,0]=b[:,1]
    b[-1,:]=b[-2,:]
    b[:,-1]=b[:,-2]
    
    
    regmax = np.zeros((h,w), dtype='uint8') #will contain the output bitmap showing local maxima.
    
    for i in range(1, h+1):
        for j in range(1, w+1):
            subim = b[i-1:i+2, j-1:j+2]
            lmax = np.amax(subim)
            lmin = np.amin(subim)
            if b[i,j] == lmax and b[i,j]>lmin : 
                regmax[i-1,j-1] = 1
                        
    for i in range(1,h):
        for j in range(w):
            if regmax[i,j]==1:
                imin=i-1
                if imin<0: imin=0
                imax=i+2
                if imax>h: imax=h
                for k in range(imin, imax):
                    jmin=j-1
                    if jmin<0: jmin=0
                    jmax = j+2
                    if jmax>w: jmax = w
                    for l in range(jmin, jmax):
                        if(img[k,l]==img[i,j]):
                            regmax[k,l]=1
    
    return regmax

"""
returns a 'cornerness' image, where each pixel-value specifies the 'degree of cornerness' of the corresponding pixel in input image
The 'cornerness' image is of size (N-2, M-2) for an input image of size (N,M) (no cornerness computed for the border pixel of input.)
@param image: 2d numpy array. Input image for which cornerness needs to be computed.
@return cornerness: 2d numpy array giving a 'cornerness'-value for the input image.
"""
def cornerMetric(image):
    image = image.astype(np.float)
   
    sensitivity_factor = 0.4
    gwin = gauss_2d((5,5), 1.5)
    
    vfilt = np.array([-1,0,1], ndmin=2)
    hfilt = vfilt.T
    A = ssg.convolve2d(image, vfilt, boundary='symm', mode='same') 
    B = ssg.convolve2d(image, hfilt, boundary='symm', mode='same')
    #crop out the valid portions of the filter-response (same size for both A and B)
    A = A[1:-2, 1:-2]
    B = B[1:-2, 1:-2]
    
    #compute products of A, B, C
    C = A*B
    A = A*A
    B = B*B
    
    #filter A, B, and C
    A = ssg.convolve2d(A, gwin, boundary='symm', mode='valid')
    B = ssg.convolve2d(B, gwin, boundary='symm', mode='valid')
    C = ssg.convolve2d(C, gwin, boundary='symm', mode='valid')

ABsum = A + B
    cornerness = (A*B) - (C*C) - sensitivity_factor *(ABsum*ABsum)
    
    return cornerness

'''
compute the corner-based similarity between 2 images (how close are the numbers of corners found in the two images?).
returns an index between 0 and 1. The smaller the better.
@param refImage: 2D numpy array (reference image)
@param testImage: 2D numpy array (test image)
@return float-scalar.
'''
def cornerSimilarity(refImage, testImage):
    
    C = cornerMetric(refImage)
    C_peaks = regionalmax(C)
    
    #imshow(C_peaks)
    
    CG = cornerMetric(testImage)
    CG_peaks = regionalmax(CG)
    
    nCornersRef = np.sum(C_peaks)
    nCornersTest = np.sum(CG_peaks)
    #print 'CornerSim::', nCornersRef, nCornersTest
    
    maxCorners = max(nCornersRef, nCornersTest)
    
    qCornersDiff = np.fabs(nCornersRef - nCornersTest)/float(maxCorners)
    
    return qCornersDiff

"""
Similarity between the edge-maps of the two input images.
Ref: I. Avcibas, N. Memon, B. Sankur: "Steganalysis using image quality metrics", IEEE Trans. Image Processing, 12, 2003.
@param refImage: 2D numpy array (reference image)
@param testImage: 2D numpy array (test image)
@return float-scalar
"""
def edgeSimilarity(refImage, testImage):
    
    #bob..sobel returns filter-responses which need to be thresholded to get the edge-map
    thinning=1
    refImage=refImage.astype(np.float)
    
    #compute edge map for reference image
    refSobel_sep = bob.ip.base.sobel(refImage) #returns 3D image. 1st dim is the edge-direction. 1st component is vertical; 2nd component is hor. responses
    refSobelX = refSobel_sep[0,:,:]
    refSobelY = refSobel_sep[1,:,:]
    refEdge = edgeThinning(refSobelX[:,:], refSobelY[:,:], thinning)
    
    
    #compute edge map for test image
    testSobel_sep = bob.ip.base.sobel(testImage)
    testSobelX = testSobel_sep[0,:,:]
    testSobelY = testSobel_sep[1,:,:]

testEdge = edgeThinning(testSobelX[:,:], testSobelY[:,:], thinning)

numPx = refImage.shape[0]*refImage.shape[1]
    numRefEdgePx = np.sum(refEdge)
    numTestEdgePx = np.sum(testEdge)
    qEdgeD = np.abs(numRefEdgePx - numTestEdgePx)/float(numPx)

return qEdgeD

"""
    function to perform edge-thining in the same way as done in Matlab. Called in edgeSimilarity()
    Returns a binary edge-map (uint8 image).
    @param  bx: vertical edge-filter responses (for example, response of 1 of the two Sobel filters)
    @param  by: horizontal edge-filter responses 
    @param  thinning: [0,1]. Default:1, implies 'do edge-thinning'. If set to 0, no edge-thinning is done.
    bx and by should be of the same shape
"""
def edgeThinning(bx, by, thinning=1):
    assert(len(bx.shape)==2) and (len(by.shape)==2), "bx and by should be 2D arrays."
    assert(bx.shape[0]==by.shape[0]) and (bx.shape[1]==by.shape[1]), "bx and by should have the same shape."
    m = bx.shape[0]
    n = by.shape[1]
    e = np.zeros([m,n], dtype=np.uint8)     # will contain the resulting edge-map.
    
#     print 'bx', bx.shape
#     print 'by', by.shape
    #compute the edge-strength from the 2 directional filter-responses
    b = np.sqrt(bx*bx + by*by)
    
    #compute the threshold a la Matlab (as described in "Digital Image Processing" book by W.K. Pratt.
    scale=4
    cutoff = scale*np.mean(b)
    #thresh = np.sqrt(cutoff)
    
    myEps = np.spacing(1)*100.0 #np.spacing(1) is the same as eps in matlab.
    #compute the edge-map a la Matlab
    for r in range(m):
        for c in range(n):
            if thinning:
                if r<0 or r>(m-1) or (c-1)<0 :
                    b1=True
                else:
                    b1=(b[r,c-1] < b[r,c])
                
                if(r<0) or  r>(m-1)  or  (c+1)>(n-1) :
                    b2=True
                else:
                    b2 = (b[r,c] > b[r,c+1])
                
                if (c<0) or  c>(n-1)  or  (r-1)<0 :
                    b3=True
                else:
                    b3=(b[r,c] > b[r-1,c])
                
                if(c<1) or  c>(n-1)  or  (r+1)>(m-1) :
                    b4=True
                else:
                    b4=(b[r,c]> b[r+1, c])
                
                c1 = (b[r,c]>cutoff)
                c2 = ((bx[r,c]>= (by[r,c]-myEps)) & b1 & b2)
                c3 = ((by[r,c]>= (bx[r,c]-myEps)) & b3 & b4)

e[r,c] = c1 & (c2 | c3)
            else:
                e[r,c] = (b[r,c]>cutoff)
    
    return e

"""
@param refImage: 2D numpy array (reference image)
@param testImage: 2D numpy array (test image)
@return sme: float-scalar. Mean difference in magnitudes of spectra of the two images.
@return spe: float-scalar. Mean difference in phases of spectra of the two images.
"""
#def spectralSimilarity(fftRef, fftTest, numPx):
def spectralSimilarity(refImage, testImage):
    
    #assume that ref and test images have the same shape
    rows=refImage.shape[0]
    cols=refImage.shape[1]
    numPx = rows*cols
    fftRef = np.fft.rfft2(refImage)
    fftTest = np.fft.rfft2(testImage)
    
    refMag=np.abs(fftRef)
    testMag=np.abs(fftTest)
    absMagDiff = np.abs(refMag - testMag)
    #SME: spectral magnitude error
    sme = np.sum(absMagDiff * absMagDiff)/float(numPx)
    
    #SPE: spectral phase error
    refPhase = np.angle(fftRef)
    testPhase = np.angle(fftTest)
    absPhaseDiff = np.abs(refPhase - testPhase)
    spe = np.sum(absPhaseDiff * absPhaseDiff)/float(numPx)

return sme, spe

"""
Cosine-Similarity between the the rows of the two input images.
Ref: I. Avcibas, N. Memon, B. Sankur: "Steganalysis using image quality metrics", IEEE Trans. Image Processing, 12, 2003.
@param refImage: 2D numpy array (reference image)
@param testImage: 2D numpy array (test image)
@param diffImage: 2D numpy array. Difference between refImage and testImage. Not strictly necessary as input but precomputed here, to save computation time.
@return mas: float-scalar. Mean angle-similarity.
@return mams: float-scalar. Mean angle-magnitude similarity.
"""
def angleSimilarity(refImage, testImage, diffImage):
    mas=None
    mams=None
    
    numPx = refImage.shape[0]*refImage.shape[1]
    
    refNorm = np.linalg.norm(refImage, axis=1)
    testNorm = np.linalg.norm(testImage, axis=1)
    thetaVec = np.zeros([refImage.shape[0], 1])
    diffNorm = np.linalg.norm(diffImage, axis=1)
    magnitVec = diffNorm/255.0 #Galbally divides by sqrt(255**2)
    magnitVec = np.reshape(magnitVec, (refImage.shape[0], 1))
    
    for i in range(refImage.shape[0]):
        refR = refImage[i,:]
        testR = testImage[i,:]
        cosTheta = np.dot(refR, testR)/(refNorm[i]*testNorm[i])
        if(cosTheta < -1.0): cosTheta = -1.0
        if(cosTheta > 1.0): cosTheta = 1.0
        theta = np.arccos(cosTheta)
        thetaVec[i]=theta
    
    tmp2 = thetaVec*2.0/np.pi
    
    #MAS: mean angle similarity
    mas = 1.0 -( sum(tmp2) / float(numPx) )
    
    tmp3 = 1.0 - tmp2 
    tmp4 = 1.0 - magnitVec    
    chi = 1.0 - (tmp3 * tmp4)
    #MAMS: mean angle-magnitude similarity 
    mams = sum(chi)/float(numPx)

return (float(mas), float(mams))

'''
Returns a 2D gaussian-filter matrix
equivalent of matlab fspecial('gaussian'...)

@param shape: tuple defining the size of the desired filter in each dimension. Elements of tuple should be 2 positive odd-integers.
@param sigma: float-scalar
@return h: 2D numpy array. Contains weights for 2D gaussian filter.
'''
def gauss_2d(shape=(3, 3), sigma=0.5): 
    """ 
    2D gaussian mask - should give the same result as MATLAB's 
    fspecial('gaussian',[shape],[sigma]) 
    """ 
    m, n = [(ss-1.)/2. for ss in shape] 
    y, x = np.ogrid[-m:m+1, -n:n+1] 
    h = np.exp(-(x*x + y*y) / (2.*sigma*sigma)) 
    h[h < np.finfo(h.dtype).eps*h.max()] = 0 
    sumh = h.sum() 
    if sumh != 0: 
        h /= sumh 
    return h

'''
Matlab-like RGB to gray...
    @param: rgbImage : numpy array for the form: [3,h,w] where h is the height of the image and w is the width of the image.
    Returns a y-image in floating-point format (range [(16/255) .. (235/255)])
'''
def matlab_rgb2gray(rgbImage):
    #g1 = 0.299*rgbFrame[0,:,:] + 0.587*rgbFrame[1,:,:] + 0.114*rgbFrame[2,:,:] #standard coeffs in CCIR601
    
    #this is how it's done in matlab...
    rgbImage = rgbImage / 255.0
    C0 = 65.481/255.0
    C1 = 128.553/255.0
    C2 = 24.966/255.0
    scaleMin = 16.0/255.0
    gray = scaleMin + (C0*rgbImage[0,:,:] + C1*rgbImage[1,:,:] + C2*rgbImage[2,:,:])

return gray

'''
computes image-quality features for a set of frames comprising a video.
    @param video4d: a '4d' video (N frames, each frame representing an r-g-b image).
    returns  a set of feature-vectors, 1 vector per frame of video4d
'''
def computeVideoIQM(video4d):
    numframes = video4d.shape[0]
    
    gwin = gauss_2d((3,3), 0.5)
    
    #process first frame separately, to get the no. of iqm features
    f=0
    rgbFrame = video4d[f,:,:,:]
    grayFrame = matlab_rgb2gray(rgbFrame) #compute gray-level image for input color-frame
    iqmSet = computeQualityFeatures(grayFrame, gwin)
    numIQMs = len(iqmSet)
    #now initialize fset to store iqm features for all frames of input video.
    fset = np.zeros([numframes, numIQMs])
    bobQFeats = np.asarray(iqmSet) # computeQualityFeatures() returns a tuple
    fset[f,:] = bobQFeats
    
    for f in range(1,numframes):
        rgbFrame = video4d[f,:,:,:]
        grayFrame = matlab_rgb2gray(rgbFrame) #compute gray-level image for input color-frame
        bobQFeats = np.asarray(computeQualityFeatures(grayFrame, gwin)) # computeQualityFeatures() returns a tuple
        fset[f,:] = bobQFeats
    
    return fset

class Algorithm:

def __init__(self):
        pass

def process(self, inputs, outputs):

video = inputs["video"].data.value
        iqmFeats = computeVideoIQM(video)

outputs["iqm_feature_set"].write({
            'value': iqmFeats
        })

return True

xxxxxxxxxx
 
import numpy as np
import scipy as sp
import scipy.signal as ssg
import scipy.ndimage.filters as snf
​
import bob
import bob.io.base
import bob.io.video
import bob.ip.color
import bob.ip.base
​
'''
functions for computing image-quality features
'''
​
"""
Main function to be called, to extract a set of image quality-features computed for the input image
@param image: 2d numpy array. Should contain input image of size [M,N] (i.e. M rows x N cols).
@return featSet: a tuple of float-scalars, each representing one image-quality measure.
"""
def computeQualityFeatures(image, gwin):
    
    smoothed = ssg.convolve2d(image, gwin, boundary='symm', mode='same') 
​
    return imageQualityMeasures(image, smoothed)
​
​
"""
actually computes similarities between the two input images, but also returns some descriptors of the reference-image that are independent of any other image.
Returns a tuple of 18 values, each of which is a float-scalar. 
The quality measures computed in this function correspond to the Image-quality features discussed in Galbally et al., 2014.
"""
def imageQualityMeasures(refImage, testImage):
    assert len(refImage.shape)==2, "refImage should be a 2D array"
    assert len(testImage.shape)==2, "testImage should be a 2D array"
    assert (refImage.shape[0] == testImage.shape[0]), "The two images should have the same width"
    assert (refImage.shape[1] == testImage.shape[1]), "The two images should have the same height"
    
    diffImg = refImage.astype(np.float) - testImage.astype(np.float) 
    diffSq = np.square(diffImg)
    sumDiffSq = np.sum(diffSq)
    absDiffImg = np.absolute(diffImg)
    
    refSq = np.square(refImage.astype(np.float))
    sumRefSq = np.sum(refSq)
    
    numPx = refImage.shape[0]*refImage.shape[1] #number of pixels in each image
    maxPxVal = 255.0;
    
    #1 MSE
    mse = float(sumDiffSq)/float(numPx)
    
    #2 PSNR
    psnr = np.inf
    if mse>0:
        psnr = 10.0*np.log10(maxPxVal*maxPxVal/mse)
    
    #3 AD: Average difference
    ad = float(np.sum(diffImg))/float(numPx)
    
    #4 SC: structural content
    testSq = np.square(testImage.astype(np.float))
    sumTestSq = np.sum(testSq)
    sc=np.inf
    if sumTestSq>0:
        sc = float(sumRefSq)/float(sumTestSq)
    
    #5 NK: normalized cross-correlation
    imgProd = refImage * testImage # element-wise product
    nk = float(np.sum(imgProd))/float(sumRefSq)
    
    #6 MD: Maximum difference
    md = float(np.amax(absDiffImg))
    
    #7 LMSE: Laplacian MSE
    #scipy implementation of laplacian is different from Matlab's version, especially at the image-borders
    # To significant differences between scipy...laplace and Matlab's del2() are:
    #    a. Matlab del2() divides the convolution result by 4, so the ratio (scipy.laplace() result)/(del2()-result) is 4
    #    b. Matlab does a different kind of processing at the boundaries, so the results at the boundaries are different in the 2 calls.
    #In Galbally's Matlab code, there is a factor of 4, which I have dropped (no difference in result), 
    #because this is implicit in scipy.ndimage.filters.laplace()
    op = snf.laplace(refImage, mode='reflect') #mode can be 'wrap', 'reflect', 'nearest', 'mirror', or ['constant' with a specified value]
    opSq = np.square(op)
    sum_opSq = np.sum(opSq)
    tmp1 = (op - (snf.laplace(testImage, mode='reflect')))
    num_op = np.square(tmp1)
    lmse = float(np.sum(num_op))/float(sum_opSq)
    
    #8 NAE: normalized abs. error
    sumRef = np.sum(np.absolute(refImage))
    nae = float(np.sum(absDiffImg))/float(sumRef)
    
    #9 SNRv: SNR in db
    snrv = 10.0*np.log10(float(sumRefSq)/float(sumDiffSq))
    
    #10 RAMDv: R-averaged max diff (r=10)
    #implementation below is same as what Galbally does in Matlab
    r=10
    sorted = np.sort(diffImg.flatten())[::-1] #the [::-1] flips the sorted vector, so that it is in descending order
    topsum = np.sum(sorted[0:r])
    ramdv = np.sqrt(float(topsum)/float(r))
    
    #11,12: MAS: Mean Angle Similarity, MAMS: Mean Angle-Magnitude Similarity
    mas, mams = angleSimilarity(refImage, testImage, diffImg) 
    
    fftRef = np.fft.fft2(refImage)
    fftTest = np.fft.fft2(testImage)
      
    #13, 14: SME: spectral magnitude error; SPE: spectral phase error
    sme, spe = spectralSimilarity(refImage, testImage) 
    
    #15 TED: Total edge difference
    #ted = edgeSimilarity(refImage, testImage)
​
    #16 TCD: Total corner difference
    #tcd = cornerSimilarity(refImage, testImage)
    
    #17, 18: GME: gradient-magnitude error; GPE: gradient phase error  
    gme, gpe = gradientSimilarity(refImage, testImage)
    
    #19 SSIM
    ssim, _ = SSIM(refImage, testImage)
    
    #20 VIF
    vif = VIFp(refImage, testImage)
    
    #21,22,23,24,25: , RRED, BIQI, JQI, NIQE: these parameters are not implemented in this code.
    
    #26 HLFI: high-low frequency index (implemented as done by Galbally in Matlab).
    hlfi=HLFreq_Index(fftRef, refImage.shape[1])
    
    return (mse, psnr, ad, sc, nk, md, lmse, nae, snrv, ramdv, mas, mams, sme, gme, gpe, ssim, vif, hlfi)
​
​
"""
SSIM: Structural Similarity index between two gray-level images. The dynamic range is assumed to be 0..255.
Ref:Z. Wang, A.C. Bovik, H.R. Sheikh and E.P. Simoncelli: 
    "Image Quality Assessment: From error measurement to Structural Similarity"
    IEEE Trans. on Image Processing, 13(1), 01/2004
    @param refImage: 2D numpy array (reference image)
    @param testImage: 2D numpy array (test image)
    Both input images should have the same dimensions. This is assumed, and not verified in this function
    @return ssim: float-scalar. The mean structural similarity between the 2 input images.
    @return ssim_map: the SSIM index map of the test image (this map is smaller than the test image).
"""
def SSIM(refImage, testImage):
    M=refImage.shape[0]
    N=refImage.shape[1]
        
    winSz=11 #window size for gaussian filter
    winSgm = 1.5 # sigma for gaussian filter
    
    #input image should be at least 11x11 in size.
    if(M<winSz) or (N<winSz):
        ssim_index = -np.inf
        ssim_map = -np.inf
        
        return ssim_index, ssim_map
    
    #construct the gaussian filter
    gwin = gauss_2d((winSz, winSz), winSgm)
    
    K1 = 0.01 # constants taken from the initial matlab implementation provided by Bovik's lab.
    K2 = 0.03
    L = 255 #dynamic range.
    
    C1 = (K1*L)*(K1*L)
    C2 = (K2*L)*(K2*L)
    
    #ssg is scipy.signal
    mu1 = ssg.convolve2d(refImage, gwin, mode='valid')
    mu2 = ssg.convolve2d(testImage, gwin, mode='valid')
    
    mu1Sq = mu1*mu1
    mu2Sq = mu2*mu2
    mu1_mu2 = mu1*mu2
    
    sigma1_sq = ssg.convolve2d((refImage*refImage), gwin, mode='valid') - mu1Sq
    sigma2_sq = ssg.convolve2d((testImage*testImage), gwin, mode='valid') - mu1Sq
    sigma12 = ssg.convolve2d((refImage*testImage), gwin, mode='valid') - mu1_mu2
    
    assert (C1>0 and C2 > 0), "Conditions for computing ssim with this code are not met. Set the Ks and L to values > 0."
    num1 = (2.0*mu1_mu2 + C1) *(2.0*sigma12 + C2) 
    den1 = (mu1Sq + mu2Sq+C1)*(sigma1_sq + sigma2_sq +C2)
    ssim_map = num1/den1
​
    ssim = np.average(ssim_map)
    
    return ssim, ssim_map
​
​
"""
VIF: Visual Information Fidelity measure.
Ref: H.R. Sheikh and A.C. Bovik: "Image Information and Visual Quality", IEEE Trans. Image Processing.
Adapted from Galbally's matlab code, which was provided by Bovik et al's LIVE lab.
    @param refImage: 2D numpy array (reference image)
    @param testImage: 2D numpy array (test image)
    Both input images should have the same dimensions. This is assumed, and not verified in this function
    @return vifp: float-scalar. Measure of visual information fidelity between the 2 input images
"""
def VIFp(refImage, testImage):
    sigma_nsq = 2.0
    num=0
    den=0
    
    #sc is scale, taking values (1,2,3,4)
    for sc in range(1,5):
        N=(2**(4-sc+1))+1
        win = gauss_2d((N,N), (float(N)/5.0))
        
        #ssg is scipy.signal
        
        if sc > 1 :
            refImage = ssg.convolve2d(refImage, win, mode='valid')
            testImage = ssg.convolve2d(testImage, win, mode='valid')
            refImage = refImage[::2, ::2]           #downsample by factor 2 in each direction
            testImage = testImage[::2, ::2]
        
        mu1 = ssg.convolve2d(refImage, win, mode='valid')
        mu2 = ssg.convolve2d(testImage, win, mode='valid')
        mu1Sq = mu1*mu1
        mu2Sq = mu2*mu2
        mu1_mu2 = mu1*mu2
        
        sigma1_sq = ssg.convolve2d((refImage*refImage), win, mode='valid') - mu1Sq
        sigma2_sq = ssg.convolve2d((testImage*testImage), win, mode='valid') - mu2Sq
        sigma12 = ssg.convolve2d((refImage*testImage), win, mode='valid') - mu1_mu2
        
        sigma1_sq[sigma1_sq<0]=0        #set negative filter responses to 0.
        sigma2_sq[sigma2_sq<0]=0
        
        g = sigma12 / (sigma1_sq + 1e-10)
        sv_sq = sigma2_sq - g*sigma12;
        
        g[(sigma1_sq < 1e-10)]=0
        sv_sq[sigma1_sq < 1e-10] = sigma2_sq[sigma1_sq < 1e-10]
        sigma1_sq[sigma1_sq<1e-10]=0
        
        g[(sigma2_sq < 1e-10)]=0
        sv_sq[sigma2_sq < 1e-10] =0
        
        sv_sq[g<0]=sigma2_sq[g<0]
        g[g<0]=0
        sv_sq[sv_sq <= 1e-10] = 1e-10     #sic. As implemented in the original matlab version...
        
        m1 = g*g*sigma1_sq
        m2 = sv_sq+sigma_nsq
        m3 = np.log10(1 + m1/m2)
        
        m4 = np.log10(1 + (sigma1_sq/sigma_nsq))
        
        num += np.sum(m3)
        den += np.sum(m4)
    
    vifp = num/den
    return vifp
​
​
"""
HLFI: relative difference between high- and low-frequency energy in image.
Ref: I. Avcibas, N. Memon, B. Sankur: "Steganalysis using image quality metrics", IEEE Trans. Image Processing, 12, 2003.
@param imgFFT: 2D numpy array of complex numbers, representing Fourier transform of test image.
@param ncols: int. Number of columns in image.
@return float-scalar.
"""
def HLFreq_Index(imgFFT, ncols):
​
    N= ncols
    colHalf = int(round(N/2))  #round it up
    freqSel = 0.15
    
    freqCol = round(freqSel*N)
    lowFreqColHalf = int(round(freqCol/2.0)) 
    
    fftRes = imgFFT 
    fftMag = np.abs(fftRes)
    totalEnergy = np.sum(fftMag)
​
    
    lowIdx = colHalf-lowFreqColHalf
    hiIdx = colHalf + lowFreqColHalf
​
    LowFreqMag = fftMag[:, lowIdx:hiIdx]
    lowFreqMagTotal = np.sum(LowFreqMag)
    
    fftMag[:, lowIdx:hiIdx]=0
    highFreqMagTotal = np.sum(fftMag)
​
    
    highLowFreqIQ = np.abs(lowFreqMagTotal - highFreqMagTotal)/float(totalEnergy)
    
    return highLowFreqIQ
    
        
     
​
"""
Image similarity based on gradient. Computes the mean phase and magnitude difference of gradient between input reference and test images.
Ref: I. Avcibas, N. Memon, B. Sankur: "Steganalysis using image quality metrics", IEEE Trans. Image Processing, 12, 2003.
    @param refImage: 2D numpy array (reference image)
    @param testImage: 2D numpy array (test image)
    Both input images should have the same dimensions. This is assumed, and not verified in this function.
    @return difGradMag: float-scalar. Mean difference in gradient-magnitude.
    @return difGradPhase: float-scalar. Mean difference in gradient-phase.
"""
def gradientSimilarity(refImage, testImage):
    
    numPx = refImage.shape[0]*refImage.shape[1] # we assume that testImage is of the same shape as refImage
    
    #compute gradient (a la matlab) for reference image
    refGrad=np.gradient(refImage,5,5) #5: spacing of 5 pixels between 2 sites of grad. evaluation.
    
    refReal = refGrad[0]
    refImag = refGrad[1]
    refGradComplex = refReal + 1j*refImag
    
    refMag=np.abs(refGradComplex)
    refPhase = np.arctan2(refImag, refReal)
​
    #compute gradient for test image
    testGrad=np.gradient(testImage,5) #5: spacing of 5 pixels between 2 sites of grad. evaluation. It applies to both dims.
    testReal = testGrad[0]
    testImag = testGrad[1]
    testGradComplex = testReal + 1j*testImag
    
    testMag=np.abs(testGradComplex)
    testPhase = np.arctan2(testImag, testReal)    
    
    absPhaseDiff = np.abs(refPhase - testPhase) 
    difGradPhase = (np.sum(absPhaseDiff))/float(numPx)
    
    absMagDiff = np.abs(refMag - testMag) 
    difGradMag = float(np.sum(absMagDiff))/float(numPx)
    
    return difGradMag, difGradPhase
​
​
"""
find local maxima using 3x3 mask.
Used in cornerSimilarity()
Should produce results very similar to matlabs imregionalmax()
@param img: 2d numpy array. Image-like, containing a 'cornerness'-index for every pixel.
@return regmax: 2d numpy array. Binary image showing corners (which are regions of local maxima in input image).
"""
def regionalmax(img):
    h = img.shape[0]
    w = img.shape[1]
    
    #extend input image borders by repeating border-values
    b = np.zeros((img.shape[0]+2,img.shape[1]+2))
    b[1:-1,1:-1] = img
    b[0,:]=b[1,:]
    b[:,0]=b[:,1]
    b[-1,:]=b[-2,:]
    b[:,-1]=b[:,-2]
    
    
    regmax = np.zeros((h,w), dtype='uint8') #will contain the output bitmap showing local maxima.
    
    for i in range(1, h+1):
        for j in range(1, w+1):
            subim = b[i-1:i+2, j-1:j+2]
            lmax = np.amax(subim)
            lmin = np.amin(subim)
            if b[i,j] == lmax and b[i,j]>lmin : 
                regmax[i-1,j-1] = 1
                        
    for i in range(1,h):
        for j in range(w):
            if regmax[i,j]==1:
                imin=i-1
                if imin<0: imin=0
                imax=i+2
                if imax>h: imax=h
                for k in range(imin, imax):
                    jmin=j-1
                    if jmin<0: jmin=0
                    jmax = j+2
                    if jmax>w: jmax = w
                    for l in range(jmin, jmax):
                        if(img[k,l]==img[i,j]):
                            regmax[k,l]=1
    
    return regmax
​
​
​
"""
returns a 'cornerness' image, where each pixel-value specifies the 'degree of cornerness' of the corresponding pixel in input image
The 'cornerness' image is of size (N-2, M-2) for an input image of size (N,M) (no cornerness computed for the border pixel of input.)
@param image: 2d numpy array. Input image for which cornerness needs to be computed.
@return cornerness: 2d numpy array giving a 'cornerness'-value for the input image.
"""
def cornerMetric(image):
    image = image.astype(np.float)
   
    sensitivity_factor = 0.4
    gwin = gauss_2d((5,5), 1.5)
    
    vfilt = np.array([-1,0,1], ndmin=2)
    hfilt = vfilt.T
    A = ssg.convolve2d(image, vfilt, boundary='symm', mode='same') 
    B = ssg.convolve2d(image, hfilt, boundary='symm', mode='same')
    #crop out the valid portions of the filter-response (same size for both A and B)
    A = A[1:-2, 1:-2]
    B = B[1:-2, 1:-2]
    
    #compute products of A, B, C
    C = A*B
    A = A*A
    B = B*B
    
    #filter A, B, and C
    A = ssg.convolve2d(A, gwin, boundary='symm', mode='valid')
    B = ssg.convolve2d(B, gwin, boundary='symm', mode='valid')
    C = ssg.convolve2d(C, gwin, boundary='symm', mode='valid')
​
    ABsum = A + B
    cornerness = (A*B) - (C*C) - sensitivity_factor *(ABsum*ABsum)
    
    return cornerness
​
​
'''
compute the corner-based similarity between 2 images (how close are the numbers of corners found in the two images?).
returns an index between 0 and 1. The smaller the better.
@param refImage: 2D numpy array (reference image)
@param testImage: 2D numpy array (test image)
@return float-scalar.
'''
def cornerSimilarity(refImage, testImage):
    
    C = cornerMetric(refImage)
    C_peaks = regionalmax(C)
    
    #imshow(C_peaks)
    
    CG = cornerMetric(testImage)
    CG_peaks = regionalmax(CG)
    
    nCornersRef = np.sum(C_peaks)
    nCornersTest = np.sum(CG_peaks)
    #print 'CornerSim::', nCornersRef, nCornersTest
    
    maxCorners = max(nCornersRef, nCornersTest)
    
    qCornersDiff = np.fabs(nCornersRef - nCornersTest)/float(maxCorners)
    
    return qCornersDiff
​
​
"""
Similarity between the edge-maps of the two input images.
Ref: I. Avcibas, N. Memon, B. Sankur: "Steganalysis using image quality metrics", IEEE Trans. Image Processing, 12, 2003.
@param refImage: 2D numpy array (reference image)
@param testImage: 2D numpy array (test image)
@return float-scalar
"""
def edgeSimilarity(refImage, testImage):
    
    #bob..sobel returns filter-responses which need to be thresholded to get the edge-map
    thinning=1
    refImage=refImage.astype(np.float)
    
    #compute edge map for reference image
    refSobel_sep = bob.ip.base.sobel(refImage) #returns 3D image. 1st dim is the edge-direction. 1st component is vertical; 2nd component is hor. responses
    refSobelX = refSobel_sep[0,:,:]
    refSobelY = refSobel_sep[1,:,:]
    refEdge = edgeThinning(refSobelX[:,:], refSobelY[:,:], thinning)
    
    
    #compute edge map for test image
    testSobel_sep = bob.ip.base.sobel(testImage)
    testSobelX = testSobel_sep[0,:,:]
    testSobelY = testSobel_sep[1,:,:]
​
    testEdge = edgeThinning(testSobelX[:,:], testSobelY[:,:], thinning)
​
    numPx = refImage.shape[0]*refImage.shape[1]
    numRefEdgePx = np.sum(refEdge)
    numTestEdgePx = np.sum(testEdge)
    qEdgeD = np.abs(numRefEdgePx - numTestEdgePx)/float(numPx)
​
    return qEdgeD
​
​
"""
    function to perform edge-thining in the same way as done in Matlab. Called in edgeSimilarity()
    Returns a binary edge-map (uint8 image).
    @param  bx: vertical edge-filter responses (for example, response of 1 of the two Sobel filters)
    @param  by: horizontal edge-filter responses 
    @param  thinning: [0,1]. Default:1, implies 'do edge-thinning'. If set to 0, no edge-thinning is done.
    bx and by should be of the same shape
"""
def edgeThinning(bx, by, thinning=1):
    assert(len(bx.shape)==2) and (len(by.shape)==2), "bx and by should be 2D arrays."
    assert(bx.shape[0]==by.shape[0]) and (bx.shape[1]==by.shape[1]), "bx and by should have the same shape."
    m = bx.shape[0]
    n = by.shape[1]
    e = np.zeros([m,n], dtype=np.uint8)     # will contain the resulting edge-map.
    
#     print 'bx', bx.shape
#     print 'by', by.shape
    #compute the edge-strength from the 2 directional filter-responses
    b = np.sqrt(bx*bx + by*by)
    
    #compute the threshold a la Matlab (as described in "Digital Image Processing" book by W.K. Pratt.
    scale=4
    cutoff = scale*np.mean(b)
    #thresh = np.sqrt(cutoff)
    
    myEps = np.spacing(1)*100.0 #np.spacing(1) is the same as eps in matlab.
    #compute the edge-map a la Matlab
    for r in range(m):
        for c in range(n):
            if thinning:
                if r<0 or r>(m-1) or (c-1)<0 :
                    b1=True
                else:
                    b1=(b[r,c-1] < b[r,c])
                
                if(r<0) or  r>(m-1)  or  (c+1)>(n-1) :
                    b2=True
                else:
                    b2 = (b[r,c] > b[r,c+1])
                
                if (c<0) or  c>(n-1)  or  (r-1)<0 :
                    b3=True
                else:
                    b3=(b[r,c] > b[r-1,c])
                
                if(c<1) or  c>(n-1)  or  (r+1)>(m-1) :
                    b4=True
                else:
                    b4=(b[r,c]> b[r+1, c])
                
                c1 = (b[r,c]>cutoff)
                c2 = ((bx[r,c]>= (by[r,c]-myEps)) & b1 & b2)
                c3 = ((by[r,c]>= (bx[r,c]-myEps)) & b3 & b4)
​
                e[r,c] = c1 & (c2 | c3)
            else:
                e[r,c] = (b[r,c]>cutoff)
    
    return e
​
​
"""
@param refImage: 2D numpy array (reference image)
@param testImage: 2D numpy array (test image)
@return sme: float-scalar. Mean difference in magnitudes of spectra of the two images.
@return spe: float-scalar. Mean difference in phases of spectra of the two images.
"""
#def spectralSimilarity(fftRef, fftTest, numPx):
def spectralSimilarity(refImage, testImage):
    
    #assume that ref and test images have the same shape
    rows=refImage.shape[0]
    cols=refImage.shape[1]
    numPx = rows*cols
    fftRef = np.fft.rfft2(refImage)
    fftTest = np.fft.rfft2(testImage)
    
    refMag=np.abs(fftRef)
    testMag=np.abs(fftTest)
    absMagDiff = np.abs(refMag - testMag)
    #SME: spectral magnitude error
    sme = np.sum(absMagDiff * absMagDiff)/float(numPx)
    
    #SPE: spectral phase error
    refPhase = np.angle(fftRef)
    testPhase = np.angle(fftTest)
    absPhaseDiff = np.abs(refPhase - testPhase)
    spe = np.sum(absPhaseDiff * absPhaseDiff)/float(numPx)
​
    return sme, spe
​
​
"""
Cosine-Similarity between the the rows of the two input images.
Ref: I. Avcibas, N. Memon, B. Sankur: "Steganalysis using image quality metrics", IEEE Trans. Image Processing, 12, 2003.
@param refImage: 2D numpy array (reference image)
@param testImage: 2D numpy array (test image)
@param diffImage: 2D numpy array. Difference between refImage and testImage. Not strictly necessary as input but precomputed here, to save computation time.
@return mas: float-scalar. Mean angle-similarity.
@return mams: float-scalar. Mean angle-magnitude similarity.
"""
def angleSimilarity(refImage, testImage, diffImage):
    mas=None
    mams=None
    
    numPx = refImage.shape[0]*refImage.shape[1]
    
    refNorm = np.linalg.norm(refImage, axis=1)
    testNorm = np.linalg.norm(testImage, axis=1)
    thetaVec = np.zeros([refImage.shape[0], 1])
    diffNorm = np.linalg.norm(diffImage, axis=1)
    magnitVec = diffNorm/255.0 #Galbally divides by sqrt(255**2)
    magnitVec = np.reshape(magnitVec, (refImage.shape[0], 1))
    
    for i in range(refImage.shape[0]):
        refR = refImage[i,:]
        testR = testImage[i,:]
        cosTheta = np.dot(refR, testR)/(refNorm[i]*testNorm[i])
        if(cosTheta < -1.0): cosTheta = -1.0
        if(cosTheta > 1.0): cosTheta = 1.0
        theta = np.arccos(cosTheta)
        thetaVec[i]=theta
    
    tmp2 = thetaVec*2.0/np.pi
    
    #MAS: mean angle similarity
    mas = 1.0 -( sum(tmp2) / float(numPx) )
    
    tmp3 = 1.0 - tmp2 
    tmp4 = 1.0 - magnitVec    
    chi = 1.0 - (tmp3 * tmp4)
    #MAMS: mean angle-magnitude similarity 
    mams = sum(chi)/float(numPx)
​
    return (float(mas), float(mams))
​
​
​
'''
Returns a 2D gaussian-filter matrix
equivalent of matlab fspecial('gaussian'...)
​
@param shape: tuple defining the size of the desired filter in each dimension. Elements of tuple should be 2 positive odd-integers.
@param sigma: float-scalar
@return h: 2D numpy array. Contains weights for 2D gaussian filter.
'''
def gauss_2d(shape=(3, 3), sigma=0.5): 
    """ 
    2D gaussian mask - should give the same result as MATLAB's 
    fspecial('gaussian',[shape],[sigma]) 
    """ 
    m, n = [(ss-1.)/2. for ss in shape] 
    y, x = np.ogrid[-m:m+1, -n:n+1] 
    h = np.exp(-(x*x + y*y) / (2.*sigma*sigma)) 
    h[h < np.finfo(h.dtype).eps*h.max()] = 0 
    sumh = h.sum() 
    if sumh != 0: 
        h /= sumh 
    return h 
​
​
'''
Matlab-like RGB to gray...
    @param: rgbImage : numpy array for the form: [3,h,w] where h is the height of the image and w is the width of the image.
    Returns a y-image in floating-point format (range [(16/255) .. (235/255)])
'''
def matlab_rgb2gray(rgbImage):
    #g1 = 0.299*rgbFrame[0,:,:] + 0.587*rgbFrame[1,:,:] + 0.114*rgbFrame[2,:,:] #standard coeffs in CCIR601
    
    #this is how it's done in matlab...
    rgbImage = rgbImage / 255.0
    C0 = 65.481/255.0
    C1 = 128.553/255.0
    C2 = 24.966/255.0
    scaleMin = 16.0/255.0
    gray = scaleMin + (C0*rgbImage[0,:,:] + C1*rgbImage[1,:,:] + C2*rgbImage[2,:,:])    
​
    return gray
​
​
'''
computes image-quality features for a set of frames comprising a video.
    @param video4d: a '4d' video (N frames, each frame representing an r-g-b image).
    returns  a set of feature-vectors, 1 vector per frame of video4d
'''
def computeVideoIQM(video4d):
    numframes = video4d.shape[0]
    
    gwin = gauss_2d((3,3), 0.5)
    
    #process first frame separately, to get the no. of iqm features
    f=0
    rgbFrame = video4d[f,:,:,:]
    grayFrame = matlab_rgb2gray(rgbFrame) #compute gray-level image for input color-frame
    iqmSet = computeQualityFeatures(grayFrame, gwin)
    numIQMs = len(iqmSet)
    #now initialize fset to store iqm features for all frames of input video.
    fset = np.zeros([numframes, numIQMs])
    bobQFeats = np.asarray(iqmSet) # computeQualityFeatures() returns a tuple
    fset[f,:] = bobQFeats
    
    for f in range(1,numframes):
        rgbFrame = video4d[f,:,:,:]
        grayFrame = matlab_rgb2gray(rgbFrame) #compute gray-level image for input color-frame
        bobQFeats = np.asarray(computeQualityFeatures(grayFrame, gwin)) # computeQualityFeatures() returns a tuple
        fset[f,:] = bobQFeats
    
    return fset
​
​
​
class Algorithm:
​
    def __init__(self):
        pass
​
    def process(self, inputs, outputs):
​
        video = inputs["video"].data.value
        iqmFeats = computeVideoIQM(video)
        
​
        outputs["iqm_feature_set"].write({
            'value': iqmFeats
        })  
​
        return True
​

The code for this algorithm in Python
The ruler at 80 columns indicate suggested POSIX line breaks (for readability).
The editor will automatically enlarge to accomodate the entirety of your input
Use keyboard shortcuts for search/replace and faster editing. For example, use Ctrl-F (PC) or Cmd-F (Mac) to search through this box

Given an input video consisting of N frames, this algorithm returns a Nxd floating-point array containing one d-dimensional feature-vector for each frame. The feature vector consists of various image-quality measures computed for the frame.

This algorithm relies on the Bob library.

Experiments

Attestation:

Privacy:

Status:

			Updated	Name	Databases/Protocols	Analyzers
				sbhatta/sbhatta/iqm-face-antispoofing-test/2/replay2-antispoofing-iqm-lda	replay/2@grandtest	sbhatta/iqm_spoof_eer_analyzer/9

This table shows the number of times this algorithm has been successfully run using the given environment. Note this does not provide sufficient information to evaluate if the algorithm will run when submitted to different conditions.

algorithms sbhatta iqm_features 1

Endpoint Groups 1

Group: main

Experiments

algorithms

sbhatta

iqm_features

1