modules/contrib/src/logpolar_bsm.cpp

   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) 2012, Willow Garage Inc., all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The names of the copyright holders may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 /*******************************************************************************************
  43
  44 The LogPolar Blind Spot Model code has been contributed by Fabio Solari and Manuela Chessa.
  45
  46 More details can be found in:
  47
  48 M. Chessa, S. P. Sabatini, F. Solari and F. Tatti (2011)
  49 A Quantitative Comparison of Speed and Reliability for Log-Polar Mapping Techniques,
  50 Computer Vision Systems - 8th International Conference,
  51 ICVS 2011, Sophia Antipolis, France, September 20-22, 2011
  52 (http://dx.doi.org/10.1007/978-3-642-23968-7_5)
  53
  54 ********************************************************************************************/
  55
  56 #include "precomp.hpp"
  57
  58 #include <cmath>
  59 #include <vector>
  60
  61 namespace cv
  62 {
  63
  64 //------------------------------------interp-------------------------------------------
  65 LogPolar_Interp::LogPolar_Interp(int w, int h, Point2i center, int _R, double _ro0, int _interp, int full, int _s, int sp)
  66 {
  67     if ( (center.x!=w/2 || center.y!=h/2) && full==0) full=1;
  68
  69     if (center.x<0) center.x=0;
  70     if (center.y<0) center.y=0;
  71     if (center.x>=w) center.x=w-1;
  72     if (center.y>=h) center.y=h-1;
  73
  74     if (full){
  75         int rtmp;
  76
  77         if (center.x<=w/2 && center.y>=h/2)
  78             rtmp=(int)std::sqrt((float)center.y*center.y + (float)(w-center.x)*(w-center.x));
  79         else if (center.x>=w/2 && center.y>=h/2)
  80             rtmp=(int)std::sqrt((float)center.y*center.y + (float)center.x*center.x);
  81         else if (center.x>=w/2 && center.y<=h/2)
  82             rtmp=(int)std::sqrt((float)(h-center.y)*(h-center.y) + (float)center.x*center.x);
  83         else //if (center.x<=w/2 && center.y<=h/2)
  84             rtmp=(int)std::sqrt((float)(h-center.y)*(h-center.y) + (float)(w-center.x)*(w-center.x));
  85
  86         M=2*rtmp; N=2*rtmp;
  87
  88         top = M/2 - center.y;
  89         bottom = M/2 - (h-center.y);
  90         left = M/2 - center.x;
  91         right = M/2 - (w - center.x);
  92
  93     }else{
  94         top=bottom=left=right=0;
  95         M=w; N=h;
  96     }
  97
  98     if (sp){
  99         int jc=M/2-1, ic=N/2-1;
 100         int _romax=std::min(ic, jc);
 101         double _a=std::exp(std::log((double)(_romax/2-1)/(double)ro0)/(double)R);
 102         S=(int) floor(2*CV_PI/(_a-1)+0.5);
 103     }
 104
 105     interp=_interp;
 106
 107     create_map(M, N, _R, _s, _ro0);
 108 }
 109
 110 void LogPolar_Interp::create_map(int _M, int _n, int _R, int _s, double _ro0)
 111 {
 112     M=_M;
 113     N=_n;
 114     R=_R;
 115     S=_s;
 116     ro0=_ro0;
 117
 118     int jc=N/2-1, ic=M/2-1;
 119     romax=std::min(ic, jc);
 120     a=std::exp(std::log((double)romax/(double)ro0)/(double)R);
 121     q=((double)S)/(2*CV_PI);
 122
 123     Rsri = Mat::zeros(S,R,CV_32FC1);
 124     Csri = Mat::zeros(S,R,CV_32FC1);
 125     ETAyx = Mat::zeros(N,M,CV_32FC1);
 126     CSIyx = Mat::zeros(N,M,CV_32FC1);
 127
 128     for(int v=0; v<S; v++)
 129     {
 130         for(int u=0; u<R; u++)
 131         {
 132             Rsri.at<float>(v,u)=(float)(ro0*std::pow(a,u)*sin(v/q)+jc);
 133             Csri.at<float>(v,u)=(float)(ro0*std::pow(a,u)*cos(v/q)+ic);
 134         }
 135     }
 136
 137     for(int j=0; j<N; j++)
 138     {
 139         for(int i=0; i<M; i++)
 140         {
 141             double theta;
 142             if(i>=ic)
 143                 theta=atan((double)(j-jc)/(double)(i-ic));
 144             else
 145                 theta=atan((double)(j-jc)/(double)(i-ic))+CV_PI;
 146
 147             if(theta<0)
 148                 theta+=2*CV_PI;
 149
 150             ETAyx.at<float>(j,i)=(float)(q*theta);
 151
 152             double ro2=(j-jc)*(j-jc)+(i-ic)*(i-ic);
 153             CSIyx.at<float>(j,i)=(float)(0.5*std::log(ro2/(ro0*ro0))/std::log(a));
 154         }
 155     }
 156 }
 157
 158 const Mat LogPolar_Interp::to_cortical(const Mat &source)
 159 {
 160     Mat out(S,R,CV_8UC1,Scalar(0));
 161
 162     Mat source_border;
 163     copyMakeBorder(source,source_border,top,bottom,left,right,BORDER_CONSTANT,Scalar(0));
 164
 165     remap(source_border,out,Csri,Rsri,interp);
 166
 167     return out;
 168 }
 169
 170
 171 const Mat LogPolar_Interp::to_cartesian(const Mat &source)
 172 {
 173     Mat out(N,M,CV_8UC1,Scalar(0));
 174
 175     Mat source_border;
 176
 177     if (interp==INTER_NEAREST || interp==INTER_LINEAR){
 178         copyMakeBorder(source,source_border,0,1,0,0,BORDER_CONSTANT,Scalar(0));
 179         Mat rowS0 = source_border.row(S);
 180         source_border.row(0).copyTo(rowS0);
 181     } else if (interp==INTER_CUBIC){
 182         copyMakeBorder(source,source_border,0,2,0,0,BORDER_CONSTANT,Scalar(0));
 183         Mat rowS0 = source_border.row(S);
 184         Mat rowS1 = source_border.row(S+1);
 185         source_border.row(0).copyTo(rowS0);
 186         source_border.row(1).copyTo(rowS1);
 187     } else if (interp==INTER_LANCZOS4){
 188         copyMakeBorder(source,source_border,0,4,0,0,BORDER_CONSTANT,Scalar(0));
 189         Mat rowS0 = source_border.row(S);
 190         Mat rowS1 = source_border.row(S+1);
 191         Mat rowS2 = source_border.row(S+2);
 192         Mat rowS3 = source_border.row(S+3);
 193         source_border.row(0).copyTo(rowS0);
 194         source_border.row(1).copyTo(rowS1);
 195         source_border.row(2).copyTo(rowS2);
 196         source_border.row(3).copyTo(rowS3);
 197     }
 198     remap(source_border,out,CSIyx,ETAyx,interp);
 199
 200     Mat out_cropped=out(Range(top,N-1-bottom),Range(left,M-1-right));
 201
 202     return out_cropped;
 203 }
 204
 205 LogPolar_Interp::~LogPolar_Interp()
 206 {
 207 }
 208
 209 //------------------------------------overlapping----------------------------------
 210
 211 LogPolar_Overlapping::LogPolar_Overlapping(int w, int h, Point2i center, int _R, double _ro0, int full, int _s, int sp)
 212 {
 213     if ( (center.x!=w/2 || center.y!=h/2) && full==0) full=1;
 214
 215     if (center.x<0) center.x=0;
 216     if (center.y<0) center.y=0;
 217     if (center.x>=w) center.x=w-1;
 218     if (center.y>=h) center.y=h-1;
 219
 220     if (full){
 221         int rtmp;
 222
 223         if (center.x<=w/2 && center.y>=h/2)
 224             rtmp=(int)std::sqrt((float)center.y*center.y + (float)(w-center.x)*(w-center.x));
 225         else if (center.x>=w/2 && center.y>=h/2)
 226             rtmp=(int)std::sqrt((float)center.y*center.y + (float)center.x*center.x);
 227         else if (center.x>=w/2 && center.y<=h/2)
 228             rtmp=(int)std::sqrt((float)(h-center.y)*(h-center.y) + (float)center.x*center.x);
 229         else //if (center.x<=w/2 && center.y<=h/2)
 230             rtmp=(int)std::sqrt((float)(h-center.y)*(h-center.y) + (float)(w-center.x)*(w-center.x));
 231
 232         M=2*rtmp; N=2*rtmp;
 233
 234         top = M/2 - center.y;
 235         bottom = M/2 - (h-center.y);
 236         left = M/2 - center.x;
 237         right = M/2 - (w - center.x);
 238
 239     }else{
 240         top=bottom=left=right=0;
 241         M=w; N=h;
 242     }
 243
 244
 245     if (sp){
 246         int jc=M/2-1, ic=N/2-1;
 247         int _romax=std::min(ic, jc);
 248         double _a=std::exp(std::log((double)(_romax/2-1)/(double)ro0)/(double)R);
 249         S=(int) floor(2*CV_PI/(_a-1)+0.5);
 250     }
 251
 252     create_map(M, N, _R, _s, _ro0);
 253 }
 254
 255 void LogPolar_Overlapping::create_map(int _M, int _n, int _R, int _s, double _ro0)
 256 {
 257     M=_M;
 258     N=_n;
 259     R=_R;
 260     S=_s;
 261     ro0=_ro0;
 262
 263     int jc=N/2-1, ic=M/2-1;
 264     romax=std::min(ic, jc);
 265     a=std::exp(std::log((double)romax/(double)ro0)/(double)R);
 266     q=((double)S)/(2*CV_PI);
 267     ind1=0;
 268
 269     Rsri=Mat::zeros(S,R,CV_32FC1);
 270     Csri=Mat::zeros(S,R,CV_32FC1);
 271     ETAyx=Mat::zeros(N,M,CV_32FC1);
 272     CSIyx=Mat::zeros(N,M,CV_32FC1);
 273     Rsr.resize(R*S);
 274     Csr.resize(R*S);
 275     Wsr.resize(R);
 276     w_ker_2D.resize(R*S);
 277
 278     for(int v=0; v<S; v++)
 279     {
 280         for(int u=0; u<R; u++)
 281         {
 282             Rsri.at<float>(v,u)=(float)(ro0*std::pow(a,u)*sin(v/q)+jc);
 283             Csri.at<float>(v,u)=(float)(ro0*std::pow(a,u)*cos(v/q)+ic);
 284             Rsr[v*R+u]=(int)floor(Rsri.at<float>(v,u));
 285             Csr[v*R+u]=(int)floor(Csri.at<float>(v,u));
 286         }
 287     }
 288
 289     bool done=false;
 290
 291     for(int i=0; i<R; i++)
 292     {
 293         Wsr[i]=ro0*(a-1)*std::pow(a,i-1);
 294         if((Wsr[i]>1)&&(done==false))
 295         {
 296             ind1=i;
 297             done =true;
 298         }
 299     }
 300
 301     for(int j=0; j<N; j++)
 302     {
 303         for(int i=0; i<M; i++)//mdf
 304         {
 305             double theta;
 306             if(i>=ic)
 307                 theta=atan((double)(j-jc)/(double)(i-ic));
 308             else
 309                 theta=atan((double)(j-jc)/(double)(i-ic))+CV_PI;
 310
 311             if(theta<0)
 312                 theta+=2*CV_PI;
 313
 314             ETAyx.at<float>(j,i)=(float)(q*theta);
 315
 316             double ro2=(j-jc)*(j-jc)+(i-ic)*(i-ic);
 317             CSIyx.at<float>(j,i)=(float)(0.5*std::log(ro2/(ro0*ro0))/std::log(a));
 318         }
 319     }
 320
 321     for(int v=0; v<S; v++)
 322         for(int u=ind1; u<R; u++)
 323         {
 324             //double sigma=Wsr[u]/2.0;
 325             double sigma=Wsr[u]/3.0;//modf
 326             int w=(int) floor(3*sigma+0.5);
 327             w_ker_2D[v*R+u].w=w;
 328             w_ker_2D[v*R+u].weights.resize((2*w+1)*(2*w+1));
 329             double dx=Csri.at<float>(v,u)-Csr[v*R+u];
 330             double dy=Rsri.at<float>(v,u)-Rsr[v*R+u];
 331             double tot=0;
 332             for(int j=0; j<2*w+1; j++)
 333                 for(int i=0; i<2*w+1; i++)
 334                 {
 335                     (w_ker_2D[v*R+u].weights)[j*(2*w+1)+i]=std::exp(-(std::pow(i-w-dx, 2)+std::pow(j-w-dy, 2))/(2*sigma*sigma));
 336                     tot+=(w_ker_2D[v*R+u].weights)[j*(2*w+1)+i];
 337                 }
 338             for(int j=0; j<(2*w+1); j++)
 339                 for(int i=0; i<(2*w+1); i++)
 340                     (w_ker_2D[v*R+u].weights)[j*(2*w+1)+i]/=tot;
 341         }
 342 }
 343
 344 const Mat LogPolar_Overlapping::to_cortical(const Mat &source)
 345 {
 346     Mat out(S,R,CV_8UC1,Scalar(0));
 347
 348     Mat source_border;
 349     copyMakeBorder(source,source_border,top,bottom,left,right,BORDER_CONSTANT,Scalar(0));
 350
 351     remap(source_border,out,Csri,Rsri,INTER_LINEAR);
 352
 353     int wm=w_ker_2D[R-1].w;
 354     std::vector<int> IMG((M+2*wm+1)*(N+2*wm+1), 0);
 355
 356     for(int j=0; j<N; j++)
 357         for(int i=0; i<M; i++)
 358             IMG[(M+2*wm+1)*(j+wm)+i+wm]=source_border.at<uchar>(j,i);
 359
 360     for(int v=0; v<S; v++)
 361         for(int u=ind1; u<R; u++)
 362         {
 363             int w=w_ker_2D[v*R+u].w;
 364             double tmp=0;
 365             for(int rf=0; rf<(2*w+1); rf++)
 366             {
 367                 for(int cf=0; cf<(2*w+1); cf++)
 368                 {
 369                     double weight=(w_ker_2D[v*R+u]).weights[rf*(2*w+1)+cf];
 370                     tmp+=IMG[(M+2*wm+1)*((rf-w)+Rsr[v*R+u]+wm)+((cf-w)+Csr[v*R+u]+wm)]*weight;
 371                 }
 372             }
 373             out.at<uchar>(v,u)=(uchar) floor(tmp+0.5);
 374         }
 375
 376     return out;
 377 }
 378
 379 const Mat LogPolar_Overlapping::to_cartesian(const Mat &source)
 380 {
 381     Mat out(N,M,CV_8UC1,Scalar(0));
 382
 383     Mat source_border;
 384     copyMakeBorder(source,source_border,0,1,0,0,BORDER_CONSTANT,Scalar(0));
 385     Mat rowS = source_border.row(S);
 386     source_border.row(0).copyTo(rowS);
 387     remap(source_border,out,CSIyx,ETAyx,INTER_LINEAR);
 388
 389     int wm=w_ker_2D[R-1].w;
 390
 391     std::vector<double> IMG((N+2*wm+1)*(M+2*wm+1), 0.);
 392     std::vector<double> NOR((N+2*wm+1)*(M+2*wm+1), 0.);
 393
 394     for(int v=0; v<S; v++)
 395         for(int u=ind1; u<R; u++)
 396         {
 397             int w=w_ker_2D[v*R+u].w;
 398             for(int j=0; j<(2*w+1); j++)
 399             {
 400                 for(int i=0; i<(2*w+1); i++)
 401                 {
 402                     int ind=(M+2*wm+1)*((j-w)+Rsr[v*R+u]+wm)+(i-w)+Csr[v*R+u]+wm;
 403                     IMG[ind]+=((w_ker_2D[v*R+u]).weights[j*(2*w+1)+i])*source.at<uchar>(v, u);
 404                     NOR[ind]+=((w_ker_2D[v*R+u]).weights[j*(2*w+1)+i]);
 405                 }
 406             }
 407         }
 408
 409     for(int i=0; i<((N+2*wm+1)*(M+2*wm+1)); i++)
 410         IMG[i]/=NOR[i];
 411
 412     //int xc=M/2-1, yc=N/2-1;
 413
 414     for(int j=wm; j<N+wm; j++)
 415         for(int i=wm; i<M+wm; i++)
 416         {
 417             /*if(NOR[(M+2*wm+1)*j+i]>0)
 418                 ret[M*(j-wm)+i-wm]=(int) floor(IMG[(M+2*wm+1)*j+i]+0.5);*/
 419             //int ro=(int)floor(std::sqrt((double)((j-wm-yc)*(j-wm-yc)+(i-wm-xc)*(i-wm-xc))));
 420             int csi=(int) floor(CSIyx.at<float>(j-wm,i-wm));
 421
 422             if((csi>=(ind1-(w_ker_2D[ind1]).w))&&(csi<R))
 423                 out.at<uchar>(j-wm,i-wm)=(uchar) floor(IMG[(M+2*wm+1)*j+i]+0.5);
 424         }
 425
 426     Mat out_cropped=out(Range(top,N-1-bottom),Range(left,M-1-right));
 427     return out_cropped;
 428 }
 429
 430 LogPolar_Overlapping::~LogPolar_Overlapping()
 431 {
 432 }
 433
 434 //----------------------------------------adjacent---------------------------------------
 435
 436 LogPolar_Adjacent::LogPolar_Adjacent(int w, int h, Point2i center, int _R, double _ro0, double smin, int full, int _s, int sp)
 437 {
 438     if ( (center.x!=w/2 || center.y!=h/2) && full==0) full=1;
 439
 440     if (center.x<0) center.x=0;
 441     if (center.y<0) center.y=0;
 442     if (center.x>=w) center.x=w-1;
 443     if (center.y>=h) center.y=h-1;
 444
 445     if (full){
 446         int rtmp;
 447
 448         if (center.x<=w/2 && center.y>=h/2)
 449             rtmp=(int)std::sqrt((float)center.y*center.y + (float)(w-center.x)*(w-center.x));
 450         else if (center.x>=w/2 && center.y>=h/2)
 451             rtmp=(int)std::sqrt((float)center.y*center.y + (float)center.x*center.x);
 452         else if (center.x>=w/2 && center.y<=h/2)
 453             rtmp=(int)std::sqrt((float)(h-center.y)*(h-center.y) + (float)center.x*center.x);
 454         else //if (center.x<=w/2 && center.y<=h/2)
 455             rtmp=(int)std::sqrt((float)(h-center.y)*(h-center.y) + (float)(w-center.x)*(w-center.x));
 456
 457         M=2*rtmp; N=2*rtmp;
 458
 459         top = M/2 - center.y;
 460         bottom = M/2 - (h-center.y);
 461         left = M/2 - center.x;
 462         right = M/2 - (w - center.x);
 463
 464     }else{
 465         top=bottom=left=right=0;
 466         M=w; N=h;
 467     }
 468
 469     if (sp){
 470         int jc=M/2-1, ic=N/2-1;
 471         int _romax=std::min(ic, jc);
 472         double _a=std::exp(std::log((double)(_romax/2-1)/(double)ro0)/(double)R);
 473         S=(int) floor(2*CV_PI/(_a-1)+0.5);
 474     }
 475
 476     create_map(M, N, _R, _s, _ro0, smin);
 477 }
 478
 479
 480 void LogPolar_Adjacent::create_map(int _M, int _n, int _R, int _s, double _ro0, double smin)
 481 {
 482     M=_M;
 483     N=_n;
 484     R=_R;
 485     S=_s;
 486     ro0=_ro0;
 487     romax=std::min(M/2.0, N/2.0);
 488
 489     a=std::exp(std::log(romax/ro0)/(double)R);
 490     q=S/(2*CV_PI);
 491
 492     A.resize(R*S);
 493     L.resize(M*N);
 494
 495     for(int i=0; i<R*S; i++)
 496         A[i]=0;
 497
 498     double xc=M/2.0, yc=N/2.0;
 499
 500     for(int j=0; j<N; j++)
 501         for(int i=0; i<M; i++)
 502         {
 503             double x=i+0.5-xc, y=j+0.5-yc;
 504             subdivide_recursively(x, y, i, j, 1, smin);
 505         }
 506 }
 507
 508
 509 void LogPolar_Adjacent::subdivide_recursively(double x, double y, int i, int j, double length, double smin)
 510 {
 511     if(length<=smin)
 512     {
 513         int u, v;
 514         if(get_uv(x, y, u, v))
 515         {
 516             pixel p;
 517             p.u=u;
 518             p.v=v;
 519             p.a=length*length;
 520             L[M*j+i].push_back(p);
 521             A[v*R+u]+=length*length;
 522         }
 523     }
 524
 525     if(length>smin)
 526     {
 527         double xs[4], ys[4];
 528         int us[4], vs[4];
 529
 530         xs[0]=xs[3]=x+length/4.0;
 531         xs[1]=xs[2]=x-length/4.0;
 532         ys[1]=ys[0]=y+length/4.0;
 533         ys[2]=ys[3]=y-length/4.0;
 534
 535         for(int z=0; z<4; z++)
 536             get_uv(xs[z], ys[z], us[z], vs[z]);
 537
 538         bool c=true;
 539         for(int w=1; w<4; w++)
 540         {
 541             if(us[w]!=us[w-1])
 542                 c=false;
 543             if(vs[w]!=vs[w-1])
 544                 c=false;
 545         }
 546
 547         if(c)
 548         {
 549             if(us[0]!=-1)
 550             {
 551                 pixel p;
 552                 p.u=us[0];
 553                 p.v=vs[0];
 554                 p.a=length*length;
 555                 L[M*j+i].push_back(p);
 556                 A[vs[0]*R+us[0]]+=length*length;
 557             }
 558         }
 559
 560         else
 561         {
 562             for(int z=0; z<4; z++)
 563                 if(us[z]!=-1)
 564                     subdivide_recursively(xs[z], ys[z], i, j, length/2.0, smin);
 565         }
 566     }
 567 }
 568
 569
 570 const Mat LogPolar_Adjacent::to_cortical(const Mat &source)
 571 {
 572     Mat source_border;
 573     copyMakeBorder(source,source_border,top,bottom,left,right,BORDER_CONSTANT,Scalar(0));
 574
 575     std::vector<double> map(R*S, 0.);
 576
 577     for(int j=0; j<N; j++)
 578         for(int i=0; i<M; i++)
 579         {
 580             for(size_t z=0; z<(L[M*j+i]).size(); z++)
 581             {
 582                 map[R*((L[M*j+i])[z].v)+((L[M*j+i])[z].u)]+=((L[M*j+i])[z].a)*(source_border.at<uchar>(j,i));
 583             }
 584         }
 585
 586     for(int i=0; i<R*S; i++)
 587         map[i]/=A[i];
 588
 589     Mat out(S,R,CV_8UC1,Scalar(0));
 590
 591     for(int i=0; i<S; i++)
 592         for(int j=0;j<R;j++)
 593             out.at<uchar>(i,j)=(uchar) floor(map[i*R+j]+0.5);
 594
 595     return out;
 596 }
 597
 598 const Mat LogPolar_Adjacent::to_cartesian(const Mat &source)
 599 {
 600     std::vector<double> map(M*N, 0.);
 601
 602     for(int j=0; j<N; j++)
 603         for(int i=0; i<M; i++)
 604         {
 605             for(size_t z=0; z<(L[M*j+i]).size(); z++)
 606             {
 607                 map[M*j+i]+=(L[M*j+i])[z].a*source.at<uchar>((L[M*j+i])[z].v,(L[M*j+i])[z].u);
 608             }
 609         }
 610
 611     Mat out(N,M,CV_8UC1,Scalar(0));
 612
 613     for(int i=0; i<N; i++)
 614         for(int j=0; j<M; j++)
 615             out.at<uchar>(i,j)=(uchar) floor(map[i*M+j]+0.5);
 616
 617     Mat out_cropped=out(Range(top,N-1-bottom),Range(left,M-1-right));
 618     return out_cropped;
 619 }
 620
 621
 622 bool LogPolar_Adjacent::get_uv(double x, double y, int&u, int&v)
 623 {
 624     double ro=std::sqrt(x*x+y*y), theta;
 625     if(x>0)
 626         theta=atan(y/x);
 627     else
 628         theta=atan(y/x)+CV_PI;
 629
 630     if(ro<ro0||ro>romax)
 631     {
 632         u=-1;
 633         v=-1;
 634         return false;
 635     }
 636     else
 637     {
 638         u= (int) floor(std::log(ro/ro0)/std::log(a));
 639         if(theta>=0)
 640             v= (int) floor(q*theta);
 641         else
 642             v= (int) floor(q*(theta+2*CV_PI));
 643         return true;
 644     }
 645 }
 646
 647 LogPolar_Adjacent::~LogPolar_Adjacent()
 648 {
 649 }
 650
 651 }