Use Faster RCNN and ResNet codes for object detection and image classification with your own training data

I have recently uploaded two repositories to GitHub, both based on publicly available codes for state-of-the-art (1) object detection and (2) image classification. I would like to leave a few notes here, though.

(1) Faster RCNN for object detection (GitHub Link).

You can use your own PASCAL VOC formatted data to train an object detector. Check out how to alter the network parameters as shown in the example files located in:

person_detection_voc2012/py-faster-rcnn/models/pascal_voc/ZF/faster_rcnn_alt_opt/*.pt

In particular, you want to change the following settings in stage1_fast_rcnn_train.pt and stage2_fast_rcnn_train.pt:

num_class:2 # in our example, person detection only has two classes: person vs background
In cls_score -- num_output:2
In bbox_pred -- num_output:8 # this value is 4*num_class

Also in stage1_rpn_train.pt and stage2_rpn_train.pt:

num_class:2

Finally, in fast_rcnn_test.pt:

In cls_score -- num_output:2
In bbox_pred -- num_output:8 # this value is 4*num_class

Additionally, you need to modify lib/datasets/pascal_voc.py:

self._classes = ('__background__', # always index 0
                 'person')

And then recompile from python prompt:

importpy_compile
py_compile.compile(r'pascal_voc.py')

You can then follow instructions from this page to train your model.

(2) Fine tuning ResNet for image classification (GitHub Link).

This one is simple to use, and you may check this out before attempting to fine tune a ResNet model.

Example scripts can be found in: finetune-resnet-flower/caffe/examples/flower463/

Network parameters can be found in: finetune-resnet-flower/caffe/models/resnet_flower463/

Note that the parameters in solver50.prototxt may not be optimal (at least for my task at hand). For better performance (of course, slower training), you can try to increase stepsize as shown below:

test_iter: 2000
test_interval: 1000
base_lr: 0.001
lr_policy: "step"
gamma: 0.1
stepsize: 100000
display: 500
max_iter: 1000000
momentum: 0.9
weight_decay: 0.0005

Also, set the batch size appropriately to reflect the graphic memory capability of your system.

Annotated feature convolution code from Deformable Parts Model

The following code is an excerpt from fconv_var_dim.cc that comes with the following paper:

A Discriminatively Trained, Multiscale, Deformable Part Model.
P. Felzenszwalb, D. McAllester, D. Ramaman.
Proceedings of the IEEE CVPR 2008.

The original code was downloaded from: http://people.cs.uchicago.edu/~pff/latent/

I only added some comments to make it easier to read. The comments illustrate how the data structures are designed (should also help understanding the multi-threaded and SSE-optimized implementations).

#include "mex.h"
#include <math.h>
#include <string.h>

/*
 * This code is used for computing filter responses.  It computes the
 * response of a set of filters with a feature map.  
 *
 * Basic version, relatively slow but very compatible.
 */

/*

The dimensions for A, B and C are as follows. For example:

A = A1*A2*32 single
B <2x1 cell> for each cell B{i} = B1i*B2i*32 single

Usually A is larger than B (as A is the feature map and B is the filter).

C <1x2 cell> for i-th cell is the response for the i-th filter.

C{i} = (A1-B1i+1)*(A2-B2i+1)*32 double

*/

struct thread_data {
  float *A;
  float *B;
  double *C;
  mxArray *mxC;
  const mwSize *A_dims;
  const mwSize *B_dims;
  mwSize C_dims[2];
};

