modules/objdetect/src/featurepyramid.cpp

   1 #include "precomp.hpp"
   2 #include "opencv2/imgproc/imgproc_c.h"
   3 #include "_latentsvm.h"
   4 #include "_lsvm_resizeimg.h"
   5
   6 #ifndef max
   7 #define max(a,b)            (((a) > (b)) ? (a) : (b))
   8 #endif
   9
  10 #ifndef min
  11 #define min(a,b)            (((a) < (b)) ? (a) : (b))
  12 #endif
  13
  14 /*
  15 // Getting feature map for the selected subimage
  16 //
  17 // API
  18 // int getFeatureMaps(const IplImage * image, const int k, featureMap **map);
  19 // INPUT
  20 // image             - selected subimage
  21 // k                 - size of cells
  22 // OUTPUT
  23 // map               - feature map
  24 // RESULT
  25 // Error status
  26 */
  27 int getFeatureMaps(const IplImage* image, const int k, CvLSVMFeatureMap **map)
  28 {
  29     int sizeX, sizeY;
  30     int p, px, stringSize;
  31     int height, width, numChannels;
  32     int i, j, kk, c, ii, jj, d;
  33     float  * datadx, * datady;
  34
  35     int   ch;
  36     float magnitude, x, y, tx, ty;
  37
  38     IplImage * dx, * dy;
  39     int *nearest;
  40     float *w, a_x, b_x;
  41
  42     float kernel[3] = {-1.f, 0.f, 1.f};
  43     CvMat kernel_dx = cvMat(1, 3, CV_32F, kernel);
  44     CvMat kernel_dy = cvMat(3, 1, CV_32F, kernel);
  45
  46     float * r;
  47     int   * alfa;
  48
  49     float boundary_x[NUM_SECTOR + 1];
  50     float boundary_y[NUM_SECTOR + 1];
  51     float max, dotProd;
  52     int   maxi;
  53
  54     height = image->height;
  55     width  = image->width ;
  56
  57     numChannels = image->nChannels;
  58
  59     dx    = cvCreateImage(cvSize(image->width, image->height),
  60                           IPL_DEPTH_32F, 3);
  61     dy    = cvCreateImage(cvSize(image->width, image->height),
  62                           IPL_DEPTH_32F, 3);
  63
  64     sizeX = width  / k;
  65     sizeY = height / k;
  66     px    = 3 * NUM_SECTOR;
  67     p     = px;
  68     stringSize = sizeX * p;
  69     allocFeatureMapObject(map, sizeX, sizeY, p);
  70
  71     cvFilter2D(image, dx, &kernel_dx, cvPoint(-1, 0));
  72     cvFilter2D(image, dy, &kernel_dy, cvPoint(0, -1));
  73
  74     float arg_vector;
  75     for(i = 0; i <= NUM_SECTOR; i++)
  76     {
  77         arg_vector    = ( (float) i ) * ( (float)(PI) / (float)(NUM_SECTOR) );
  78         boundary_x[i] = cosf(arg_vector);
  79         boundary_y[i] = sinf(arg_vector);
  80     }/*for(i = 0; i <= NUM_SECTOR; i++) */
  81
  82     r    = (float *)malloc( sizeof(float) * (width * height));
  83     alfa = (int   *)malloc( sizeof(int  ) * (width * height * 2));
  84
  85     for(j = 1; j < height - 1; j++)
  86     {
  87         datadx = (float*)(dx->imageData + dx->widthStep * j);
  88         datady = (float*)(dy->imageData + dy->widthStep * j);
  89         for(i = 1; i < width - 1; i++)
  90         {
  91             c = 0;
  92             x = (datadx[i * numChannels + c]);
  93             y = (datady[i * numChannels + c]);
  94
  95             r[j * width + i] =sqrtf(x * x + y * y);
  96             for(ch = 1; ch < numChannels; ch++)
  97             {
  98                 tx = (datadx[i * numChannels + ch]);
  99                 ty = (datady[i * numChannels + ch]);
 100                 magnitude = sqrtf(tx * tx + ty * ty);
 101                 if(magnitude > r[j * width + i])
 102                 {
 103                     r[j * width + i] = magnitude;
 104                     c = ch;
 105                     x = tx;
 106                     y = ty;
 107                 }
 108             }/*for(ch = 1; ch < numChannels; ch++)*/
 109
 110             max  = boundary_x[0] * x + boundary_y[0] * y;
 111             maxi = 0;
 112             for (kk = 0; kk < NUM_SECTOR; kk++)
 113             {
 114                 dotProd = boundary_x[kk] * x + boundary_y[kk] * y;
 115                 if (dotProd > max)
 116                 {
 117                     max  = dotProd;
 118                     maxi = kk;
 119                 }
 120                 else
 121                 {
 122                     if (-dotProd > max)
 123                     {
 124                         max  = -dotProd;
 125                         maxi = kk + NUM_SECTOR;
 126                     }
 127                 }
 128             }
 129             alfa[j * width * 2 + i * 2    ] = maxi % NUM_SECTOR;
 130             alfa[j * width * 2 + i * 2 + 1] = maxi;
 131         }/*for(i = 0; i < width; i++)*/
 132     }/*for(j = 0; j < height; j++)*/
 133
 134     nearest = (int  *)malloc(sizeof(int  ) *  k);
 135     w       = (float*)malloc(sizeof(float) * (k * 2));
 136
 137     for(i = 0; i < k / 2; i++)
 138     {
 139         nearest[i] = -1;
 140     }/*for(i = 0; i < k / 2; i++)*/
 141     for(i = k / 2; i < k; i++)
 142     {
 143         nearest[i] = 1;
 144     }/*for(i = k / 2; i < k; i++)*/
 145
 146     for(j = 0; j < k / 2; j++)
 147     {
 148         b_x = k / 2 + j + 0.5f;
 149         a_x = k / 2 - j - 0.5f;
 150         w[j * 2    ] = 1.0f/a_x * ((a_x * b_x) / ( a_x + b_x));
 151         w[j * 2 + 1] = 1.0f/b_x * ((a_x * b_x) / ( a_x + b_x));
 152     }/*for(j = 0; j < k / 2; j++)*/
 153     for(j = k / 2; j < k; j++)
 154     {
 155         a_x = j - k / 2 + 0.5f;
 156         b_x =-j + k / 2 - 0.5f + k;
 157         w[j * 2    ] = 1.0f/a_x * ((a_x * b_x) / ( a_x + b_x));
 158         w[j * 2 + 1] = 1.0f/b_x * ((a_x * b_x) / ( a_x + b_x));
 159     }/*for(j = k / 2; j < k; j++)*/
 160
 161
 162     for(i = 0; i < sizeY; i++)
 163     {
 164       for(j = 0; j < sizeX; j++)
 165       {
 166         for(ii = 0; ii < k; ii++)
 167         {
 168           for(jj = 0; jj < k; jj++)
 169           {
 170             if ((i * k + ii > 0) &&
 171                 (i * k + ii < height - 1) &&
 172                 (j * k + jj > 0) &&
 173                 (j * k + jj < width  - 1))
 174             {
 175               d = (k * i + ii) * width + (j * k + jj);
 176               (*map)->map[ i * stringSize + j * (*map)->numFeatures + alfa[d * 2    ]] +=
 177                   r[d] * w[ii * 2] * w[jj * 2];
 178               (*map)->map[ i * stringSize + j * (*map)->numFeatures + alfa[d * 2 + 1] + NUM_SECTOR] +=
 179                   r[d] * w[ii * 2] * w[jj * 2];
 180               if ((i + nearest[ii] >= 0) &&
 181                   (i + nearest[ii] <= sizeY - 1))
 182               {
 183                 (*map)->map[(i + nearest[ii]) * stringSize + j * (*map)->numFeatures + alfa[d * 2    ]             ] +=
 184                   r[d] * w[ii * 2 + 1] * w[jj * 2 ];
 185                 (*map)->map[(i + nearest[ii]) * stringSize + j * (*map)->numFeatures + alfa[d * 2 + 1] + NUM_SECTOR] +=
 186                   r[d] * w[ii * 2 + 1] * w[jj * 2 ];
 187               }
 188               if ((j + nearest[jj] >= 0) &&
 189                   (j + nearest[jj] <= sizeX - 1))
 190               {
 191                 (*map)->map[i * stringSize + (j + nearest[jj]) * (*map)->numFeatures + alfa[d * 2    ]             ] +=
 192                   r[d] * w[ii * 2] * w[jj * 2 + 1];
 193                 (*map)->map[i * stringSize + (j + nearest[jj]) * (*map)->numFeatures + alfa[d * 2 + 1] + NUM_SECTOR] +=
 194                   r[d] * w[ii * 2] * w[jj * 2 + 1];
 195               }
 196               if ((i + nearest[ii] >= 0) &&
 197                   (i + nearest[ii] <= sizeY - 1) &&
 198                   (j + nearest[jj] >= 0) &&
 199                   (j + nearest[jj] <= sizeX - 1))
 200               {
 201                 (*map)->map[(i + nearest[ii]) * stringSize + (j + nearest[jj]) * (*map)->numFeatures + alfa[d * 2    ]             ] +=
 202                   r[d] * w[ii * 2 + 1] * w[jj * 2 + 1];
 203                 (*map)->map[(i + nearest[ii]) * stringSize + (j + nearest[jj]) * (*map)->numFeatures + alfa[d * 2 + 1] + NUM_SECTOR] +=
 204                   r[d] * w[ii * 2 + 1] * w[jj * 2 + 1];
 205               }
 206             }
 207           }/*for(jj = 0; jj < k; jj++)*/
 208         }/*for(ii = 0; ii < k; ii++)*/
 209       }/*for(j = 1; j < sizeX - 1; j++)*/
 210     }/*for(i = 1; i < sizeY - 1; i++)*/
 211
 212     cvReleaseImage(&dx);
 213     cvReleaseImage(&dy);
 214
 215
 216     free(w);
 217     free(nearest);
 218
 219     free(r);
 220     free(alfa);
 221
 222     return LATENT_SVM_OK;
 223 }
 224
 225 /*
 226 // Feature map Normalization and Truncation
 227 //
 228 // API
 229 // int normalizeAndTruncate(featureMap *map, const float alfa);
 230 // INPUT
 231 // map               - feature map
 232 // alfa              - truncation threshold
 233 // OUTPUT
 234 // map               - truncated and normalized feature map
 235 // RESULT
 236 // Error status
 237 */
 238 int normalizeAndTruncate(CvLSVMFeatureMap *map, const float alfa)
 239 {
 240     int i,j, ii;
 241     int sizeX, sizeY, p, pos, pp, xp, pos1, pos2;
 242     float * partOfNorm; // norm of C(i, j)
 243     float * newData;
 244     float   valOfNorm;
 245
 246     sizeX     = map->sizeX;
 247     sizeY     = map->sizeY;
 248     partOfNorm = (float *)malloc (sizeof(float) * (sizeX * sizeY));
 249
 250     p  = NUM_SECTOR;
 251     xp = NUM_SECTOR * 3;
 252     pp = NUM_SECTOR * 12;
 253
 254     for(i = 0; i < sizeX * sizeY; i++)
 255     {
 256         valOfNorm = 0.0f;
 257         pos = i * map->numFeatures;
 258         for(j = 0; j < p; j++)
 259         {
 260             valOfNorm += map->map[pos + j] * map->map[pos + j];
 261         }/*for(j = 0; j < p; j++)*/
 262         partOfNorm[i] = valOfNorm;
 263     }/*for(i = 0; i < sizeX * sizeY; i++)*/
 264
 265     sizeX -= 2;
 266     sizeY -= 2;
 267
 268     newData = (float *)malloc (sizeof(float) * (sizeX * sizeY * pp));
 269     //normalization
 270     for(i = 1; i <= sizeY; i++)
 271     {
 272         for(j = 1; j <= sizeX; j++)
 273         {
 274             valOfNorm = sqrtf(
 275                 partOfNorm[(i    )*(sizeX + 2) + (j    )] +
 276                 partOfNorm[(i    )*(sizeX + 2) + (j + 1)] +
 277                 partOfNorm[(i + 1)*(sizeX + 2) + (j    )] +
 278                 partOfNorm[(i + 1)*(sizeX + 2) + (j + 1)]) + FLT_EPSILON;
 279             pos1 = (i  ) * (sizeX + 2) * xp + (j  ) * xp;
 280             pos2 = (i-1) * (sizeX    ) * pp + (j-1) * pp;
 281             for(ii = 0; ii < p; ii++)
 282             {
 283                 newData[pos2 + ii        ] = map->map[pos1 + ii    ] / valOfNorm;
 284             }/*for(ii = 0; ii < p; ii++)*/
 285             for(ii = 0; ii < 2 * p; ii++)
 286             {
 287                 newData[pos2 + ii + p * 4] = map->map[pos1 + ii + p] / valOfNorm;
 288             }/*for(ii = 0; ii < 2 * p; ii++)*/
 289             valOfNorm = sqrtf(
 290                 partOfNorm[(i    )*(sizeX + 2) + (j    )] +
 291                 partOfNorm[(i    )*(sizeX + 2) + (j + 1)] +
 292                 partOfNorm[(i - 1)*(sizeX + 2) + (j    )] +
 293                 partOfNorm[(i - 1)*(sizeX + 2) + (j + 1)]) + FLT_EPSILON;
 294             for(ii = 0; ii < p; ii++)
 295             {
 296                 newData[pos2 + ii + p    ] = map->map[pos1 + ii    ] / valOfNorm;
 297             }/*for(ii = 0; ii < p; ii++)*/
 298             for(ii = 0; ii < 2 * p; ii++)
 299             {
 300                 newData[pos2 + ii + p * 6] = map->map[pos1 + ii + p] / valOfNorm;
 301             }/*for(ii = 0; ii < 2 * p; ii++)*/
 302             valOfNorm = sqrtf(
 303                 partOfNorm[(i    )*(sizeX + 2) + (j    )] +
 304                 partOfNorm[(i    )*(sizeX + 2) + (j - 1)] +
 305                 partOfNorm[(i + 1)*(sizeX + 2) + (j    )] +
 306                 partOfNorm[(i + 1)*(sizeX + 2) + (j - 1)]) + FLT_EPSILON;
 307             for(ii = 0; ii < p; ii++)
 308             {
 309                 newData[pos2 + ii + p * 2] = map->map[pos1 + ii    ] / valOfNorm;
 310             }/*for(ii = 0; ii < p; ii++)*/
 311             for(ii = 0; ii < 2 * p; ii++)
 312             {
 313                 newData[pos2 + ii + p * 8] = map->map[pos1 + ii + p] / valOfNorm;
 314             }/*for(ii = 0; ii < 2 * p; ii++)*/
 315             valOfNorm = sqrtf(
 316                 partOfNorm[(i    )*(sizeX + 2) + (j    )] +
 317                 partOfNorm[(i    )*(sizeX + 2) + (j - 1)] +
 318                 partOfNorm[(i - 1)*(sizeX + 2) + (j    )] +
 319                 partOfNorm[(i - 1)*(sizeX + 2) + (j - 1)]) + FLT_EPSILON;
 320             for(ii = 0; ii < p; ii++)
 321             {
 322                 newData[pos2 + ii + p * 3 ] = map->map[pos1 + ii    ] / valOfNorm;
 323             }/*for(ii = 0; ii < p; ii++)*/
 324             for(ii = 0; ii < 2 * p; ii++)
 325             {
 326                 newData[pos2 + ii + p * 10] = map->map[pos1 + ii + p] / valOfNorm;
 327             }/*for(ii = 0; ii < 2 * p; ii++)*/
 328         }/*for(j = 1; j <= sizeX; j++)*/
 329     }/*for(i = 1; i <= sizeY; i++)*/
 330     //truncation
 331     for(i = 0; i < sizeX * sizeY * pp; i++)
 332     {
 333         if(newData [i] > alfa) newData [i] = alfa;
 334     }/*for(i = 0; i < sizeX * sizeY * pp; i++)*/
 335     //swap data
 336
 337     map->numFeatures  = pp;
 338     map->sizeX = sizeX;
 339     map->sizeY = sizeY;
 340
 341     free (map->map);
 342     free (partOfNorm);
 343
 344     map->map = newData;
 345
 346     return LATENT_SVM_OK;
 347 }
 348
 349 /*
 350 // Feature map reduction
 351 // In each cell we reduce dimension of the feature vector
 352 // according to original paper special procedure
 353 //
 354 // API
 355 // int PCAFeatureMaps(featureMap *map)
 356 // INPUT
 357 // map               - feature map
 358 // OUTPUT
 359 // map               - feature map
 360 // RESULT
 361 // Error status
 362 */
 363 int PCAFeatureMaps(CvLSVMFeatureMap *map)
 364 {
 365     int i,j, ii, jj, k;
 366     int sizeX, sizeY, p,  pp, xp, yp, pos1, pos2;
 367     float * newData;
 368     float val;
 369     float nx, ny;
 370
 371     sizeX = map->sizeX;
 372     sizeY = map->sizeY;
 373     p     = map->numFeatures;
 374     pp    = NUM_SECTOR * 3 + 4;
 375     yp    = 4;
 376     xp    = NUM_SECTOR;
 377
 378     nx    = 1.0f / sqrtf((float)(xp * 2));
 379     ny    = 1.0f / sqrtf((float)(yp    ));
 380
 381     newData = (float *)malloc (sizeof(float) * (sizeX * sizeY * pp));
 382
 383     for(i = 0; i < sizeY; i++)
 384     {
 385         for(j = 0; j < sizeX; j++)
 386         {
 387             pos1 = ((i)*sizeX + j)*p;
 388             pos2 = ((i)*sizeX + j)*pp;
 389             k = 0;
 390             for(jj = 0; jj < xp * 2; jj++)
 391             {
 392                 val = 0;
 393                 for(ii = 0; ii < yp; ii++)
 394                 {
 395                     val += map->map[pos1 + yp * xp + ii * xp * 2 + jj];
 396                 }/*for(ii = 0; ii < yp; ii++)*/
 397                 newData[pos2 + k] = val * ny;
 398                 k++;
 399             }/*for(jj = 0; jj < xp * 2; jj++)*/
 400             for(jj = 0; jj < xp; jj++)
 401             {
 402                 val = 0;
 403                 for(ii = 0; ii < yp; ii++)
 404                 {
 405                     val += map->map[pos1 + ii * xp + jj];
 406                 }/*for(ii = 0; ii < yp; ii++)*/
 407                 newData[pos2 + k] = val * ny;
 408                 k++;
 409             }/*for(jj = 0; jj < xp; jj++)*/
 410             for(ii = 0; ii < yp; ii++)
 411             {
 412                 val = 0;
 413                 for(jj = 0; jj < 2 * xp; jj++)
 414                 {
 415                     val += map->map[pos1 + yp * xp + ii * xp * 2 + jj];
 416                 }/*for(jj = 0; jj < xp; jj++)*/
 417                 newData[pos2 + k] = val * nx;
 418                 k++;
 419             } /*for(ii = 0; ii < yp; ii++)*/
 420         }/*for(j = 0; j < sizeX; j++)*/
 421     }/*for(i = 0; i < sizeY; i++)*/
 422     //swap data
 423
 424     map->numFeatures = pp;
 425
 426     free (map->map);
 427
 428     map->map = newData;
 429
 430     return LATENT_SVM_OK;
 431 }
 432
 433
 434 static int getPathOfFeaturePyramid(IplImage * image,
 435                             float step, int numStep, int startIndex,
 436                             int sideLength, CvLSVMFeaturePyramid **maps)
 437 {
 438     CvLSVMFeatureMap *map;
 439     IplImage *scaleTmp;
 440     float scale;
 441     int   i;
 442
 443     for(i = 0; i < numStep; i++)
 444     {
 445         scale = 1.0f / powf(step, (float)i);
 446         scaleTmp = resize_opencv (image, scale);
 447         getFeatureMaps(scaleTmp, sideLength, &map);
 448         normalizeAndTruncate(map, VAL_OF_TRUNCATE);
 449         PCAFeatureMaps(map);
 450         (*maps)->pyramid[startIndex + i] = map;
 451         cvReleaseImage(&scaleTmp);
 452     }/*for(i = 0; i < numStep; i++)*/
 453     return LATENT_SVM_OK;
 454 }
 455
 456 /*
 457 // Getting feature pyramid
 458 //
 459 // API
 460 // int getFeaturePyramid(IplImage * image, const filterObject **all_F,
 461                       const int n_f,
 462                       const int lambda, const int k,
 463                       const int startX, const int startY,
 464                       const int W, const int H, featurePyramid **maps);
 465 // INPUT
 466 // image             - image
 467 // OUTPUT
 468 // maps              - feature maps for all levels
 469 // RESULT
 470 // Error status
 471 */
 472 int getFeaturePyramid(IplImage * image, CvLSVMFeaturePyramid **maps)
 473 {
 474     IplImage *imgResize;
 475     float step;
 476     int   numStep;
 477     int   maxNumCells;
 478     int   W, H;
 479
 480     if(image->depth == IPL_DEPTH_32F)
 481     {
 482         imgResize = image;
 483     }
 484     else
 485     {
 486         imgResize = cvCreateImage(cvSize(image->width , image->height) ,
 487                                   IPL_DEPTH_32F , 3);
 488         cvConvert(image, imgResize);
 489     }
 490
 491     W = imgResize->width;
 492     H = imgResize->height;
 493
 494     step = powf(2.0f, 1.0f / ((float)LAMBDA));
 495     maxNumCells = W / SIDE_LENGTH;
 496     if( maxNumCells > H / SIDE_LENGTH )
 497     {
 498         maxNumCells = H / SIDE_LENGTH;
 499     }
 500     numStep = (int)(logf((float) maxNumCells / (5.0f)) / logf( step )) + 1;
 501
 502     allocFeaturePyramidObject(maps, numStep + LAMBDA);
 503
 504     getPathOfFeaturePyramid(imgResize, step   , LAMBDA, 0,
 505                             SIDE_LENGTH / 2, maps);
 506     getPathOfFeaturePyramid(imgResize, step, numStep, LAMBDA,
 507                             SIDE_LENGTH    , maps);
 508
 509     if(image->depth != IPL_DEPTH_32F)
 510     {
 511         cvReleaseImage(&imgResize);
 512     }
 513
 514     return LATENT_SVM_OK;
 515 }