Merge branch 'optimizations' into 'master'

Optimizations

See merge request !6

MANIFEST.in

 include README.rst bootstrap-buildout.py buildout.cfg develop.cfg requirements.txt version.txt
 recursive-include doc *.py *.rst *.png *.ico
-recursive-include bob/ip/qualitymeasure/data *.mp4 *.jpg *.h5 *.hdf5
+recursive-include bob/ip/qualitymeasure/data *.mp4 *.jpg *.h5 *.hdf5 *.ppm

bob/ip/qualitymeasure/__init__.py

 #!/usr/bin/env python
 # vim: set fileencoding=utf-8 :
 
+from ._library import remove_highlights
 from .galbally_iqm_features import compute_quality_features
 from .msu_iqa_features import compute_msu_iqa_features
 
 
+
 def get_config():
     """
     Returns a string containing the configuration information.

bob/ip/qualitymeasure/main.cpp

+// include directly and indirectly dependent libraries
+#ifdef NO_IMPORT_ARRAY
+#undef NO_IMPORT_ARRAY
+#endif
+
+
+#include <bob.blitz/cppapi.h>
+#include <bob.blitz/cleanup.h>
+#include <bob.extension/documentation.h>
+
+
+// declare C++ function
+void remove_highlights(   blitz::Array<float ,3> &img,
+                          blitz::Array<float ,3> &diff,
+                          blitz::Array<float ,3> &sfi,
+                          blitz::Array<float ,3> &residue,
+                          float epsilon,
+                          bool  skip_diffuse,
+                          bool  check_nan_inf);
+
+// use the documentation classes to document the function
+static bob::extension::FunctionDoc remove_highlights_doc = bob::extension::FunctionDoc(
+  "remove_highlights",
+  "This function implements a specular highlight removal algorithm.",
+  "This function implements a specular highlight removal algorithm which, by \
+     using an iterative process, separates the specular component from the \
+     diffuse component of the picture. It returns a color incorect specular free \
+     image, the diffuse component and the specular residue, respectively. It is \
+     based on the original code by Robby T. Tan \
+     reference: \
+     separating reflection components of textured surfaces using a single image \
+     by Robby T. Tan, Katsushi Ikeuchi, \
+     IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), \
+     27(2), pp.179-193, February, 2005"
+)
+.add_prototype("image", "specular_free_image, diffuse_image, specular_residue, epsilon")
+.add_parameter("image", "array_like (3D, float32)", "The image to process")
+.add_return("specular_free_image", "array_like (3D, float32)", "Specular free image with altered colors.")
+.add_return("diffuse_image", "array_like (3D, float32)", "Diffuse component image, free of specularity.")
+.add_return("specular_residue", "array_like (3D, float32)", "Specular residue of the image.")
+.add_return("epsilon", "float32", "Parameter that specifies the number of iterations.")
+;
+
+
+// declare the function
+static PyObject* PyRemoveHighlights(PyObject*, PyObject* args, PyObject* kwargs) {
+
+  BOB_TRY
+
+  static const char* const_kwlist[] =
+    {"array", "startEps", "skip_diffuse", "check_nan_inf", 0};
+  static char** kwlist = const_cast<char**>(const_kwlist);
+
+  PyBlitzArrayObject* array;
+  double  epsilon  = 0.5f;
+  bool    skip_diffuse  = false;
+  bool    check_nan_inf = false;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&|dii", kwlist,
+                                  &PyBlitzArray_Converter, &array,
+                                  &epsilon, &skip_diffuse, &check_nan_inf))
+    return 0;
+
+  // check that the array has the expected properties
+  if (array->type_num != NPY_FLOAT32|| array->ndim != 3){
+    PyErr_Format(PyExc_TypeError,
+                "remove_highlights : Only 3D arrays of type float32 are allowed");
+    return 0;
+  }
+
+  // extract the actual blitz array from the Python type
+  blitz::Array<float ,3> img = *PyBlitzArrayCxx_AsBlitz<float , 3>(array);
+
+  // results
+  int dim_x = img.shape()[2];
+  int dim_y = img.shape()[1];
+
+  blitz::Array<float ,3> diffuse_img(3, dim_y, dim_x);
+  blitz::Array<float ,3> speckle_free_img(3, dim_y, dim_x);
+  blitz::Array<float ,3> speckle_img(3, dim_y, dim_x);
+
+  diffuse_img       = 0;
+  speckle_free_img  = 0;
+  speckle_img       = 0;
+
+  // call the C++ function
+  remove_highlights(img, diffuse_img, speckle_free_img, speckle_img,
+                    (float)epsilon, skip_diffuse, check_nan_inf);
+
+  // convert the blitz array back to numpy and return it
+  PyObject *ret_tuple = PyTuple_New(3);
+  PyTuple_SetItem(ret_tuple, 0, PyBlitzArrayCxx_AsNumpy(speckle_free_img));
+  PyTuple_SetItem(ret_tuple, 1, PyBlitzArrayCxx_AsNumpy(diffuse_img));
+  PyTuple_SetItem(ret_tuple, 2, PyBlitzArrayCxx_AsNumpy(speckle_img));
+
+  return ret_tuple;
+
+  // handle exceptions that occurred in this function
+  BOB_CATCH_FUNCTION("remove_highlights", 0)
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+/////// Python module declaration ////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+// module-wide methods
+static PyMethodDef module_methods[] = {
+  {
+    "remove_highlights",
+    (PyCFunction)PyRemoveHighlights,
+    METH_VARARGS|METH_KEYWORDS,
+    remove_highlights_doc.doc()
+  },
+  {NULL}  // Sentinel
+};
+
+// module documentation
+PyDoc_STRVAR(module_docstr, "Tan Specular Highlights");
+
+// module definition
+#if PY_VERSION_HEX >= 0x03000000
+static PyModuleDef module_definition = {
+  PyModuleDef_HEAD_INIT,
+  BOB_EXT_MODULE_NAME,
+  module_docstr,
+  -1,
+  module_methods,
+  0, 0, 0, 0
+};
+#endif
+
+// create the module
+static PyObject* create_module (void) {
+
+# if PY_VERSION_HEX >= 0x03000000
+  PyObject* module = PyModule_Create(&module_definition);
+  auto module_ = make_xsafe(module);
+  const char* ret = "O";
+# else
+  PyObject* module = Py_InitModule3(BOB_EXT_MODULE_NAME, module_methods, module_docstr);
+  const char* ret = "N";
+# endif
+  if (!module) return 0;
+
+  if (PyModule_AddStringConstant(module, "__version__", BOB_EXT_MODULE_VERSION) < 0) return 0;
+
+  /* imports bob.blitz C-API + dependencies */
+  if (import_bob_blitz() < 0) return 0;
+
+  return Py_BuildValue(ret, module);
+}
+
+PyMODINIT_FUNC BOB_EXT_ENTRY_NAME (void) {
+# if PY_VERSION_HEX >= 0x03000000
+  return
+# endif
+    create_module();
+}

bob/ip/qualitymeasure/msu_iqa_features.py

@@ -9,10 +9,10 @@ import scipy.signal as ssg
 import bob.ip.base
 import bob.ip.color
 from . import galbally_iqm_features as iqm
-from . import tan_specular_highlights as tsh
 
-''' Utility functions '''
+from . import remove_highlights
 
+''' Utility functions '''
 
 def matlab_rgb2hsv(rgbImage):
     # first normalize the range of values to 0-1
@@ -26,7 +26,7 @@ def matlab_rgb2hsv(rgbImage):
     h = hsv[0, :, :]
     s = hsv[1, :, :]
     v = hsv[2, :, :]
-#
+
     return (h, s, v)
 
 
@@ -155,8 +155,18 @@ def compute_iqa_specularity_features(rgbImage, startEps=0.05):
     """
 
     # separate the specular and diffuse components of input color image.
