1 // This is based on Rainer Lienhart contribution. Below is the original copyright:
3 /*M///////////////////////////////////////////////////////////////////////////////////////
5 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
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.
12 // University of Augsburg License Agreement
13 // For Open Source MultiMedia Computing (MMC) Library
15 // Copyright (C) 2007, University of Augsburg, Germany, all rights reserved.
16 // Third party copyrights are property of their respective owners.
18 // Redistribution and use in source and binary forms, with or without modification,
19 // are permitted provided that the following conditions are met:
21 // * Redistribution's of source code must retain the above copyright notice,
22 // this list of conditions and the following disclaimer.
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.
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.
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.
44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
45 // Author: Rainer Lienhart
46 // email: Rainer.Lienhart@informatik.uni-augsburg.de
47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
49 // Please cite the following two papers:
50 // 1. Shechtman, E., Irani, M.:
51 // Matching local self-similarities across images and videos.
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.
63 SelfSimDescriptor::SelfSimDescriptor()
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;
72 SelfSimDescriptor::SelfSimDescriptor(int _ssize, int _lsize,
73 int _startDistanceBucket,
74 int _numberOfDistanceBuckets, int _numberOfAngles)
78 startDistanceBucket = _startDistanceBucket;
79 numberOfDistanceBuckets = _numberOfDistanceBuckets;
80 numberOfAngles = _numberOfAngles;
83 SelfSimDescriptor::SelfSimDescriptor(const SelfSimDescriptor& ss)
85 smallSize = ss.smallSize;
86 largeSize = ss.largeSize;
87 startDistanceBucket = ss.startDistanceBucket;
88 numberOfDistanceBuckets = ss.numberOfDistanceBuckets;
89 numberOfAngles = ss.numberOfAngles;
92 SelfSimDescriptor::~SelfSimDescriptor()
96 SelfSimDescriptor& SelfSimDescriptor::operator = (const SelfSimDescriptor& ss)
100 smallSize = ss.smallSize;
101 largeSize = ss.largeSize;
102 startDistanceBucket = ss.startDistanceBucket;
103 numberOfDistanceBuckets = ss.numberOfDistanceBuckets;
104 numberOfAngles = ss.numberOfAngles;
109 size_t SelfSimDescriptor::getDescriptorSize() const
111 return numberOfAngles*(numberOfDistanceBuckets - startDistanceBucket);
114 Size SelfSimDescriptor::getGridSize( Size imgSize, Size winStride ) const
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);
123 // TODO: optimized with SSE2
124 void SelfSimDescriptor::SSD(const Mat& img, Point pt, Mat& ssd) const
126 int x, y, dx, dy, r0 = largeSize/2, r1 = smallSize/2;
128 for( y = -r0; y <= r0; y++ )
130 float* sptr = ssd.ptr<float>(y+r0) + r0;
131 for( x = -r0; x <= r0; x++ )
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++ )
139 int t = src0[dx] - src1[dx];
142 sptr[x] = (float)sum;
148 void SelfSimDescriptor::compute(const Mat& img, vector<float>& descriptors, Size winStride,
149 const vector<Point>& locations) const
151 CV_Assert( img.depth() == CV_8U );
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);
165 int nthreads = cvGetNumThreads();
166 #pragma omp parallel for num_threads(nthreads)
168 for( i = 0; i < nwindows; i++ )
171 float* feature0 = &descriptors[fsize*i];
172 float* feature = &tempFeature[0];
175 if( !locations.empty() )
178 if( pt.x < border || pt.x >= img.cols - border ||
179 pt.y < border || pt.y >= img.rows - border )
181 for( j = 0; j < fsize; j++ )
187 pt = Point((i % gridSize.width)*winStride.width + border,
188 (i / gridSize.width)*winStride.height + border);
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));
198 for( j = 0; j <= fsize; j++ )
199 feature[j] = FLT_MAX;
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++ )
207 const schar *mappingMaskPtr = mappingMask.ptr<schar>(y);
208 const float *ssdPtr = ssd.ptr<float>(y);
209 for( x = 0 ; x < ssd.cols; x++ )
211 int index = mappingMaskPtr[x];
212 feature[index] = std::max(feature[index], ssdPtr[x]);
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);
224 void SelfSimDescriptor::computeLogPolarMapping(Mat& mappingMask) const
226 mappingMask.create(largeSize, largeSize, CV_8S);
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]
234 int radius = largeSize/2, angleBucketSize = 360 / numberOfAngles;
235 int fsize = (int)getDescriptorSize();
236 double inv_log10m = (double)numberOfDistanceBuckets/log10((double)radius);
238 for (int y=-radius ; y<=radius ; y++)
240 schar* mrow = mappingMask.ptr<schar>(y+radius);
241 for (int x=-radius ; x<=radius ; x++)
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 )
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;
254 mrow[x + radius] = saturate_cast<schar>(index);