initial commit

README.rst

+.. vim: set fileencoding=utf-8 :
+.. Sat  3 Dec 20:18:15 2016 CET
+
+.. image:: http://img.shields.io/badge/docs-stable-yellow.png
+   :target: http://pythonhosted.org/bob.ip.qualitymeasure/index.html
+.. image:: http://img.shields.io/badge/docs-latest-orange.png
+   :target: https://www.idiap.ch/software/bob/docs/latest/bob/bob.ip.qualitymeasure/master/index.html
+.. image:: https://gitlab.idiap.ch/bob/bob.ip.qualitymeasure/badges/master/build.svg
+   :target: https://gitlab.idiap.ch/bob/bob.ip.qualitymeasure/commits/master
+.. image:: https://img.shields.io/badge/gitlab-project-0000c0.svg
+   :target: https://gitlab.idiap.ch/bob/bob.ip.qualitymeasure
+.. image:: http://img.shields.io/pypi/v/bob.ip.qualitymeasure.png
+   :target: https://pypi.python.org/pypi/bob.ip.qualitymeasure
+.. image:: http://img.shields.io/pypi/dm/bob.ip.qualitymeasure.png
+   :target: https://pypi.python.org/pypi/bob.ip.qualitymeasure
+
+
+==================================================
+ Bob's library of image-quality feature-extractors
+==================================================
+
+This package is part of the signal-processing and machine learning toolbox
+Bob_. It provides functions for extracting image-quality features proposed
+for PAD experiments by different research group.
+
+
+Installation
+------------
+
+Follow our `installation`_ instructions. Then, using the Python interpreter
+provided by the distribution, bootstrap and buildout this package::
+
+  $ python bootstrap-buildout.py
+  $ ./bin/buildout
+
+
+Contact
+-------
+
+For questions or reporting issues to this software package, contact our
+development `mailing list`_.
+
+
+.. Place your references here:
+.. _bob: https://www.idiap.ch/software/bob
+.. _installation: https://gitlab.idiap.ch/bob/bob/wikis/Installation
+.. _mailing list: https://groups.google.com/forum/?fromgroups#!forum/bob-devel

bob/__init__.py

+# see https://docs.python.org/3/library/pkgutil.html
+from pkgutil import extend_path
+__path__ = extend_path(__path__, __name__)

bob/ip/__init__.py

+# see https://docs.python.org/3/library/pkgutil.html
+from pkgutil import extend_path
+__path__ = extend_path(__path__, __name__)

bob/ip/qualitymeasure/__init__.py

+#!/usr/bin/env python
+# vim: set fileencoding=utf-8 :
+
+
+
+def get_config():
+    """
+    Returns a string containing the configuration information.
+
+    """
+    import bob.extension
+    return bob.extension.get_config(__name__)
+
+
+
+# gets sphinx autodoc done right - don't remove it
+__all__ = [_ for _ in dir() if not _.startswith('_')]
+

bob/ip/qualitymeasure/compute_iqm_features.py

