modules/objdetect/src/featurepyramid.cpp

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