|
| 1 | +// This file is part of OpenCV project. |
| 2 | +// It is subject to the license terms in the LICENSE file found in the top-level directory |
| 3 | +// of this distribution and at http://opencv.org/license.html.#include "precomp.hpp" |
| 4 | +#include "precomp.hpp" |
| 5 | +#include "opencv2/imgproc.hpp" |
| 6 | +#include "opencv2/ximgproc/color_match.hpp" |
| 7 | + |
| 8 | +using namespace std; |
| 9 | + |
| 10 | +namespace cv { namespace ximgproc { |
| 11 | + |
| 12 | +void createQuaternionImage(InputArray _img, OutputArray _qimg) |
| 13 | +{ |
| 14 | + int type = _img.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type); |
| 15 | + CV_CheckType(depth, depth == CV_8U || depth == CV_32F || depth == CV_64F, "Depth must be CV_8U, CV_32F or CV_64F"); |
| 16 | + CV_Assert(_img.dims() == 2 && cn == 3); |
| 17 | + vector<Mat> qplane(4); |
| 18 | + vector<Mat> plane; |
| 19 | + split(_img, plane); |
| 20 | + qplane[0] = Mat::zeros(_img.size(), CV_64FC1); |
| 21 | + for (int i = 0; i < cn; i++) |
| 22 | + plane[i].convertTo(qplane[3-i], CV_64F); |
| 23 | + merge(qplane, _qimg); |
| 24 | +} |
| 25 | + |
| 26 | +void qconj(InputArray _img, OutputArray _qimg) |
| 27 | +{ |
| 28 | + int type = _img.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type); |
| 29 | + CV_CheckType(depth, depth == CV_32F || depth == CV_64F, "Depth must be CV_32F or CV_64F"); |
| 30 | + CV_Assert(_img.dims() == 2 && cn == 4); |
| 31 | + vector<Mat> qplane(4), plane; |
| 32 | + split(_img, plane); |
| 33 | + qplane[0] = plane[0]; |
| 34 | + qplane[1] = -plane[1]; |
| 35 | + qplane[2] = -plane[2]; |
| 36 | + qplane[3] = -plane[3]; |
| 37 | + merge(qplane, _qimg); |
| 38 | +} |
| 39 | + |
| 40 | +void qunitary(InputArray _img, OutputArray _qimg) |
| 41 | +{ |
| 42 | + int type = _img.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type); |
| 43 | + CV_Assert((depth == CV_64F) && _img.dims() == 2 && cn == 4); |
| 44 | + _img.copyTo(_qimg); |
| 45 | + Mat qimg = _qimg.getMat(); |
| 46 | + qimg.forEach<Vec4d>([](Vec4d &p, const int * /*position*/) -> void { |
| 47 | + double d = p[0] * p[0] + p[1] * p[1] + p[2] * p[2] + p[3] * p[3]; |
| 48 | + d = 1 / sqrt(d); |
| 49 | + p *= d; |
| 50 | + }); |
| 51 | +} |
| 52 | + |
| 53 | +void qdft(InputArray _img, OutputArray _qimg, int flags, bool sideLeft) |
| 54 | +{ |
| 55 | + CV_INSTRUMENT_REGION(); |
| 56 | + |
| 57 | + int type = _img.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type); |
| 58 | + CV_Assert(depth == CV_64F && _img.dims() == 2 && cn == 4); |
| 59 | + float c; |
| 60 | + if (sideLeft) |
| 61 | + c = 1; // Left qdft |
| 62 | + else |
| 63 | + c = -1; // right qdft |
| 64 | + |
| 65 | + vector<Mat> q; |
| 66 | + Mat img; |
| 67 | + img = _img.getMat(); |
| 68 | + |
| 69 | + CV_Assert(getOptimalDFTSize(img.rows) == img.rows && getOptimalDFTSize(img.cols) == img.cols); |
| 70 | + |
| 71 | + split(img, q); |
| 72 | + Mat c1r; |
| 73 | + Mat c1i; // Imaginary part of c1 =x' |
| 74 | + Mat c2r; // Real part of c2 =y' |
| 75 | + Mat c2i; // Imaginary part of c2=z' |
| 76 | + c1r = q[0].clone(); |
| 77 | + c1i = (q[1] + q[2] + q[3]) / sqrt(3); |
| 78 | + c2r = (q[2] - q[3]) / sqrt(2); |
| 79 | + c2i = c * (q[3] + q[2] - 2 * q[1]) / sqrt(6); |
| 80 | + vector<Mat> vc1 = { c1r,c1i }, vc2 = { c2r,c2i }; |
| 81 | + Mat c1, c2, C1, C2; |
| 82 | + merge(vc1, c1); |
| 83 | + merge(vc2, c2); |
| 84 | + if (flags& DFT_INVERSE) |
| 85 | + { |
| 86 | + dft(c1, C1, DFT_COMPLEX_OUTPUT | DFT_INVERSE|DFT_SCALE); |
| 87 | + dft(c2, C2, DFT_COMPLEX_OUTPUT | DFT_INVERSE | DFT_SCALE); |
| 88 | + } |
| 89 | + else |
| 90 | + { |
| 91 | + dft(c1, C1, DFT_COMPLEX_OUTPUT); |
| 92 | + dft(c2, C2, DFT_COMPLEX_OUTPUT); |
| 93 | + } |
| 94 | + split(C1, vc1); |
| 95 | + split(C2, vc2); |
| 96 | + vector<Mat> qdft(4); |
| 97 | + qdft[0] = vc1[0].clone(); |
| 98 | + qdft[1] = vc1[1] / sqrt(3) - c * 2 * vc2[1] / sqrt(6); |
| 99 | + qdft[2] = vc1[1] / sqrt(3) + vc2[0] / sqrt(2) + c * vc2[1] / sqrt(6); |
| 100 | + qdft[3] = vc1[1] / sqrt(3) - vc2[0] / sqrt(2) + c * vc2[1] / sqrt(6); |
| 101 | + Mat dst0; |
| 102 | + merge(qdft, dst0); |
| 103 | + dst0.copyTo(_qimg); |
| 104 | +} |
| 105 | + |
| 106 | + |
| 107 | +void qmultiply(InputArray src1, InputArray src2, OutputArray dst) |
| 108 | +{ |
| 109 | + int type = src1.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type); |
| 110 | + CV_Assert(depth == CV_64F && src1.dims() == 2 && cn == 4); |
| 111 | + type = src2.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type); |
| 112 | + CV_Assert(depth == CV_64F && src2.dims() == 2 && cn == 4); |
| 113 | + vector<Mat> q3(4); |
| 114 | + if (src1.rows() == src2.rows() && src1.cols() == src2.cols()) |
| 115 | + { |
| 116 | + vector<Mat> q1, q2; |
| 117 | + split(src1, q1); |
| 118 | + split(src2, q2); |
| 119 | + q3[0] = q1[0].mul(q2[0]) - q1[1].mul(q2[1]) - q1[2].mul(q2[2]) - q1[3].mul(q2[3]); |
| 120 | + q3[1] = q1[0].mul(q2[1]) + q1[1].mul(q2[0]) + q1[2].mul(q2[3]) - q1[3].mul(q2[2]); |
| 121 | + q3[2] = q1[0].mul(q2[2]) - q1[1].mul(q2[3]) + q1[2].mul(q2[0]) + q1[3].mul(q2[1]); |
| 122 | + q3[3] = q1[0].mul(q2[3]) + q1[1].mul(q2[2]) - q1[2].mul(q2[1]) + q1[3].mul(q2[0]); |
| 123 | + } |
| 124 | + else if (src1.rows() == 1 && src1.cols() == 1) |
| 125 | + { |
| 126 | + vector<Mat> q2; |
| 127 | + Vec4d q1 = src1.getMat().at<Vec4d>(0, 0); |
| 128 | + split(src2, q2); |
| 129 | + q3[0] = q1[0] * q2[0] - q1[1] * q2[1] - q1[2] * q2[2] - q1[3] * q2[3]; |
| 130 | + q3[1] = q1[0] * q2[1] + q1[1] * q2[0] + q1[2] * q2[3] - q1[3] * q2[2]; |
| 131 | + q3[2] = q1[0] * q2[2] - q1[1] * q2[3] + q1[2] * q2[0] + q1[3] * q2[1]; |
| 132 | + q3[3] = q1[0] * q2[3] + q1[1] * q2[2] - q1[2] * q2[1] + q1[3] * q2[0]; |
| 133 | + } |
| 134 | + else if (src2.rows() == 1 && src2.cols() == 1) |
| 135 | + { |
| 136 | + vector<Mat> q1; |
| 137 | + split(src1, q1); |
| 138 | + Vec4d q2 = src2.getMat().at<Vec4d>(0, 0); |
| 139 | + q3[0] = q1[0] * q2[0] - q1[1] * q2[1] - q1[2] * q2[2] - q1[3] * q2[3]; |
| 140 | + q3[1] = q1[0] * q2[1] + q1[1] * q2[0] + q1[2] * q2[3] - q1[3] * q2[2]; |
| 141 | + q3[2] = q1[0] * q2[2] - q1[1] * q2[3] + q1[2] * q2[0] + q1[3] * q2[1]; |
| 142 | + q3[3] = q1[0] * q2[3] + q1[1] * q2[2] - q1[2] * q2[1] + q1[3] * q2[0]; |
| 143 | + } |
| 144 | + else |
| 145 | + CV_Assert(src1.rows() == src2.rows() && src1.cols() == src2.cols()); |
| 146 | + merge(q3, dst); |
| 147 | + |
| 148 | +} |
| 149 | + |
| 150 | +void colorMatchTemplate(InputArray _image, InputArray _templ, OutputArray _result) |
| 151 | +{ |
| 152 | + CV_INSTRUMENT_REGION(); |
| 153 | + Mat image = _image.getMat(), imageF; |
| 154 | + CV_Assert(image.channels() == 3); |
| 155 | + Mat colorTemplate = _templ.getMat(), colorTemplateF; |
| 156 | + CV_Assert(colorTemplate.channels() == 3); |
| 157 | + int rr = max(getOptimalDFTSize(image.rows), getOptimalDFTSize(colorTemplate.rows)); |
| 158 | + int cc = max(getOptimalDFTSize(image.cols), getOptimalDFTSize(colorTemplate.cols)); |
| 159 | + Mat logo(rr, cc, CV_64FC3, Scalar::all(0)); |
| 160 | + Mat img = Mat(rr, cc, CV_64FC3, Scalar::all(0)); |
| 161 | + colorTemplate.convertTo(colorTemplateF, CV_64F, 1 / 256.), |
| 162 | + colorTemplateF.copyTo(logo(Rect(0, 0, colorTemplate.cols, colorTemplate.rows))); |
| 163 | + image.convertTo(imageF, CV_64F, 1 / 256.); |
| 164 | + imageF.copyTo(img(Rect(0, 0, image.cols, image.rows))); |
| 165 | + Mat qimg, qlogo; |
| 166 | + Mat qimgFFT, qimgIFFT, qlogoFFT; |
| 167 | + // Create quaternion image |
| 168 | + createQuaternionImage(img, qimg); |
| 169 | + createQuaternionImage(logo, qlogo); |
| 170 | + // quaternion fourier transform |
| 171 | + qdft(qimg, qimgFFT, 0, true); |
| 172 | + qdft(qimg, qimgIFFT, DFT_INVERSE, true); |
| 173 | + qdft(qlogo, qlogoFFT, 0, false); |
| 174 | + double sqrtnn = sqrt(static_cast<int>(qimgFFT.rows*qimgFFT.cols)); |
| 175 | + qimgFFT /= sqrtnn; |
| 176 | + qimgIFFT *= sqrtnn; |
| 177 | + qlogoFFT /= sqrtnn; |
| 178 | + Mat mu(1, 1, CV_64FC4, Scalar(0, 1, 1, 1) / sqrt(3.)); |
| 179 | + Mat qtmp, qlogopara, qlogoortho; |
| 180 | + qmultiply(mu, qlogoFFT, qtmp); |
| 181 | + qmultiply(qtmp, mu, qtmp); |
| 182 | + subtract(qlogoFFT, qtmp, qlogopara); |
| 183 | + qlogopara = qlogopara / 2; |
| 184 | + subtract(qlogoFFT, qlogopara, qlogoortho); |
| 185 | + Mat qcross1, qcross2, cqf, cqfi; |
| 186 | + qconj(qimgFFT, cqf); |
| 187 | + qconj(qimgIFFT, cqfi); |
| 188 | + qmultiply(cqf, qlogopara, qcross1); |
| 189 | + qmultiply(cqfi, qlogoortho, qcross2); |
| 190 | + Mat pwsp = qcross1 + qcross2; |
| 191 | + Mat crossCorr, pwspUnitary; |
| 192 | + qunitary(pwsp, pwspUnitary); |
| 193 | + qdft(pwspUnitary, crossCorr, DFT_INVERSE, false); |
| 194 | + vector<Mat> p; |
| 195 | + split(crossCorr, p); |
| 196 | + Mat imgcorr = (p[0].mul(p[0]) + p[1].mul(p[1]) + p[2].mul(p[2]) + p[3].mul(p[3])); |
| 197 | + sqrt(imgcorr, _result); |
| 198 | +} |
| 199 | +} |
| 200 | +} |
0 commit comments