-    speckleFreeImg, diffuseImg, speckleImg = tsh.remove_highlights(
-        rgbImage.astype(float), startEps, verboseFlag=False)
+    # original python version
+    #speckleFreeImg, diffuseImg, speckleImg = tsh.remove_highlights(
+    #    rgbImage.astype(float), startEps, verboseFlag=False)
+
+    speckleFreeImg, diffuseImg, speckleImg = remove_highlights(
+            rgbImage.astype(np.float32), startEps)
+
+    speckleImg[np.where(np.isinf(speckleImg))] = 0
+    speckleImg[np.where(np.isnan(speckleImg))] = 0
+    speckleImg[np.where(speckleImg < 0)] = 0
+    speckleImg[np.where(speckleImg > 255)] = 255
+
     # speckleImg contains the specular-component
 
     if len(speckleImg.shape) == 3:
@@ -302,7 +312,6 @@ def marzilianoBlur(image):
 
     return blurMetric
 
-
 def calmoment(channel, regionMask=None):
     """ returns the first 3 statistical moments (mean, standard-dev., skewness)
     and 2 other first-order statistical measures of input image

bob/ip/qualitymeasure/script/remove_highlights.py

+#!/usr/bin/env python
+# vim: set fileencoding=utf-8 :
+
+
+import sys
+import argparse
+import time
+
+import bob.io.base
+import bob.io.image
+import numpy as np
+
+from bob.ip.qualitymeasure import remove_highlights
+
+def main(command_line_parameters=None):
+    """Remove the specular component of the input image and write result to
+    a file.
+    """
+
+    argParser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter)
+    argParser.add_argument(
+        '-i',
+        '--input',
+        dest='inpImg',
+        default=None,
+        help='filename of image to be processed (including complete path). ')
+
+    argParser.add_argument(
+        '-o',
+        '--output',
+        dest='outImg',
+        default=None,
+        help='filename of specular free image.')
+
+    argParser.add_argument(
+        '-r',
+        '--residue',
+        dest='resImg',
+        default=None,
+        help='filename of specular residue image.')
+
+    argParser.add_argument(
+        '-e',
+        '--epsilon',
+        dest='epsilon',
+        default=0.5,
+        help='value of epsilon parameter.')
+
+    args = argParser.parse_args(command_line_parameters)
+
+    if not args.inpImg:
+        argParser.error('Specify parameter --input')
+    if not args.outImg:
+        argParser.error('Specify parameter --output')
+
+    # 1. open input image
+    print("Opening file: %s" % args.inpImg)
+    img = bob.io.image.load(args.inpImg)
+
+    # 2. compute
+    print("Extracting diffuse component...")
+
+    sfi, diff, residue = remove_highlights(img.astype(np.float32), float(args.epsilon))
+
+    # 1. save output image
+    print("Saving output file: %s" % args.outImg)
+    bob.io.base.save(diff.astype('uint8'), args.outImg)
+
+    if args.resImg:
+        print("Saving residue file: %s" % args.resImg)
+        bob.io.base.save(residue.astype('uint8'), args.resImg)
+
+    print('Done')
+
+
+if __name__ == '__main__':
+    main(sys.argv[1:])

bob/ip/qualitymeasure/tan_specular_highlights.cpp

+/**
+ *  Original version of the specular highlights removal code by Robby T. Tan
+ *  reference:
+ *  "separating reflection components of textured surfaces using a single image"
+ *  by Robby T. Tan, Katsushi Ikeuchi,
+ *  IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
+ *  27(2), pp.179-193, February, 2005
+ *
+ *  This is the original implementation based on the C++ code provided by Prof.
+ *  Robby Tan but using Blitz++ arrays:
+ *       http://tanrobby.github.io/code.html#
+ *       http://tanrobby.github.io/code/highlight.zip
+ *
+ *  It also implements a modification used in the previous python version that
+ *  ignores pixels marked G_DIFFUSE, leading to a smaller number of
+ *  iterations per epsilon value, and a flag that ensure no nan and inf values
+ *  are produced.
+ */
+
+#include <blitz/array.h>
+
+#define SPECULARX		10
+#define SPECULARY		11
+#define DIFFUSE			12
+#define BOUNDARY		13
+#define NOISE       14
+#define CAMERA_DARK 15
+
+void specular_free_image( blitz::Array<float ,3> &src,
+                          blitz::Array<int,2>   &src_i,
+					                blitz::Array<float ,3> &sfi,
+                          bool check_nan_inf);
+
+void iteration(           blitz::Array<float ,3> &src,
+                          blitz::Array<int,2>   &src_i,
+                          blitz::Array<float ,3> &sfi,
+				                  float  epsilon,
+                          bool skip_diffuse,
+                          bool check_nan_inf);
+
+int init(                 blitz::Array<float ,3> &src,
+                          blitz::Array<int,2>   &src_i,
+                          blitz::Array<float ,3> &sfi,
+				                  float  epsilon,
+                          bool skip_diffuse,
+                          bool check_nan_inf);
+
+int reset_labels(         blitz::Array<int,2>   &src_i);
+
+// the main function to remove highlights from a single image
+void remove_highlights(   blitz::Array<float ,3> &img,
+                          blitz::Array<float ,3> &diff,
+                          blitz::Array<float ,3> &sfi,
+                          blitz::Array<float ,3> &residue,
+                          float  epsilon,
+                          bool skip_diffuse,
+                          bool check_nan_inf)
+{
+  // flags
+  int dim_x = img.shape()[2];
+  int dim_y = img.shape()[1];
+
+  blitz::Array<int,2> img_i(dim_y, dim_x);
+
+  //SPECULAR-FREE IMAGE
+
+  specular_free_image(img, img_i, sfi, check_nan_inf);
+
+  //ITERATIVE PART
+  float  step =0.01f;
+
+  // copy source
+  diff = img;
+
+  while( epsilon >= 0.0 )
+  {
+    // run the main iteration
+    //printf("*");
+    iteration(diff, img_i, sfi, epsilon, check_nan_inf, skip_diffuse);
+    epsilon -= step;
+    //printf(": %f\n",epsilon);
+  }
+
+  // compute residue
+  residue = img - diff;
+}
+
+// utilities
+
+inline float  tot(float  r, float  g, float  b)
+{
+  return r + g + b;
+}
+
+inline float  max(float  r, float  g, float  b)
+{
+  float  max_ = r;
+  if(g>max_) max_ = g;
+  if(b>max_) max_ = b;
+  return max_;
+}
+
+inline float  r_chroma(float  r, float  g, float  b)
+{
+  float  tot_ = tot(r,g,b);
+  if (tot_!=0)  return r/tot_;
+  else          return 0;
+}
+
+inline float  g_chroma(float  r, float  g, float  b)
+{
+  float  tot_ = tot(r,g,b);
+  if (tot_!=0)  return g/tot_;
+  else          return 0;
+}
+
+inline float  b_chroma(float  r, float  g, float  b)
+{
+  float  tot_ = tot(r,g,b);
+  if (tot_!=0)  return b/tot_;
+  else          return 0;
+}
+
+inline float  max_chroma(float  r, float  g, float  b)
+{
+  float  tot_ = tot(r,g,b);
+  if (tot_!=0)  return max(r,g,b)/tot_;
+  else          return 0;
+}
+
+// remove the specular component from source
+
+void specular_free_image( blitz::Array<float ,3> &src,
+                          blitz::Array<int,2>   &src_i,
+					                blitz::Array<float ,3> &sfi,
+                          bool check_nan_inf)
+{
+  float  Lambda=0.6f;
+  float  camDark=10.0f;	// for pixels that are too dark
+  float  lambdaConst = 3.0f * Lambda - 1.0f;
+
+  //SPECULAR FREE IMAGE
+  int dim_x = src.shape()[2];
+  int dim_y = src.shape()[1];
+
+  int y,x;
+  for(y = 0; y < dim_y; y++){
+    for(x = 0; x < dim_x; x++){
+      //get the rgb values
+      float  r = src(0,y,x);
+      float  g = src(1,y,x);
+      float  b = src(2,y,x);
+
+      //copy the rgb to the sfi
+      sfi(0,y,x) = r;
+      sfi(1,y,x) = g;
+      sfi(2,y,x) = b;
+
+      // init
+      src_i(y,x) = 0;
+
+      //check for camera dark and achromatic pixels
+      if(((r < camDark) &&
+          (g < camDark) &&
+          (b < camDark)))
+      {
+	       src_i(y,x) = CAMERA_DARK;
+	       continue;
+      }
+
+      //perform the specular-to-diffuse mechanism
+      float  c    = max_chroma(r,g,b);
+      float  numr = max(r,g,b) * (3.0f * c - 1.0f);
+      float  denm = c * lambdaConst;
+
+      float  dI   = numr / denm;
+      if(denm == 0 && check_nan_inf) dI = 0;
+
+      float  sI = (tot(r,g,b) - dI)/3.0f;
+
+      float  dr,dg,db;
+      dr = (r - sI);
+      dg = (g - sI);
+      db = (b - sI);
+
+      if(dr<0) dr=0;
+      if(dg<0) dg=0;
+      if(db<0) db=0;
+
+      if(dr>255) dr=255;
+      if(dg>255) dg=255;
+      if(db>255) db=255;
+
+      sfi(0,y,x) = dr;
+      sfi(1,y,x) = dg;
+      sfi(2,y,x) = db;
+
+    }
+  }
+
+}
+
+// to apply specular to diffuse equation or mechanism
+
+inline int specular_2_diffuse(int y, int x, blitz::Array<float ,3> &iro,
+                              blitz::Array<int,2> &iro_i,
+                              float  maxChroma,
+                              bool check_nan_inf)
+{
+  float  c = max_chroma(iro(0,y,x), iro(1,y,x), iro(2,y,x));
+  float  m = max(iro(0,y,x), iro(1,y,x), iro(2,y,x));
+  float  t = tot(iro(0,y,x), iro(1,y,x), iro(2,y,x));
+  float  numr = (m*(3.0f*c - 1.0f));
+  float  denm = (c*(3.0f*maxChroma - 1.0f));
+
+  if(check_nan_inf && abs(denm) < 0.000000001)
+  {
+    iro_i(y,x)=NOISE;
+    return 1;
+  }
+
+  float  dI = numr / denm;
+
+  float  sI = (t - dI)/3.0f;
+
+  float  nr = (iro(0,y,x) - sI);
+  float  ng = (iro(1,y,x) - sI);
+  float  nb = (iro(2,y,x) - sI);
+
+  if(nr<=0 || ng<=0 || nb<=0)
+  {
+    iro_i(y,x)=NOISE;
+    return 1;
+  }
+  else
+  {
+    iro(0,y,x) = nr;
+    iro(1,y,x) = ng;
+    iro(2,y,x) = nb;
+
+    return 0;
+  }
+}
+
+// specular reduction mechanism
+
+void iteration( blitz::Array<float ,3> &src,
+                blitz::Array<int,2>   &src_i,
+                blitz::Array<float ,3> &sfi,
+				        float  epsilon,
+                bool skip_diffuse,
+                bool check_nan_inf)
+{
+  int x,y;
+  int dim_x = src.shape()[2];
+  int dim_y = src.shape()[1];
+
+  float  thR = 0.1f, thG = 0.1f;
+
+  // to have the initial labels
+  int count = init(src,src_i,sfi,epsilon,skip_diffuse,check_nan_inf);
+  int pcount;
+
+  while(1)
+  {
+    for(y=0;y<dim_y-1;y++)
+    {
+	    for(x=0;x<dim_x-1;x++)
+      {
+
+        if(src_i(y,x)==CAMERA_DARK) continue;
+
+        //get the rgb values
+        float  r = src(0,y,x);
+        float  g = src(1,y,x);
+        float  b = src(2,y,x);