An Image Processing example of Using CUDA

The example shows the steps of  capturing video frames from camera and building a background. Several OpenCV functions are used.

In host part, main.cpp:

int *bgy;
int *bgy_CU;
BYTE *imY_CU;

dev = findCudaDevice(argc, (const char **)argv);
checkCudaErrors( cudaGetDeviceProperties(&deviceProp, dev) );
cvNamedWindow("video", 0);
    cvNamedWindow("BACKGROUND", 0);
    cvNamedWindow("BACKGROUND_GPU", 0);
    CvCapture *cap = cvCaptureFromCAM(0);
    IplImage* pFrame = 0; 
    IplImage *pFrame2 = cvCreateImage(cvSize(320, 240), IPL_DEPTH_8U, 3);
    imageY      = cvCreateImage(cvGetSize(pFrame2), 8, 1);
    imageBGY      = cvCreateImage(cvGetSize(pFrame2), 8, 1);

    bgy = new int [IMAGE_SIZE];    
    checkCudaErrors( cudaMalloc((void **)&bgy_CU,  IMAGE_SIZE*sizeof(int)  ) );
    checkCudaErrors( cudaMalloc((void **)&imY_CU,  IMAGE_SIZE*sizeof(BYTE)  ) );
    checkCudaErrors( cudaMalloc((void **)&imBY_CU, IMAGE_SIZE*sizeof(BYTE)  ) );

    while(pFrame = cvQueryFrame( cap ))
    {
        cvResize(pFrame, pFrame2, 1);
        cvShowImage("video", pFrame2);
        cvCvtColor(pFrame2, imageY, CV_BGR2GRAY);
        for (int i = 0; i < IMAGE_SIZE; i++)                      
        {
            imageY->imageData[i] = (unsigned char)pFrame2->imageData[i*3+1];
        }    
        imY = (BYTE*)imageY->imageData;
        imBY = (BYTE*)imageBGY->imageData;

        checkCudaErrors( cudaMemcpy(imY_CU, imY, IMAGE_SIZE*sizeof(BYTE), cudaMemcpyHostToDevice) );

        sdkCreateTimer(&hTimer);
        for(int iter = 0; iter < 10; iter++)
        {
            if(iter == 0){
                checkCudaErrors( cudaDeviceSynchronize() );
                sdkResetTimer(&hTimer);
                sdkStartTimer(&hTimer);
            }

            EstBackground_CUDA(imBY_CU, imY_CU, bgy_CU, IMAGE_SIZE);
            checkCudaErrors( cudaMemcpy(imBY, imBY_CU, IMAGE_SIZE * sizeof(BYTE), cudaMemcpyDeviceToHost) );
            cvShowImage("BACKGROUND_GPU", imageBGY);
        }
        checkCudaErrors( cudaDeviceSynchronize() );
        sdkStopTimer(&hTimer);
        double dAvgSecs = 1.0e-3 * (double)sdkGetTimerValue(&hTimer) / (double)10;
        shrLog("GPU time (average) : %.5f sec, \n\n", dAvgSecs * 1.0e6);
  
        if( cvWaitKey(1) >= 0 )
              break;
    }

    delete []bgy;
    cudaFree(bgy_CU);
    cudaFree(imY_CU);
    cudaFree(imBY_CU);

    cvDestroyWindow("video");
    cvDestroyWindow("BACKGROUND");
    cvDestroyWindow("BACKGROUND_GPU");
    cvReleaseImage(&pFrame2);
    cvReleaseCapture(&cap);

In foo.h:

extern "C" void EstBackground_CUDA(unsigned char *imBY_CU, unsigned char *imY_CU, int *bgy_CU, int IMAGE_SIZE);

In foo.cu:

__global__ void EstBackground_CUDAKernel(unsigned char *imBY_CU, unsigned char *imY_CU, int *bgy_CU, int IMAGE_SIZE)
{
    int  dist, sdist;
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    
    //cuPrintf("%d, %d, %d \n", threadIdx.x , blockIdx.x , blockDim.x);
    
    while (i < IMAGE_SIZE)
    {
        imBY_CU[i] = ...// run algorithms
        i += blockDim.x * gridDim.x;
    }    
    //__syncthreads(); // this is sync thread inside one block, no use here.
}
extern "C"  void EstBackground_CUDA(unsigned char *imBY_CU, unsigned char *imY_CU, int *bgy_CU, int IMAGE_SIZE)
{
    //cudaPrintfInit();
    EstBackground_CUDAKernel<<<256, 256>>>(imBY_CU, imY_CU, bgy_CU, IMAGE_SIZE);
    //cudaPrintfDisplay(stdout, true);
    //cudaPrintfEnd();
}

Print to console inside CUDA Kernel

Need to download “cuPrintf.cu” and “cuPrintf.cuh” at .
I am running “main.cpp” at Host and “foo.cu” at Device.
In main.cpp, call

foo( ... parameters ...);

In foo.h:

extern "C" void foo(... parameters ...);
extern "C" void test_print();

In foo.cu

__global__ void test_print()
{
    int tid;
    tid = blockIdx.x * blockDim.x + threadIdx.x;
    cuPrintf("%d\n", tid);
}
extern "C"  void foo(... parameters ...){
    cudaPrintfInit();
    test_print<<<32,8>>>();
    cudaPrintfDisplay(stdout, true);
    cudaPrintfEnd();
}

Note: This will affect the performance. In my application, the execute time raised from about 500ms to 1000ms.

HOG

Histogram of Oriented Gradients.

• Use derivative mask [1 0 –1], [1 0 –1]^T in one or both directions of the image. The gradient magnitude will be used.
• A Cell (ex. 6x6 pixels) and a Block (ex. 3x3 cells) are defined as following figure.
• Create oriented histogram bin: the bin can be spread over 0-180 degree if using signed values, and 0-360 degree if using unsigned values. For example, using a 20 degree wide bins, we get 9 bins. Every pixel (use gradient value) inside the cell is used as a weighted vote for the bins.
• Local normalization in block. The normalization factor could be: , where v is the non-normalized vector containing all histograms in a given block, ||v||2 is the 2-norm: , e is a small constant.
• The HOG descriptor is the vector of the components of the normalized cell histograms from all of the block region […,…,…]. The blocks are usually have 1/4 or 1/2 overlaps. 

HOG is widely used descriptor due to the speed of computations, inherent robustness to slight object variation/deformations, and the ability to capture a coarse spatial layout of features.

Update: 4/17/2012
Some implementation details:
Suppose the image size is 320x240, cell size is 8x8, so we have 40x30 blocks.  Suppose 18 orientations are used, we define two vectors for x,y directions:

double uu[9] = {1.0000, 

        0.9397, 

        0.7660, 

        0.500, 

        0.1736, 

        -0.1736, 

        -0.5000, 

        -0.7660, 

        -0.9397};

double vv[9] = {0.0000, 

        0.3420, 

        0.6428, 

        0.8660, 

        0.9848, 

        0.9848, 

        0.8660, 

        0.6428, 

        0.3420};

To check which orientation the current pixel locates at:

loop dx, dy within image size

double best_dot = 0;

int best_o = 0;

for (int o = 0; o < 9; o++) {

    double dot = uu[o]*dx + vv[o]*dy;

    if (dot > best_dot) {

      best_dot = dot;

      best_o = o;

    } else if (-dot > best_dot) {

      best_dot = -dot;

      best_o = o+9;

    }

}

Create the histograms for 18 orientations.  The pixel values in each orientation are added up.

Computer energy in each cell by summing over orientations

Output features. The features could be:

contrast-sensitive features (size 18): it is the pixel energy in each orientation

contrast-insensitive features (size 9): the energy in each orientation (no sign)