+#!/usr/bin/env python
+# vim: set fileencoding=utf-8 :
+
+'''
+Created on 13 Oct 2015
+
+@author: sbhatta
+'''
+
+import os, sys
+import argparse
+
+import bob.io.base
+import bob.io.video
+import bob.ip.color
+import numpy as np
+import galbally_iqm_features as iqm
+import antispoofing.utils.db as bobdb
+
+
+#'''
+#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 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
+#    #scaleMax = 235.0/255.0
+#    gray = scaleMin + (C0*rgbImage[0,:,:] + C1*rgbImage[1,:,:] + C2*rgbImage[2,:,:])    
+#
+#    return gray
+
+
+# """
+# loads a video, and returns a feature-vector for each frame of video
+# """
+# def computeIQM_1video(vidPath):
+#     inputVideo = bob.io.video.reader(vidPath)
+#     vin = inputVideo.load()
+#     numframes = vin.shape[0]
+#     fset = np.zeros([numframes, 21])
+#     for f in range(numframes):
+#         rgbFrame = vin[f,:,:,:]
+#         grayFrame = matlab_rgb2gray(rgbFrame) #compute gray-level image for input color-frame
+#         bobQFeats = np.asarray(iqm.computeQualityFeatures(grayFrame)) # computeQualityFeatures() returns a tuple
+#         fset[f,:] = bobQFeats
+#     
+#     return fset
+#
+
+
+'''
+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]
+    
+    #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 = iqm.compute_quality_features(grayFrame)
+    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(iqm.compute_quality_features(grayFrame)) # computeQualityFeatures() returns a tuple
+        bobQFeats = np.asarray(iqm.compute_quality_features(rgbFrame)) # computeQualityFeatures() returns a tuple
+        fset[f,:] = bobQFeats
+    
+    return fset
+
+
+'''
+loads a video, and returns a feature-vector for each frame of video
+'''
+def computeIQM_1video(vidPath):
+    inputVideo = bob.io.video.reader(vidPath)
+    vin = inputVideo.load()
+    return computeVideoIQM(vin)
+    
+
+
+def main(command_line_parameters=None):
+    
+    #code for parsing command line args.
+    argParser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    argParser.add_argument('-f', '--print_num_files', action='store_true', dest='printNumFiles',
+      default=False, help='Option to print no. of files that will be processed. (Default: %(default)s)')
+    
+    argParser.add_argument('-i', '--input_videofile', dest='inpFile', default = None,
+       help='filename of video to be processed (including complete path). Video expected in .mov format.')
+    
+    argParser.add_argument('-o', '--output_featurefile', dest='outFile', default = None,
+       help='filename where computed features will be stored. Output file will be in hdf5 format.')
+    
+    args = argParser.parse_args(command_line_parameters)
+    #make sure the user specifies a folder where feature-files exist
+    if not args.inpFile: argParser.error('Specify parameter --input_videofile')
+    if not args.outFile: argParser.error('Specify parameter --output_featurefile')
+
+    #1. load video file
+    infile = args.inpFile #k.make_path(videoRoot, '.mov')
+    #2. compute features, 1 vector per frame of input video.
+    bobIqmFeats = computeIQM_1video(infile)
+    #3. save features in file 
+    outfile = args.outFile #k.make_path(featRoot, '.h5')
+    ohf = bob.io.base.HDF5File(outfile, 'w')
+    ohf.set('bobiqm', bobIqmFeats)
+    del ohf
+    
+
+if __name__ == '__main__':
+    main(sys.argv[1:])

bob/ip/qualitymeasure/compute_msuiqa_features.py

+'''
+Created on 13 Oct 2015
+
+@author: sbhatta
+'''
+
+import os, sys
+import argparse
+
+import bob.io.base
+import bob.io.image #under the hood: loads Bob plugin for image file
+
+import bob.io.video
+import bob.ip.color
+import numpy as np
+import msu_iqa_features as iqa
+#import MSU_MaskedIQAFeats as iqa
+import antispoofing.utils.db as bobdb
+
+
+'''
+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 computeVideoIQA(video4d, validFrames):
+    numframes = video4d.shape[0]
+    
+    #process first frame separately, to get the no. of iqm features
+    
+    numValidFrames = np.sum(validFrames)
+    k=0
+    while validFrames[k] == 0: k+=1
+    print 'first valid frame: ', k
+    rgbFrame = video4d[k,:,:,:]
+    
+    iqmSet = iqa.computeMsuIQAFeatures(rgbFrame)
+
+    numIQMs = len(iqmSet)
+    #now initialize fset to store iqm features for all frames of input video.
+    fset = np.zeros([numValidFrames, numIQMs])
+    msuQFeats = np.asarray(iqmSet) # computeQualityFeatures() returns a tuple
+    fset[0,:] = msuQFeats
+    print 'fset shape:', fset.shape
+    j=1
+    for f in range(k+1,numframes):
+        if validFrames[f]==1:
+            rgbFrame = video4d[f,:,:,:]
+            #grayFrame = matlab_rgb2gray(rgbFrame) #compute gray-level image for input color-frame
+            msuQFeats = np.asarray(iqa.computeMsuIQAFeatures(rgbFrame)) # computeQualityFeatures() returns a tuple
+            fset[j,:] = msuQFeats
+            #print j, f
+            j += 1
+            
+  
+    return fset
+
+
+'''
+loads a video, and returns a feature-vector for each frame of video
+'''
+def computeIQA_1video(videoFile, frameQualFile):
+    inputVideo = bob.io.video.reader(videoFile)
+    
+    #load input video
+    vin = inputVideo.load()
+    numFrames = vin.shape[0]
+    
+    
+    if frameQualFile is not None:
+        f = bob.io.base.HDF5File(frameQualFile) #read only
+        validFrames = (f.read('/frameQuality')).flatten() #reads list of frame-quality indicators
+        validFrames[np.where(validFrames <> 1)]=0
+    else:
+        validFrames = np.ones(numFrames)
+    #print validFrames
+#     print type(validFrames)
+    numValidFrames = np.sum(validFrames)
+    
+    print 'valid frames:', numValidFrames, 'of', numFrames
+    
+    #bob.io.base.save(vin[0,:,:,:].astype('uint8'), '/idiap/temp/sbhatta/msudb_colorImg.png')
+    
+    import time
+    startTime = time.time()
+    fset = computeVideoIQA(vin, validFrames)
+    print("Time for one video --- %s seconds ---" % (time.time() - startTime))
+    
+    return fset
+    
+
+'''
+'''
+def parse_arguments(arguments):
+        #code for parsing command line args.
+    argParser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+
+#     # verbose
+    argParser.add_argument('-v', '--verbose', dest='verbose', metavar='INT', type=int, choices=[0, 1, 2], default=1,
+      help='Prints (hopefully helpful) messages  (Default: %(default)s)')
+    
+    argParser.add_argument('-db', '--dbase_path', dest='dbPath', default = None, #'/idiap/user/sbhatta/work/Antispoofing/ImageQualityMeasures',
+       help='path where database videos exist.')
+
+    argParser.add_argument('-op', '--output_path', dest='outPath', default = None,
+       help='path where face-files will be stored.')
+    
+    argParser.add_argument('-nf', '--numFiles', action='store_true', dest='numFiles',
+      default=False, help='Option to print no. of files that will be processed. (Default: %(default)s)')
+    
+    argParser.add_argument('-f', '--singleFile', dest='singleFile', default=None, 
+      help='filename (including complete path) of video to be used to test this code: %(default)s)')
+    
+    argParser.add_argument('-ve', '--video_ext', dest='vidExtn', default=None, choices = ['avi', 'mov', 'mp4'],
+      help='filename (including complete path) of video to be used to test this code: %(default)s)')
+    
+    
+    bobdb.Database.create_parser(argParser, implements_any_of='video')
+    args = argParser.parse_args(arguments)
+    
+    database = args.cls(args)
+   
+    if args.singleFile is None:
+        #make sure the user specifies a folder where feature-files exist
+        if not args.dbPath: argParser.error('Specify parameter --dbase_path')
+    else:
+        folder = os.path.dirname(args.singleFile)
+        filename = os.path.basename(args.singleFile)
+        args.dbPath = folder
+        args.singleFile = filename
+            
+    if not args.outPath: argParser.error('Specify parameter --output_path')
+        
+    return (args, database)
+
+
+'''
+'''
+def main(arguments):
+    
+    args, database = parse_arguments(arguments)
+    
+    inpDir = args.dbPath
+    outDir = args.outPath
+    assert os.path.exists(inpDir), "Input database folder %s does not exist" %inpDir
+    if args.verbose>0: print 'Loading data from',inpDir
+    
+    if args.singleFile is None:
+            
+        tr_realFiles, tr_attackFiles = database.get_train_data()
+        dv_realFiles, dv_attackFiles = database.get_devel_data()
+        ts_realFiles, ts_attackFiles = database.get_test_data()
+        allFiles = tr_realFiles + dv_realFiles + ts_realFiles + tr_attackFiles + dv_attackFiles + ts_attackFiles
+        del tr_realFiles, tr_attackFiles, dv_realFiles, dv_attackFiles, ts_realFiles, ts_attackFiles
+        
+        numFiles = len(allFiles)
+        if args.numFiles:
+            print 'Number of files to be processed:',numFiles
+            print 'exiting'
+            return
+        
+    #     print numFiles
+    #     numFiles = 1        #test
+        
+        # if we are on a grid environment, just find what I have to process.
+        fileSet = allFiles[0:numFiles]
+        if os.environ.has_key('SGE_TASK_ID'):
+            pos = int(os.environ['SGE_TASK_ID']) - 1
+            
+            if pos >= numFiles:
+                raise RuntimeError, "Grid request for job %d on a setup with %d jobs" % (pos, numFiles)
+            fileSet = [allFiles[pos]] # objects = [objects[pos]]
+
+        print 'processing', len(fileSet), ' files'
+        k1=0
+        for k in fileSet:
+            #1. load video file
+            print 'filenum:', k1
+    #         infile = k.make_path(videoRoot, '.avi')
+    #         outfile = k.make_path(featRoot, '.h5')
+            print k
+            if args.vidExtn is None:
+                inVidFile = k.videofile(inpDir)  #k.make_path(inpDir, '.avi')
+            else:
+                inVidFile = k.make_path(inpDir, ('.' + args.vidExtn))
+            inFrameQFile = None #k.make_path(inpDir, '_frameQ.h5')
+            outFeatFile = k.make_path(outDir, '.h5')
+            head, tail = os.path.split(outFeatFile)
+            if not os.path.exists(head): os.makedirs(head)      #create output folder, if it doesn't exist
+            
+            print inFrameQFile
+            print outFeatFile
+            
+    #         if True: #not os.path.isfile(outFeatFile):
+            msuIQAFeats = computeIQA_1video(inVidFile, inFrameQFile)
+    
+            #4. save features in file 
+            ohf = bob.io.base.HDF5File(outFeatFile, 'w')
+            ohf.set('msuiqa', msuIQAFeats)
+            del ohf
+            
+    #         assert 0>1, 'stop'
+            k1 += 1    
+    else:
+        # test feature-computation with a single file specified as input
+        filePart = os.path.splitext(args.singleFile)[0]
+        inVidFile = os.path.join(args.dbPath, filePart)+ '.avi'
+        inFrameQFile = os.path.join(args.dbPath, filePart)+ '_frameQ.h5'
+        
+        outFeatFile = os.path.join(outDir, filePart)+ '.h5'
+        head, tail = os.path.split(outFeatFile)
+        if not os.path.exists(head): os.makedirs(head)      #create output folder, if it doesn't exist
+            
+        print 'single file:', inVidFile
+        print inFrameQFile
+        print outFeatFile
+        
+        msuIQAFeats = computeIQA_1video(inVidFile, inFrameQFile)
+        #4. save features in file 
+        ohf = bob.io.base.HDF5File(outFeatFile, 'w')
+        ohf.set('msuiqa', msuIQAFeats)
+        del ohf
+
+# special fn to extract first frame from video-file and store it as hdf5
+def extractTestFrame():
+    videoFile = '/idiap/home/sbhatta/work/git/refactoring/bob.db.msu_mfsd_mod/bob/db/msu_mfsd_mod/test_images/real/real_client022_android_SD_scene01.mp4'
+    inputVideo = bob.io.video.reader(videoFile)
+    
+    #load input video
+    vin = inputVideo.load()
+    numFrames = vin.shape[0]
+    outframe = vin[0]
+    outfile = '/idiap/home/sbhatta/work/git/refactoring/bob.db.msu_mfsd_mod/bob/db/msu_mfsd_mod/test_images/real_client022_android_SD_scene01_frame0_correct.hdf5'
+    ohf = bob.io.base.HDF5File(outfile, 'w')
+    ohf.set('color_frame', outframe)
+    del ohf
+
+if __name__ == '__main__':
+#    extractTestFrame()
+    main(sys.argv[1:])

