diff --git a/README.rst b/README.rst
new file mode 100644
index 0000000000000000000000000000000000000000..3897ec13d4a968c57ed2de819680259f2ef391cd
--- /dev/null
+++ b/README.rst
@@ -0,0 +1,47 @@
+.. 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
diff --git a/bob/__init__.py b/bob/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..2ab1e28b150f0549def9963e9e87de3fdd6b2579
--- /dev/null
+++ b/bob/__init__.py
@@ -0,0 +1,3 @@
+# see https://docs.python.org/3/library/pkgutil.html
+from pkgutil import extend_path
+__path__ = extend_path(__path__, __name__)
diff --git a/bob/ip/__init__.py b/bob/ip/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..2ab1e28b150f0549def9963e9e87de3fdd6b2579
--- /dev/null
+++ b/bob/ip/__init__.py
@@ -0,0 +1,3 @@
+# see https://docs.python.org/3/library/pkgutil.html
+from pkgutil import extend_path
+__path__ = extend_path(__path__, __name__)
diff --git a/bob/ip/qualitymeasure/__init__.py b/bob/ip/qualitymeasure/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..d5d61f05021867a5aac356dc8a567c2a4a8d4347
--- /dev/null
+++ b/bob/ip/qualitymeasure/__init__.py
@@ -0,0 +1,18 @@
+#!/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('_')]
+
diff --git a/bob/ip/qualitymeasure/compute_iqm_features.py b/bob/ip/qualitymeasure/compute_iqm_features.py
new file mode 100644
index 0000000000000000000000000000000000000000..92904d71a42af1628159a1c827f0d9ee81b1cdef
--- /dev/null
+++ b/bob/ip/qualitymeasure/compute_iqm_features.py
@@ -0,0 +1,128 @@
+#!/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:])
diff --git a/bob/ip/qualitymeasure/compute_msuiqa_features.py b/bob/ip/qualitymeasure/compute_msuiqa_features.py
new file mode 100644
index 0000000000000000000000000000000000000000..6addad0837d2853562d861da37fd6159a81f4e28
--- /dev/null
+++ b/bob/ip/qualitymeasure/compute_msuiqa_features.py
@@ -0,0 +1,240 @@
+'''
+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:])
diff --git a/bob/ip/qualitymeasure/galbally_iqm_features.py b/bob/ip/qualitymeasure/galbally_iqm_features.py
new file mode 100644
index 0000000000000000000000000000000000000000..a1aaa7bcc5d9a147e29ab914fffd2a481e91d9de
--- /dev/null
+++ b/bob/ip/qualitymeasure/galbally_iqm_features.py
@@ -0,0 +1,807 @@
+#!/usr/bin/env python
+# vim: set fileencoding=utf-8 :
+
+'''
+Created on 25 Sep 2015
+
+@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 bob.ip.base
+
+
+
+"""
+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 compute_quality_features(image):
+    """Extract a set of image quality-features computed for the input image.
+    :param image: 2d or 3d numpy array. Should represent input image of shape [M,N] (M rows x N cols).
+    If 2D, image should contain a gray-image of shape [M,N].
+    If 3d, image should have a shape [3,M,N], and should contain an RGB-image.
+
+    :return featSet: a tuple of float-scalars, each representing one image-quality measure.
+    This function returns a subset of the image-quality features (for face anti-spoofing) that have been 
+    described by Galbally et al. in their paper:
+    "Image Quality Assessment for Fake Biometric Detection: Application to Iris, Fingerprint, and Face Recognition",
+     IEEE Trans. on Image Processing Vol 23(2), 2014.
+    """
+    
+    if len(image.shape) == 3:
+        if(image.shape[0]==3): 
+            gray_image = matlab_rgb2gray(rgbFrame) #compute gray-level image for input color-frame
+        else:
+            print('error. Wrong kind of input image')
+    else:
+        if len(image.shape) == 2:
+            gray_image = image
+        else:
+            print('error -- wrong kind of input')
+
+    gwin = gauss_2d((3,3), 0.5)
+    if gray_image: 
+        smoothed = ssg.convolve2d(gray_image, gwin, boundary='symm', mode='same') 
+    
+        """
+        Some of the image-quality measures computed here require a reference image.
+        For these measures, we use the input-image itself as a reference-image, and we compute
+        the quality-measure of a smoothed version of the input-image. The assumption in this
+        approach is that smoothing degrades a spoof-image more than it does a genuine image
+        (see Galbally's paper referenced above).
+        """
+        featSet = image_quality_measures(gray_image, smoothed)
+    
+        return featSet
+
+    else:
+        return None
+
+
+
+"""
+actually computes various measures of similarity 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 image_quality_measures(refImage, testImage):
+    """Compute image-quality measures for testImage and return a tuple of quality-measures.
+       Some of the quality-measures require a reference-image, but others are 'no-reference' measures.
+       :input refImage: 2d numpy array. Should represent input 8-bit gray-level image of size [M,N].
+       :input testImage: 2d numpy array. Should represent input 8-bit gray-level image of size [M,N]..
+       :return 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.
+    """
+    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
+    mse00 = float(sumDiffSq)/float(numPx)
+    
+    #2 PSNR
+    psnr01 = np.inf
+    if mse00>0:
+        psnr01 = 10.0*np.log10(maxPxVal*maxPxVal/mse00)
+    
+    #3 AD: Average difference
+    ad02 = float(np.sum(diffImg))/float(numPx)
+    
+    #4 SC: structural content
+    testSq = np.square(testImage.astype(np.float))
+    sumTestSq = np.sum(testSq)
+    sc03=np.inf
+    if sumTestSq>0:
+        sc03 = float(sumRefSq)/float(sumTestSq)
+    
+    #5 NK: normalized cross-correlation
+    imgProd = refImage * testImage # element-wise product
+    nk04 = float(np.sum(imgProd))/float(sumRefSq)
+    
+    #6 MD: Maximum difference
+    md05 = 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)
+    lmse06 = float(np.sum(num_op))/float(sum_opSq)
+    
+    #8 NAE: normalized abs. error
+    sumRef = np.sum(np.absolute(refImage))
+    nae07 = float(np.sum(absDiffImg))/float(sumRef)
+    
+    #9 SNRv: SNR in db
+    snrv08 = 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])
+    ramdv09 = np.sqrt(float(topsum)/float(r))
+    
+    #11,12: MAS: Mean Angle Similarity, MAMS: Mean Angle-Magnitude Similarity
+    mas10, mams11 = angle_similarity(refImage, testImage, diffImg) 
+    
+    fftRef = np.fft.fft2(refImage)
+#     fftTest = np.fft.fft2(testImage)
+      
+    #13, 14: SME: spectral magnitude error; SPE: spectral phase error
+    sme12, spe13 = spectral_similarity(refImage, testImage) #spectralSimilarity(fftRef, fftTest, numPx)
+    
+    #15 TED: Total edge difference
+    ted14 = edge_similarity(refImage, testImage)
+
+    #16 TCD: Total corner difference
+    tcd15 = corner_similarity(refImage, testImage)
+    
+    #17, 18: GME: gradient-magnitude error; GPE: gradient phase error  
+    gme16, gpe17 = gradient_similarity(refImage, testImage)
+    
+    #19 SSIM
+    ssim18, _ = ssim(refImage, testImage)
+    
+    #20 VIF
+    vif19 = vif(refImage, testImage)
+    
+    #21,22,23,24,25: RRED, BIQI, JQI, NIQE: these parameters are not computed here.
+    
+    #26 HLFI: high-low frequency index (implemented as done by Galbally in Matlab).
+    hlfi25=high_low_freq_index(fftRef, refImage.shape[1])
+    
+    return (mse00, psnr01, ad02, sc03, nk04, md05, lmse06, nae07, snrv08, ramdv09, mas10, mams11, sme12, gme16, gpe17, ssim18, vif19, hlfi25)
+
+
+"""
+Matlab-like RGB to gray...
+"""
+def matlab_rgb2gray(rgbImage): 
+    '''converts color rgbImage to gray to produce exactly the same result as Matlab would.
+    Inputs:
+    rgbimage: numpy array of shape [3, height, width]
+    Return:
+    numpy array of shape [height, width] containing a gray-image with floating-point pixel values, in the range[(16.0/255) .. (235.0/255)]
+    '''
+    #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
+
+
+"""
+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):
+    """Compute and return SSIM between two images.
+    @param refImage: 2D numpy array (reference image)
+    @param testImage: 2D numpy array (test image)
+    Returns ssim and ssim_map
+     @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).
+    """
+    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)
+    #refImage=refImage.astype(np.float)
+    #testImage=testImage.astype(np.float)
+    
+    #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 vif(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
+        #print N, sc
+        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 high_low_freq_index(imgFFT, ncols):
+
+    N= ncols
+    colHalf = int(round(N/2)) # (N/2) + (N % 2) #round it up
+    freqSel = 0.15
+    
+    freqCol = round(freqSel*N)
+    lowFreqColHalf = int(round(freqCol/2.0)) 
+    
+    fftRes = imgFFT #np.fft.fft2(image)
+    fftMag = np.abs(fftRes)
+    totalEnergy = np.sum(fftMag)
+    #print totalEnergy
+    
+    lowIdx = colHalf-lowFreqColHalf
+    hiIdx = colHalf + lowFreqColHalf
+    #print lowIdx, hiIdx
+    LowFreqMag = fftMag[:, lowIdx:hiIdx]
+    lowFreqMagTotal = np.sum(LowFreqMag)
+    
+    fftMag[:, lowIdx:hiIdx]=0
+    highFreqMagTotal = np.sum(fftMag)
+    #print 'partial freq. sums:', lowFreqMagTotal, highFreqMagTotal
+    
+    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 gradient_similarity(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
+
+'''
+'''
+def testRegionalMax():
+    A = 10*np.ones([10,10])
+    A[1:4, 1:4] = 22
+    A[5:8, 5:8] = 33
+    A[1,6]=44
+    A[2,7]=45
+    A[3,8]=44
+    rm = regionalmax(A)
+    print A
+    print rm
+
+"""
+find local maxima using 3x3 mask.
+Used in corner_similarity()
+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
+    
+
+# 
+# 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()
+
+'''
+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 corner_similarity(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 edge_similarity(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 = edge_thinning(refSobelX[:,:], refSobelY[:,:], thinning)
+    
+    
+    #compute edge map for test image
+    testSobel_sep = bob.ip.base.sobel(testImage)
+    testSobelX = testSobel_sep[0,:,:]
+    testSobelY = testSobel_sep[1,:,:]
+#     testSobelX = ssg.convolve2d(testImage, x_mask, boundary='symm', mode='valid')
+#     testSobelY = ssg.convolve2d(testImage, y_mask, boundary='symm', mode='valid')
+    testEdge = edge_thinning(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 edge_similarity()
+    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 edge_thinning(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 spectral_similarity(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 angle_similarity(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
+    
+    #the following (commented out) code should be more efficient than the for-loop above, but seems to be buggy.
+#     rowDP= np.diag(np.dot(refImage,testImage.T))
+#     normRefrows= np.linalg.norm(refImage,axis=1)
+#     normTestrows= np.linalg.norm(testImage,axis=1)
+#     
+#     cosThetaVec = rowDP/(normRefrows * normTestrows)
+#     cosThetaVec = np.nan_to_num(cosThetaVec) #nan occurs when one of the norms is 0, ie. vector is all 0s. In that case set cosTheta to 0
+#     
+#     thetaVec = np.arccos(cosThetaVec)
+#     
+    
+    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'...)
+
+Works correctly.
+@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 
+
+
+'''
+returns a smoothed version of input gray-level image
+smoothing is done using a fixed 3x3 gaussian filter (sigma=0.5)
+'''
+def gaussianSmooth(image):
+    # perform Gaussian filtering
+    gaussian = bob.ip.base.Gaussian(sigma = (0.5, 0.5), radius = (1, 1)) #radius=1 creates a 3x3 filter.
+    return gaussian(image)
+
+
+
+# '''
+# '''
+# def testQualityMeasures(image, smoothed):
+#     frameFeatSet = imageQualityMeasures(image, smoothed)
+#     
+#     print frameFeatSet
+#     
+#     
+# if __name__ == '__main__':
+#     image, smoothed = testGaussianFilter()
+#     
+#     #testQualityMeasures(image, smoothed)
+#     diffImg = image - smoothed
+#     #test MAMS
+#     mas, mams = angle_similarity2(image, smoothed, diffImg)
+#     print 'mas', mas
+#     print 'mams', mams
+
+
+    
+    
+    
diff --git a/bob/ip/qualitymeasure/msu_iqa_features.py b/bob/ip/qualitymeasure/msu_iqa_features.py
new file mode 100644
index 0000000000000000000000000000000000000000..fe2ceff60f7c63e921fc8aef5e0110b574d3cefa
--- /dev/null
+++ b/bob/ip/qualitymeasure/msu_iqa_features.py
@@ -0,0 +1,437 @@
+'''
+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
+
+
+
diff --git a/bootstrap-buildout.py b/bootstrap-buildout.py
new file mode 100644
index 0000000000000000000000000000000000000000..a4599211f741c468cd37a29861d1c7f2c3a641d1
--- /dev/null
+++ b/bootstrap-buildout.py
@@ -0,0 +1,210 @@
+##############################################################################
+#
+# 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)
diff --git a/buildout.cfg b/buildout.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..2a4e052ffb9a5cc941f834a7582c9ac887e1e5e6
--- /dev/null
+++ b/buildout.cfg
@@ -0,0 +1,14 @@
+; 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
diff --git a/develop.cfg b/develop.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..7ceb8035b7d9664ddb837f6c13e384e6644721a2
--- /dev/null
+++ b/develop.cfg
@@ -0,0 +1,43 @@
+; 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
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b56f837b44aad515e9cbd80d02e5f9c6a86c3602
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,8 @@
+setuptools
+docopt
+bob.extension
+bob.io.base
+bob.io.image
+bob.io.video
+bob.ip.draw
+bob.ip.color
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..fafa2f8f5adfc385843c5c40846287a19c52b400
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,54 @@
+#!/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',
+    ],
+)
diff --git a/version.txt b/version.txt
new file mode 100644
index 0000000000000000000000000000000000000000..41432f00d9ce57fadd55cc7dd27b391ddf5ca0b9
--- /dev/null
+++ b/version.txt
@@ -0,0 +1 @@
+1.0.0a0