1 /*M/////////////////////////////////////////////////////////////////////////////////////// 2 // 3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4 // 5 // By downloading, copying, installing or using the software you agree to this license. 6 // If you do not agree to this license, do not download, install, 7 // copy or use the software. 8 // 9 // 10 // License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. 14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved. 15 // Copyright (C) 2013, OpenCV Foundation, 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 the copyright holders 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 Intel Corporation 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 #include "precomp.hpp" 45 46 /****************************************************************************************\ 47 * PCA * 48 \****************************************************************************************/ 49 50 namespace cv 51 { 52 53 PCA::PCA() {} 54 55 PCA::PCA(InputArray data, InputArray _mean, int flags, int maxComponents) 56 { 57 operator()(data, _mean, flags, maxComponents); 58 } 59 60 PCA::PCA(InputArray data, InputArray _mean, int flags, double retainedVariance) 61 { 62 operator()(data, _mean, flags, retainedVariance); 63 } 64 65 PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, int maxComponents) 66 { 67 Mat data = _data.getMat(), _mean = __mean.getMat(); 68 int covar_flags = CV_COVAR_SCALE; 69 int len, in_count; 70 Size mean_sz; 71 72 CV_Assert( data.channels() == 1 ); 73 if( flags & CV_PCA_DATA_AS_COL ) 74 { 75 len = data.rows; 76 in_count = data.cols; 77 covar_flags |= CV_COVAR_COLS; 78 mean_sz = Size(1, len); 79 } 80 else 81 { 82 len = data.cols; 83 in_count = data.rows; 84 covar_flags |= CV_COVAR_ROWS; 85 mean_sz = Size(len, 1); 86 } 87 88 int count = std::min(len, in_count), out_count = count; 89 if( maxComponents > 0 ) 90 out_count = std::min(count, maxComponents); 91 92 // "scrambled" way to compute PCA (when cols(A)>rows(A)): 93 // B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y 94 if( len <= in_count ) 95 covar_flags |= CV_COVAR_NORMAL; 96 97 int ctype = std::max(CV_32F, data.depth()); 98 mean.create( mean_sz, ctype ); 99 100 Mat covar( count, count, ctype ); 101 102 if( !_mean.empty() ) 103 { 104 CV_Assert( _mean.size() == mean_sz ); 105 _mean.convertTo(mean, ctype); 106 covar_flags |= CV_COVAR_USE_AVG; 107 } 108 109 calcCovarMatrix( data, covar, mean, covar_flags, ctype ); 110 eigen( covar, eigenvalues, eigenvectors ); 111 112 if( !(covar_flags & CV_COVAR_NORMAL) ) 113 { 114 // CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A 115 // CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A' 116 Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols); 117 if( data.type() != ctype || tmp_mean.data == mean.data ) 118 { 119 data.convertTo( tmp_data, ctype ); 120 subtract( tmp_data, tmp_mean, tmp_data ); 121 } 122 else 123 { 124 subtract( data, tmp_mean, tmp_mean ); 125 tmp_data = tmp_mean; 126 } 127 128 Mat evects1(count, len, ctype); 129 gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1, 130 (flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0); 131 eigenvectors = evects1; 132 133 // normalize eigenvectors 134 int i; 135 for( i = 0; i < out_count; i++ ) 136 { 137 Mat vec = eigenvectors.row(i); 138 normalize(vec, vec); 139 } 140 } 141 142 if( count > out_count ) 143 { 144 // use clone() to physically copy the data and thus deallocate the original matrices 145 eigenvalues = eigenvalues.rowRange(0,out_count).clone(); 146 eigenvectors = eigenvectors.rowRange(0,out_count).clone(); 147 } 148 return *this; 149 } 150 151 void PCA::write(FileStorage& fs ) const 152 { 153 CV_Assert( fs.isOpened() ); 154 155 fs << "name" << "PCA"; 156 fs << "vectors" << eigenvectors; 157 fs << "values" << eigenvalues; 158 fs << "mean" << mean; 159 } 160 161 void PCA::read(const FileNode& fs) 162 { 163 CV_Assert( !fs.empty() ); 164 String name = (String)fs["name"]; 165 CV_Assert( name == "PCA" ); 166 167 cv::read(fs["vectors"], eigenvectors); 168 cv::read(fs["values"], eigenvalues); 169 cv::read(fs["mean"], mean); 170 } 171 172 template <typename T> 173 int computeCumulativeEnergy(const Mat& eigenvalues, double retainedVariance) 174 { 175 CV_DbgAssert( eigenvalues.type() == DataType<T>::type ); 176 177 Mat g(eigenvalues.size(), DataType<T>::type); 178 179 for(int ig = 0; ig < g.rows; ig++) 180 { 181 g.at<T>(ig, 0) = 0; 182 for(int im = 0; im <= ig; im++) 183 { 184 g.at<T>(ig,0) += eigenvalues.at<T>(im,0); 185 } 186 } 187 188 int L; 189 190 for(L = 0; L < eigenvalues.rows; L++) 191 { 192 double energy = g.at<T>(L, 0) / g.at<T>(g.rows - 1, 0); 193 if(energy > retainedVariance) 194 break; 195 } 196 197 L = std::max(2, L); 198 199 return L; 200 } 201 202 PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, double retainedVariance) 203 { 204 Mat data = _data.getMat(), _mean = __mean.getMat(); 205 int covar_flags = CV_COVAR_SCALE; 206 int len, in_count; 207 Size mean_sz; 208 209 CV_Assert( data.channels() == 1 ); 210 if( flags & CV_PCA_DATA_AS_COL ) 211 { 212 len = data.rows; 213 in_count = data.cols; 214 covar_flags |= CV_COVAR_COLS; 215 mean_sz = Size(1, len); 216 } 217 else 218 { 219 len = data.cols; 220 in_count = data.rows; 221 covar_flags |= CV_COVAR_ROWS; 222 mean_sz = Size(len, 1); 223 } 224 225 CV_Assert( retainedVariance > 0 && retainedVariance <= 1 ); 226 227 int count = std::min(len, in_count); 228 229 // "scrambled" way to compute PCA (when cols(A)>rows(A)): 230 // B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y 231 if( len <= in_count ) 232 covar_flags |= CV_COVAR_NORMAL; 233 234 int ctype = std::max(CV_32F, data.depth()); 235 mean.create( mean_sz, ctype ); 236 237 Mat covar( count, count, ctype ); 238 239 if( !_mean.empty() ) 240 { 241 CV_Assert( _mean.size() == mean_sz ); 242 _mean.convertTo(mean, ctype); 243 } 244 245 calcCovarMatrix( data, covar, mean, covar_flags, ctype ); 246 eigen( covar, eigenvalues, eigenvectors ); 247 248 if( !(covar_flags & CV_COVAR_NORMAL) ) 249 { 250 // CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A 251 // CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A' 252 Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols); 253 if( data.type() != ctype || tmp_mean.data == mean.data ) 254 { 255 data.convertTo( tmp_data, ctype ); 256 subtract( tmp_data, tmp_mean, tmp_data ); 257 } 258 else 259 { 260 subtract( data, tmp_mean, tmp_mean ); 261 tmp_data = tmp_mean; 262 } 263 264 Mat evects1(count, len, ctype); 265 gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1, 266 (flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0); 267 eigenvectors = evects1; 268 269 // normalize all eigenvectors 270 int i; 271 for( i = 0; i < eigenvectors.rows; i++ ) 272 { 273 Mat vec = eigenvectors.row(i); 274 normalize(vec, vec); 275 } 276 } 277 278 // compute the cumulative energy content for each eigenvector 279 int L; 280 if (ctype == CV_32F) 281 L = computeCumulativeEnergy<float>(eigenvalues, retainedVariance); 282 else 283 L = computeCumulativeEnergy<double>(eigenvalues, retainedVariance); 284 285 // use clone() to physically copy the data and thus deallocate the original matrices 286 eigenvalues = eigenvalues.rowRange(0,L).clone(); 287 eigenvectors = eigenvectors.rowRange(0,L).clone(); 288 289 return *this; 290 } 291 292 void PCA::project(InputArray _data, OutputArray result) const 293 { 294 Mat data = _data.getMat(); 295 CV_Assert( !mean.empty() && !eigenvectors.empty() && 296 ((mean.rows == 1 && mean.cols == data.cols) || (mean.cols == 1 && mean.rows == data.rows))); 297 Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols); 298 int ctype = mean.type(); 299 if( data.type() != ctype || tmp_mean.data == mean.data ) 300 { 301 data.convertTo( tmp_data, ctype ); 302 subtract( tmp_data, tmp_mean, tmp_data ); 303 } 304 else 305 { 306 subtract( data, tmp_mean, tmp_mean ); 307 tmp_data = tmp_mean; 308 } 309 if( mean.rows == 1 ) 310 gemm( tmp_data, eigenvectors, 1, Mat(), 0, result, GEMM_2_T ); 311 else 312 gemm( eigenvectors, tmp_data, 1, Mat(), 0, result, 0 ); 313 } 314 315 Mat PCA::project(InputArray data) const 316 { 317 Mat result; 318 project(data, result); 319 return result; 320 } 321 322 void PCA::backProject(InputArray _data, OutputArray result) const 323 { 324 Mat data = _data.getMat(); 325 CV_Assert( !mean.empty() && !eigenvectors.empty() && 326 ((mean.rows == 1 && eigenvectors.rows == data.cols) || 327 (mean.cols == 1 && eigenvectors.rows == data.rows))); 328 329 Mat tmp_data, tmp_mean; 330 data.convertTo(tmp_data, mean.type()); 331 if( mean.rows == 1 ) 332 { 333 tmp_mean = repeat(mean, data.rows, 1); 334 gemm( tmp_data, eigenvectors, 1, tmp_mean, 1, result, 0 ); 335 } 336 else 337 { 338 tmp_mean = repeat(mean, 1, data.cols); 339 gemm( eigenvectors, tmp_data, 1, tmp_mean, 1, result, GEMM_1_T ); 340 } 341 } 342 343 Mat PCA::backProject(InputArray data) const 344 { 345 Mat result; 346 backProject(data, result); 347 return result; 348 } 349 350 } 351 352 void cv::PCACompute(InputArray data, InputOutputArray mean, 353 OutputArray eigenvectors, int maxComponents) 354 { 355 PCA pca; 356 pca(data, mean, 0, maxComponents); 357 pca.mean.copyTo(mean); 358 pca.eigenvectors.copyTo(eigenvectors); 359 } 360 361 void cv::PCACompute(InputArray data, InputOutputArray mean, 362 OutputArray eigenvectors, double retainedVariance) 363 { 364 PCA pca; 365 pca(data, mean, 0, retainedVariance); 366 pca.mean.copyTo(mean); 367 pca.eigenvectors.copyTo(eigenvectors); 368 } 369 370 void cv::PCAProject(InputArray data, InputArray mean, 371 InputArray eigenvectors, OutputArray result) 372 { 373 PCA pca; 374 pca.mean = mean.getMat(); 375 pca.eigenvectors = eigenvectors.getMat(); 376 pca.project(data, result); 377 } 378 379 void cv::PCABackProject(InputArray data, InputArray mean, 380 InputArray eigenvectors, OutputArray result) 381 { 382 PCA pca; 383 pca.mean = mean.getMat(); 384 pca.eigenvectors = eigenvectors.getMat(); 385 pca.backProject(data, result); 386 } 387
