Category Archives: Programming

In-focus Region Detection Using Harmonic Variance of DCT Band-pass Filtering

Detecting in-focus regions in a low depth-of-field (DoF) image has important applications in scene understanding, object-based coding, image quality assessment and depth estimation because such regions may indicate semantically meaningful objects.

Here I show a challenging example,  a spider web. Note that some parts of the defocus regions still demonstrate strong patch variances.

spider-web-macro-photography

In our paper, we show that by simply using DCT band-pass filtering and harmonic variance, we could get very robust in-focus measurement, as shown in the following picture. The idea is simple. Defocus regions only cover some low frequency bands because of lens blur, while in-focus region cover the full spectrum bands. By using DCT band-pass filtering, we could extract those frequency responses. If we just sum those responses all together, the final estimate could be well biased due to very high responses of just some low frequency channels. Harmonic mean solves this problem by penalizing those outliers and could generate very robust outputs.

hv

In short, this algorithm works by
1) apply 63 DCT band-pass filtering (for 8×8 patches) to the input image, excluding the first DC/average filtering
2) compute patch variance for each 63 filtered images. Now at every pixel location (x, y), we have 63 variance measurements
3) compute harmonic mean of these 63 variance as harmonic variance
4) optionally, we could apply sigmoid function to the harmonic variance image to restrict the outputs to be [0, 1].

If you want to estimate the actual noise level of the input image, please refer to our paper for kurtosis analysis.

import numpy as np
import scipy.misc as misc
import scipy.ndimage as ndimage
import scipy.stats as stats
import matplotlib.pyplot as plt

def Sigmoid(x, ratio) :
    return 1 / (1 + np.exp(-ratio * x) )

def Tanh(x, ratio) :
    return 2 * Sigmoid(2 * x, ratio) - 1

def GenerateDCT8x8Base() :

    N = 8
    x = np.array( np.arange(0, N ), dtype=np.float32, ndmin=2 )

    C = np.cos( np.transpose(2*x+1) * x * np.pi / (2.0*N) ) * np.sqrt(2.0/N)
    C[:, 0] = C[:, 0] / np.sqrt(2.0)

    plt.close('all')
    plt.figure()

    filters = []
    for i in range( N ) :
        for j in range( N ) :

            X = np.transpose(np.array(C[:, i], ndmin=2)) * np.transpose( C[:, j])
            im = misc.imresize(X, (128, 128), interp='nearest')
            plt.subplot(N, N, i * N + j + 1)
            plt.imshow(im, cmap=plt.cm.gray)
            plt.axis('off')

            filters.append( X )

    del filters[0] 

    return filters

def ComputeHarmonicVariance(im, kernels):

    size = list(im.shape)
    size.append(len(kernels))

    k_size = kernels[0].shape
    k_avg = np.ones( k_size, np.float32) / (k_size[0] * k_size[1])

    outs = np.empty( size, dtype=np.float32)
    for i in range(len(kernels)) :
        data = ndimage.convolve(im, kernels[i])
        data2 = data * data

        EX2 = ndimage.convolve(data2, k_avg)
        E2X = np.square(ndimage.convolve(data, k_avg))

        outs[:,:,i] = EX2 - E2X

    variance = np.empty( im.shape, dtype=np.float32 )
    for y in range(size[0]):
        for x in range(size[1]):
            data = outs[y, x, :]
            variance[y, x] = stats.hmean(outs[y,x,:]+1e-8)

    return variance            

im = misc.imread('Spider-Web-Macro-Photography.jpg').astype(np.float32)
im = im[:,:,0] # extract Red channel for demonstration
kernels = GenerateDCT8x8Base()
hv = ComputeHarmonicVariance(im, kernels)
plt.figure(), plt.imshow(Tanh(hv, 0.01), cmap=plt.cm.gray), plt.title('hv t 0.01')
Advertisement

Multiply Strings

Given two numbers represented as strings, return multiplication of the numbers as a string.
Note: The numbers can be arbitrarily large and are non-negative.

string stringMultip(string &A, string &B) {

	string res(A.size() + B.size(), '0');
	
	string::reverse_iterator rit1, rit2, rit3, rit4;
	int num_a, num_b, num_c, num_carry, num_res, i;
	
	for(rit1 = A.rbegin(), rit3 = res.rbegin(); rit1 != A.rend(); rit1 ++, rit3 ++) {
		
		num_a = (*rit1) - '0';
		num_carry = 0;
		
		for(rit2 = B.rbegin(), rit4 = rit3; rit2 != B.rend(); rit2 ++, rit4 ++) {
			
			num_b = (*rit2) - '0';
			num_c = (*rit4) - '0';
			
			num_res = num_a * num_b + num_c + num_carry;
			num_carry = num_res / 10;
			(*rit4) = '0' + (num_res % 10);
		}

		// put the carry back to the result
		(*rit4) = '0' + num_carry;
	}

	// remove leading '0's
	for(i = 0; i < res.size(); i ++)
		if(res[i] != '0') 
			break;

	return res.erase(0, i);
}

The following provides a simple test case.

#include <iostream>
#include <string>

using namespace std;