// convolve A and B (A is feature map, B is filter)
void process(void *thread_arg) {
  thread_data *args = (thread_data *)thread_arg;
  float *A = args->A;
  float *B = args->B;
  double *C = args->C;
  const mwSize *A_dims = args->A_dims;
  const mwSize *B_dims = args->B_dims;
  const mwSize *C_dims = args->C_dims;
  int num_features = args->A_dims[2];

  // for each feature dimension (e.g., 32 dim HOG)
  for (int f = 0; f < num_features; f++) {
    double *dst = C;
    // calibrate the pointer to the base at the f-th dimension
    float *A_src = A + f*A_dims[0]*A_dims[1];      
    float *B_src = B + f*B_dims[0]*B_dims[1];
    // for each pixel (x,y) on the response map
    for (int x = 0; x < C_dims[1]; x++) {
      for (int y = 0; y < C_dims[0]; y++) {
        double val = 0;
        // operate along the second dimension (column) of B (filter)
        // xp is the pointer/offset for the second dimension of B
        for (int xp = 0; xp < B_dims[1]; xp++) {
          // calibrate the pointer to the current column
          // for A it's (x+xp)*(number of elements per column) + (row offset)
          float *A_off = A_src + (x+xp)*A_dims[0] + y;
          // for B it's xp*(number of elements per column)
          float *B_off = B_src + xp*B_dims[0];
          // depending on the number of elements per row
          // if number of elements per row <= 20 then use the loop-less code below
          // the more general default case which uses for-loop follows
          // as for most filters number of elements per row <= 20
          // this would acclerate the process.
          switch(B_dims[0]) {
            case 20: val += A_off[19] * B_off[19];
            case 19: val += A_off[18] * B_off[18];
            case 18: val += A_off[17] * B_off[17];
            case 17: val += A_off[16] * B_off[16];
            case 16: val += A_off[15] * B_off[15];
            case 15: val += A_off[14] * B_off[14];
            case 14: val += A_off[13] * B_off[13];
            case 13: val += A_off[12] * B_off[12];
            case 12: val += A_off[11] * B_off[11];
            case 11: val += A_off[10] * B_off[10];
            case 10: val += A_off[9] * B_off[9];
            case 9: val += A_off[8] * B_off[8];
            case 8: val += A_off[7] * B_off[7];
            case 7: val += A_off[6] * B_off[6];
            case 6: val += A_off[5] * B_off[5];
            case 5: val += A_off[4] * B_off[4];
            case 4: val += A_off[3] * B_off[3];
            case 3: val += A_off[2] * B_off[2];
            case 2: val += A_off[1] * B_off[1];
            case 1: val += A_off[0] * B_off[0];
              break;
            default:            
              for (int yp = 0; yp < B_dims[0]; yp++) {
                val += *(A_off++) * *(B_off++);
              }
         }
       }
       *(dst++) += val;
      }
    }
  }
}

// matlab entry point
// C = fconv(A, cell of B, start, end);
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
  if (nrhs != 4)
    mexErrMsgTxt("Wrong number of inputs"); 
  if (nlhs != 1)
    mexErrMsgTxt("Wrong number of outputs");

  // get A
  const mxArray *mxA = prhs[0];
  if (mxGetNumberOfDimensions(mxA) != 3 || 
      mxGetClassID(mxA) != mxSINGLE_CLASS)
    mexErrMsgTxt("Invalid input: A");

  // get B and start/end
  const mxArray *cellB = prhs[1];
  mwSize num_bs = mxGetNumberOfElements(cellB);  
  int start = (int)mxGetScalar(prhs[2]) - 1;
  int end = (int)mxGetScalar(prhs[3]) - 1;
  if (start < 0 || end >= num_bs || start > end)
    mexErrMsgTxt("Invalid input: start/end");
  int len = end-start+1;

  // output cell
  plhs[0] = mxCreateCellMatrix(1, len);

  // do convolutions
  thread_data td;
  const mwSize *A_dims = mxGetDimensions(mxA);
  float *A = (float *)mxGetPr(mxA);
  for (int i = 0; i < len; i++) {
    const mxArray *mxB = mxGetCell(cellB, i+start);
    td.A_dims = A_dims;
    td.A = A;
    td.B_dims = mxGetDimensions(mxB);
    td.B = (float *)mxGetPr(mxB);
    if (mxGetNumberOfDimensions(mxB) != 3 ||
        mxGetClassID(mxB) != mxSINGLE_CLASS ||
        td.A_dims[2] != td.B_dims[2])
      mexErrMsgTxt("Invalid input: B");

    // compute size of output
    int height = td.A_dims[0] - td.B_dims[0] + 1;
    int width = td.A_dims[1] - td.B_dims[1] + 1;
    if (height < 1 || width < 1)
      mexErrMsgTxt("Invalid input: B should be smaller than A");
    td.C_dims[0] = height;
    td.C_dims[1] = width;
    td.mxC = mxCreateNumericArray(2, td.C_dims, mxDOUBLE_CLASS, mxREAL);
    td.C = (double *)mxGetPr(td.mxC);
    process((void *)&td);
    mxSetCell(plhs[0], i, td.mxC);
  }
}