bob/ip/qualitymeasure/msu_iqa_features.py

+'''
+Created on 9 Feb 2016
+
+@author: sbhatta
+'''
+
+
+
+#import re
+#import os
+import math
+
+import numpy as np
+import scipy as sp
+import scipy.signal as ssg
+import scipy.ndimage.filters as snf
+import IQMFeatures as iqm
+import bob.ip.base
+import bob.ip.color
+
+########## Utility functions ###########
+'''
+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 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
+    #scaleMax = 235.0/255.0
+    gray = scaleMin + (C0*rgbImage[0,:,:] + C1*rgbImage[1,:,:] + C2*rgbImage[2,:,:])    
+
+    return gray
+
+'''
+'''
+def matlab_rgb2hsv(rgbImage):
+    # first normalize the range of values to 0-1
+   
+    isUint8 = True
+    if isUint8: rgbImage = rgbImage.astype(np.float64)/255.0
+        
+    hsv = np.zeros_like(rgbImage)
+    bob.ip.color.rgb_to_hsv(rgbImage, hsv)
+    h = hsv[0,:,:]
+    s = hsv[1,:,:]
+    v = hsv[2,:,:]
+#     
+    return (h, s, v)
+
+
+def imshow(image):
+    import matplotlib
+    from matplotlib import pyplot as plt
+    if len(image.shape)==3:
+        #imshow() expects color image in a slightly different format, so first rearrange the 3d data for imshow...
+        outImg = image.tolist()
+        print len(outImg)
+        result = np.dstack((outImg[0], outImg[1]))
+        outImg = np.dstack((result, outImg[2]))
+        plt.imshow((outImg*255.0).astype(np.uint8)) #[:,:,1], cmap=mpl.cm.gray)
+         
+    else:
+        if(len(image.shape)==2):
+            #display gray image.
+            plt.imshow(image.astype(np.uint8), cmap=matplotlib.cm.gray)
+             
+    plt.show()
+    
+
+########### End of Utilities ##############
+
+########### Auxilliary functions  ##############
+
+"""
+"""
+def sobelEdgeMap(image, orientation='both'):
+    
+    #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
+    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,:,:]
+    if orientation is 'horizontal':
+        refEdge = iqm.edgeThinning(refSobelX[:,:], refSobelX[:,:], thinning)
+    else:
+        if orientation is 'vertical':
+            refEdge = iqm.edgeThinning(refSobelY[:,:], refSobelY[:,:], thinning)
+        else:
+            refEdge = iqm.edgeThinning(refSobelX[:,:], refSobelY[:,:], thinning)
+            
+
+    return refEdge
+
+########### End of Aux. functions ##############
+
+'''
+'''
+#def computeMsuIQAFeatures(rgbImage, printFV=False):
+def computeMsuIQAFeatures(rgbImage):
+    assert len(rgbImage.shape)==3, 'computeMsuIQAFeatures():: image should be a 3D array (containing a rgb image)'
+#     hsv = np.zeros_like(rgbImage)
+#     bob.ip.color.rgb_to_hsv(rgbImage, hsv)
+#     h = hsv[0,:,:]
+#     s = hsv[1,:,:]
+#     v = hsv[2,:,:]
+    h,s,v = matlab_rgb2hsv(rgbImage) #defined above. Calls Bob's rgb_to_hsv() after rescaling the input image.
+    
+    #print "computeMsuIQAFeatures():: check bob.ip.color.rgb_to_hsv conversion"
+    grayImage = np.zeros_like(h, dtype='uint8')
+    bob.ip.color.rgb_to_gray(rgbImage, grayImage)
+    
+    blurFeat = blurriness(grayImage)
+#     print 'blurriness:', blurFeat
+    
+    pinaBlur = marzilianoBlur(grayImage)
+    pinaBlur /= 30.0
+#     print 'pinaBlur:',pinaBlur
+    
+    colorHist, totNumColors = calColorHist(rgbImage)
+    totNumColors /= 2000.0 #as done in Matlab code provided by MSU.
+#     print "color hist shape:", colorHist.shape
+#     print colorHist[0:11]
+#     print 'totNumColors', totNumColors
+    
+    # calculate mean, deviation and skewness of each channel
+    # use histogram shifting for the hue channel
+    #print h.shape
+    momentFeatsH = calmoment_shift(h)
+    #print 'H-moments:', momentFeatsH
+        
+    momentFeats = momentFeatsH.copy()
+    momentFeatsS = calmoment(s)
+    #print 'S-moments:', momentFeatsS
+    momentFeats = np.hstack((momentFeats, momentFeatsS))
+    momentFeatsV = calmoment(v)
+    #print 'V-moments:', momentFeatsV
+    momentFeats = np.hstack((momentFeats, momentFeatsV))
+    
+    fv = momentFeats.copy()
+    #print 'moment features:', fv
+    
+    fv = np.hstack((fv, colorHist))
+    fv = np.hstack((fv, totNumColors))
+    fv = np.hstack((fv, blurFeat))
+    fv = np.hstack((fv, pinaBlur))
+
+    return fv
+
+
+"""
+Implements the method proposed by Marziliano et al. for determining the average width of vertical edges, as a measure of blurredness in an image.
+This function is a Python version of 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. 
+"""
+def marzilianoBlur(image):
+    assert len(image.shape)==2, 'marzilianoBlur():: input image should be a 2D array (gray level image)'
+       
+    edgeMap = sobelEdgeMap(image, 'vertical')        # compute vertical edge-map of image using sobel 
+    
+    #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().
+#     edgeMap = bob.io.base.load('/idiap/temp/sbhatta/msudb_faceEdgeMap.png')
+#     imshow(edgeMap)
+    
+    blurImg = image
+    C = blurImg.shape[1]    #number of cols in image
+    (row, col) = edgeMap.nonzero()  # row, col contain the indices of the pixels that comprise edge-map.
+    
+    blurMetric = 0
+    if len(row) > 0:
+
+        #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:', 1+col
+        #print 'lexsort_row:', 1+row
+        #This was only used for debugging (to compare with Matlab code). In fact it is not necessary, so it is commented out.
+    
+        
+        edgeWidths = np.zeros_like(row, dtype=int)
+        
+        firstRow = row[0]
+    #     print 'firstRow:',firstRow
+        for i in range(len(row)):
+            rEdge = row[i]
+            cEdge = col[i]
+    #         if rEdge == firstRow: print "edgePoints:", (i, rEdge, cEdge)
+            
+            cStart = 0          # instead of setting them to 'inf' as in MSU's Matlab version
+            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
+                    cStart = j
+                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
+                    cStart = j
+            
+            #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
+                    cEnd = j
+                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
+                    cEnd = j
+            
+            edgeWidths[i] = cEnd - cStart
+            
+            #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!'
+            
+        # Final metric computation
+        blurMetric = np.mean(edgeWidths)
+            
+    #compute histogram of edgeWidths ...(later)
+                                # binnum = 100;
+                                # t = ((1:binnum) - .5) .* C ./ binnum;
+                                # whist = hist(width_array, t) ./ length(width_array);
+
+    return blurMetric
+
+
+"""
+    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
+"""
+def calmoment( channel, regionMask=None ):
+    assert len(channel.shape) == 2, 'calmoment():: channel should be a 2D array (a single color-channel)'    
+
+    t = np.arange(0.05, 1.05, 0.05) + 0.025                 # t = 0.05:0.05:1;
+#     t = np.arange(0.05, 1.05, 0.05) + 0.025                 # t = 0.05:0.05:1;
+#     np.insert(t, 0, -np.inf)
+#     t[-1]= np.inf
+#     print type(t)
+#     print t
+
+    nPix = np.prod(channel.shape)                   # pixnum = length(channel(:));
+    m = np.mean(channel)                            # m = mean(channel(:));
+    d = np.std(channel)                             # d = sqrt(sum((channel(:) - m) .^ 2) / pixnum);
+    s = np.sum(np.power( ((channel - m)/d), 3))/nPix  # s = sum(((channel(:) - m) ./ d) .^ 3) / pixnum;
+    #print 't:', t
+    myHH = np.histogram(channel, t)[0]
+    #print myHH
+    hh = myHH.astype(float)/nPix              # hh = hist(channel(:),t) / pixnum;
+    #print 'numPix:', nPix
+    #print 'histogram:',hh
+    
+    #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])
+    s0 = np.sum(hh[0:2])
+    #print s0
+    H = np.hstack((H,s0))
+    s1 = np.sum(hh[-2:])
+    #print s1
+    H = np.hstack((H, s1) )
+    #print 'calmoment:',H.shape
+    
+    return H
+
+
+'''
+'''
+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;
+                                                    # tag = find(channel>1);
+    channel[np.where(channel>1.0)] -= 1.0           # channel(tag) = channel(tag) - 1;
+
+#     t = np.arange(0.05, 1.05, 0.05)                 # t = 0.05:0.05:1;
+#     nPix = np.prod(channel.shape)                   # pixnum = length(channel(:));
+                                                    #     m = mean(channel(:));
+                                                    #     d = sqrt(sum((channel(:) - m) .^ 2) / pixnum);
+                                                    #     s = sum(((channel(:) - m) ./ d) .^ 3) / pixnum;
+                                                    #     hh = hist(channel(:),t) / pixnum;
+                                                    #     
+                                                    #     H = [m d s sum(hh(1:2)) sum(hh(end-1:end))];
+    
+    H = calmoment(channel)
+    
+    return H
+
+
+
+"""
+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 
+"""
+def calColorHist(image, m=100):
+    #1. compute color histogram of image (normalized, if specified)
+    numBins = 32
+    maxval=255
+    #print "calColorHist():: ", image.shape
+    cHist = rgbhist(image, maxval, numBins, 1)
+
+    #2. determine top 100 colors of image from histogram
+    #cHist.sort() cHist = cHist[::-1]
+    y = sorted(cHist, reverse=True) # [Y, I] = sort(H,'descend');
+    cHist=y[0:m]                    # H = Y(1:m)';
+
+    c = np.cumsum(y)                # C = cumsum(Y);
+#     print 'cumsum shape:', c.shape
+#     thresholdedC = np.where(c>0.999)
+#    # print thresholdedC.shape
+#     print 'thresholdedC:', thresholdedC[0][0] #:200]
+    numClrs = np.where(c>0.99)[0][0]    # clrnum = find(C>.99,1,'first') - 1;
+    
+    cHist = np.array(cHist)
+    return cHist, numClrs
+
+
+'''
+computes 3d color histogram of image
+'''
+def rgbhist(image, maxval, nBins, normType=0):
+    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'
+    H = np.zeros((nBins, nBins, nBins), dtype=np.uint32)  # zeros([nBins nBins nBins]);
+
+#     testImage = image[:,0:3,0:3].copy()
+#     print testImage.shape
+#     print image.shape
+#     print testImage[0,:,:]
+#     print ''
+#     print testImage[1,:,:]
+#     print ''
+#     print testImage[2,:,:]
+#     print ''
+#     print testImage.reshape(3, 9, order='C').T
+#     
+#     assert(0>1), "Stop!"
+    
+    decimator = (maxval+1)/nBins
+    numPix = image.shape[1]*image.shape[2]
+    im = image.reshape(3, numPix).copy() # im = reshape(I,[size(I,1)*size(I,2) 3]);
+    im=im.T
+    
+    for i in range(0, numPix):          # for i=1:size(I,1)*size(I,2)
+        p = (im[i,:]).astype(float)     # p = double(im(i,:));
+        p = np.floor(p/decimator)       # p = floor(p/(maxval/nBins))+1;
+        H[p[0], p[1], p[2]] += 1        # H(p(1),p(2),p(3)) = H(p(1),p(2),p(3)) + 1;
+                                        # end
+    #totalNBins = np.prod(H.shape)
+    #H = H.reshape(totalNBins, 1, order='F') same as H = reshape(H, nBins**3, 1)
+    H = H.ravel()        #H = H(:);
+#     print 'H type:',type(H[0])
+#     print H.shape
+    # Un-Normalized histogram
+    
+    if normType ==1: H =  H.astype(np.float32) / np.sum(H)    # l1 normalization
+#     else:
+#         if normType==2:
+#             H = normc(H);            # l2 normalization
+    
+    return H
+
+
+"""
+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
+"""
+def blurriness(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
+    khor = kver.T
+    
+    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');
+    
+    #DFver = np.abs(np.diff(image.astype('int32'), axis=0)) # abs(img(2:end,:) - img(1:end-1,:));
+    #DFhor = np.abs(np.diff(image.astype('int32'), axis=1)) #abs(img(:,2:end) - img(:,1:end-1));
+    # 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.
+    DFver = np.diff(image.astype('int16'), axis=0)
+    DFver[np.where(DFver<0)]=0
+
+    DFhor = np.diff(image.astype('int16'), axis=1)
+    DFhor[np.where(DFhor<0)]=0
+
+    DBver = np.abs(np.diff(Bver, axis=0)) # abs(Bver(2:end,:) - Bver(1:end-1,:));
+    DBhor = np.abs(np.diff(Bhor, axis=1)) #abs(Bhor(:,2:end) - Bhor(:,1:end-1));
+
+    Vver = DFver.astype(float) - DBver.astype(float)
+    Vhor = DFhor.astype(float) - DBhor.astype(float)
+    Vver[Vver<0]=0 #Vver(find(Vver<0)) = 0;
+    Vhor[Vhor<0]=0 #Vhor(find(Vhor<0)) = 0;
+
+    SFver = np.sum(DFver)
+    SFhor = np.sum(DFhor) #sum(DFhor(:));
+
+    SVver = np.sum(Vver) #sum(Vver(:));
+    SVhor = np.sum(Vhor) #sum(Vhor(:));
+
+    BFver = (SFver - SVver) / SFver;
+    BFhor = (SFhor - SVhor) / SFhor;
+
+    blurF = max(BFver, BFhor) #max([BFver BFhor]);
+    
+    return blurF
+
+
+

bootstrap-buildout.py

+##############################################################################
+#
+# Copyright (c) 2006 Zope Foundation and Contributors.
+# All Rights Reserved.
+#
+# This software is subject to the provisions of the Zope Public License,
+# Version 2.1 (ZPL).  A copy of the ZPL should accompany this distribution.
+# THIS SOFTWARE IS PROVIDED "AS IS" AND ANY AND ALL EXPRESS OR IMPLIED
+# WARRANTIES ARE DISCLAIMED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF TITLE, MERCHANTABILITY, AGAINST INFRINGEMENT, AND FITNESS
+# FOR A PARTICULAR PURPOSE.
+#
+##############################################################################
+"""Bootstrap a buildout-based project
+
+Simply run this script in a directory containing a buildout.cfg.
+The script accepts buildout command-line options, so you can
+use the -c option to specify an alternate configuration file.
+"""
+
+import os
+import shutil
+import sys
+import tempfile
+
+from optparse import OptionParser
+
+__version__ = '2015-07-01'
+# See zc.buildout's changelog if this version is up to date.
+
+tmpeggs = tempfile.mkdtemp(prefix='bootstrap-')
+
+usage = '''\
+[DESIRED PYTHON FOR BUILDOUT] bootstrap.py [options]
+
+Bootstraps a buildout-based project.
+
+Simply run this script in a directory containing a buildout.cfg, using the
+Python that you want bin/buildout to use.
+
+Note that by using --find-links to point to local resources, you can keep
+this script from going over the network.
+'''
+
+parser = OptionParser(usage=usage)
+parser.add_option("--version",
+                  action="store_true", default=False,
+                  help=("Return bootstrap.py version."))
+parser.add_option("-t", "--accept-buildout-test-releases",
+                  dest='accept_buildout_test_releases',
+                  action="store_true", default=False,
+                  help=("Normally, if you do not specify a --version, the "
+                        "bootstrap script and buildout gets the newest "
+                        "*final* versions of zc.buildout and its recipes and "
+                        "extensions for you.  If you use this flag, "
+                        "bootstrap and buildout will get the newest releases "
+                        "even if they are alphas or betas."))
+parser.add_option("-c", "--config-file",
+                  help=("Specify the path to the buildout configuration "
+                        "file to be used."))
+parser.add_option("-f", "--find-links",
+                  help=("Specify a URL to search for buildout releases"))
+parser.add_option("--allow-site-packages",
+                  action="store_true", default=False,
+                  help=("Let bootstrap.py use existing site packages"))
+parser.add_option("--buildout-version",
+                  help="Use a specific zc.buildout version")
+parser.add_option("--setuptools-version",
+                  help="Use a specific setuptools version")
+parser.add_option("--setuptools-to-dir",
+                  help=("Allow for re-use of existing directory of "
+                        "setuptools versions"))
+
+options, args = parser.parse_args()
+if options.version:
+    print("bootstrap.py version %s" % __version__)
+    sys.exit(0)
+
+
+######################################################################
+# load/install setuptools
+
+try:
+    from urllib.request import urlopen
+except ImportError:
+    from urllib2 import urlopen
+
+ez = {}
+if os.path.exists('ez_setup.py'):
+    exec(open('ez_setup.py').read(), ez)
+else:
+    exec(urlopen('https://bootstrap.pypa.io/ez_setup.py').read(), ez)
+
+if not options.allow_site_packages:
+    # ez_setup imports site, which adds site packages
+    # this will remove them from the path to ensure that incompatible versions
+    # of setuptools are not in the path
+    import site
+    # inside a virtualenv, there is no 'getsitepackages'.
+    # We can't remove these reliably
+    if hasattr(site, 'getsitepackages'):
+        for sitepackage_path in site.getsitepackages():
+            # Strip all site-packages directories from sys.path that
+            # are not sys.prefix; this is because on Windows
+            # sys.prefix is a site-package directory.
+            if sitepackage_path != sys.prefix:
+                sys.path[:] = [x for x in sys.path
+                               if sitepackage_path not in x]
+
+setup_args = dict(to_dir=tmpeggs, download_delay=0)
+
+if options.setuptools_version is not None:
+    setup_args['version'] = options.setuptools_version
+if options.setuptools_to_dir is not None:
+    setup_args['to_dir'] = options.setuptools_to_dir
+
+ez['use_setuptools'](**setup_args)
+import setuptools
+import pkg_resources
+
+# This does not (always?) update the default working set.  We will
+# do it.
+for path in sys.path:
+    if path not in pkg_resources.working_set.entries:
+        pkg_resources.working_set.add_entry(path)
+
+######################################################################
+# Install buildout
+
+ws = pkg_resources.working_set
+
+setuptools_path = ws.find(
+    pkg_resources.Requirement.parse('setuptools')).location
+
+# Fix sys.path here as easy_install.pth added before PYTHONPATH
+cmd = [sys.executable, '-c',
+       'import sys; sys.path[0:0] = [%r]; ' % setuptools_path +
+       'from setuptools.command.easy_install import main; main()',
+       '-mZqNxd', tmpeggs]
+
+find_links = os.environ.get(
+    'bootstrap-testing-find-links',
+    options.find_links or
+    ('http://downloads.buildout.org/'
+     if options.accept_buildout_test_releases else None)
+    )
+if find_links:
+    cmd.extend(['-f', find_links])
+
+requirement = 'zc.buildout'
+version = options.buildout_version
+if version is None and not options.accept_buildout_test_releases:
+    # Figure out the most recent final version of zc.buildout.
+    import setuptools.package_index
+    _final_parts = '*final-', '*final'
+
+    def _final_version(parsed_version):
+        try:
+            return not parsed_version.is_prerelease
+        except AttributeError:
+            # Older setuptools
+            for part in parsed_version:
+                if (part[:1] == '*') and (part not in _final_parts):
+                    return False
+            return True
+
+    index = setuptools.package_index.PackageIndex(
+        search_path=[setuptools_path])
+    if find_links:
+        index.add_find_links((find_links,))
+    req = pkg_resources.Requirement.parse(requirement)
+    if index.obtain(req) is not None:
+        best = []
+        bestv = None
+        for dist in index[req.project_name]:
+            distv = dist.parsed_version
+            if _final_version(distv):
+                if bestv is None or distv > bestv:
+                    best = [dist]
+                    bestv = distv
+                elif distv == bestv:
+                    best.append(dist)
+        if best:
+            best.sort()
+            version = best[-1].version
+if version:
+    requirement = '=='.join((requirement, version))
+cmd.append(requirement)
+
+import subprocess
+if subprocess.call(cmd) != 0:
+    raise Exception(
+        "Failed to execute command:\n%s" % repr(cmd)[1:-1])
+
+######################################################################
+# Import and run buildout
+
+ws.add_entry(tmpeggs)
+ws.require(requirement)
+import zc.buildout.buildout
+
+if not [a for a in args if '=' not in a]:
+    args.append('bootstrap')
+
+# if -c was provided, we push it back into args for buildout' main function
+if options.config_file is not None:
+    args[0:0] = ['-c', options.config_file]
+
+zc.buildout.buildout.main(args)
+shutil.rmtree(tmpeggs)

buildout.cfg

+; vim: set fileencoding=utf-8 :
+; Tue 16 Aug 17:30:19 CEST 2016
+
+[buildout]
+parts = scripts
+develop = .
+eggs = bob.ip.qualitymeasure
+extensions = bob.buildout
+newest = false
+verbose = true
+
+[scripts]
+recipe = bob.buildout:scripts
+dependent-scripts = true

develop.cfg

+; vim: set fileencoding=utf-8 :
+; Andre Anjos <andre.anjos@idiap.ch>
+; Mon 16 Apr 08:29:18 2012 CEST
+
+[buildout]
+parts = scripts
+eggs = bob.ip.qualitymeasure
+extensions = bob.buildout
+             mr.developer
+auto-checkout = *
+develop = src/bob.extension
+          src/bob.blitz
+          src/bob.core
+          src/bob.io.base
+          src/bob.io.image
+          src/bob.io.video
+          src/bob.math
+          src/bob.sp
+          src/bob.ip.base
+          src/bob.ip.color
+          .
+
+; options for bob.buildout extension
+debug = true
+verbose = true
+newest = false
+
+[sources]
+bob.extension = git https://gitlab.idiap.ch/bob/bob.extension
+bob.blitz = git https://gitlab.idiap.ch/bob/bob.blitz
+bob.core = git https://gitlab.idiap.ch/bob/bob.core
+bob.io.base = git https://gitlab.idiap.ch/bob/bob.io.base
+bob.io.image = git https://gitlab.idiap.ch/bob/bob.io.image
+bob.io.video = git https://gitlab.idiap.ch/bob/bob.io.video
+bob.math = git https://gitlab.idiap.ch/bob/bob.math
+bob.sp = git https://gitlab.idiap.ch/bob/bob.sp
+bob.ip.base = git https://gitlab.idiap.ch/bob/bob.ip.base
+bob.ip.color = git https://gitlab.idiap.ch/bob/bob.ip.color
+
+
+[scripts]
+recipe = bob.buildout:scripts
+dependent-scripts = true

requirements.txt

+setuptools
+docopt
+bob.extension
+bob.io.base
+bob.io.image
+bob.io.video
+bob.ip.draw
+bob.ip.color

setup.py

+#!/usr/bin/env python
+# vim: set fileencoding=utf-8 :
+# Andre Anjos <andre.anjos@idiap.ch>
+# Tue 16 Feb 2016 15:26:26 CET
+
+from setuptools import setup, find_packages, dist
+dist.Distribution(dict(setup_requires=['bob.extension']))
+
+from bob.extension.utils import load_requirements
+requirements = load_requirements()
+
+version = open("version.txt").read().rstrip()
+
+setup(
+
+    name='bob.ip.qualitymeeasure',
+    version=version,
+    description='Image-quality feature-extractors for PAD applications',
+    url='http://gitlab.idiap.ch/bob/bob.ip.qualitymeasure',
+    license='BSD',
+    author='Andre Anjos',
+    author_email='andre.anjos@idiap.ch',
+    keywords='bob, mage-quality, face',
+    maintainer="Sushil Bhattacharjee",
+    maintainer_email="sbhatta@idiap.ch",
+    long_description=open('README.rst').read(),
+
+    # This line is required for any distutils based packaging.
+    packages=find_packages(),
+    include_package_data=True,
+
+    install_requires = requirements,
+
+    entry_points={
+      # scripts should be declared using this entry:
+      'console_scripts': [
+        'detect_landmarks.py = bob.ip.facelandmarks.script.detect_landmarks:main',
+      ],
+    },
+
+    # Classifiers are important if you plan to distribute this package through
+    # PyPI. You can find the complete list of classifiers that are valid and
+    # useful here (http://pypi.python.org/pypi?%3Aaction=list_classifiers).
+    classifiers = [
+      'Framework :: Bob',
+      'Development Status :: 4 - Beta',
+      'Intended Audience :: Developers',
+      'License :: OSI Approved :: BSD License',
+      'Natural Language :: English',
+      'Programming Language :: Python',
+      'Programming Language :: Python :: 3',
+      'Topic :: Software Development :: Libraries :: Python Modules',
+    ],
+)

version.txt

+1.0.0a0