From c8e206c2abca4c975442860f61cdf2500764036d Mon Sep 17 00:00:00 2001 From: Vadim Pisarevsky Date: Sun, 18 Mar 2012 22:29:13 +0000 Subject: [PATCH] added LogPolar Blind Spot Model (thanks to Fabio Solari for the contribution) --- .../include/opencv2/contrib/contrib.hpp | 214 +++++- modules/contrib/src/logpolar_bsm.cpp | 652 ++++++++++++++++++ samples/cpp/logpolar_bsm.cpp | 82 +++ 3 files changed, 944 insertions(+), 4 deletions(-) create mode 100644 modules/contrib/src/logpolar_bsm.cpp create mode 100644 samples/cpp/logpolar_bsm.cpp diff --git a/modules/contrib/include/opencv2/contrib/contrib.hpp b/modules/contrib/include/opencv2/contrib/contrib.hpp index eaab0a023f..f27c0c1801 100644 --- a/modules/contrib/include/opencv2/contrib/contrib.hpp +++ b/modules/contrib/include/opencv2/contrib/contrib.hpp @@ -44,6 +44,7 @@ #define __OPENCV_CONTRIB_HPP__ #include "opencv2/core/core.hpp" +#include "opencv2/imgproc/imgproc.hpp" #include "opencv2/features2d/features2d.hpp" #include "opencv2/objdetect/objdetect.hpp" @@ -633,14 +634,219 @@ namespace cv TRANSLATION = 2, RIGID_BODY_MOTION = 4 }; - CV_EXPORTS bool RGBDOdometry( cv::Mat& Rt, const Mat& initRt, - const cv::Mat& image0, const cv::Mat& depth0, const cv::Mat& mask0, - const cv::Mat& image1, const cv::Mat& depth1, const cv::Mat& mask1, - const cv::Mat& cameraMatrix, float minDepth, float maxDepth, float maxDepthDiff, + CV_EXPORTS bool RGBDOdometry( Mat& Rt, const Mat& initRt, + const Mat& image0, const Mat& depth0, const Mat& mask0, + const Mat& image1, const Mat& depth1, const Mat& mask1, + const Mat& cameraMatrix, float minDepth, float maxDepth, float maxDepthDiff, const std::vector& iterCounts, const std::vector& minGradientMagnitudes, int transformType=RIGID_BODY_MOTION ); + + /** + *Bilinear interpolation technique. + * + *The value of a desired cortical pixel is obtained through a bilinear interpolation of the values + *of the four nearest neighbouring Cartesian pixels to the center of the RF. + *The same principle is applied to the inverse transformation. + * + *More details can be found in http://dx.doi.org/10.1007/978-3-642-23968-7_5 + */ + class CV_EXPORTS LogPolar_Interp + { + public: + + LogPolar_Interp() {} + + /** + *Constructor + *\param w the width of the input image + *\param h the height of the input image + *\param center the transformation center: where the output precision is maximal + *\param R the number of rings of the cortical image (default value 70 pixel) + *\param ro0 the radius of the blind spot (default value 3 pixel) + *\param full \a 1 (default value) means that the retinal image (the inverse transform) is computed within the circumscribing circle. + * \a 0 means that the retinal image is computed within the inscribed circle. + *\param S the number of sectors of the cortical image (default value 70 pixel). + * Its value is usually internally computed to obtain a pixel aspect ratio equals to 1. + *\param sp \a 1 (default value) means that the parameter \a S is internally computed. + * \a 0 means that the parameter \a S is provided by the user. + */ + LogPolar_Interp(int w, int h, Point2i center, int R=70, double ro0=3.0, + int interp=INTER_LINEAR, int full=1, int S=117, int sp=1); + /** + *Transformation from Cartesian image to cortical (log-polar) image. + *\param source the Cartesian image + *\return the transformed image (cortical image) + */ + const Mat to_cortical(const Mat &source); + /** + *Transformation from cortical image to retinal (inverse log-polar) image. + *\param source the cortical image + *\return the transformed image (retinal image) + */ + const Mat to_cartesian(const Mat &source); + /** + *Destructor + */ + ~LogPolar_Interp(); + + protected: + + Mat Rsri; + Mat Csri; + + int S, R, M, N; + int top, bottom,left,right; + double ro0, romax, a, q; + int interp; + + Mat ETAyx; + Mat CSIyx; + + void create_map(int M, int N, int R, int S, double ro0); + }; + + /** + *Overlapping circular receptive fields technique + * + *The Cartesian plane is divided in two regions: the fovea and the periphery. + *The fovea (oversampling) is handled by using the bilinear interpolation technique described above, whereas in + *the periphery we use the overlapping Gaussian circular RFs. + * + *More details can be found in http://dx.doi.org/10.1007/978-3-642-23968-7_5 + */ + class CV_EXPORTS LogPolar_Overlapping + { + public: + LogPolar_Overlapping() {} + + /** + *Constructor + *\param w the width of the input image + *\param h the height of the input image + *\param center the transformation center: where the output precision is maximal + *\param R the number of rings of the cortical image (default value 70 pixel) + *\param ro0 the radius of the blind spot (default value 3 pixel) + *\param full \a 1 (default value) means that the retinal image (the inverse transform) is computed within the circumscribing circle. + * \a 0 means that the retinal image is computed within the inscribed circle. + *\param S the number of sectors of the cortical image (default value 70 pixel). + * Its value is usually internally computed to obtain a pixel aspect ratio equals to 1. + *\param sp \a 1 (default value) means that the parameter \a S is internally computed. + * \a 0 means that the parameter \a S is provided by the user. + */ + LogPolar_Overlapping(int w, int h, Point2i center, int R=70, + double ro0=3.0, int full=1, int S=117, int sp=1); + /** + *Transformation from Cartesian image to cortical (log-polar) image. + *\param source the Cartesian image + *\return the transformed image (cortical image) + */ + const Mat to_cortical(const Mat &source); + /** + *Transformation from cortical image to retinal (inverse log-polar) image. + *\param source the cortical image + *\return the transformed image (retinal image) + */ + const Mat to_cartesian(const Mat &source); + /** + *Destructor + */ + ~LogPolar_Overlapping(); + + protected: + + Mat Rsri; + Mat Csri; + vector Rsr; + vector Csr; + vector Wsr; + + int S, R, M, N, ind1; + int top, bottom,left,right; + double ro0, romax, a, q; + + struct kernel + { + kernel() { w = 0; } + vector weights; + int w; + }; + + Mat ETAyx; + Mat CSIyx; + vector w_ker_2D; + + void create_map(int M, int N, int R, int S, double ro0); + }; + + /** + * Adjacent receptive fields technique + * + *All the Cartesian pixels, whose coordinates in the cortical domain share the same integer part, are assigned to the same RF. + *The precision of the boundaries of the RF can be improved by breaking each pixel into subpixels and assigning each of them to the correct RF. + *This technique is implemented from: Traver, V., Pla, F.: Log-polar mapping template design: From task-level requirements + *to geometry parameters. Image Vision Comput. 26(10) (2008) 1354-1370 + * + *More details can be found in http://dx.doi.org/10.1007/978-3-642-23968-7_5 + */ + class CV_EXPORTS LogPolar_Adjacent + { + public: + LogPolar_Adjacent() {} + + /** + *Constructor + *\param w the width of the input image + *\param h the height of the input image + *\param center the transformation center: where the output precision is maximal + *\param R the number of rings of the cortical image (default value 70 pixel) + *\param ro0 the radius of the blind spot (default value 3 pixel) + *\param smin the size of the subpixel (default value 0.25 pixel) + *\param full \a 1 (default value) means that the retinal image (the inverse transform) is computed within the circumscribing circle. + * \a 0 means that the retinal image is computed within the inscribed circle. + *\param S the number of sectors of the cortical image (default value 70 pixel). + * Its value is usually internally computed to obtain a pixel aspect ratio equals to 1. + *\param sp \a 1 (default value) means that the parameter \a S is internally computed. + * \a 0 means that the parameter \a S is provided by the user. + */ + LogPolar_Adjacent(int w, int h, Point2i center, int R=70, double ro0=3.0, double smin=0.25, int full=1, int S=117, int sp=1); + /** + *Transformation from Cartesian image to cortical (log-polar) image. + *\param source the Cartesian image + *\return the transformed image (cortical image) + */ + const Mat to_cortical(const Mat &source); + /** + *Transformation from cortical image to retinal (inverse log-polar) image. + *\param source the cortical image + *\return the transformed image (retinal image) + */ + const Mat to_cartesian(const Mat &source); + /** + *Destructor + */ + ~LogPolar_Adjacent(); + + protected: + struct pixel + { + pixel() { u = v = 0; a = 0.; } + int u; + int v; + double a; + }; + int S, R, M, N; + int top, bottom,left,right; + double ro0, romax, a, q; + vector > L; + vector A; + + void subdivide_recursively(double x, double y, int i, int j, double length, double smin); + bool get_uv(double x, double y, int&u, int&v); + void create_map(int M, int N, int R, int S, double ro0, double smin); + }; } + #include "opencv2/contrib/retina.hpp" #endif diff --git a/modules/contrib/src/logpolar_bsm.cpp b/modules/contrib/src/logpolar_bsm.cpp new file mode 100644 index 0000000000..b84e1a6df0 --- /dev/null +++ b/modules/contrib/src/logpolar_bsm.cpp @@ -0,0 +1,652 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2012, Willow Garage Inc., all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistribution's of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistribution's in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * The names of the copyright holders may not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall the Intel Corporation or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +/************************************************************************************* + +The LogPolar Blind Spot Model code has been contributed by Fabio Solari. + +More details can be found in: + +M. Chessa, S. P. Sabatini, F. Solari and F. Tatti (2011) +A Quantitative Comparison of Speed and Reliability for Log-Polar Mapping Techniques, +Computer Vision Systems - 8th International Conference, +ICVS 2011, Sophia Antipolis, France, September 20-22, 2011 +(http://dx.doi.org/10.1007/978-3-642-23968-7_5) + +***************************************************************************************/ + +#include "precomp.hpp" + +#include +#include + +namespace cv +{ + +//------------------------------------interp------------------------------------------- +LogPolar_Interp::LogPolar_Interp(int w, int h, Point2i center, int R, double ro0, int interp, int full, int S, int sp) +{ + if ( (center.x!=w/2 || center.y!=h/2) && full==0) full=1; + + if (center.x<0) center.x=0; + if (center.y<0) center.y=0; + if (center.x>=w) center.x=w-1; + if (center.y>=h) center.y=h-1; + + if (full){ + int rtmp; + + if (center.x<=w/2 && center.y>=h/2) + rtmp=(int)sqrt((float)center.y*center.y + (float)(w-center.x)*(w-center.x)); + if (center.x>=w/2 && center.y>=h/2) + rtmp=(int)sqrt((float)center.y*center.y + (float)center.x*center.x); + if (center.x>=w/2 && center.y<=h/2) + rtmp=(int)sqrt((float)(h-center.y)*(h-center.y) + (float)center.x*center.x); + if (center.x<=w/2 && center.y<=h/2) + rtmp=(int)sqrt((float)(h-center.y)*(h-center.y) + (float)(w-center.x)*(w-center.x)); + + M=2*rtmp; N=2*rtmp; + + top = M/2 - center.y; + bottom = M/2 - (h-center.y); + left = M/2 - center.x; + right = M/2 - (w - center.x); + + }else{ + top=bottom=left=right=0; + M=w; N=h; + } + + if (sp){ + int jc=M/2-1, ic=N/2-1; + int romax=min(ic, jc); + double a=exp(log((double)(romax/2-1)/(double)ro0)/(double)R); + S=(int) floor(2*M_PI/(a-1)+0.5); + } + + this->interp=interp; + + create_map(M, N, R, S, ro0); +} + +void LogPolar_Interp::create_map(int M, int N, int R, int S, double ro0) +{ + this->M=M; + this->N=N; + this->R=R; + this->S=S; + this->ro0=ro0; + + int jc=N/2-1, ic=M/2-1; + romax=min(ic, jc); + a=exp(log((double)romax/(double)ro0)/(double)R); + q=((double)S)/(2*M_PI); + + Rsri = Mat::zeros(S,R,CV_32FC1); + Csri = Mat::zeros(S,R,CV_32FC1); + ETAyx = Mat::zeros(N,M,CV_32FC1); + CSIyx = Mat::zeros(N,M,CV_32FC1); + + for(int v=0; v(v,u)=(float)(ro0*pow(a,u)*sin(v/q)+jc); + Csri.at(v,u)=(float)(ro0*pow(a,u)*cos(v/q)+ic); + } + } + + for(int j=0; j=ic) + theta=atan((double)(j-jc)/(double)(i-ic)); + else + theta=atan((double)(j-jc)/(double)(i-ic))+M_PI; + + if(theta<0) + theta+=2*M_PI; + + ETAyx.at(j,i)=(float)(q*theta); + + double ro2=(j-jc)*(j-jc)+(i-ic)*(i-ic); + CSIyx.at(j,i)=(float)(0.5*log(ro2/(ro0*ro0))/log(a)); + } + } +} + +const Mat LogPolar_Interp::to_cortical(const Mat &source) +{ + Mat out(S,R,CV_8UC1,Scalar(0)); + + Mat source_border; + copyMakeBorder(source,source_border,top,bottom,left,right,BORDER_CONSTANT,Scalar(0)); + + remap(source_border,out,Csri,Rsri,interp); + + return out; +} + + +const Mat LogPolar_Interp::to_cartesian(const Mat &source) +{ + Mat out(N,M,CV_8UC1,Scalar(0)); + + Mat source_border; + + if (interp==INTER_NEAREST || interp==INTER_LINEAR){ + copyMakeBorder(source,source_border,0,1,0,0,BORDER_CONSTANT,Scalar(0)); + Mat rowS0 = source_border.row(S); + source_border.row(0).copyTo(rowS0); + } else if (interp==INTER_CUBIC){ + copyMakeBorder(source,source_border,0,2,0,0,BORDER_CONSTANT,Scalar(0)); + Mat rowS0 = source_border.row(S); + Mat rowS1 = source_border.row(S+1); + source_border.row(0).copyTo(rowS0); + source_border.row(1).copyTo(rowS1); + } else if (interp==INTER_LANCZOS4){ + copyMakeBorder(source,source_border,0,4,0,0,BORDER_CONSTANT,Scalar(0)); + Mat rowS0 = source_border.row(S); + Mat rowS1 = source_border.row(S+1); + Mat rowS2 = source_border.row(S+2); + Mat rowS3 = source_border.row(S+3); + source_border.row(0).copyTo(rowS0); + source_border.row(1).copyTo(rowS1); + source_border.row(2).copyTo(rowS2); + source_border.row(3).copyTo(rowS3); + } + remap(source_border,out,CSIyx,ETAyx,interp); + + Mat out_cropped=out(Range(top,N-1-bottom),Range(left,M-1-right)); + + return out_cropped; +} + +LogPolar_Interp::~LogPolar_Interp() +{ +} + +//------------------------------------overlapping---------------------------------- + +LogPolar_Overlapping::LogPolar_Overlapping(int w, int h, Point2i center, int R, double ro0, int full, int S, int sp) +{ + if ( (center.x!=w/2 || center.y!=h/2) && full==0) full=1; + + if (center.x<0) center.x=0; + if (center.y<0) center.y=0; + if (center.x>=w) center.x=w-1; + if (center.y>=h) center.y=h-1; + + if (full){ + int rtmp; + + if (center.x<=w/2 && center.y>=h/2) + rtmp=(int)sqrt((float)center.y*center.y + (float)(w-center.x)*(w-center.x)); + if (center.x>=w/2 && center.y>=h/2) + rtmp=(int)sqrt((float)center.y*center.y + (float)center.x*center.x); + if (center.x>=w/2 && center.y<=h/2) + rtmp=(int)sqrt((float)(h-center.y)*(h-center.y) + (float)center.x*center.x); + if (center.x<=w/2 && center.y<=h/2) + rtmp=(int)sqrt((float)(h-center.y)*(h-center.y) + (float)(w-center.x)*(w-center.x)); + + M=2*rtmp; N=2*rtmp; + + top = M/2 - center.y; + bottom = M/2 - (h-center.y); + left = M/2 - center.x; + right = M/2 - (w - center.x); + + }else{ + top=bottom=left=right=0; + M=w; N=h; + } + + + if (sp){ + int jc=M/2-1, ic=N/2-1; + int romax=min(ic, jc); + double a=exp(log((double)(romax/2-1)/(double)ro0)/(double)R); + S=(int) floor(2*M_PI/(a-1)+0.5); + } + + create_map(M, N, R, S, ro0); +} + +void LogPolar_Overlapping::create_map(int M, int N, int R, int S, double ro0) +{ + this->M=M; + this->N=N; + this->R=R; + this->S=S; + this->ro0=ro0; + + int jc=N/2-1, ic=M/2-1; + romax=min(ic, jc); + a=exp(log((double)romax/(double)ro0)/(double)R); + q=((double)S)/(2*M_PI); + ind1=0; + + Rsri=Mat::zeros(S,R,CV_32FC1); + Csri=Mat::zeros(S,R,CV_32FC1); + ETAyx=Mat::zeros(N,M,CV_32FC1); + CSIyx=Mat::zeros(N,M,CV_32FC1); + Rsr.resize(R*S); + Csr.resize(R*S); + Wsr.resize(R); + w_ker_2D.resize(R*S); + + for(int v=0; v(v,u)=(float)(ro0*pow(a,u)*sin(v/q)+jc); + Csri.at(v,u)=(float)(ro0*pow(a,u)*cos(v/q)+ic); + Rsr[v*R+u]=(int)floor(Rsri.at(v,u)); + Csr[v*R+u]=(int)floor(Csri.at(v,u)); + } + } + + bool done=false; + + for(int i=0; i1)&&(done==false)) + { + ind1=i; + done =true; + } + } + + for(int j=0; j=ic) + theta=atan((double)(j-jc)/(double)(i-ic)); + else + theta=atan((double)(j-jc)/(double)(i-ic))+M_PI; + + if(theta<0) + theta+=2*M_PI; + + ETAyx.at(j,i)=(float)(q*theta); + + double ro2=(j-jc)*(j-jc)+(i-ic)*(i-ic); + CSIyx.at(j,i)=(float)(0.5*log(ro2/(ro0*ro0))/log(a)); + } + } + + for(int v=0; v(v,u)-Csr[v*R+u]; + double dy=Rsri.at(v,u)-Rsr[v*R+u]; + double tot=0; + for(int j=0; j<2*w+1; j++) + for(int i=0; i<2*w+1; i++) + { + (w_ker_2D[v*R+u].weights)[j*(2*w+1)+i]=exp(-(pow(i-w-dx, 2)+pow(j-w-dy, 2))/(2*sigma*sigma)); + tot+=(w_ker_2D[v*R+u].weights)[j*(2*w+1)+i]; + } + for(int j=0; j<(2*w+1); j++) + for(int i=0; i<(2*w+1); i++) + (w_ker_2D[v*R+u].weights)[j*(2*w+1)+i]/=tot; + } +} + +const Mat LogPolar_Overlapping::to_cortical(const Mat &source) +{ + Mat out(S,R,CV_8UC1,Scalar(0)); + + Mat source_border; + copyMakeBorder(source,source_border,top,bottom,left,right,BORDER_CONSTANT,Scalar(0)); + + remap(source_border,out,Csri,Rsri,INTER_LINEAR); + + int wm=w_ker_2D[R-1].w; + vector IMG((M+2*wm+1)*(N+2*wm+1), 0); + + for(int j=0; j(j,i); + + for(int v=0; v(v,u)=(uchar) floor(tmp+0.5); + } + + return out; +} + +const Mat LogPolar_Overlapping::to_cartesian(const Mat &source) +{ + Mat out(N,M,CV_8UC1,Scalar(0)); + + Mat source_border; + copyMakeBorder(source,source_border,0,1,0,0,BORDER_CONSTANT,Scalar(0)); + Mat rowS = source_border.row(S); + source_border.row(0).copyTo(rowS); + remap(source_border,out,CSIyx,ETAyx,INTER_LINEAR); + + int wm=w_ker_2D[R-1].w; + + vector IMG((N+2*wm+1)*(M+2*wm+1), 0.); + vector NOR((N+2*wm+1)*(M+2*wm+1), 0.); + + for(int v=0; v(v, u); + NOR[ind]+=((w_ker_2D[v*R+u]).weights[j*(2*w+1)+i]); + } + } + } + + for(int i=0; i<((N+2*wm+1)*(M+2*wm+1)); i++) + IMG[i]/=NOR[i]; + + //int xc=M/2-1, yc=N/2-1; + + for(int j=wm; j0) + ret[M*(j-wm)+i-wm]=(int) floor(IMG[(M+2*wm+1)*j+i]+0.5);*/ + //int ro=(int)floor(sqrt((double)((j-wm-yc)*(j-wm-yc)+(i-wm-xc)*(i-wm-xc)))); + int csi=(int) floor(CSIyx.at(j-wm,i-wm)); + + if((csi>=(ind1-(w_ker_2D[ind1]).w))&&(csi(j-wm,i-wm)=(uchar) floor(IMG[(M+2*wm+1)*j+i]+0.5); + } + + Mat out_cropped=out(Range(top,N-1-bottom),Range(left,M-1-right)); + return out_cropped; +} + +LogPolar_Overlapping::~LogPolar_Overlapping() +{ +} + +//----------------------------------------adjacent--------------------------------------- + +LogPolar_Adjacent::LogPolar_Adjacent(int w, int h, Point2i center, int R, double ro0, double smin, int full, int S, int sp) +{ + if ( (center.x!=w/2 || center.y!=h/2) && full==0) full=1; + + if (center.x<0) center.x=0; + if (center.y<0) center.y=0; + if (center.x>=w) center.x=w-1; + if (center.y>=h) center.y=h-1; + + if (full){ + int rtmp; + + if (center.x<=w/2 && center.y>=h/2) + rtmp=(int)sqrt((float)center.y*center.y + (float)(w-center.x)*(w-center.x)); + if (center.x>=w/2 && center.y>=h/2) + rtmp=(int)sqrt((float)center.y*center.y + (float)center.x*center.x); + if (center.x>=w/2 && center.y<=h/2) + rtmp=(int)sqrt((float)(h-center.y)*(h-center.y) + (float)center.x*center.x); + if (center.x<=w/2 && center.y<=h/2) + rtmp=(int)sqrt((float)(h-center.y)*(h-center.y) + (float)(w-center.x)*(w-center.x)); + + M=2*rtmp; N=2*rtmp; + + top = M/2 - center.y; + bottom = M/2 - (h-center.y); + left = M/2 - center.x; + right = M/2 - (w - center.x); + + }else{ + top=bottom=left=right=0; + M=w; N=h; + } + + if (sp){ + int jc=M/2-1, ic=N/2-1; + int romax=min(ic, jc); + double a=exp(log((double)(romax/2-1)/(double)ro0)/(double)R); + S=(int) floor(2*M_PI/(a-1)+0.5); + } + + create_map(M, N, R, S, ro0, smin); +} + + +void LogPolar_Adjacent::create_map(int M, int N, int R, int S, double ro0, double smin) +{ + LogPolar_Adjacent::M=M; + LogPolar_Adjacent::N=N; + LogPolar_Adjacent::R=R; + LogPolar_Adjacent::S=S; + LogPolar_Adjacent::ro0=ro0; + romax=min(M/2.0, N/2.0); + + a=exp(log(romax/ro0)/(double)R); + q=S/(2*M_PI); + + A.resize(R*S); + L.resize(M*N); + + for(int i=0; ismin) + { + double xs[4], ys[4]; + int us[4], vs[4]; + + xs[0]=xs[3]=x+length/4.0; + xs[1]=xs[2]=x-length/4.0; + ys[1]=ys[0]=y+length/4.0; + ys[2]=ys[3]=y-length/4.0; + + for(int z=0; z<4; z++) + get_uv(xs[z], ys[z], us[z], vs[z]); + + bool c=true; + for(int w=1; w<4; w++) + { + if(us[w]!=us[w-1]) + c=false; + if(vs[w]!=vs[w-1]) + c=false; + } + + if(c) + { + if(us[0]!=-1) + { + pixel p; + p.u=us[0]; + p.v=vs[0]; + p.a=length*length; + L[M*j+i].push_back(p); + A[vs[0]*R+us[0]]+=length*length; + } + } + + else + { + for(int z=0; z<4; z++) + if(us[z]!=-1) + subdivide_recursively(xs[z], ys[z], i, j, length/2.0, smin); + } + } +} + + +const Mat LogPolar_Adjacent::to_cortical(const Mat &source) +{ + Mat source_border; + copyMakeBorder(source,source_border,top,bottom,left,right,BORDER_CONSTANT,Scalar(0)); + + vector map(R*S, 0.); + + for(int j=0; j(j,i)); + } + } + + for(int i=0; i(i,j)=(uchar) floor(map[i*R+j]+0.5); + + return out; +} + +const Mat LogPolar_Adjacent::to_cartesian(const Mat &source) +{ + vector map(M*N, 0.); + + for(int j=0; j((L[M*j+i])[z].v,(L[M*j+i])[z].u); + } + } + + Mat out(N,M,CV_8UC1,Scalar(0)); + + for(int i=0; i(i,j)=(uchar) floor(map[i*M+j]+0.5); + + Mat out_cropped=out(Range(top,N-1-bottom),Range(left,M-1-right)); + return out_cropped; +} + + +bool LogPolar_Adjacent::get_uv(double x, double y, int&u, int&v) +{ + double ro=sqrt(x*x+y*y), theta; + if(x>0) + theta=atan(y/x); + else + theta=atan(y/x)+M_PI; + + if(roromax) + { + u=-1; + v=-1; + return false; + } + else + { + u= (int) floor(log(ro/ro0)/log(a)); + if(theta>=0) + v= (int) floor(q*theta); + else + v= (int) floor(q*(theta+2*M_PI)); + return true; + } +} + +LogPolar_Adjacent::~LogPolar_Adjacent() +{ +} + +} + diff --git a/samples/cpp/logpolar_bsm.cpp b/samples/cpp/logpolar_bsm.cpp new file mode 100644 index 0000000000..619959bbe0 --- /dev/null +++ b/samples/cpp/logpolar_bsm.cpp @@ -0,0 +1,82 @@ +/*Authors +* Manuela Chessa, Fabio Solari, Fabio Tatti, Silvio P. Sabatini +* +* manuela.chessa@unige.it, fabio.solari@unige.it +* +* PSPC-lab - University of Genoa +*/ + +#include "opencv2/opencv.hpp" +#include +#include + +using namespace cv; +using namespace std; + +void help() +{ + cout << "LogPolar Blind Spot Model sample.\nShortcuts:" + "\n\tn for nearest pixel technique" + "\n\tb for bilinear interpolation technique" + "\n\to for overlapping circular receptive fields" + "\n\ta for adjacent receptive fields" + "\n\tq or ESC quit\n"; +} + +int main(int argc, char** argv) +{ + Mat img = imread(argc > 1 ? argv[1] : "lena.jpg",1); // open the image + if(img.empty()) // check if we succeeded + { + cout << "can not load image\n"; + return 0; + } + help(); + + Size s=img.size(); + int w=s.width, h=s.height; + int ro0=3; //radius of the blind spot + int R=120; //number of rings + + //Creation of the four different objects that implement the four log-polar transformations + //Off-line computation + Point2i center(w/2,h/2); + LogPolar_Interp nearest(w, h, center, R, ro0, INTER_NEAREST); + LogPolar_Interp bilin(w,h, center,R,ro0); + LogPolar_Overlapping overlap(w,h,center,R,ro0); + LogPolar_Adjacent adj(w,h,center,R,ro0,0.25); + + namedWindow("Cartesian",1); + namedWindow("retinal",1); + namedWindow("cortical",1); + int wk='n'; + Mat Cortical, Retinal; + + //On-line computation + for(;;) + { + if(wk=='n'){ + Cortical=nearest.to_cortical(img); + Retinal=nearest.to_cartesian(Cortical); + }else if (wk=='b'){ + Cortical=bilin.to_cortical(img); + Retinal=bilin.to_cartesian(Cortical); + }else if (wk=='o'){ + Cortical=overlap.to_cortical(img); + Retinal=overlap.to_cartesian(Cortical); + }else if (wk=='a'){ + Cortical=adj.to_cortical(img); + Retinal=adj.to_cartesian(Cortical); + } + + imshow("Cartesian", img); + imshow("cortical", Cortical); + imshow("retinal", Retinal); + + int c=waitKey(15); + if (c>0) wk=c; + if(wk =='q' || (wk & 255) == 27) break; + } + + return 0; +}