added modified c++ algorithm

bob/ip/qualitymeasure/__init__.py

@@ -5,6 +5,7 @@ from .galbally_iqm_features import compute_quality_features
 from .msu_iqa_features import compute_msu_iqa_features
 
 from ._library import remove_highlights
+from ._library_orig import remove_highlights_orig
 
 
 def get_config():

bob/ip/qualitymeasure/main.cpp

@@ -9,13 +9,12 @@
 #include <bob.extension/documentation.h>
 
 
-// declare C++ functions
-void remove_highlights_orig(  blitz::Array<float ,3> &img,
-                              blitz::Array<float ,3> &diff,
-                              blitz::Array<float ,3> &sfi,
-                              blitz::Array<float ,3> &residue,
-                              float  epsilon);
-
+// 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);
 
 // declare the function
 static PyObject* PyRemoveHighlights(PyObject*, PyObject* args, PyObject* kwargs) {
@@ -26,7 +25,7 @@ static PyObject* PyRemoveHighlights(PyObject*, PyObject* args, PyObject* kwargs)
   static char** kwlist = const_cast<char**>(const_kwlist);
 
   PyBlitzArrayObject* array;
-  double epsilon = 0.5f;
+  double epsilon  = 0.5f;
 
   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&|d", kwlist,
                                   &PyBlitzArray_Converter, &array,
@@ -55,7 +54,7 @@ static PyObject* PyRemoveHighlights(PyObject*, PyObject* args, PyObject* kwargs)
   speckle_img       = 0;
 
   // call the C++ function
-  remove_highlights_orig(img, diffuse_img, speckle_free_img, speckle_img, (float)epsilon);
+  remove_highlights(img, diffuse_img, speckle_free_img, speckle_img, (float)epsilon);
 
   // convert the blitz array back to numpy and return it
   PyObject *ret_tuple = PyTuple_New(3);

bob/ip/qualitymeasure/main_orig.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_orig(  blitz::Array<float ,3> &img,
+                              blitz::Array<float ,3> &diff,
+                              blitz::Array<float ,3> &sfi,
+                              blitz::Array<float ,3> &residue,
+                              float  epsilon);
+
+// declare the function
+static PyObject* PyRemoveHighlightsOrig(PyObject*, PyObject* args, PyObject* kwargs) {
+
+  BOB_TRY
+
+  static const char* const_kwlist[] = {"array", "startEps", 0};
+  static char** kwlist = const_cast<char**>(const_kwlist);
+
+  PyBlitzArrayObject* array;
+  double epsilon  = 0.5f;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&|d", kwlist,
+                                  &PyBlitzArray_Converter, &array,
+                                  &epsilon)) 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_orig(img, diffuse_img, speckle_free_img, speckle_img, (float)epsilon);
+
+  // 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_orig", 0)
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+/////// Python module declaration ////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+// module-wide methods
+static PyMethodDef module_methods[] = {
+  {
+    "remove_highlights_orig",
+    (PyCFunction)PyRemoveHighlightsOrig,
+    METH_VARARGS|METH_KEYWORDS,
+    "remove_highlights [doc]"
+  },
+  {NULL}  // Sentinel
+};
+
+// module documentation
+PyDoc_STRVAR(module_docstr, "Exemplary Python Bindings");
+
+// 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/script/remove_highlights.py

@@ -13,7 +13,7 @@ import argparse
 import bob.io.base
 import numpy as np
 