texture features (size 4): check the neighbor 4 pixel energy 

So the output features size is 40x30x31. (31=18+9+4)

Self-Organizing Map (SOM)

SOM is a data visualization technology that reduces the dimensions of data through the use of self-organizing neural network, to help us to understand the high dimensional data.

Initialize the map with random weight vectors.
for t = 0 ~ 1
    select a sample randomly from the set of training data
    every node is examined and find the best match unit   --- (a)
    choose neighbors and scale neighbors                          --- (b)
    increase t
end

(a): go through all the weight vector and calculate the distance of each weight to the sample.

(b): different methods to choose neighbors, such like Gaussian, or within a radius R. The new value is: current value * (1 –t) + sample vector * t.

Disadvantages: Need a value for each dimension of each sample. It is very computationally expensive.

Data mining (2)

The data classification process includes two steps: Learning and Classification. The class label of each training item is provided, is is called supervised learning. It contrasts with unsupervised learning (clustering), in which the class label of each training item is not known, and the number of set of classes to be learned may not be known in advance either.

Evaluate the classification and prediction method:
Accuracy, speed, robustness (even given noisy data or missing data), scalability (can applied to large amount of data), interpretability

Backpropagation (BP) is a neural network learning algorithm. The advantages of neural networks include the high tolerance of noisy data as well as the ability to classify patterns on which they have not been trained.

A multiplayer neural network includes input layer, hidden layer(s), and output layer. We call it two-layer neural network if there are only there three layers (input layer is not counted because it serves only to pass the input values to the next layer). If it contains two hidden layers, it is called a three-layer neural network.

Before training begins, we have to decide: number of units in input layer, number of hidden layers, number of units in each hidden layer, and the number of units in output layer. Normalize the input data will speed up the learning process.

SVM uses a nonlinear mapping to transform the original training data into a higher dimension. Within this new dimension, it searches the linear optimal separating hyperplane (decision boundary- separate the items of one class from another). With an appropriate nonlinear mapping to a sufficiently high dimension, data from two classes can always to be separated by a hyperplane. SVM finds this hyperplane using support vectors (some “essential” training items) and margins (defined by the support vectors). SVM searches for hyperplane with the largest margin, that is the maximum marginal hyperplane (MMH). The complexity of the learned classifier is decided by the number of support vectors rather than the dimensionality of the data. Hence, SVM is less sensitive to overfitting than other method. The support vectors are essential or critical training items, they lie closet to the decision boundary (MMH). If all the other training items are removed and repeat the training process, we get the same separating hyperplane. For nonlinar SVM, we can get it by extending the approach for linear SVM: first, transform the original input data into a higher dimensional space using a nonlinear mapping. 2nd, search for a linear separating hyperplane in the new spance. For example, a 3D input vector X={x1,x2,x3} is mapped to a 6D space Z using the mappings \phi_1(X)=x1, \phi_1(X)=x2… \phi_4(X)=(x1)^2, \phi_5(X)=x1x2, \phi_6(X)=x1x3. A decision hyperplane in the new spance is d(X)=WZ+b. Instead of computing the dot product on the transformed data items, it turns out that is is mathematically equivalent to apply a kernel function K(Xi, Xj)=\phi(Xi).\phi(Xj) – In other word, every \phi(Xi).\phi(Xj)  appears in the training algorithm, we can replace it with a kernel function  K(Xi, Xj). The the calculations are made in the original input space, which is much lower dimensionality.

Cross-validation

In k-fold cross validation, the initial data are randomly partitioned into k mutually exclusive subsets or 'folds', D1, D2, ...Dk, each of approximately equal size. Training and testing is performed k times. In iteration i, partition Di is reserved as the test set, and the remaining partitions are collectively used to train the model. So each sample is used the same number of times for training and once for testing. For classification, the accuracy estimate is the overall number of correct classification from the k iterations, divided by the total number of items in the initial data.
In general, 10-fold cross-validation is recommended for estimating accuracy due to its relatively low bias and variance.

