如何做图像的对数极坐标变换

whatman

大家好！
将两副图像进行对数极坐标变换，然后再做fft2变换，用cart2pol做完不知道该怎么存，要将对数极坐标下的点存成二维矩阵，且分辨率不知道该怎么表示出来，是不是不能用cart2pol做啊！
请指点！

xjzuo

自己看一下"基于图像对数极坐标变换的多分辨率相关匹配算法"这篇文章吧.

happy

下面三个函数完成极对数坐标转换

1. function [rout,g,b] = imlogpolar(varargin)
   
2. %IMLOGPOLAR Compute logarithmic polar transformation of image.
   
3. %   B = IMLOGPOLAR(A,NRHO,NTHETA,METHOD) computes the logarithmic
   
4. %   polar transformation of image A, generating a log polar image
   
5. %   of size NRHO by NTHETA.  METHOD describes the interpolation
   
6. %   method.  METHOD is a string that can have one of these values:
   
7. %
   
8. %        'nearest'  (default) nearest neighbor interpolation
   
9. %
   
10. %        'bilinear' bilinear interpolation
    
11. %
    
12. %        'bicubic'  bicubic interpolation
    
13. %
    
14. %   If you omit the METHOD argument, IMLOGPOLAR uses the default
    
15. %   method of 'nearest'.
    
16. %
    
17. %   B = IMLOGPOLAR(A,NRHO,NTHETA,METHOD,CTR) assumes that the 2x1
    
18. %   vector CTR contains the coordinates of the origin in image A.  
    
19. %   If CTR is not supplied, the default is CTR = [(m+1)/2,(n+1)/2],
    
20. %   where A has n rows and m columns.
    
21. %
    
22. %   B = IMLOGPOLAR(A,NRHO,NTHETA,METHOD,CTR,SHAPE) where SHAPE is a
    
23. %   string that can have one of these values:
    
24. %
    
25. %        'full' - returns log polar transformation containing ALL
    
26. %                 pixels from image A (the circumscribed circle
    
27. %                 centered at CTR)
    
28. %
    
29. %        'valid' - returns log polar transformation containing only
    
30. %                 pixels from the largest inscribed circle in image A
    
31. %                 centered at CTR.
    
32. %
    
33. %   If you omit the SHAPE argument, IMLOGPOLAR uses the default shape
    
34. %   of 'valid'.  If you specify the shape 'full', invalid values on the
    
35. %   periphery of B are set to NaN.
    
36. %
    
37. %   Class Support
    
38. %   -------------
    
39. %   The input image can be of class uint8 or double. The output
    
40. %   image is of the same class as the input image.
    
41. %
    
42. %   Example
    
43. %   -------
    
44. %        I = imread('ic.tif');
    
45. %        J = imlogpolar(I,64,64,'bilinear');
    
46. %        imshow(I), figure, imshow(J)
    
47. %
    
48. %   See also IMCROP, IMRESIZE, IMROTATE.
    
49. 
    
50. %   Nathan D. Cahill 8-16-01, modified from:
    
51. %   Clay M. Thompson 8-4-92
    
52. %   Copyright 1993-1998 The MathWorks, Inc.  All Rights Reserved.
    
53. %   $Revision: 5.10 $  $Date: 1997/11/24 15:35:33 $
    
54. 
    
55. % Grandfathered:
    
56. %   Without output arguments, IMLOGPOLAR(...) displays the transformed
    
57. %   image in the current axis.  
    
58. 
    
59. [Image,rows,cols,Nrho,Ntheta,Method,Center,Shape,ClassIn] = parse_inputs(varargin{:});
    
60. 
    
61. threeD = (ndims(Image)==3); % Determine if input includes a 3-D array
    
62. 
    
63. if threeD,
    
64.    [r,g,b] = transformImage(Image,rows,cols,Nrho,Ntheta,Method,Center,Shape);
    
65.    if nargout==0,
    
66.       imshow(r,g,b);
    
67.       return;
    
68.    elseif nargout==1,
    
69.       if strcmp(ClassIn,'uint8');
    
70.          rout = repmat(uint8(0),[size(r),3]);
    
71.          rout(:,:,1) = uint8(round(r*255));
    
72.          rout(:,:,2) = uint8(round(g*255));
    
73.          rout(:,:,3) = uint8(round(b*255));
    
74.       else
    
75.          rout = zeros([size(r),3]);
    
76.          rout(:,:,1) = r;
    
77.          rout(:,:,2) = g;
    
78.          rout(:,:,3) = b;
    
79.       end
    
80.    else % nargout==3
    
81.       if strcmp(ClassIn,'uint8')
    
82.          rout = uint8(round(r*255));
    
83.          g = uint8(round(g*255));
    
84.          b = uint8(round(b*255));
    
85.       else
    
86.          rout = r;        % g,b are already defined correctly above
    
87.       end
    
88.    end
    
89. else
    
90.    r = transformImage(Image,rows,cols,Nrho,Ntheta,Method,Center,Shape);
    
91.    if nargout==0,
    
92.       imshow(r);
    