Improved HOG feature in Deformable Parts Model (CVPR 2008)

The following code is an excerpt from features.cc that comes with the following paper:

A Discriminatively Trained, Multiscale, Deformable Part Model.
P. Felzenszwalb, D. McAllester, D. Ramaman.
Proceedings of the IEEE CVPR 2008.

The original code was downloaded from: http://people.cs.uchicago.edu/~pff/latent/

I only added some comments to make it easier to read. This improved HOG feature performs a little bit better than the original HOG implementation in OpenCV. It is easy to understand, tweak, and very popular in the community.

#include <math.h>
#include "mex.h"

// small value, used to avoid division by zero
#define eps 0.0001

// unit vectors used to compute gradient orientation
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};

static inline double min(double x, double y) { return (x <= y ? x : y); }
static inline double max(double x, double y) { return (x <= y ? y : x); }

static inline int min(int x, int y) { return (x <= y ? x : y); }
static inline int max(int x, int y) { return (x <= y ? y : x); }

// main function:
// takes a double color image and a bin size 
// returns HOG features
mxArray *process(const mxArray *mximage, const mxArray *mxsbin) {
  double *im = (double *)mxGetPr(mximage);
  const int *dims = mxGetDimensions(mximage);
  if (mxGetNumberOfDimensions(mximage) != 3 ||
      dims[2] != 3 || // must be a 3-channel image
      mxGetClassID(mximage) != mxDOUBLE_CLASS)
    mexErrMsgTxt("Invalid input");

  int sbin = (int)mxGetScalar(mxsbin); // pixel size for the HOG cells

  // memory for caching orientation histograms & their norms
  int blocks[2];
  blocks[0] = (int)round((double)dims[0]/(double)sbin); // how many cells along the cloumn
  blocks[1] = (int)round((double)dims[1]/(double)sbin); // how many cells along the row
  double *hist = (double *)mxCalloc(blocks[0]*blocks[1]*18, sizeof(double)); // level histogram into 18 orientation bins
  double *norm = (double *)mxCalloc(blocks[0]*blocks[1], sizeof(double)); // store the norm for gradient in each cell

  // memory for HOG features
  int out[3];
  out[0] = max(blocks[0]-2, 0);
  out[1] = max(blocks[1]-2, 0);
  out[2] = 27+4+1; 
  mxArray *mxfeat = mxCreateNumericArray(3, out, mxDOUBLE_CLASS, mxREAL);
  double *feat = (double *)mxGetPr(mxfeat);
  
  int visible[2];
  visible[0] = blocks[0]*sbin; // number of pixels along cloumn
  visible[1] = blocks[1]*sbin; // number of pixels along row
  
  for (int x = 1; x < visible[1]-1; x++) { // NOTE: in this loop row comes first, and x and y are iterated from 1 to visible[]-1,
    for (int y = 1; y < visible[0]-1; y++) { // in this way both x and y cooridnates leave out 1 pixel at the edges

      // first color channel
      double *s = im + min(x, dims[1]-2)*dims[0] + min(y, dims[0]-2); // NOTE: matlab matrices are stored in cloumns, to get access to rows 
      // we need to multiply by dims[0] (number of cells along the cloumn)
      double dy = *(s+1) - *(s-1);
      double dx = *(s+dims[0]) - *(s-dims[0]);
      double v = dx*dx + dy*dy;

      // second color channel
      s += dims[0]*dims[1]; // moves to the next channel
      double dy2 = *(s+1) - *(s-1);
      double dx2 = *(s+dims[0]) - *(s-dims[0]);
      double v2 = dx2*dx2 + dy2*dy2;

      // third color channel
      s += dims[0]*dims[1]; // moves to the next channel
      double dy3 = *(s+1) - *(s-1);
      double dx3 = *(s+dims[0]) - *(s-dims[0]);
      double v3 = dx3*dx3 + dy3*dy3;

      // pick channel with strongest gradient
      if (v2 > v) {
        v = v2;
        dx = dx2;
        dy = dy2;
      } 
      if (v3 > v) {
        v = v3;
        dx = dx3;
        dy = dy3;
      }

      // snap to one of 18 orientations
      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;
        }
      }
      
      // add to 4 histograms around pixel using linear interpolation
      double xp = ((double)x+0.5)/(double)sbin - 0.5;
      double yp = ((double)y+0.5)/(double)sbin - 0.5;
      int ixp = (int)floor(xp);
      int iyp = (int)floor(yp);
      double vx0 = xp-ixp;
      double vy0 = yp-iyp;
      double vx1 = 1.0-vx0;
      double vy1 = 1.0-vy0;
      v = sqrt(v);

      if (ixp >= 0 && iyp >= 0) {
        *(hist + ixp*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
         vx1*vy1*v;
      }

      if (ixp+1 < blocks[1] && iyp >= 0) {
        *(hist + (ixp+1)*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
         vx0*vy1*v;
      }

      if (ixp >= 0 && iyp+1 < blocks[0]) {
        *(hist + ixp*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
         vx1*vy0*v;
      }

      if (ixp+1 < blocks[1] && iyp+1 < blocks[0]) {
        *(hist + (ixp+1)*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
         vx0*vy0*v;
      }
    }
  }

  // compute energy in each block by summing over orientations
  for (int o = 0; o < 9; o++) {
    double *src1 = hist + o*blocks[0]*blocks[1];
    double *src2 = hist + (o+9)*blocks[0]*blocks[1];
    double *dst = norm;
    double *end = norm + blocks[1]*blocks[0];
    while (dst < end) {
      *(dst++) += (*src1 + *src2) * (*src1 + *src2);
      src1++;
      src2++;
    }
  }

  // compute features
  for (int x = 0; x < out[1]; x++) {
    for (int y = 0; y < out[0]; y++) {
      double *dst = feat + x*out[0] + y;      
      double *src, *p, n1, n2, n3, n4;

      p = norm + (x+1)*blocks[0] + y+1;
      n1 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
      p = norm + (x+1)*blocks[0] + y;
      n2 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
      p = norm + x*blocks[0] + y+1;
      n3 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
      p = norm + x*blocks[0] + y;      
      n4 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);

      double t1 = 0;
      double t2 = 0;
      double t3 = 0;
      double t4 = 0;

      // contrast-sensitive features
      src = hist + (x+1)*blocks[0] + (y+1);
      for (int o = 0; o < 18; o++) {
        double h1 = min(*src * n1, 0.2);
        double h2 = min(*src * n2, 0.2);
        double h3 = min(*src * n3, 0.2);
        double h4 = min(*src * n4, 0.2);
        *dst = 0.5 * (h1 + h2 + h3 + h4);
        t1 += h1;
        t2 += h2;
        t3 += h3;
        t4 += h4;
        dst += out[0]*out[1];
        src += blocks[0]*blocks[1];
      }

      // contrast-insensitive features
      src = hist + (x+1)*blocks[0] + (y+1);
      for (int o = 0; o < 9; o++) {
        double sum = *src + *(src + 9*blocks[0]*blocks[1]);
        double h1 = min(sum * n1, 0.2);
        double h2 = min(sum * n2, 0.2);
        double h3 = min(sum * n3, 0.2);
        double h4 = min(sum * n4, 0.2);
        *dst = 0.5 * (h1 + h2 + h3 + h4);
        dst += out[0]*out[1];
        src += blocks[0]*blocks[1];
      }

      // texture features
      *dst = 0.2357 * t1;
      dst += out[0]*out[1];
      *dst = 0.2357 * t2;
      dst += out[0]*out[1];
      *dst = 0.2357 * t3;
      dst += out[0]*out[1];
      *dst = 0.2357 * t4;

      // truncation feature
      dst += out[0]*out[1];
      *dst = 0;
    }
  }

  mxFree(hist);
  mxFree(norm);
  return mxfeat;
}

// matlab entry point
// F = features(image, bin)
// image should be color with double values
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
  if (nrhs != 2)
    mexErrMsgTxt("Wrong number of inputs"); 
  if (nlhs != 1)
    mexErrMsgTxt("Wrong number of outputs");
  plhs[0] = process(prhs[0], prhs[1]);
}