Data mining (1)

• Data mining is more clear to be called “knowledge mining from data”. Knowledge extraction, data/pattern analysis.
• The knowledge extraction/discovery is a sequence of:
  1. Data cleaning (remove noise and inconsistent data)
  2. Data integration (from multi sources)
  3. Data selection
  4. Data transformation ( to specific form)
  5. Data mining (extract data pattern)
  6. Patter evaluation
  7. Knowledge presentation
• Classification is the process of finding a model/function that describes and distinguishes data classes/concepts, for the purpose of being able too use the model to predict the class of objects whose class label is unknown. The derived model is based on the analysis of a set of training data whose class label is known.
• Cluster analysis. Unlike classification and prediction, which analyze class-labeled data, clustering analyzes data objects without knowing a known class label.
• Data preparing, such as data normalization.
  

SVM basic

Linear Classification:
with all the data of x, x -> f() ->y, you want to separate them with: . There are many solutions there (w and b), which one is better? -- the one with the maximum margin (Margin Width ).

Support vectors are those data points that the margin pushed up against (on the border of the margin). The maximum margin implies that only the support vectors are important, the other training examples are ignorable.



The margin width
maximum is the same as minimum (dot production)
To solve w and b:
Minimize subject to
The solution is the quadratic optimization problem:

There are computation of inner products x_i^T x_j between all pairs of training points. 
The solution has the form:

Each non-zero \alpha_i indicates the corresponding x_i is a support vector. The classifying function is:

There are inner product between the test point x and the support vector point x_i.

Soft margin is to minimize
Linear classifier is a separating hyperplane, the support vectors (the most important training points) define the hyperplane

Non-linear SVM: mapping data to a higher-dimensional feature space where the training data is separable.
Linear classifier relies on dot production between vectors, the non-linear relies on kernel function. Kernel function is some function that corresponds to an inner product in feature space. Kernel function examples:
Linear:
Polynomial:
Gaussian(radial-basis function network), is the distance between closest points with different classification.
Non-linear SVM locates a separating hyperplane in the feature space and classify points in that space. It doesn't need to represent the space explicitly, simply by defining a kernel function. The kernel function plays the role of the dot product in feature space.

Weakness of SVM
sensitive to noise.
It only consider two classes (binary class). Didn't we use it to classify multiply categories? If you have m categories, you have m SVM leans. SVM1 leans "output ==1" vs "output != 1". ......SVM m leans "output == m" vs "output != m". In prediction part, predict the new input with each SVM and find out the best probability.
--most contents are from: http://www.cs.cmu.edu/~awm/tutorials

libsvm

1. Collect Training Data: (with your own codes to create TrainingData)
   

fprintf(fp, "3 "); // label: 1, 2 ...
for (j = 0; j < 12; j++)
fprintf(fp, "%d:%f ", j, mel_cep[j]);
fprintf(fp, "\n");

2. Scale training data:
   

svm-scale.exe -l -1 -u 1 -s range1 TrainingData > TrainingData.scale 

3. Get parameters and Train data with parameters
   

python grid.py TrainingData.scale

Got the best parameters c, g (2.0, 2.0)

svm-train.exe -c 2 -g 2 TrainingData.scale

The output is TrainingData.scale.model
4. Test/Run classificaiton: same way to get test data:
   

//fprintf(fp, "1 "); // if you know the label, can get the accuracy output with this
for (j = 0; j < 12; j++)
fprintf(fp, "%d:%f ", j, mel_cep[j]);
fprintf(fp, "\n");

Scale the test data with the same scaling factor with training data:

svm-scale.exe -r range1 test_libsvm.t > test_libsvm.t.scale

Run the classification/Prediction: (write your own codes here)

svm-predict.exe test_libsvm.t.scale TrainingData.scale.model test_libsvm.t.predict

