From 5f16fde64a9f48f291a53095e78932dfbff1106b Mon Sep 17 00:00:00 2001
From: David Geissbuhler <david.geissbuhler@idiap.ch>
Date: Wed, 28 Jun 2017 16:22:25 +0200
Subject: [PATCH] added modified c++ algorithm
---
bob/ip/qualitymeasure/__init__.py | 1 +
bob/ip/qualitymeasure/main.cpp | 17 +-
bob/ip/qualitymeasure/main_orig.cpp | 128 +++++
.../script/remove_highlights.py | 4 +-
.../tan_specular_highlights.cpp | 457 ++++++++++++++++++
setup.py | 16 +-
6 files changed, 611 insertions(+), 12 deletions(-)
create mode 100644 bob/ip/qualitymeasure/main_orig.cpp
create mode 100644 bob/ip/qualitymeasure/tan_specular_highlights.cpp
diff --git a/bob/ip/qualitymeasure/__init__.py b/bob/ip/qualitymeasure/__init__.py
index f9dbfca..3b35a74 100644
--- a/bob/ip/qualitymeasure/__init__.py
+++ b/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():
diff --git a/bob/ip/qualitymeasure/main.cpp b/bob/ip/qualitymeasure/main.cpp
index fbe072a..6dd3b12 100644
--- a/bob/ip/qualitymeasure/main.cpp
+++ b/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);
diff --git a/bob/ip/qualitymeasure/main_orig.cpp b/bob/ip/qualitymeasure/main_orig.cpp
new file mode 100644
index 0000000..f964ccd
--- /dev/null
+++ b/bob/ip/qualitymeasure/main_orig.cpp
@@ -0,0 +1,128 @@
+// 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();
+}
diff --git a/bob/ip/qualitymeasure/script/remove_highlights.py b/bob/ip/qualitymeasure/script/remove_highlights.py
index a314e07..17610af 100644
--- a/bob/ip/qualitymeasure/script/remove_highlights.py
+++ b/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)
diff --git a/bob/ip/qualitymeasure/tan_specular_highlights.cpp b/bob/ip/qualitymeasure/tan_specular_highlights.cpp
new file mode 100644
index 0000000..8054f7a
--- /dev/null
+++ b/bob/ip/qualitymeasure/tan_specular_highlights.cpp
@@ -0,0 +1,457 @@
+/**
+ * @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));
+ float sfi_tot_x = tot(sfi(0,y,x+1), sfi(1,y,x+1), sfi(2,y,x+1));
+ float sfi_tot_y = tot(sfi(0,y+1,x), sfi(1,y+1,x), sfi(2,y+1,x));
+
+ float dlog_src_x = log( fabs( src_tot_x - src_tot_0 ) );
+ float dlog_src_y = log( fabs( src_tot_y - src_tot_0 ) );
+
+ float dlog_sfi_x = log( fabs( sfi_tot_x - sfi_tot_0 ) );
+ float dlog_sfi_y = log( fabs( sfi_tot_y - sfi_tot_0 ) );
+
+ float dlogx = (dlog_src_x - dlog_sfi_x);
+ float dlogy = (dlog_src_y - dlog_sfi_y);
+
+ dlogx=fabs(dlogx);
+ dlogy=fabs(dlogy);
+
+ // specular in the x direction
+ if(dlogx > epsilon || isinf(dlog_src_x) || isinf(dlog_sfi_x))
+ {
+ src_i(y,x) = SPECULARX;
+ count++;
+ continue; // go to the next pixel
+ }
+
+ //specular in the y direction
+ if(dlogy > epsilon || isinf(dlog_src_y) || isinf(dlog_sfi_y))
+ {
+ src_i(y,x)= SPECULARY;
+ count++;
+ continue;
+ }
+
+ src_i(y,x) = DIFFUSE;
+
+ }
+ }
+
+ return count; // return the number of specular pixels
+}
+
+// to reset the label of the pixels
+
+int reset_labels(blitz::Array<int,2> &src_i)
+{
+ int dim_x = src_i.shape()[1];
+ int dim_y = src_i.shape()[0];
+
+ for(int y=0;y<dim_y;y++)
+ {
+ for(int x=0;x<dim_x;x++)
+ {
+ if(src_i(y,x)==CAMERA_DARK) continue;
+ src_i(y,x)=0;
+ }
+ }
+ return 0;
+}
diff --git a/setup.py b/setup.py
index 8371b84..e2a2b74 100644
--- a/setup.py
+++ b/setup.py
@@ -58,7 +58,7 @@ setup(
# list of files compiled into this extension
[
# the pure C++ code
- "bob/ip/qualitymeasure/tan_specular_highlights_orig.cpp",
+ "bob/ip/qualitymeasure/tan_specular_highlights.cpp",
# the Python bindings
"bob/ip/qualitymeasure/main.cpp",
],
@@ -66,6 +66,20 @@ setup(
version = version,
bob_packages = bob_packages,
),
+
+ # The second extension contains the actual C++ code and the Python bindings
+ Extension("bob.ip.qualitymeasure._library_orig",
+ # list of files compiled into this extension
+ [
+ # the pure C++ code
+ "bob/ip/qualitymeasure/tan_specular_highlights_orig.cpp",
+ # the Python bindings
+ "bob/ip/qualitymeasure/main_orig.cpp",
+ ],
+ # additional parameters, see Extension documentation
+ version = version,
+ bob_packages = bob_packages,
+ ),
],
# Important! We need to tell setuptools that we want the extension to be
--
GitLab