phasecorr.cpp

00001 /*********************************************************************
00002  * Software License Agreement (BSD License)
00003  *
00004  *  Copyright (c) 2008-2011, William Lucas
00005  *  All rights reserved.
00006  *
00007  *  Redistribution and use in source and binary forms, with or without
00008  *  modification, are permitted provided that the following conditions
00009  *  are met:
00010  *
00011  *   * Redistributions of source code must retain the above copyright
00012  *     notice, this list of conditions and the following disclaimer.
00013  *   * Redistributions in binary form must reproduce the above
00014  *     copyright notice, this list of conditions and the following
00015  *     disclaimer in the documentation and/or other materials provided
00016  *     with the distribution.
00017  *   * Neither the name of the Willow Garage nor the names of its
00018  *     contributors may be used to endorse or promote products derived
00019  *     from this software without specific prior written permission.
00020  *
00021  *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00022  *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00023  *  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
00024  *  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
00025  *  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
00026  *  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
00027  *  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
00028  *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
00029  *  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00030  *  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
00031  *  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
00032  *  POSSIBILITY OF SUCH DAMAGE.
00033  *********************************************************************/
00034 
00035 #include "precomp.hpp"
00036 #include <vector>
00037 
00038 namespace cv
00039 {
00040 
00041 static void magSpectrums( InputArray _src, OutputArray _dst)
00042 {
00043     Mat src = _src.getMat();
00044     int depth = src.depth(), cn = src.channels(), type = src.type();
00045     int rows = src.rows, cols = src.cols;
00046     int j, k;
00047 
00048     CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 );
00049 
00050     if(src.depth() == CV_32F)
00051         _dst.create( src.rows, src.cols, CV_32FC1 );
00052     else
00053         _dst.create( src.rows, src.cols, CV_64FC1 );
00054 
00055     Mat dst = _dst.getMat();
00056     dst.setTo(0);//Mat elements are not equal to zero by default!
00057 
00058     bool is_1d = (rows == 1 || (cols == 1 && src.isContinuous() && dst.isContinuous()));
00059 
00060     if( is_1d )
00061         cols = cols + rows - 1, rows = 1;
00062 
00063     int ncols = cols*cn;
00064     int j0 = cn == 1;
00065     int j1 = ncols - (cols % 2 == 0 && cn == 1);
00066 
00067     if( depth == CV_32F )
00068     {
00069         const float* dataSrc = src.ptr<float>();
00070         float* dataDst = dst.ptr<float>();
00071 
00072         size_t stepSrc = src.step/sizeof(dataSrc[0]);
00073         size_t stepDst = dst.step/sizeof(dataDst[0]);
00074 
00075         if( !is_1d && cn == 1 )
00076         {
00077             for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
00078             {
00079                 if( k == 1 )
00080                     dataSrc += cols - 1, dataDst += cols - 1;
00081                 dataDst[0] = dataSrc[0]*dataSrc[0];
00082                 if( rows % 2 == 0 )
00083                     dataDst[(rows-1)*stepDst] = dataSrc[(rows-1)*stepSrc]*dataSrc[(rows-1)*stepSrc];
00084 
00085                 for( j = 1; j <= rows - 2; j += 2 )
00086                 {
00087                     dataDst[j*stepDst] = (float)std::sqrt((double)dataSrc[j*stepSrc]*dataSrc[j*stepSrc] +
00088                                                           (double)dataSrc[(j+1)*stepSrc]*dataSrc[(j+1)*stepSrc]);
00089                 }
00090 
00091                 if( k == 1 )
00092                     dataSrc -= cols - 1, dataDst -= cols - 1;
00093             }
00094         }
00095 
00096         for( ; rows--; dataSrc += stepSrc, dataDst += stepDst )
00097         {
00098             if( is_1d && cn == 1 )
00099             {
00100                 dataDst[0] = dataSrc[0]*dataSrc[0];
00101                 if( cols % 2 == 0 )
00102                     dataDst[j1] = dataSrc[j1]*dataSrc[j1];
00103             }
00104 
00105             for( j = j0; j < j1; j += 2 )
00106             {
00107                 dataDst[j] = (float)std::sqrt((double)dataSrc[j]*dataSrc[j] + (double)dataSrc[j+1]*dataSrc[j+1]);
00108             }
00109         }
00110     }
00111     else
00112     {
00113         const double* dataSrc = src.ptr<double>();
00114         double* dataDst = dst.ptr<double>();
00115 
00116         size_t stepSrc = src.step/sizeof(dataSrc[0]);
00117         size_t stepDst = dst.step/sizeof(dataDst[0]);
00118 
00119         if( !is_1d && cn == 1 )
00120         {
00121             for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
00122             {
00123                 if( k == 1 )
00124                     dataSrc += cols - 1, dataDst += cols - 1;
00125                 dataDst[0] = dataSrc[0]*dataSrc[0];
00126                 if( rows % 2 == 0 )
00127                     dataDst[(rows-1)*stepDst] = dataSrc[(rows-1)*stepSrc]*dataSrc[(rows-1)*stepSrc];
00128 
00129                 for( j = 1; j <= rows - 2; j += 2 )
00130                 {
00131                     dataDst[j*stepDst] = std::sqrt(dataSrc[j*stepSrc]*dataSrc[j*stepSrc] +
00132                                                    dataSrc[(j+1)*stepSrc]*dataSrc[(j+1)*stepSrc]);
00133                 }
00134 
00135                 if( k == 1 )
00136                     dataSrc -= cols - 1, dataDst -= cols - 1;
00137             }
00138         }
00139 
00140         for( ; rows--; dataSrc += stepSrc, dataDst += stepDst )
00141         {
00142             if( is_1d && cn == 1 )
00143             {
00144                 dataDst[0] = dataSrc[0]*dataSrc[0];
00145                 if( cols % 2 == 0 )
00146                     dataDst[j1] = dataSrc[j1]*dataSrc[j1];
00147             }
00148 
00149             for( j = j0; j < j1; j += 2 )
00150             {
00151                 dataDst[j] = std::sqrt(dataSrc[j]*dataSrc[j] + dataSrc[j+1]*dataSrc[j+1]);
00152             }
00153         }
00154     }
00155 }
00156 
00157 static void divSpectrums( InputArray _srcA, InputArray _srcB, OutputArray _dst, int flags, bool conjB)
00158 {
00159     Mat srcA = _srcA.getMat(), srcB = _srcB.getMat();
00160     int depth = srcA.depth(), cn = srcA.channels(), type = srcA.type();
00161     int rows = srcA.rows, cols = srcA.cols;
00162     int j, k;
00163 
00164     CV_Assert( type == srcB.type() && srcA.size() == srcB.size() );
00165     CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 );
00166 
00167     _dst.create( srcA.rows, srcA.cols, type );
00168     Mat dst = _dst.getMat();
00169 
00170     bool is_1d = (flags & DFT_ROWS) || (rows == 1 || (cols == 1 &&
00171              srcA.isContinuous() && srcB.isContinuous() && dst.isContinuous()));
00172 
00173     if( is_1d && !(flags & DFT_ROWS) )
00174         cols = cols + rows - 1, rows = 1;
00175 
00176     int ncols = cols*cn;
00177     int j0 = cn == 1;
00178     int j1 = ncols - (cols % 2 == 0 && cn == 1);
00179 
00180     if( depth == CV_32F )
00181     {
00182         const float* dataA = srcA.ptr<float>();
00183         const float* dataB = srcB.ptr<float>();
00184         float* dataC = dst.ptr<float>();
00185         float eps = FLT_EPSILON; // prevent div0 problems
00186 
00187         size_t stepA = srcA.step/sizeof(dataA[0]);
00188         size_t stepB = srcB.step/sizeof(dataB[0]);
00189         size_t stepC = dst.step/sizeof(dataC[0]);
00190 
00191         if( !is_1d && cn == 1 )
00192         {
00193             for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
00194             {
00195                 if( k == 1 )
00196                     dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
00197                 dataC[0] = dataA[0] / (dataB[0] + eps);
00198                 if( rows % 2 == 0 )
00199                     dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / (dataB[(rows-1)*stepB] + eps);
00200                 if( !conjB )
00201                     for( j = 1; j <= rows - 2; j += 2 )
00202                     {
00203                         double denom = (double)dataB[j*stepB]*dataB[j*stepB] +
00204                                        (double)dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + (double)eps;
00205 
00206                         double re = (double)dataA[j*stepA]*dataB[j*stepB] +
00207                                     (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
00208 
00209                         double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] -
00210                                     (double)dataA[j*stepA]*dataB[(j+1)*stepB];
00211 
00212                         dataC[j*stepC] = (float)(re / denom);
00213                         dataC[(j+1)*stepC] = (float)(im / denom);
00214                     }
00215                 else
00216                     for( j = 1; j <= rows - 2; j += 2 )
00217                     {
00218 
00219                         double denom = (double)dataB[j*stepB]*dataB[j*stepB] +
00220                                        (double)dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + (double)eps;
00221 
00222                         double re = (double)dataA[j*stepA]*dataB[j*stepB] -
00223                                     (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
00224 
00225                         double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] +
00226                                     (double)dataA[j*stepA]*dataB[(j+1)*stepB];
00227 
00228                         dataC[j*stepC] = (float)(re / denom);
00229                         dataC[(j+1)*stepC] = (float)(im / denom);
00230                     }
00231                 if( k == 1 )
00232                     dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1;
00233             }
00234         }
00235 
00236         for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC )
00237         {
00238             if( is_1d && cn == 1 )
00239             {
00240                 dataC[0] = dataA[0] / (dataB[0] + eps);
00241                 if( cols % 2 == 0 )
00242                     dataC[j1] = dataA[j1] / (dataB[j1] + eps);
00243             }
00244 
00245             if( !conjB )
00246                 for( j = j0; j < j1; j += 2 )
00247                 {
00248                     double denom = (double)(dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps);
00249                     double re = (double)(dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1]);
00250                     double im = (double)(dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1]);
00251                     dataC[j] = (float)(re / denom);
00252                     dataC[j+1] = (float)(im / denom);
00253                 }
00254             else
00255                 for( j = j0; j < j1; j += 2 )
00256                 {
00257                     double denom = (double)(dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps);
00258                     double re = (double)(dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1]);
00259                     double im = (double)(dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1]);
00260                     dataC[j] = (float)(re / denom);
00261                     dataC[j+1] = (float)(im / denom);
00262                 }
00263         }
00264     }
00265     else
00266     {
00267         const double* dataA = srcA.ptr<double>();
00268         const double* dataB = srcB.ptr<double>();
00269         double* dataC = dst.ptr<double>();
00270         double eps = DBL_EPSILON; // prevent div0 problems
00271 
00272         size_t stepA = srcA.step/sizeof(dataA[0]);
00273         size_t stepB = srcB.step/sizeof(dataB[0]);
00274         size_t stepC = dst.step/sizeof(dataC[0]);
00275 
00276         if( !is_1d && cn == 1 )
00277         {
00278             for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
00279             {
00280                 if( k == 1 )
00281                     dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
00282                 dataC[0] = dataA[0] / (dataB[0] + eps);
00283                 if( rows % 2 == 0 )
00284                     dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / (dataB[(rows-1)*stepB] + eps);
00285                 if( !conjB )
00286                     for( j = 1; j <= rows - 2; j += 2 )
00287                     {
00288                         double denom = dataB[j*stepB]*dataB[j*stepB] +
00289                                        dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + eps;
00290 
00291                         double re = dataA[j*stepA]*dataB[j*stepB] +
00292                                     dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
00293 
00294                         double im = dataA[(j+1)*stepA]*dataB[j*stepB] -
00295                                     dataA[j*stepA]*dataB[(j+1)*stepB];
00296 
00297                         dataC[j*stepC] = re / denom;
00298                         dataC[(j+1)*stepC] = im / denom;
00299                     }
00300                 else
00301                     for( j = 1; j <= rows - 2; j += 2 )
00302                     {
00303                         double denom = dataB[j*stepB]*dataB[j*stepB] +
00304                                        dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + eps;
00305 
00306                         double re = dataA[j*stepA]*dataB[j*stepB] -
00307                                     dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
00308 
00309                         double im = dataA[(j+1)*stepA]*dataB[j*stepB] +
00310                                     dataA[j*stepA]*dataB[(j+1)*stepB];
00311 
00312                         dataC[j*stepC] = re / denom;
00313                         dataC[(j+1)*stepC] = im / denom;
00314                     }
00315                 if( k == 1 )
00316                     dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1;
00317             }
00318         }
00319 
00320         for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC )
00321         {
00322             if( is_1d && cn == 1 )
00323             {
00324                 dataC[0] = dataA[0] / (dataB[0] + eps);
00325                 if( cols % 2 == 0 )
00326                     dataC[j1] = dataA[j1] / (dataB[j1] + eps);
00327             }
00328 
00329             if( !conjB )
00330                 for( j = j0; j < j1; j += 2 )
00331                 {
00332                     double denom = dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps;
00333                     double re = dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1];
00334                     double im = dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1];
00335                     dataC[j] = re / denom;
00336                     dataC[j+1] = im / denom;
00337                 }
00338             else
00339                 for( j = j0; j < j1; j += 2 )
00340                 {
00341                     double denom = dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps;
00342                     double re = dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1];
00343                     double im = dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1];
00344                     dataC[j] = re / denom;
00345                     dataC[j+1] = im / denom;
00346                 }
00347         }
00348     }
00349 
00350 }
00351 
00352 static void fftShift(InputOutputArray _out)
00353 {
00354     Mat out = _out.getMat();
00355 
00356     if(out.rows == 1 && out.cols == 1)
00357     {
00358         // trivially shifted.
00359         return;
00360     }
00361 
00362     std::vector<Mat> planes;
00363     split(out, planes);
00364 
00365     int xMid = out.cols >> 1;
00366     int yMid = out.rows >> 1;
00367 
00368     bool is_1d = xMid == 0 || yMid == 0;
00369 
00370     if(is_1d)
00371     {
00372         xMid = xMid + yMid;
00373 
00374         for(size_t i = 0; i < planes.size(); i++)
00375         {
00376             Mat tmp;
00377             Mat half0(planes[i], Rect(0, 0, xMid, 1));
00378             Mat half1(planes[i], Rect(xMid, 0, xMid, 1));
00379 
00380             half0.copyTo(tmp);
00381             half1.copyTo(half0);
00382             tmp.copyTo(half1);
00383         }
00384     }
00385     else
00386     {
00387         for(size_t i = 0; i < planes.size(); i++)
00388         {
00389             // perform quadrant swaps...
00390             Mat tmp;
00391             Mat q0(planes[i], Rect(0,    0,    xMid, yMid));
00392             Mat q1(planes[i], Rect(xMid, 0,    xMid, yMid));
00393             Mat q2(planes[i], Rect(0,    yMid, xMid, yMid));
00394             Mat q3(planes[i], Rect(xMid, yMid, xMid, yMid));
00395 
00396             q0.copyTo(tmp);
00397             q3.copyTo(q0);
00398             tmp.copyTo(q3);
00399 
00400             q1.copyTo(tmp);
00401             q2.copyTo(q1);
00402             tmp.copyTo(q2);
00403         }
00404     }
00405 
00406     merge(planes, out);
00407 }
00408 
00409 static Point2d weightedCentroid(InputArray _src, cv::Point  peakLocation, cv::Size weightBoxSize, double* response)
00410 {
00411     Mat src = _src.getMat();
00412 
00413     int type = src.type();
00414     CV_Assert( type == CV_32FC1 || type == CV_64FC1 );
00415 
00416     int minr = peakLocation.y - (weightBoxSize.height >> 1);
00417     int maxr = peakLocation.y + (weightBoxSize.height >> 1);
00418     int minc = peakLocation.x - (weightBoxSize.width  >> 1);
00419     int maxc = peakLocation.x + (weightBoxSize.width  >> 1);
00420 
00421     Point2d centroid;
00422     double sumIntensity = 0.0;
00423 
00424     // clamp the values to min and max if needed.
00425     if(minr < 0)
00426     {
00427         minr = 0;
00428     }
00429 
00430     if(minc < 0)
00431     {
00432         minc = 0;
00433     }
00434 
00435     if(maxr > src.rows - 1)
00436     {
00437         maxr = src.rows - 1;
00438     }
00439 
00440     if(maxc > src.cols - 1)
00441     {
00442         maxc = src.cols - 1;
00443     }
00444 
00445     if(type == CV_32FC1)
00446     {
00447         const float* dataIn = src.ptr<float>();
00448         dataIn += minr*src.cols;
00449         for(int y = minr; y <= maxr; y++)
00450         {
00451             for(int x = minc; x <= maxc; x++)
00452             {
00453                 centroid.x   += (double)x*dataIn[x];
00454                 centroid.y   += (double)y*dataIn[x];
00455                 sumIntensity += (double)dataIn[x];
00456             }
00457 
00458             dataIn += src.cols;
00459         }
00460     }
00461     else
00462     {
00463         const double* dataIn = src.ptr<double>();
00464         dataIn += minr*src.cols;
00465         for(int y = minr; y <= maxr; y++)
00466         {
00467             for(int x = minc; x <= maxc; x++)
00468             {
00469                 centroid.x   += (double)x*dataIn[x];
00470                 centroid.y   += (double)y*dataIn[x];
00471                 sumIntensity += dataIn[x];
00472             }
00473 
00474             dataIn += src.cols;
00475         }
00476     }
00477 
00478     if(response)
00479         *response = sumIntensity;
00480 
00481     sumIntensity += DBL_EPSILON; // prevent div0 problems...
00482 
00483     centroid.x /= sumIntensity;
00484     centroid.y /= sumIntensity;
00485 
00486     return centroid;
00487 }
00488 
00489 }
00490 
00491 cv::Point2d cv::phaseCorrelate(InputArray _src1, InputArray _src2, InputArray _window, double* response)
00492 {
00493     Mat src1 = _src1.getMat();
00494     Mat src2 = _src2.getMat();
00495     Mat window = _window.getMat();
00496 
00497     CV_Assert( src1.type() == src2.type());
00498     CV_Assert( src1.type() == CV_32FC1 || src1.type() == CV_64FC1 );
00499     CV_Assert( src1.size == src2.size);
00500 
00501     if(!window.empty())
00502     {
00503         CV_Assert( src1.type() == window.type());
00504         CV_Assert( src1.size == window.size);
00505     }
00506 
00507     int M = getOptimalDFTSize(src1.rows);
00508     int N = getOptimalDFTSize(src1.cols);
00509 
00510     Mat padded1, padded2, paddedWin;
00511 
00512     if(M != src1.rows || N != src1.cols)
00513     {
00514         copyMakeBorder(src1, padded1, 0, M - src1.rows, 0, N - src1.cols, BORDER_CONSTANT, Scalar::all(0));
00515         copyMakeBorder(src2, padded2, 0, M - src2.rows, 0, N - src2.cols, BORDER_CONSTANT, Scalar::all(0));
00516 
00517         if(!window.empty())
00518         {
00519             copyMakeBorder(window, paddedWin, 0, M - window.rows, 0, N - window.cols, BORDER_CONSTANT, Scalar::all(0));
00520         }
00521     }
00522     else
00523     {
00524         padded1 = src1;
00525         padded2 = src2;
00526         paddedWin = window;
00527     }
00528 
00529     Mat FFT1, FFT2, P, Pm, C;
00530 
00531     // perform window multiplication if available
00532     if(!paddedWin.empty())
00533     {
00534         // apply window to both images before proceeding...
00535         multiply(paddedWin, padded1, padded1);
00536         multiply(paddedWin, padded2, padded2);
00537     }
00538 
00539     // execute phase correlation equation
00540     // Reference: http://en.wikipedia.org/wiki/Phase_correlation
00541     dft(padded1, FFT1, DFT_REAL_OUTPUT);
00542     dft(padded2, FFT2, DFT_REAL_OUTPUT);
00543 
00544     mulSpectrums(FFT1, FFT2, P, 0, true);
00545 
00546     magSpectrums(P, Pm);
00547     divSpectrums(P, Pm, C, 0, false); // FF* / |FF*| (phase correlation equation completed here...)
00548 
00549     idft(C, C); // gives us the nice peak shift location...
00550 
00551     fftShift(C); // shift the energy to the center of the frame.
00552 
00553     // locate the highest peak
00554     Point peakLoc;
00555     minMaxLoc(C, NULL, NULL, NULL, &peakLoc);
00556 
00557     // get the phase shift with sub-pixel accuracy, 5x5 window seems about right here...
00558     Point2d t;
00559     t = weightedCentroid(C, peakLoc, Size(5, 5), response);
00560 
00561     // max response is M*N (not exactly, might be slightly larger due to rounding errors)
00562     if(response)
00563         *response /= M*N;
00564 
00565     // adjust shift relative to image center...
00566     Point2d center((double)padded1.cols / 2.0, (double)padded1.rows / 2.0);
00567 
00568     return (center - t);
00569 }
00570 
00571 
00572 void cv::createHanningWindow(OutputArray _dst, cv::Size winSize, int type)
00573 {
00574     CV_Assert( type == CV_32FC1 || type == CV_64FC1 );
00575 
00576     _dst.create(winSize, type);
00577     Mat dst = _dst.getMat();
00578 
00579     int rows = dst.rows, cols = dst.cols;
00580 
00581     AutoBuffer<double> _wc(cols);
00582     double * const wc = (double *)_wc;
00583 
00584     double coeff0 = 2.0 * CV_PI / (double)(cols - 1), coeff1 = 2.0f * CV_PI / (double)(rows - 1);
00585     for(int j = 0; j < cols; j++)
00586         wc[j] = 0.5 * (1.0 - cos(coeff0 * j));
00587 
00588     if(dst.depth() == CV_32F)
00589     {
00590         for(int i = 0; i < rows; i++)
00591         {
00592             float* dstData = dst.ptr<float>(i);
00593             double wr = 0.5 * (1.0 - cos(coeff1 * i));
00594             for(int j = 0; j < cols; j++)
00595                 dstData[j] = (float)(wr * wc[j]);
00596         }
00597     }
00598     else
00599     {
00600         for(int i = 0; i < rows; i++)
00601         {
00602             double* dstData = dst.ptr<double>(i);
00603             double wr = 0.5 * (1.0 - cos(coeff1 * i));
00604             for(int j = 0; j < cols; j++)
00605                 dstData[j] = wr * wc[j];
00606         }
00607     }
00608 
00609     // perform batch sqrt for SSE performance gains
00610     cv::sqrt(dst, dst);
00611 }
00612