===================
To use the cross validation to check the accuracy:

svm-train.exe -v 5 TrainingData
svm-train.exe -v 5 TrainingData.scale
svm-train.exe -c 2 -g 2 -v 5 TrainingData.scale

without normalize (scale), Cross Validation Accuracy = 87%
with scale (-1 1), Cross Validation Accuracy = 92%
with parameter selection –c 2 –g 2, Cross Validation Accuracy = 95%

To get the probability information, in training part, use:
svm-train.exe -c 2 -g 2 -b 1 TrainingData.scale
In test/run part, make "int predict_probability=1;" and the program will call "predict_label = svm_predict_probability(model,x,prob_estimates);" , which writes the probability value in the output file.

[update:8/29/2011]:
To run step 3 successfully, make sure that you have that grid.py in the same directory, and in grid.py line 22:
make sure gnuplot installed in correct location, and svmtrain_exe is in correct directory too (may copy to a 'window' folder in current directory).
In file TrainingData.scale.model, take a look at 'label'. It could be 6 2 3 4 5 1, instead of 1 2 3 4 5 6

3 image features for content based image retrieval

All are from MPEG-7.
Scalable Color Descriptor: a color histogram in HSV color space. For example, Hue is divided into 8 fuzzy areas: (0) Red to Orange, (1) Orange, (2) Yellow, (3) Green, (4) Cyan, (5) Blue, (6) Magenta and (7) Blue to Red. S is divided into 2 fuzzy regions, The first area, in combination with the position of the pixel in channel V, is used to define if the color is clear enough to be ranked in one of the categories which are described in H histogram, or if it is a shade of white or gray color.
Color Layout:
Edge Histogram Descriptor: the filters for edges from different angles:

• vertical: [1 –1; 1 –1]
• horizontal: [1 1; –1 –1]
• 45 diagonal: [sqrt(2) 0; 0 -sqrt(2)]
• 135 diagonal: [0 sqrt(2); -sqrt(2) 0]
• non-directional: [2 –2; –2 2]

the respective edge magnitudes are ‘m_xx’, each image block is derived to 4 sub blocks, the average gray level of each sub block is a_x(i,j), where x is 1,2,3,4, then m_xx = | sum_0^3 a_x(i,j) * f_x(k) |, where xx is vertical, horizontal, etc., and its max value is ‘max’. We get the edge histogram with:

if (max < TEdge) EdgeHist[0]++; // no edge
else                            // different direction
if (m_nd > T0) EdgeHist[1]++;
if (m_h  > T1) EdgeHist[2]++;
if (m_v  > T1) EdgeHist[3]++;
if (m_45 > T2) EdgeHist[4]++;
if (m_135> T2) EdgeHist[5]++;

Minimum Mean Square Error Filter (Wiener Filter)

or least square error filter.
It comes from inverse filter. The inverse filter doesn't handle noise. This method consider images and noise as random variables, and the objective is to find an estimate f' of the uncorrupted image f such that the mean square error between them is minimized
if the noise is zero, the Wiener filter reduces to the inverse filter.

The details and equations are at Digital Image Processing 3rd Ed. by Gonzalez, Page 353

mean shift

Get "Mass Center"

for(int i=0;i< height;i++)
for(int j=0;j width;j++)
M00+=I(i,j);

for(int i=0;i height;i++)
for(int j=0;j width;j++)
{
M10+=i*I(i,j);
M01+=j*I(i,j);
}

Mass Center is：

Xc=M10/M00;
Yc=M01/M00

Four steps for Mean Shift algorithm:
1.initialize the size and position of window.
2.calculate Mass Center of the window.
3.adjust Mass Center of the window.
4.repeat 2 and 3 till the center qualified.

OpenCV function: cvMeanShift
int cvMeanShift(IplImage* imgprob,CvRect windowIn,
CvTermCriteria criteria,CvConnectedComp* out);

