c++版 mulSpectrums

来源：互联网发布：qq会员免费开通软件编辑：程序博客网时间：2024/06/06 21:39
但是没有测试效果：//---------------------------- #define VAL(buf, elem) (((T*)((char*)data ## buf + (step ## buf * (elem))))[0]) #define MUL_SPECTRUMS_COL(A, B, C) \     VAL(C, 0) = VAL(A, 0) * VAL(B, 0); \     for (size_t j = 1; j <= rows - 2; j += 2) \     { \         double a_re = VAL(A, j), a_im = VAL(A, j + 1); \         double b_re = VAL(B, j), b_im = VAL(B, j + 1); \         if (conjB) b_im = -b_im; \         double c_re = a_re * b_re - a_im * b_im; \         double c_im = a_re * b_im + a_im * b_re; \         VAL(C, j) = (T)c_re; VAL(C, j + 1) = (T)c_im; \     } \     if ((rows & 1) == 0) \         VAL(C, rows-1) = VAL(A, rows-1) * VAL(B, rows-1) template <typename T, bool conjB> static inline void mulSpectrums_processCol_noinplace(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows) {     MUL_SPECTRUMS_COL(A, B, C); } template <typename T, bool conjB> static inline void mulSpectrums_processCol_inplaceA(const T* dataB, T* dataAC, size_t stepB, size_t stepAC, size_t rows) {     MUL_SPECTRUMS_COL(AC, B, AC); } template <typename T, bool conjB, bool inplaceA> static inline void mulSpectrums_processCol(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows) {     if (inplaceA)         mulSpectrums_processCol_inplaceA<T, conjB>(dataB, dataC, stepB, stepC, rows);     else         mulSpectrums_processCol_noinplace<T, conjB>(dataA, dataB, dataC, stepA, stepB, stepC, rows); } #undef MUL_SPECTRUMS_COL #undef VAL #define VAL(buf, elem) (data ## buf[(elem)]) #define MUL_SPECTRUMS_ROW(A, B, C) \     for (size_t j = j0; j < j1; j += 2) \     { \         double a_re = VAL(A, j), a_im = VAL(A, j + 1); \         double b_re = VAL(B, j), b_im = VAL(B, j + 1); \         if (conjB) b_im = -b_im; \         double c_re = a_re * b_re - a_im * b_im; \         double c_im = a_re * b_im + a_im * b_re; \         VAL(C, j) = (T)c_re; VAL(C, j + 1) = (T)c_im; \     } template <typename T, bool conjB> static inline void mulSpectrums_processRow_noinplace(const T* dataA, const T* dataB, T* dataC, size_t j0, size_t j1) {     MUL_SPECTRUMS_ROW(A, B, C); } template <typename T, bool conjB> static inline void mulSpectrums_processRow_inplaceA(const T* dataB, T* dataAC, size_t j0, size_t j1) {     MUL_SPECTRUMS_ROW(AC, B, AC); } template <typename T, bool conjB, bool inplaceA> static inline void mulSpectrums_processRow(const T* dataA, const T* dataB, T* dataC, size_t j0, size_t j1) {     if (inplaceA)         mulSpectrums_processRow_inplaceA<T, conjB>(dataB, dataC, j0, j1);     else         mulSpectrums_processRow_noinplace<T, conjB>(dataA, dataB, dataC, j0, j1); } #undef MUL_SPECTRUMS_ROW #undef VAL template <typename T, bool conjB, bool inplaceA> static inline void mulSpectrums_processCols(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows, size_t cols) {     mulSpectrums_processCol<T, conjB, inplaceA>(dataA, dataB, dataC, stepA, stepB, stepC, rows);     if ((cols & 1) == 0)     {         mulSpectrums_processCol<T, conjB, inplaceA>(dataA + cols - 1, dataB + cols - 1, dataC + cols - 1, stepA, stepB, stepC, rows);     } } template <typename T, bool conjB, bool inplaceA> static inline void mulSpectrums_processRows(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows, size_t cols, size_t j0, size_t j1, bool is_1d_CN1) {     while (rows-- > 0)     {         if (is_1d_CN1)             dataC[0] = dataA[0]*dataB[0];         mulSpectrums_processRow<T, conjB, inplaceA>(dataA, dataB, dataC, j0, j1);         if (is_1d_CN1 && (cols & 1) == 0)             dataC[j1] = dataA[j1]*dataB[j1];         dataA = (const T*)(((char*)dataA) + stepA);         dataB = (const T*)(((char*)dataB) + stepB);         dataC =       (T*)(((char*)dataC) + stepC);     } } template <typename T, bool conjB, bool inplaceA> static inline void mulSpectrums_Impl_(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows, size_t cols, size_t j0, size_t j1, bool is_1d, bool isCN1) {     if (!is_1d && isCN1)     {         mulSpectrums_processCols<T, conjB, inplaceA>(dataA, dataB, dataC, stepA, stepB, stepC, rows, cols);     }     mulSpectrums_processRows<T, conjB, inplaceA>(dataA, dataB, dataC, stepA, stepB, stepC, rows, cols, j0, j1, is_1d && isCN1); } template <typename T, bool conjB> static inline void mulSpectrums_Impl(const T* dataA, const T* dataB, T* dataC, size_t stepA, size_t stepB, size_t stepC, size_t rows, size_t cols, size_t j0, size_t j1, bool is_1d, bool isCN1) {     if (dataA == dataC)         mulSpectrums_Impl_<T, conjB, true>(dataA, dataB, dataC, stepA, stepB, stepC, rows, cols, j0, j1, is_1d, isCN1);     else         mulSpectrums_Impl_<T, conjB, false>(dataA, dataB, dataC, stepA, stepB, stepC, rows, cols, j0, j1, is_1d, isCN1); } void mulSpectrums( InputArray _srcA, InputArray _srcB,                        OutputArray _dst, int flags, bool conjB ) { //    CV_INSTRUMENT_REGION() //    CV_OCL_RUN(_dst.isUMat() && _srcA.dims() <= 2 && _srcB.dims() <= 2, //            ocl_mulSpectrums(_srcA, _srcB, _dst, flags, conjB))     Mat srcA = _srcA.getMat(), srcB = _srcB.getMat();     int depth = srcA.depth(), cn = srcA.channels(), type = srcA.type();     size_t rows = srcA.rows, cols = srcA.cols;     CV_Assert( type == srcB.type() && srcA.size() == srcB.size() );     CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 );     _dst.create( srcA.rows, srcA.cols, type );     Mat dst = _dst.getMat();     // correct inplace support     // Case 'dst.data == srcA.data' is handled by implementation,     // because it is used frequently (filter2D, matchTemplate)     if (dst.data == srcB.data)         srcB = srcB.clone(); // workaround for B only     bool is_1d = (flags & DFT_ROWS)         || (rows == 1)         || (cols == 1 && srcA.isContinuous() && srcB.isContinuous() && dst.isContinuous());     if( is_1d && !(flags & DFT_ROWS) )         cols = cols + rows - 1, rows = 1;     bool isCN1 = cn == 1;     size_t j0 = isCN1 ? 1 : 0;     size_t j1 = cols*cn - (((cols & 1) == 0 && cn == 1) ? 1 : 0);     if (depth == CV_32F)     {         const float* dataA = srcA.ptr<float>();         const float* dataB = srcB.ptr<float>();         float* dataC = dst.ptr<float>();         if (!conjB)             mulSpectrums_Impl<float, false>(dataA, dataB, dataC, srcA.step, srcB.step, dst.step, rows, cols, j0, j1, is_1d, isCN1);         else             mulSpectrums_Impl<float, true>(dataA, dataB, dataC, srcA.step, srcB.step, dst.step, rows, cols, j0, j1, is_1d, isCN1);     }     else     {         const double* dataA = srcA.ptr<double>();         const double* dataB = srcB.ptr<double>();         double* dataC = dst.ptr<double>();         if (!conjB)             mulSpectrums_Impl<double, false>(dataA, dataB, dataC, srcA.step, srcB.step, dst.step, rows, cols, j0, j1, is_1d, isCN1);         else             mulSpectrums_Impl<double, true>(dataA, dataB, dataC, srcA.step, srcB.step, dst.step, rows, cols, j0, j1, is_1d, isCN1);     } } //-------------------------int main(){    bool HOG = true;    bool FIXEDWINDOW = false;    bool MULTISCALE = true;    int w=10;    int h=10;    fftwf_complex * r   =            (fftwf_complex*) fftwf_malloc(sizeof (fftwf_complex) *w*h );    for (size_t i = 0; i < h; i++) {        for (size_t j = 0; j < w; j++) {                r[i*w+j][0]= i*j*0.1+1;                r[i*w+j][1]= i*j*0.2+1;        }    }    fftwf_complex * r2   =            (fftwf_complex*) fftwf_malloc(sizeof (fftwf_complex) *w*h );    for (size_t i = 0; i < h; i++) {        for (size_t j = 0; j < w; j++) {                r2[i*w+j][0]= i*j*0.5+2;                r2[i*w+j][1]= i*j*0.9+2;        }    }      SimpleTracker  tracker(HOG, FIXEDWINDOW, MULTISCALE );fftwf_complex * res=tracker.mulSpectrums(r,r2,w,h,0);for (size_t i = 0; i < h; i++) {    for (size_t j = 0; j < w; j++) {           printf("%f+%fi ", r2[i*w+j][0],r2[i*w+j][1]);    }     printf("\n");}
0 0