#include "precomp.hpp" #include "opencv2/face/mace.hpp" namespace cv { namespace face { // //! Rearrange the quadrants of Fourier image //! so that the origin is at the image center // static void shiftDFT(const Mat &src, Mat &dst) { Size size = src.size(); if (dst.empty() || (dst.size().width != size.width || dst.size().height != size.height)) { dst.create(src.size(), src.type()); } int cx = size.width/2; int cy = size.height/2; // image center Mat q1 = src(Rect(0, 0, cx,cy)); Mat q2 = src(Rect(cx,0, cx,cy)); Mat q3 = src(Rect(cx,cy,cx,cy)); Mat q4 = src(Rect(0, cy,cx,cy)); Mat d1 = dst(Rect(0, 0, cx,cy)); Mat d2 = dst(Rect(cx,0, cx,cy)); Mat d3 = dst(Rect(cx,cy,cx,cy)); Mat d4 = dst(Rect(0, cy,cx,cy)); if (src.data != dst.data){ q3.copyTo(d1); q4.copyTo(d2); q1.copyTo(d3); q2.copyTo(d4); } else { Mat tmp; q3.copyTo(tmp); q1.copyTo(d3); tmp.copyTo(d1); q4.copyTo(tmp); q2.copyTo(d4); tmp.copyTo(d2); } } // Computes 64-bit "cyclic redundancy check" sum, as specified in ECMA-182 static uint64 crc64( const uchar* data, size_t size, uint64 crc0=0 ) { static uint64 table[256]; static bool initialized = false; if( !initialized ) { for( int i = 0; i < 256; i++ ) { uint64 c = i; for( int j = 0; j < 8; j++ ) c = ((c & 1) ? CV_BIG_UINT(0xc96c5795d7870f42) : 0) ^ (c >> 1); table[i] = c; } initialized = true; } uint64 crc = ~crc0; for( size_t idx = 0; idx < size; idx++ ) crc = table[(uchar)crc ^ data[idx]] ^ (crc >> 8); return ~crc; } struct MACEImpl CV_FINAL : MACE { Mat_<Vec2d> maceFilter; // filled from compute() Mat convFilter; // optional random convolution (cancellable) int IMGSIZE; // images will get resized to this double threshold; // minimal "sameness" threshold from the train images MACEImpl(int siz) : IMGSIZE(siz), threshold(DBL_MAX) {} void salt(const String &passphrase) CV_OVERRIDE { theRNG().state = ((int64)crc64((uchar*)passphrase.c_str(), passphrase.size())); convFilter.create(IMGSIZE, IMGSIZE, CV_64F); randn(convFilter, 0, 1.0/(IMGSIZE*IMGSIZE)); } Mat dftImage(Mat img) const { Mat gray; resize(img, gray, Size(IMGSIZE,IMGSIZE)) ; if (gray.channels() > 1) cvtColor(gray, gray, COLOR_BGR2GRAY); equalizeHist(gray, gray); gray.convertTo(gray, CV_64F); if (! convFilter.empty()) { // optional, but unfortunately, it has to happen after resize/equalize ops. filter2D(gray, gray, CV_64F, convFilter); } Mat input[2] = {gray, Mat(gray.size(), gray.type(), 0.0)}; Mat complexInput; merge(input, 2, complexInput); Mat_<Vec2d> dftImg(IMGSIZE*2, IMGSIZE*2, 0.0); complexInput.copyTo(dftImg(Rect(0,0,IMGSIZE,IMGSIZE))); dft(dftImg, dftImg); return std::move(dftImg); } // compute the mace filter: `h = D(-1) * X * (X(+) * D(-1) * X)(-1) * C` void compute(std::vector<Mat> images) { return compute(images, false); } void compute(std::vector<Mat> images, bool isdft) { int size = (int)images.size(); int IMGSIZE_2X = IMGSIZE * 2; int TOTALPIXEL = IMGSIZE_2X * IMGSIZE_2X; Mat_<double> D(TOTALPIXEL, 1, 0.0); Mat_<Vec2d> S(TOTALPIXEL, size, Vec2d(0,0)); Mat_<Vec2d> SPLUS(size, TOTALPIXEL, Vec2d(0,0)); for (int i=0; i<size; i++) { Mat_<Vec2d> dftImg = isdft ? images[i] : dftImage(images[i]); for (int l=0; l<IMGSIZE_2X; l++) { for (int m=0; m<IMGSIZE_2X; m++) { int j = l * IMGSIZE_2X + m; Vec2d s = dftImg(l, m); S(j, i) = s; SPLUS(i, j) = Vec2d(s[0], -s[1]); D(j, 0) += (s[0]*s[0]) + (s[1]*s[1]); } } } Mat_<double> DSQ; cv::sqrt(D, DSQ); Mat_<double> DINV = TOTALPIXEL * size / DSQ; Mat_<Vec2d> DINV_S(TOTALPIXEL, size); Mat_<Vec2d> SPLUS_DINV(size, TOTALPIXEL); for (int l=0; l<size; l++) { for (int m=0; m<TOTALPIXEL; m++) { DINV_S(m, l) = S(m,l) * DINV(m,0); SPLUS_DINV(l, m) = SPLUS(l,m) * DINV(m,0); } } Mat_<Vec2d> SPLUS_DINV_S = SPLUS_DINV * S; Mat_<double> SPLUS_DINV_S_INV_1(2*size, 2*size, 0.0); for (int l=0; l<size; l++) { for (int m=0; m<size; m++) { Vec2d s = SPLUS_DINV_S(l, m); SPLUS_DINV_S_INV_1(l, m) = s[0]; SPLUS_DINV_S_INV_1(l+size, m+size) = s[0]; SPLUS_DINV_S_INV_1(l, m+size) = s[1]; SPLUS_DINV_S_INV_1(l+size, m) = -s[1]; } } invert(SPLUS_DINV_S_INV_1, SPLUS_DINV_S_INV_1); Mat_<Vec2d> SPLUS_DINV_S_INV(size, size); for (int l=0; l<size; l++) { for (int m=0; m<size; m++) { SPLUS_DINV_S_INV(l, m) = Vec2d(SPLUS_DINV_S_INV_1(l,m), SPLUS_DINV_S_INV_1(l,m+size)); } } Mat_<Vec2d> Hmace = DINV_S * SPLUS_DINV_S_INV; Mat_<Vec2d> C(size, 1, Vec2d(1,0)); maceFilter = Mat(Hmace * C).reshape(2,IMGSIZE_2X); } // get the lowest (worst) positive train correlation, // our lower bound threshold for the "same()" test later double computeThreshold(const std::vector<Mat> &images, bool isdft) const { double best=DBL_MAX; for (size_t i=0; i<images.size(); i++) { double d = correlate(images[i], isdft); if (d < best) { best = d; } } return best; } // convolute macefilter and dft image, // calculate the peak to sidelobe ratio // on the real part of the inverse dft double correlate(const Mat &img) const { return correlate(img, false); } double correlate(const Mat &img, bool isdft) const { if (maceFilter.empty()) return -1; // not trained. int IMGSIZE_2X = IMGSIZE * 2; Mat dftImg = isdft ? img : dftImage(img); mulSpectrums(dftImg, maceFilter, dftImg, DFT_ROWS, true); dft(dftImg, dftImg, DFT_INVERSE|DFT_SCALE, 0); Mat chn[2]; split(dftImg, chn); Mat_<double> re; shiftDFT(chn[0], re); double m1,M1; minMaxLoc(re, &m1, &M1, 0, 0); double peakCorrPlaneEnergy = M1 / sqrt(sum(re)[0]); re -= m1; double value=0; double num=0; int rad_1=int(floor((double)(45.0/64.0)*(double)IMGSIZE)); int rad_2=int(floor((double)(27.0/64.0)*(double)IMGSIZE)); // cache a few pow's and sqrts std::vector<double> r2(IMGSIZE_2X); Mat_<double> radtab(IMGSIZE_2X,IMGSIZE_2X); for (int l=0; l<IMGSIZE_2X; l++) { r2[l] = (l-IMGSIZE) * (l-IMGSIZE); } for (int l=0; l<IMGSIZE_2X; l++) { for (int m=l+1; m<IMGSIZE_2X; m++) { double rad = sqrt(r2[m] + r2[l]); radtab(l,m) = radtab(m,l) = rad; } } // mean of the sidelobe area: for (int l=0; l<IMGSIZE_2X; l++) { for (int m=0; m<IMGSIZE_2X; m++) { double rad = radtab(l,m); if (rad < rad_1) { if (rad > rad_2) { value += re(l,m); num++; } } } } value /= num; // normalize it double std2=0; for (int l=0; l<IMGSIZE_2X; l++) { for (int m=0; m<IMGSIZE_2X; m++) { double rad = radtab(l,m); if (rad < rad_1) { if (rad > rad_2) { double d = (value - re(l,m)); std2 += d * d; } } } } std2 /= num; std2 = sqrt(std2); double sca = re(IMGSIZE, IMGSIZE); double peakToSideLobeRatio = (sca - value) / std2; return 100.0 * peakToSideLobeRatio * peakCorrPlaneEnergy; } // MACE interface void train(InputArrayOfArrays input) CV_OVERRIDE { std::vector<Mat> images, dftImg; input.getMatVector(images); for (size_t i=0; i<images.size(); i++) { // cache dft images dftImg.push_back(dftImage(images[i])); } compute(dftImg, true); threshold = computeThreshold(dftImg, true); } bool same(InputArray img) const CV_OVERRIDE { return correlate(img.getMat()) >= threshold; } // cv::Algorithm: bool empty() const CV_OVERRIDE { return maceFilter.empty() || IMGSIZE == 0; } String getDefaultName () const CV_OVERRIDE { return String("MACE"); } void clear() CV_OVERRIDE { maceFilter.release(); convFilter.release(); } void write(cv::FileStorage &fs) const CV_OVERRIDE { fs << "mace" << maceFilter; fs << "conv" << convFilter; fs << "threshold" << threshold; } void read(const cv::FileNode &fn) CV_OVERRIDE { fn["mace"] >> maceFilter; fn["conv"] >> convFilter; fn["threshold"] >> threshold; IMGSIZE = maceFilter.cols/2; } }; cv::Ptr<MACE> MACE::create(int siz) { return makePtr<MACEImpl>(siz); } cv::Ptr<MACE> MACE::load(const String &filename, const String &objname) { return Algorithm::load<MACE>(filename, objname); } } /* namespace face */ } /* namespace cv */