RM

1. setup bk, this is an easy one coz the region is simple compared with others
2. Get edge inside monitored region. For each frame, white/black >a%? too many edges, need to change th. 3. Pixel based check all around one edge pixel, changed perc >%? Candidate target now, but removed or blocked?
4. check motion. if the corresponding part has motion there, it is a block.

DCT

• DCT is a Fourier-related transform similar to DFT, but using only real number.

• DCT is similar to the Fast Fourier Transform (FFT), but can approximate lines well with fewer coefficients (see following figure, the effect after IDCT and IDFT).
• The DCT concentrates most of the power on the lower frequencies.

Da


DFT:


DCT:


DCTs use only cosine functions, while DFTs use both cosine and sine functions.
The two-dimensional DCT is:





some of above are from http://en.wikipedia.org/wiki/Discrete_cosine_transform

Interpolatioin

Nearest neighbor, fast, drawback: produce undesirable artifacts, such as distortion of straight edges in images of high resolution. It produces a checkboard effect that is particularly objectionable at high factors of magnification.

cubic convolution interpolation
bilinear interpolation v(x,y) = ax + by + cxy + d

Opening, Closing

dilation,erosion
Opening smooths the contour of an object , breaks narrow isthmuses, and eliminates thin protrusions. Opening A by B is the erosion of A by B, followed by a dilation of B.
Closing also tends to smooth sections of contours but, as opposed to opening, it generally fuses narrow breaks and long thin gulfs, eliminates small holes, and fills gaps in contour. Closing of A by B is the dilation of A by B, followed by erosion of B.

One problem about 'minSize'

After minScale, 'minSize' could be zero and obversely it is not what we want. So we may use:

if (objSize_width>2 && objSize_height>2 && objSize_width>minScale && ... objSize

where '2' is 2 blocks

normalize cross correlation

c = normxcorr2(Template, imageA);
figure, surf(c), shading flat

normxcorr2 only works on grayscale images
http://www.mathworks.com/products/demos/image/cross_correlation/imreg.html

Off the grill

Finally off the grill, but go to next smaller grill again ...

AE autoexposure

The most important aspect of the exposure duration is to guarantee that the acquired image falls in a good region of the sensor’s sensitivity range. The imaging community uses a measure called exposure value (EV) to specify the relationship between the f-number, F, and exposure duration, T :
f-number or focal ration [tex]f/#=f/D[/tex], where f is the focal length, and D is the diameter of the entrance pupil.
[tex] EV = log_2 (F^2/T) = 2 log_2(F) - log_2(T) [/tex]
The exposure value becomes smaller as the exposure duration increases, and it becomes larger as the f-number grows.
Most auto exposure algorithms work this way:

Take a picture with a pre-determined exposure value (EVpre)
Convert the RGB values to brightness, B.
Derive a single number Bpre (like center-weighted mean, median, or more complicated weighted method as in matrix-metering) from the brightness picture
Based on linearity assumption and equation (1), the optimum exposure value EVopt should be the one that, the picture we take at this EVopt will give us a number close to an pre-defined ideal value Bopt, or: [tex] EV_{opt} = EV_{pre} + log_2(B_{pre}) - log_2(B_{opt}).

The ideal value Bopt for each algorithm is typically selected empirically. For the moment, however, let’s assume that Bopt is known.

Different algorithms mainly differ in how they derive the single number Bpre from the picture. Some simple algorithms include:

Mean: Bpre is the mean brightness across the whole picture.

Center-Weighted Mean: Bpre is the weighted mean of the center area and the rest area. It puts more weight on the center part than the surrounding area. There are many alternatives as how you choose the center area and how much weight for it. Here we choose center to be the center 25% area and weights [0.8 0.2] for center and surrounding respectively.

Spot: Bpre is the mean of the center 3% area.

Median: Bpre is the median brightness of the whole picture.

Green: Bpre is the mean of the green channel only.