int main ()
{
	string A("123456");
	string B("294987");
	string C("123");
	string D("234");	

	cout << stringMultip(A, B) <<endl;
	cout << stringMultip(C, D) <<endl;
	
	return 0;
}

Binary Tree Postorder Traversal using Stack

The problem:
Given a binary tree, return the postorder traversal of its nodes’ values. For example, given binary tree:
a              1
          /      \
       4           2
         \       /
e           10  3
return [10, 4, 3, 2, 1]

#include <iostream>
#include <stack>
#include <string>

using namespace std;

struct TreeNode {
    int val;
    TreeNode *left;
    TreeNode *right;
    TreeNode(int x): val(x), left(nullptr), right(nullptr) {};
};

int postorderTraversal(TreeNode *root) {

    const TreeNode *p = root;
    stack<const TreeNode*> s;

    if(p!=nullptr) {
        TreeNode x(p->val);
        s.push(&x);

        if(p->right != nullptr) s.push(p->right);
        if(p->left != nullptr) s.push(p->left);
    }

    while(!s.empty()) {

        p = s.top();
        s.pop();

        if(p->left == nullptr && p->right == nullptr)
            printf("%d\n", p->val);
        else {
            TreeNode x(p->val);
            s.push(&x);

            if(p->right != nullptr) s.push(p->right);
            if(p->left != nullptr) s.push(p->left);
        }
    }

    return 0;
}

int main() {

    TreeNode root(1);
    TreeNode node1(2);
    TreeNode node2(3);
    TreeNode node3(4);
    TreeNode node4(10);

    root.left = &node3;
    root.right = &node1;
    node1.left = &node2;
    node3.right = &node4;

    postorderTraversal(&root);

    return 0;
}

The similar code structure could be used for binary tree inorder traversal. The key idea here is to create a temporary node which holds only the value and push it into the stack. Whenever we meet a node without both childern, we print out its value.

Longest Valid Parentheses

The Problem:

Given a string containing just the characters ‘(‘ and ‘)’, find the length of the longest valid (well formed) parentheses sub-strings.

#include <iostream>
#include <string>
#include <stack>

using namespace std;

int longest_valid_parentheses(string const& s) {

    int count = 0, max_count = 0;
    stack<char> stk;

    for(auto c : s) {

        if(c == '(')
            stk.push(c);
        else {
            if( stk.empty() ) //|| stk.top() != '(')
                count = 0;
            else {
                count += 2;
                stk.pop();
            }
        }

        max_count = count > max_count ? count : max_count;
    }

    return max_count;
}

You can test the above code using the following test cases @ http://www.cpp.sh/.

int main()
{
  string test1("(()");
  string test2("(()()(");
  string test3(")()()))(((())))");

  cout << "there are " << longest_valid_parentheses(test1) << " valid parentheses in " << test1 << endl;
  cout << "there are " << longest_valid_parentheses(test2) << " valid parentheses in " << test2 << endl;
  cout << "there are " << longest_valid_parentheses(test3) << " valid parentheses in " << test3 << endl;

  return 0;
}

The output would look like as the follows:

there are 2 valid parentheses in (() 
there are 4 valid parentheses in (()()( 
there are 8 valid parentheses in )()()))(((())))

Remove all the duplicates from a sorted linked list

Given a sorted linked list, remove all the duplicates from the list. For example, given 1->2->3->3->6, return 1->2->6. Another example, given 1->1->1->2->4, return 2->4.

typedef struct node Node;

/* The function removes duplicates from a sorted list */
void removeDuplicates(Node** head)
{ 
  /* do nothing if the list is empty */
  Node* current = *head;  
  if(current == NULL)
     return;

  /* Pointer to traverse the linked list */
  Node* pre_head = (Node*) malloc(sizeof(Node));
  pre_head->next = *head;

  /* Pointer to store the next pointer of a node to be deleted*/
  Node* next_next;
  Node* pre = pre_head;
  bool deleted = false;
 
  /* Traverse the list till last node */
  while(current->next != NULL)
  {
   /* Compare current node with next node */
    if(current->data == current->next->data) {       
		/*The sequence of steps is important*/              
		next_next = current->next->next;
		free(current->next);
		current->next = next_next;
		deleted = true;
    }
    else if (deleted)  {
		pre->next = current->next;
		free(current);
		current = pre->next;
		deleted = false;
    }
	else {
		pre = current;
		current = current->next;
    }
  }

  /* remove the last duplicated node*/
  if (deleted) {
	  pre->next = NULL;
	  free(current);
  }
  
  *head = pre_head->next;
  free(pre_head);
}
 

To test the code, please use the testing program from http://www.geeksforgeeks.org/remove-duplicates-from-a-sorted-linked-list/.

Dutch national flag problem – matlab

Dutch national flag problem is a classic programming problem. It groups the input array into three sub-groups in just one pass. For detailed information, please refer to http://en.wikipedia.org/wiki/Dutch_national_flag_problem

function x = category_sort(x, low, high)

len = length(x);
p = 1; q = len;
i = 1;