93.       return;
    
94.    end
    
95.    if strcmp(ClassIn,'uint8')
    
96.       if islogical(Image)
    
97.          r = im2uint8(logical(round(r)));   
    
98.       else
    
99.          r = im2uint8(r);
    
100.       end
     
101.    end
     
102.    rout = r;
     
103. end

复制代码

1. function [A,Ar,Ac,Nrho,Ntheta,Method,Center,Shape,Class] = parse_inputs(varargin)
   
2. % Outputs:  A       the input image
   
3. %           Nrho    the desired number of rows of transformed image
   
4. %           Ntheta  the desired number of columns of transformed image
   
5. %           Method  interpolation method (nearest,bilinear,bicubic)
   
6. %           Center  origin of input image
   
7. %           Shape   output size (full,valid)
   
8. %           Class   storage class of A
   
9. 
   
10. error(nargchk(3,6,nargin));
    
11. 
    
12. A = varargin{1};
    
13. 
    
14. Ar = size(A,1);     % Ar = number of rows of the input image
    
15. Ac = size(A,2);     % Ac = number of columns of the input image
    
16. 
    
17. Nrho = varargin{2};
    
18. Ntheta = varargin{3};
    
19. Class = class(A);
    
20. 
    
21. if nargin < 4
    
22.     Method = '';
    
23. else
    
24.     Method = varargin{4};
    
25. end
    
26. if isempty(Method)
    
27.     Method = 'nearest';
    
28. end
    
29. Method = lower(Method);
    
30. if ~any(strcmp(Method,{'nearest','bilinear','bicubic'}))
    
31.     error('Method must be one of ''nearest'', ''bilinear'', or ''bicubic''.');
    
32. end
    
33. 
    
34. if nargin < 5
    
35.     Center = [];
    
36. else
    
37.     Center = varargin{5};
    
38. end
    
39. if isempty(Center)
    
40.     Center = [(Ac+1)/2 (Ar+1)/2];
    
41. end
    
42. if length(Center(:))~=2
    
43.     error('Center should be 1x2 array.');
    
44. end
    
45. if
    
46. any(Center(:)>[Ac;Ar] | Center(:)<1) % THIS LINE USED TO READ 'if
    
47. any(Center(:)>[Ar;Ac] | Center(:)<1)' but Ar and Ac should be
    
48. swapped round -- look at line 40 for whty this should be. A.I.Wilmer,
    
49. 12th Oct 2002
    
