vcs.maemo.org Git - opencv/blob - src/cvaux/cvselfsimilarity.cpp

   1 // This is based on Rainer Lienhart contribution. Below is the original copyright:
   2 //
   3 /*M///////////////////////////////////////////////////////////////////////////////////////
   4 //
   5 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   6 //
   7 //  By downloading, copying, installing or using the software you agree to this license.
   8 //  If you do not agree to this license, do not download, install,
   9 //  copy or use the software.
  10 //
  11 //
  12 //                University of Augsburg License Agreement
  13 //                For Open Source MultiMedia Computing (MMC) Library
  14 //
  15 // Copyright (C) 2007, University of Augsburg, Germany, all rights reserved.
  16 // Third party copyrights are property of their respective owners.
  17 //
  18 // Redistribution and use in source and binary forms, with or without modification,
  19 // are permitted provided that the following conditions are met:
  20 //
  21 //   * Redistribution's of source code must retain the above copyright notice,
  22 //     this list of conditions and the following disclaimer.
  23 //
  24 //   * Redistribution's in binary form must reproduce the above copyright notice,
  25 //     this list of conditions and the following disclaimer in the documentation
  26 //     and/or other materials provided with the distribution.
  27 //
  28 //   * The name of University of Augsburg, Germany may not be used to endorse or promote products
  29 //     derived from this software without specific prior written permission.
  30 //
  31 // This software is provided by the copyright holders and contributors "as is" and
  32 // any express or implied warranties, including, but not limited to, the implied
  33 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  34 // In no event shall the University of Augsburg, Germany or contributors be liable for any direct,
  35 // indirect, incidental, special, exemplary, or consequential damages
  36 // (including, but not limited to, procurement of substitute goods or services;
  37 // loss of use, data, or profits; or business interruption) however caused
  38 // and on any theory of liability, whether in contract, strict liability,
  39 // or tort (including negligence or otherwise) arising in any way out of
  40 // the use of this software, even if advised of the possibility of such damage.
  41 //
  42 //M*/
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  45 //   Author:    Rainer Lienhart
  46 //              email: Rainer.Lienhart@informatik.uni-augsburg.de
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  48
  49 // Please cite the following two papers:
  50 // 1. Shechtman, E., Irani, M.:
  51 //              Matching local self-similarities across images and videos.
  52 //              CVPR, (2007)
  53 // 2. Eva Horster, Thomas Greif, Rainer Lienhart, Malcolm Slaney.
  54 //              Comparing Local Feature Descriptors in pLSA-Based Image Models.
  55 //              30th Annual Symposium of the German Association for
  56 //      Pattern Recognition (DAGM) 2008, Munich, Germany, June 2008.
  57
  58 #include "_cvaux.h"
  59
  60 namespace cv
  61 {
  62
  63 SelfSimDescriptor::SelfSimDescriptor()
  64 {
  65     smallSize = DEFAULT_SMALL_SIZE;
  66     largeSize = DEFAULT_LARGE_SIZE;
  67     numberOfAngles = DEFAULT_NUM_ANGLES;
  68     startDistanceBucket = DEFAULT_START_DISTANCE_BUCKET;
  69     numberOfDistanceBuckets = DEFAULT_NUM_DISTANCE_BUCKETS;
  70 }
  71
  72 SelfSimDescriptor::SelfSimDescriptor(int _ssize, int _lsize,
  73                                      int _startDistanceBucket,
  74                                      int _numberOfDistanceBuckets, int _numberOfAngles)
  75 {
  76     smallSize = _ssize;
  77     largeSize = _lsize;
  78     startDistanceBucket = _startDistanceBucket;
  79     numberOfDistanceBuckets = _numberOfDistanceBuckets;
  80     numberOfAngles = _numberOfAngles;
  81 }
  82
  83 SelfSimDescriptor::SelfSimDescriptor(const SelfSimDescriptor& ss)
  84 {
  85     smallSize = ss.smallSize;
  86     largeSize = ss.largeSize;
  87     startDistanceBucket = ss.startDistanceBucket;
  88     numberOfDistanceBuckets = ss.numberOfDistanceBuckets;
  89     numberOfAngles = ss.numberOfAngles;
  90 }
  91
  92 SelfSimDescriptor::~SelfSimDescriptor()
  93 {
  94 }
  95
  96 SelfSimDescriptor& SelfSimDescriptor::operator = (const SelfSimDescriptor& ss)
  97 {
  98     if( this != &ss )
  99     {
 100         smallSize = ss.smallSize;
 101         largeSize = ss.largeSize;
 102         startDistanceBucket = ss.startDistanceBucket;
 103         numberOfDistanceBuckets = ss.numberOfDistanceBuckets;
 104         numberOfAngles = ss.numberOfAngles;
 105     }
 106     return *this;
 107 }
 108
 109 size_t SelfSimDescriptor::getDescriptorSize() const
 110 {
 111     return numberOfAngles*(numberOfDistanceBuckets - startDistanceBucket);
 112 }
 113
 114 Size SelfSimDescriptor::getGridSize( Size imgSize, Size winStride ) const
 115 {
 116     winStride.width = std::max(winStride.width, 1);
 117     winStride.height = std::max(winStride.height, 1);
 118     int border = largeSize/2 + smallSize/2;
 119     return Size(std::max(imgSize.width - border*2 + winStride.width - 1, 0)/winStride.width,
 120                 std::max(imgSize.height - border*2 + winStride.height - 1, 0)/winStride.height);
 121 }
 122
 123 // TODO: optimized with SSE2
 124 void SelfSimDescriptor::SSD(const Mat& img, Point pt, Mat& ssd) const
 125 {
 126         int x, y, dx, dy, r0 = largeSize/2, r1 = smallSize/2;
 127     int step = img.step;
 128     for( y = -r0; y <= r0; y++ )
 129     {
 130         float* sptr = ssd.ptr<float>(y+r0) + r0;
 131         for( x = -r0; x <= r0; x++ )
 132         {
 133             int sum = 0;
 134             const uchar* src0 = img.ptr<uchar>(y + pt.y - r1) + x + pt.x;
 135             const uchar* src1 = img.ptr<uchar>(pt.y - r1) + pt.x;
 136             for( dy = -r1; dy <= r1; dy++, src0 += step, src1 += step )
 137                 for( dx = -r1; dx <= r1; dx++ )
 138                 {
 139                     int t = src0[dx] - src1[dx];
 140                     sum += t*t;
 141                 }
 142             sptr[x] = (float)sum;
 143         }
 144     }
 145 }
 146
 147
 148 void SelfSimDescriptor::compute(const Mat& img, vector<float>& descriptors, Size winStride,
 149                                 const vector<Point>& locations) const
 150 {
 151     CV_Assert( img.depth() == CV_8U );
 152
 153     winStride.width = std::max(winStride.width, 1);
 154     winStride.height = std::max(winStride.height, 1);
 155     Size gridSize = getGridSize(img.size(), winStride);
 156     int i, nwindows = locations.empty() ? gridSize.width*gridSize.height : (int)locations.size();
 157     int border = largeSize/2 + smallSize/2;
 158     int fsize = (int)getDescriptorSize();
 159     vector<float> tempFeature(fsize+1);
 160     descriptors.resize(fsize*nwindows + 1);
 161     Mat ssd(largeSize, largeSize, CV_32F), mappingMask;
 162     computeLogPolarMapping(mappingMask);
 163
 164 #if 0 //def _OPENMP
 165     int nthreads = cvGetNumThreads();
 166     #pragma omp parallel for num_threads(nthreads)
 167 #endif
 168     for( i = 0; i < nwindows; i++ )
 169     {
 170         Point pt;
 171         float* feature0 = &descriptors[fsize*i];
 172         float* feature = &tempFeature[0];
 173         int x, y, j;
 174
 175         if( !locations.empty() )
 176         {
 177             pt = locations[i];
 178             if( pt.x < border || pt.x >= img.cols - border ||
 179                 pt.y < border || pt.y >= img.rows - border )
 180             {
 181                 for( j = 0; j < fsize; j++ )
 182                     feature0[j] = 0.f;
 183                 continue;
 184             }
 185         }
 186         else
 187             pt = Point((i % gridSize.width)*winStride.width + border,
 188                        (i / gridSize.width)*winStride.height + border);
 189
 190         SSD(img, pt, ssd);
 191
 192             // Determine in the local neighborhood the largest difference and use for normalization
 193             float var_noise = 1000.f;
 194         for( y = -1; y <= 1 ; y++ )
 195                     for( x = -1 ; x <= 1 ; x++ )
 196                 var_noise = std::max(var_noise, ssd.at<float>(largeSize/2+y, largeSize/2+x));
 197
 198         for( j = 0; j <= fsize; j++ )
 199             feature[j] = FLT_MAX;
 200
 201             // Derive feature vector before exp(-x) computation
 202             // Idea: for all  x,a >= 0, a=const.   we have:
 203             //       max [ exp( -x / a) ] = exp ( -min(x) / a )
 204             // Thus, determine min(ssd) and store in feature[...]
 205         for( y = 0; y < ssd.rows; y++ )
 206         {
 207                     const schar *mappingMaskPtr = mappingMask.ptr<schar>(y);
 208                     const float *ssdPtr = ssd.ptr<float>(y);
 209                     for( x = 0 ; x < ssd.cols; x++ )
 210             {
 211                 int index = mappingMaskPtr[x];
 212                 feature[index] = std::max(feature[index], ssdPtr[x]);
 213                     }
 214             }
 215
 216         var_noise = -1.f/var_noise;
 217         for( j = 0; j < fsize; j++ )
 218             feature0[j] = feature[j]*var_noise;
 219         Mat _f(1, fsize, CV_32F, feature0);
 220         cv::exp(_f, _f);
 221     }
 222 }
 223
 224 void SelfSimDescriptor::computeLogPolarMapping(Mat& mappingMask) const
 225 {
 226     mappingMask.create(largeSize, largeSize, CV_8S);
 227
 228     // What we want is
 229     //           log_m (radius) = numberOfDistanceBuckets
 230     //  <==> log_10 (radius) / log_10 (m) = numberOfDistanceBuckets
 231     //  <==> log_10 (radius) / numberOfDistanceBuckets = log_10 (m)
 232     //  <==> m = 10 ^ log_10(m) = 10 ^ [log_10 (radius) / numberOfDistanceBuckets]
 233     //
 234     int radius = largeSize/2, angleBucketSize = 360 / numberOfAngles;
 235     int fsize = (int)getDescriptorSize();
 236     double inv_log10m = (double)numberOfDistanceBuckets/log10((double)radius);
 237
 238     for (int y=-radius ; y<=radius ; y++)
 239     {
 240         schar* mrow = mappingMask.ptr<schar>(y+radius);
 241         for (int x=-radius ; x<=radius ; x++)
 242         {
 243             int index = fsize;
 244             float dist = (float)std::sqrt((float)x*x + (float)y*y);
 245             int distNo = dist > 0 ? cvRound(log10(dist)*inv_log10m) : 0;
 246             if( startDistanceBucket <= distNo && distNo < numberOfDistanceBuckets )
 247             {
 248                 float angle = std::atan2( (float)y, (float)x ) / (float)CV_PI * 180.0f;
 249                 if (angle < 0) angle += 360.0f;
 250                 int angleInt = (cvRound(angle) + angleBucketSize/2) % 360;
 251                 int angleIndex = angleInt / angleBucketSize;
 252                 index = (distNo-startDistanceBucket)*numberOfAngles + angleIndex;
 253             }
 254             mrow[x + radius] = saturate_cast<schar>(index);
 255         }
 256     }
 257 }
 258
 259 }