while i <= q
    if x(i) < low     % move the current element to the front
        tmp = x(p);
        x(p) = x(i);
        x(i) = tmp;

        p = p + 1;  
        i = i + 1;     
    elseif x(i) >= high     % move the current element to the back
        tmp = x(q);
        x(q) = x(i);
        x(i) = tmp;
                
        q = q - 1;
    else    % move to the next element in the array
        i = i + 1;
    end    
end

We can test the code by simply running the following script:

x = 3 * rand(20, 1);
x = category_sort(x, 1, 2)

Implement strtok in C

To implement the C function ‘strtok’, one doesn’t need to allocate extra memory for the input string before modifying it, since the caller function is supposed to make copy of the input string and deallocate the original input string. The function only needs to define static variables to store the starting of the string.

The caller function provides a string as the input during the first call, and ‘null’ for the following function calls. For details about how this function works, please refer to http://www.cplusplus.com/reference/cstring/strtok/.

char* sp = NULL; /* the start position of the string */

char* strtok1(char* str, const char* delimiters) {

	int i = 0;
	int len = strlen(delimiters);

	/* check in the delimiters */
	if(len == 0)
		printf("delimiters are empty\n");

	/* if the original string has nothing left */
	if(!str && !sp)
		return NULL;

	/* initialize the sp during the first call */
	if(str && !sp)
		sp = str;

	/* find the start of the substring, skip delimiters */
	char* p_start = sp;
	while(true) {
		for(i = 0; i < len; i ++) {
			if(*p_start == delimiters[i]) {
				p_start ++;
				break;
			}
		}

		if(i == len) {
		       sp = p_start;
		       break;
		}
	}

	/* return NULL if nothing left */
	if(*sp == '\0') {
		sp = NULL;
		return sp;
	}

	/* find the end of the substring, and
        replace the delimiter with null */
	while(*sp != '\0') {
		for(i = 0; i < len; i ++) {
			if(*sp == delimiters[i]) {
				*sp = '\0';
				break;
			}
		}

		sp ++;
		if (i < len)
			break;
	}

	return p_start;
}

This following is the test code from cplusplus.com.

#include <stdio.h>
#include <string.h>

int main ()
{
  char str[] ="- This, a sample string.";
  char * pch;
  printf ("Splitting string \"%s\" into tokens:\n",str);
  pch = strtok1(str," ,.-");
  while (pch != NULL)
  {
    printf ("%s\n",pch);
    pch = strtok1(NULL, " ,.-");
  }
  return 0;
}

————————-
As pointed out by Ian in the following comments [9.25.2014], there is a much simpler solution to this problem. Feel free to test it and let us know if it has any bugs.

char* strtok1(char *str, const char* delim) {
	static char* _buffer;
	if(str != NULL) _buffer = str;
	if(_buffer[0] == '\0') return NULL;

	char *ret = _buffer, *b;
	const char *d;

	for(b = _buffer; *b !='\0'; b++) {
		for(d = delim; *d != '\0'; d++) {
			if(*b == *d) {
				*b = '\0';
				_buffer = b+1;

				// skip the beginning delimiters
				if(b == ret) { 
					ret++; 
					continue; 
				}
				return ret;
			}
		}
	}

	return ret;
}

Forming a 2D window from a 1D Function

I came to this question because of the Kaiser window used for my BM3D implementation. I used Matlab to generate the Kaiser window I want for my C code so I don’t bother to write that function in my project. However, Matlab only has 1D version of the window functions. The following demonstrates two ways to generate the 2D window from its 1D function.

——————————————————————-

There are basically two ways of forming a 2D window from a 1D function. The
first is the outer product formulation; the second is the rotational
formulation. To create a outer product window (which your kaiser examples
look like they are trying to do), you can use

w1D = hamming(16); % Some 1D window
w2D = w1D(:) * w1D(:).’; % Outer product

To use a rotational formulation (which generally gives more circular
contours), you can use

L = 16;
w1D = hamming(L); % Some 1D window
M = (L-1)/2;
xx = linspace(-M,M,L);
[x,y] = meshgrid(xx);
r = sqrt( x.^2 + y.^2 );
w2D = zeros(L);
w2D(r<=M) = interp1(xx,w1D,r(r<=M));

http://www.mathworks.com/matlabcentral/newsreader/view_thread/23588

Using opencv from python on windows

http://stackoverflow.com/questions/4709301/installing-opencv-on-windows-7-for-python-2-7

The official OpenCV2.2 installer does not install the Python bindings into your Python directory. There should be a Python2.7 directory inside your OpenCV 2.2.0 installation directory. Copy the whole Lib folder from OpenCV\Python2.7\ to C:\Python27\. Then set the windows environmental variable PATH to include your OpenCV\bin directory together with C:\Python.

Alternatively use the opencv-python installers at http://www.lfd.uci.edu/~gohlke/pythonlibs/#opencv.

Topcoder Algorithm Tutorials

http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=alg_index

Topcoder has many excellent tuorials for algorithms that we meet everyday, for example, dijkstra and maxflow, which all were written by topcoder members. Usually, at the end of each tutorial, the author would often explain some applications for this specific kind of algorithm, and most of the applications are from the topcoder’s competition problems. It’s a nice place to start with if you want to improve your programming skills.