-from bob.ip.qualitymeasure import remove_highlights
+from bob.ip.qualitymeasure import remove_highlights_orig
 
 def main(command_line_parameters=None):
     """Remove the specular component of the input image and write result to
@@ -51,7 +51,7 @@ def main(command_line_parameters=None):
     # 2. compute
     print("Extracting diffuse component...")
 
-    sfi, diff, residue = remove_highlights(img.astype(np.float32), 0.5)
+    sfi, diff, residue = remove_highlights_orig(img.astype(np.float32), 0.5)
 
     # 1. save output image
     print("Saving output file: %s" % args.outImg)

bob/ip/qualitymeasure/tan_specular_highlights.cpp

+/**
+ * @author David Geissbuhler <andre.anjos@idiap.ch>
+ * @date Tue 27 Jun 15:54:00 2016
+ *
+ *  Modified 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
+ *       http://tanrobby.github.io/code.html#
+ *       http://tanrobby.github.io/code/highlight.zip
+ *
+ *  This is a modified implementation based on the C++ code provided by Prof.
+ *  Robby Tan but using Blitz++ arrays and with a modification that also
+ *  ignores pixels marked G_DIFFUSE. This leads to a smaller number of
+ *  iterations per epsilon value, while producing the similar results.
+ *
+ */
+
+#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);
+
+void iteration(           blitz::Array<float ,3> &src,
+                          blitz::Array<int,2>   &src_i,
+                          blitz::Array<float ,3> &sfi,
+				                  float  epsilon);
+
+int init(                 blitz::Array<float ,3> &src,
+                          blitz::Array<int,2>   &src_i,
+                          blitz::Array<float ,3> &sfi,
+				                  float  epsilon);
+
+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)
+{
+  // 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);
+
+  //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);
+    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)
+{
+  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;
+
+      if(denm == 0) dI = 0;
+      else dI = numr / denm;
+
+      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)
+{
+  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(fabs(denm) > 0.000000001)
+  {
+    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;
+    }
+  }
+  else
+  {
+    iro_i(y,x)=NOISE;
+    return 1;
+  }
+}
+
+// specular reduction mechanism
+
+void iteration( blitz::Array<float ,3> &src,
+                blitz::Array<int,2>   &src_i,
+                blitz::Array<float ,3> &sfi,
+				        float  epsilon)
+{
+  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);
+  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);
+
+        float  cr = r_chroma(r,g,b);		// red chroma
+        float  cg = g_chroma(r,g,b);		// green chroma
+
+        //get the rgb values
+        float  rx = src(0,y,x+1);
+        float  gx = src(1,y,x+1);
+        float  bx = src(2,y,x+1);
+
+        float  cr_next_x = r_chroma(rx,gx,bx); // red chroma
+        float  cg_next_x = g_chroma(rx,gx,bx); // green chroma
+
+        //get the rgb values
+        float  ry = src(0,y+1,x);
+        float  gy = src(1,y+1,x);
+        float  by = src(2,y+1,x);
+
+        float  cr_next_y = r_chroma(ry,gy,by); // red chroma
+        float  cg_next_y = g_chroma(ry,gy,by); // green chroma
+
+        // derivatives
+        float  drx = cr_next_x-cr;//pixel right
+        float  dgx = cg_next_x-cg;
+        float  dry = cr_next_y-cr;//pixel below
+        float  dgy = cg_next_y-cg;
+
+        if(src_i(y,x) == SPECULARX)
+        {
+          //if it is  a boundary in the x direction
+          if(fabs(drx) > thR && fabs(dgx) > thG)
+          {
+            //pixel right
+            src_i(y,x)=BOUNDARY;
+            continue;
+          }
+
+          //if it is a noise
+          if(fabs(max_chroma(r,g,b) - max_chroma(rx,gx,bx)) < 0.01)
+          {
+            src_i(y,x)=NOISE;
+            continue;
+          }
+
+          //reduce the specularity at x direction
+          if(max_chroma(r,g,b) < max_chroma(rx,gx,bx))
+          {
+            specular_2_diffuse(y,x,src,src_i,max_chroma(rx,gx,bx));
+            src_i(y,x)=DIFFUSE;
+            src_i(y,x+1)=DIFFUSE;
+          }
+          else
+          {
+            specular_2_diffuse(y,x+1,src,src_i,max_chroma(r,g,b));
+            src_i(y,x)=DIFFUSE;
+            src_i(y,x+1)=DIFFUSE;
+          }
+        }
+
+	      if(src_i(y,x) == SPECULARY)
+	      {
+	        //if it is a boundary in the y direction
+	        if(fabs(dry) > thR && fabs(dgy) > thG)
+          {
+            //pixel right
+		        src_i(y,x)=BOUNDARY;
+		        continue;
+	        }
+
+	        //if it is a noise
+	        if(fabs(max_chroma(r,g,b) - max_chroma(ry,gy,by))<0.01)
+          {
+		        src_i(y,x)=NOISE;
+		        continue;
+	        }
+
+	        //reduce the specularity in y direction
+	        if(max_chroma(r,g,b) < max_chroma(ry,gy,by))
+          {
+            specular_2_diffuse(y,x,src,src_i,max_chroma(ry,gy,by));
+            src_i(y,x)=DIFFUSE;
+            src_i(y+1,x)=DIFFUSE;
+          }
+          else
+          {
+            specular_2_diffuse(y+1,x,src,src_i,max_chroma(r,g,b));
+            src_i(y,x)=DIFFUSE;
+            src_i(y+1,x)=DIFFUSE;
+          }
+	      }
+	    }
+    }
+
+    pcount=count;
+    count = init(src,src_i,sfi,epsilon);
+
+    if(count==0)
+	    break;
+    if(pcount<=count)
+	    break;
+  }
+
+  reset_labels(src_i);
+}
+
+// to have initial labels
+
+int init(   blitz::Array<float ,3> &src,
+            blitz::Array<int,2>   &src_i,
+            blitz::Array<float ,3> &sfi,
+				    float  epsilon)
+{
+  int dim_x = src.shape()[2];
+  int dim_y = src.shape()[1];
+  int x,y;	// pixel iterators
+
+  int count=0;
+
+  for(y = 1; y < dim_y - 1; y++){
+    for(x = 1; x < dim_x - 1; x++){
+      switch(src_i(y,x))
+      {
+        case BOUNDARY:
+        case NOISE:
+        case CAMERA_DARK:
+        case DIFFUSE:
+	        continue;
+	      break;
+      }
+
+      float  src_tot_0 = tot(src(0,y,x), src(1,y,x), src(2,y,x));
+      float  src_tot_x = tot(src(0,y,x+1), src(1,y,x+1), src(2,y,x+1));
+      float  src_tot_y = tot(src(0,y+1,x), src(1,y+1,x), src(2,y+1,x));
+
+      float  sfi_tot_0 = tot(sfi(0,y,x), sfi(1,y,x), sfi(2,y,x));