|
| 1 | +#include <stdlib.h> |
| 2 | +#include <math.h> |
| 3 | +#include <mex.h> |
| 4 | +#include <vector> |
| 5 | + |
| 6 | +using namespace std; |
| 7 | + |
| 8 | +void HoG(double *pixels, double *params, int *img_size, double *dth_des, unsigned int grayscale){ |
| 9 | + |
| 10 | + const float pi = 3.1415926536; |
| 11 | + |
| 12 | + int nb_bins = (int) params[0]; |
| 13 | + double cwidth = params[1]; |
| 14 | + int block_size = (int) params[2]; |
| 15 | + int orient = (int) params[3]; |
| 16 | + double clip_val = params[4]; |
| 17 | + |
| 18 | + int img_width = img_size[1]; |
| 19 | + int img_height = img_size[0]; |
| 20 | + |
| 21 | + int hist1= 2+ceil(-0.5 + img_height/cwidth); |
| 22 | + int hist2= 2+ceil(-0.5 + img_width/cwidth); |
| 23 | + |
| 24 | + double bin_size = (1+(orient==1))*pi/nb_bins; |
| 25 | + |
| 26 | + float dx[3], dy[3], grad_or, grad_mag, temp_mag; |
| 27 | + float Xc, Yc, Oc, block_norm; |
| 28 | + int x1, x2, y1, y2, bin1, bin2; |
| 29 | + int des_indx = 0; |
| 30 | + |
| 31 | + vector<vector<vector<double> > > h(hist1, vector<vector<double> > (hist2, vector<double> (nb_bins, 0.0) ) ); |
| 32 | + vector<vector<vector<double> > > block(block_size, vector<vector<double> > (block_size, vector<double> (nb_bins, 0.0) ) ); |
| 33 | + |
| 34 | + //Calculate gradients (zero padding) |
| 35 | + |
| 36 | + for(unsigned int y=0; y<img_height; y++) { |
| 37 | + for(unsigned int x=0; x<img_width; x++) { |
| 38 | + if (grayscale == 1){ |
| 39 | + if(x==0) dx[0] = pixels[y +(x+1)*img_height]; |
| 40 | + else{ |
| 41 | + if (x==img_width-1) dx[0] = -pixels[y + (x-1)*img_height]; |
| 42 | + else dx[0] = pixels[y+(x+1)*img_height] - pixels[y + (x-1)*img_height]; |
| 43 | + } |
| 44 | + if(y==0) dy[0] = -pixels[y+1+x*img_height]; |
| 45 | + else{ |
| 46 | + if (y==img_height-1) dy[0] = pixels[y-1+x*img_height]; |
| 47 | + else dy[0] = -pixels[y+1+x*img_height] + pixels[y-1+x*img_height]; |
| 48 | + } |
| 49 | + } |
| 50 | + else{ |
| 51 | + if(x==0){ |
| 52 | + dx[0] = pixels[y +(x+1)*img_height]; |
| 53 | + dx[1] = pixels[y +(x+1)*img_height + img_height*img_width]; |
| 54 | + dx[2] = pixels[y +(x+1)*img_height + 2*img_height*img_width]; |
| 55 | + } |
| 56 | + else{ |
| 57 | + if (x==img_width-1){ |
| 58 | + dx[0] = -pixels[y + (x-1)*img_height]; |
| 59 | + dx[1] = -pixels[y + (x-1)*img_height + img_height*img_width]; |
| 60 | + dx[2] = -pixels[y + (x-1)*img_height + 2*img_height*img_width]; |
| 61 | + } |
| 62 | + else{ |
| 63 | + dx[0] = pixels[y+(x+1)*img_height] - pixels[y + (x-1)*img_height]; |
| 64 | + dx[1] = pixels[y+(x+1)*img_height + img_height*img_width] - pixels[y + (x-1)*img_height + img_height*img_width]; |
| 65 | + dx[2] = pixels[y+(x+1)*img_height + 2*img_height*img_width] - pixels[y + (x-1)*img_height + 2*img_height*img_width]; |
| 66 | + |
| 67 | + } |
| 68 | + } |
| 69 | + if(y==0){ |
| 70 | + dy[0] = -pixels[y+1+x*img_height]; |
| 71 | + dy[1] = -pixels[y+1+x*img_height + img_height*img_width]; |
| 72 | + dy[2] = -pixels[y+1+x*img_height + 2*img_height*img_width]; |
| 73 | + } |
| 74 | + else{ |
| 75 | + if (y==img_height-1){ |
| 76 | + dy[0] = pixels[y-1+x*img_height]; |
| 77 | + dy[1] = pixels[y-1+x*img_height + img_height*img_width]; |
| 78 | + dy[2] = pixels[y-1+x*img_height + 2*img_height*img_width]; |
| 79 | + } |
| 80 | + else{ |
| 81 | + dy[0] = -pixels[y+1+x*img_height] + pixels[y-1+x*img_height]; |
| 82 | + dy[1] = -pixels[y+1+x*img_height + img_height*img_width] + pixels[y-1+x*img_height + img_height*img_width]; |
| 83 | + dy[2] = -pixels[y+1+x*img_height + 2*img_height*img_width] + pixels[y-1+x*img_height + 2*img_height*img_width]; |
| 84 | + } |
| 85 | + } |
| 86 | + } |
| 87 | + |
| 88 | + grad_mag = sqrt(dx[0]*dx[0] + dy[0]*dy[0]); |
| 89 | + grad_or= atan2(dy[0], dx[0]); |
| 90 | + |
| 91 | + if (grayscale == 0){ |
| 92 | + temp_mag = grad_mag; |
| 93 | + for (unsigned int cli=1;cli<3;++cli){ |
| 94 | + temp_mag= sqrt(dx[cli]*dx[cli] + dy[cli]*dy[cli]); |
| 95 | + if (temp_mag>grad_mag){ |
| 96 | + grad_mag=temp_mag; |
| 97 | + grad_or= atan2(dy[cli], dx[cli]); |
| 98 | + } |
| 99 | + } |
| 100 | + } |
| 101 | + |
| 102 | + if (grad_or<0) grad_or+=pi + (orient==1) * pi; |
| 103 | + |
| 104 | + // trilinear interpolation |
| 105 | + |
| 106 | + bin1 = (int)floor(0.5 + grad_or/bin_size) - 1; |
| 107 | + bin2 = bin1 + 1; |
| 108 | + x1 = (int)floor(0.5+ x/cwidth); |
| 109 | + x2 = x1+1; |
| 110 | + y1 = (int)floor(0.5+ y/cwidth); |
| 111 | + y2 = y1 + 1; |
| 112 | + |
| 113 | + Xc = (x1+1-1.5)*cwidth + 0.5; |
| 114 | + Yc = (y1+1-1.5)*cwidth + 0.5; |
| 115 | + |
| 116 | + Oc = (bin1+1+1-1.5)*bin_size; |
| 117 | + |
| 118 | + if (bin2==nb_bins){ |
| 119 | + bin2=0; |
| 120 | + } |
| 121 | + if (bin1<0){ |
| 122 | + bin1=nb_bins-1; |
| 123 | + } |
| 124 | + |
| 125 | + h[y1][x1][bin1]= h[y1][x1][bin1]+grad_mag*(1-((x+1-Xc)/cwidth))*(1-((y+1-Yc)/cwidth))*(1-((grad_or-Oc)/bin_size)); |
| 126 | + h[y1][x1][bin2]= h[y1][x1][bin2]+grad_mag*(1-((x+1-Xc)/cwidth))*(1-((y+1-Yc)/cwidth))*(((grad_or-Oc)/bin_size)); |
| 127 | + h[y2][x1][bin1]= h[y2][x1][bin1]+grad_mag*(1-((x+1-Xc)/cwidth))*(((y+1-Yc)/cwidth))*(1-((grad_or-Oc)/bin_size)); |
| 128 | + h[y2][x1][bin2]= h[y2][x1][bin2]+grad_mag*(1-((x+1-Xc)/cwidth))*(((y+1-Yc)/cwidth))*(((grad_or-Oc)/bin_size)); |
| 129 | + h[y1][x2][bin1]= h[y1][x2][bin1]+grad_mag*(((x+1-Xc)/cwidth))*(1-((y+1-Yc)/cwidth))*(1-((grad_or-Oc)/bin_size)); |
| 130 | + h[y1][x2][bin2]= h[y1][x2][bin2]+grad_mag*(((x+1-Xc)/cwidth))*(1-((y+1-Yc)/cwidth))*(((grad_or-Oc)/bin_size)); |
| 131 | + h[y2][x2][bin1]= h[y2][x2][bin1]+grad_mag*(((x+1-Xc)/cwidth))*(((y+1-Yc)/cwidth))*(1-((grad_or-Oc)/bin_size)); |
| 132 | + h[y2][x2][bin2]= h[y2][x2][bin2]+grad_mag*(((x+1-Xc)/cwidth))*(((y+1-Yc)/cwidth))*(((grad_or-Oc)/bin_size)); |
| 133 | + } |
| 134 | + } |
| 135 | + |
| 136 | + |
| 137 | + |
| 138 | + //Block normalization |
| 139 | + |
| 140 | + for(unsigned int x=1; x<hist2-block_size; x++){ |
| 141 | + for (unsigned int y=1; y<hist1-block_size; y++){ |
| 142 | + |
| 143 | + block_norm=0; |
| 144 | + for (unsigned int i=0; i<block_size; i++){ |
| 145 | + for(unsigned int j=0; j<block_size; j++){ |
| 146 | + for(unsigned int k=0; k<nb_bins; k++){ |
| 147 | + block_norm+=h[y+i][x+j][k]*h[y+i][x+j][k]; |
| 148 | + } |
| 149 | + } |
| 150 | + } |
| 151 | + |
| 152 | + block_norm=sqrt(block_norm); |
| 153 | + for (unsigned int i1=0; i1<block_size; i1++){ |
| 154 | + for(unsigned int j1=0; j1<block_size; j1++){ |
| 155 | + for(unsigned int k1=0; k1<nb_bins; k1++){ |
| 156 | + if (block_norm>0){ |
| 157 | + block[i1][j1][k1]=h[y+i1][x+j1][k1]/block_norm; |
| 158 | + if (block[i1][j1][k1]>clip_val) block[i1][j1][k1]=clip_val; |
| 159 | + } |
| 160 | + } |
| 161 | + } |
| 162 | + } |
| 163 | + |
| 164 | + block_norm=0; |
| 165 | + for (unsigned int i2=0; i2<block_size; i2++){ |
| 166 | + for(unsigned int j2=0; j2<block_size; j2++){ |
| 167 | + for(unsigned int k2=0; k2<nb_bins; k2++){ |
| 168 | + block_norm+=block[i2][j2][k2]*block[i2][j2][k2]; |
| 169 | + } |
| 170 | + } |
| 171 | + } |
| 172 | + |
| 173 | + block_norm=sqrt(block_norm); |
| 174 | + for (unsigned int i3=0; i3<block_size; i3++){ |
| 175 | + for(unsigned int j3=0; j3<block_size; j3++){ |
| 176 | + for(unsigned int k3=0; k3<nb_bins; k3++){ |
| 177 | + if (block_norm>0) dth_des[des_indx]=block[i3][j3][k3]/block_norm; |
| 178 | + else dth_des[des_indx]=0.0; |
| 179 | + des_indx++; |
| 180 | + } |
| 181 | + } |
| 182 | + } |
| 183 | + } |
| 184 | + } |
| 185 | +} |
| 186 | + |
| 187 | + |
| 188 | +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { |
| 189 | + |
| 190 | + double *pixels, *dth_des, *params; |
| 191 | + int nb_bins, block_size; |
| 192 | + int img_size[2]; |
| 193 | + unsigned int grayscale = 1; |
| 194 | + |
| 195 | + if (nlhs>1) mexErrMsgTxt("Too many output arguments"); |
| 196 | + if (nrhs==0) mexErrMsgTxt("No Image -> No HoG"); |
| 197 | + |
| 198 | + pixels = (double*) mxGetPr(prhs[0]); |
| 199 | + img_size[0] = mxGetM(prhs[0]); |
| 200 | + img_size[1] = mxGetN(prhs[0]); |
| 201 | + if (mxGetNumberOfDimensions(prhs[0])==3){ |
| 202 | + img_size[1] /= 3; |
| 203 | + grayscale = 0; |
| 204 | + } |
| 205 | + |
| 206 | + if (nrhs>1){ |
| 207 | + params = mxGetPr(prhs[1]); |
| 208 | + if (params[0]<=0) mexErrMsgTxt("Number of orientation bins must be positive"); |
| 209 | + if (params[1]<=0) mexErrMsgTxt("Cell size must be positive"); |
| 210 | + if (params[2]<=0) mexErrMsgTxt("Block size must be positive"); |
| 211 | + } |
| 212 | + else { |
| 213 | + params = new double[5]; |
| 214 | + params[0]=9; |
| 215 | + params[1]=8; |
| 216 | + params[2]=2; |
| 217 | + params[3]=0; |
| 218 | + params[4]=0.2; |
| 219 | + } |
| 220 | + |
| 221 | + nb_bins = (int) params[0]; |
| 222 | + block_size = (int) params[2]; |
| 223 | + |
| 224 | + int hist1= 2+ceil(-0.5 + img_size[0]/params[1]); |
| 225 | + int hist2= 2+ceil(-0.5 + img_size[1]/params[1]); |
| 226 | + |
| 227 | + plhs[0] = mxCreateDoubleMatrix((hist1-2-(block_size-1))*(hist2-2-(block_size-1))*nb_bins*block_size*block_size, 1, mxREAL); |
| 228 | + dth_des = mxGetPr(plhs[0]); |
| 229 | + |
| 230 | + HoG(pixels, params, img_size, dth_des, grayscale); |
| 231 | + if (nrhs==1) delete[] params; |
| 232 | +} |
0 commit comments