50. num2str(['Center is
    
51. ',num2str(Center(1)),',',num2str(Center(2)) ' with size of image =
    
52. ',num2str(Ar),'x',num2str(Ac),' (rows,columns)'])
    
53.     warning('Center supplied is not within image boundaries.');
    
54. end
    
55. 
    
56. if nargin < 6
    
57.     Shape = '';
    
58. else
    
59.     Shape = varargin{6};
    
60. end
    
61. if isempty(Shape)
    
62.     Shape = 'valid';
    
63. end
    
64. Shape = lower(Shape);
    
65. if ~any(strcmp(Shape,{'full','valid'}))
    
66.     error('Shape must be one of ''full'' or ''valid''.');
    
67. end
    
68. 
    
69. if isa(A, 'uint8'),     % Convert A to Double grayscale for interpolation
    
70.    if islogical(A)
    
71.       A = double(A);
    
72.    else
    
73.       A = double(A)/255;
    
74.    end
    
75. end

复制代码

1. function [r,g,b] = transformImage(A,Ar,Ac,Nrho,Ntheta,Method,Center,Shape)
   
2. % Inputs:   A       the input image
   
3. %           Nrho    the desired number of rows of transformed image
   
4. %           Ntheta  the desired number of columns of transformed image
   
5. %           Method  interpolation method (nearest,bilinear,bicubic)
   
6. %           Center  origin of input image
   
7. %           Shape   output size (full,valid)
   
8. %           Class   storage class of A
   
9. 
   
10. global rho;
    
11. 
    
12. theta = linspace(0,2*pi,Ntheta+1); theta(end) = [];
    
13. 
    
14. switch Shape
    
15. case 'full'
    
16.     corners = [1 1;Ar 1;Ar Ac;1 Ac];
    
17.     d = max(sqrt(sum((repmat(Center(:)',4,1)-corners).^2,2)));
    
18. case 'valid'
    
19.     d = min([Ac-Center(1) Center(1)-1 Ar-Center(2) Center(2)-1]);
    
20. end
    
21. minScale = 1;
    
22. rho = logspace(log10(minScale),log10(d),Nrho)';  % default 'base 10' logspace - play with d to change the scale of the log axis
    
23. 
    
24. % convert polar coordinates to cartesian coordinates and center
    
25. xx = rho*cos(theta+pi) + Center(1);
    
26. yy = rho*sin(theta+pi) + Center(2);
    
27. 
    
28. if nargout==3
    
29.   if strcmp(Method,'nearest'), % Nearest neighbor interpolation
    
30.     r=interp2(A(:,:,1),xx,yy,'nearest');
    
31.     g=interp2(A(:,:,2),xx,yy,'nearest');
    
32.     b=interp2(A(:,:,3),xx,yy,'nearest');
    
33.   elseif strcmp(Method,'bilinear'), % Linear interpolation
    
34.     r=interp2(A(:,:,1),xx,yy,'linear');
    
35.     g=interp2(A(:,:,2),xx,yy,'linear');
    
36.     b=interp2(A(:,:,3),xx,yy,'linear');
    
37.   elseif strcmp(Method,'bicubic'), % Cubic interpolation
    
38.     r=interp2(A(:,:,1),xx,yy,'cubic');
    
39.     g=interp2(A(:,:,2),xx,yy,'cubic');
    
40.     b=interp2(A(:,:,3),xx,yy,'cubic');
    
41.   else
    
42.     error(['Unknown interpolation method: ',method]);
    
43.   end
    
44.   % any pixels outside , pad with black
    
45.   mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1);
    
46.   r(mask)=NaN;
    
47.   g(mask)=NaN;
    
48.   b(mask)=NaN;
    
49. else
    
50.   if strcmp(Method,'nearest'), % Nearest neighbor interpolation
    
51.     r=interp2(A,xx,yy,'nearest');
    
52.   elseif strcmp(Method,'bilinear'), % Linear interpolation
    
53.     r=interp2(A,xx,yy,'linear');
    
54.   elseif strcmp(Method,'bicubic'), % Cubic interpolation
    
55.     r=interp2(A,xx,yy,'cubic');
    
56.   else
    
57.     error(['Unknown interpolation method: ',method]);
    
58.   end
    
59.   % any pixels outside warp, pad with black
    
60.   mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1);
    
61.   r(mask)=NaN;
    
62. end

复制代码

whatman

呵呵，谢谢楼上两位，这个我看过，可是有点不明白，极坐标下求得一个dleta函数，而dleta位置理论值为非整数值，用find函数得到的都是整数，在matlab中怎么达到这个要求啊？
谢谢了！

happy

原帖由 whatman 于 2007-6-24 09:59 发表

                               
登录/注册后可看大图

呵呵，谢谢楼上两位，这个我看过，可是有点不明白，极坐标下求得一个dleta函数，而dleta位置理论值为非整数值，用find函数得到的都是整数，在matlab中怎么达到这个要求啊？
谢谢了！

搞清楚find函数就不成问题了

whatman

好的，谢谢楼上

zhanghongling

我也遇到同样的问题，能不能指教一下啊

whatman

我也只是用了一下上面的程序计算了一下，深入的还没有研究明白，可以共同切磋啊！

[ 本帖最后由 无水1324 于 2007-11-7 13:12 编辑 ]

zhanghongling

高手解释一下：d = max(sqrt(sum((repmat(Center(:)',4,1)-corners).^2,2)));这实现的是什么功能啊？

花如月

由内向外，把程序逐句化简

[ 本帖最后由 eight 于 2007-10-31 22:00 编辑 ]

zhanghongling

我引用上面的程序计算出来的旋转度数不是90，就是180，270，360度的，可事实上旋转角度并不是这么多，这是怎么回事啊？你的计算结果也是这样吗？以下是我的调用程序：
clear, close all, clc
f=imread('mein.bmp');
g=imrotate(f,60,'bicubic','crop');
f1=imcrop(f,[50,50,128,128]);
g1=imcrop(g,[80,80,128,128]);
R=medfilt2(f1);
b=medfilt2(g1);
figure,
subplot(221);imshow(R);
subplot(222); imshow(b);
[m,n]=size(R);
b1 = imlogpolar(R,128,128,'bicubic',[(m+1)/2,(n+1)/2],'full');
b2 = imlogpolar(b,128,128,'bicubic',[(m+1)/2,(n+1)/2],'full');
subplot(223);imshow(b1);
subplot(224);imshow(b2);
F=fft2(b1);
G=fft2(b2);
FG = conj(F).*G;
cc=FG./abs(FG);
s=ifft2(cc);
figure;
mesh(abs(s));
M = max(max(s));
disp('输出最大峰值;');
disp(M);
[i j] = find(s == M);
theta=(128-j+1)/128*360;
disp(['输出平移坐标: (' num2str(i) ',' num2str(j) ') pixels '] );
disp(['Estimated shuofang: ' num2str(10^(i-1))  ] );
disp(['Estimated rotation: ' num2str(theta) '度' ] );
是我的调用程序的问题吗？

whatman

这个程序不是多大的角度都可以的，只有小角度的才可以，并且跟你所选图像的大小有关系，调用函数中间有个分辨率的关系

zhanghongling

谢谢楼上，大角度也是可以准确求出来的。可是缩放参数该怎么求啊，经过photoshop处理过的图像求出的缩放参数总是不对，到底问题在哪儿呢？:@o

lmeagle01

楼主做的怎么样了，有结果了吗？

dorawan

我怎么调用有问题呢？？现在我要求的是对图像的频谱图映射到对数极坐标，并将谱投影到角度轴上，然后通过频谱峰值计算角度，请问该怎么做？？