声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 12867|回复: 18

[图像处理] 如何做图像的对数极坐标变换

[复制链接]
发表于 2007-6-22 09:29 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

使用道具 举报

发表于 2007-6-23 13:00 | 显示全部楼层
自己看一下"基于图像对数极坐标变换的多分辨率相关匹配算法"这篇文章吧.
发表于 2007-6-23 21:30 | 显示全部楼层
下面三个函数完成极对数坐标转换

  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. %   Nathan D. Cahill 8-16-01, modified from:
  50. %   Clay M. Thompson 8-4-92
  51. %   Copyright 1993-1998 The MathWorks, Inc.  All Rights Reserved.
  52. %   $Revision: 5.10 $  $Date: 1997/11/24 15:35:33 $

  53. % Grandfathered:
  54. %   Without output arguments, IMLOGPOLAR(...) displays the transformed
  55. %   image in the current axis.  

  56. [Image,rows,cols,Nrho,Ntheta,Method,Center,Shape,ClassIn] = parse_inputs(varargin{:});

  57. threeD = (ndims(Image)==3); % Determine if input includes a 3-D array

  58. if threeD,
  59.    [r,g,b] = transformImage(Image,rows,cols,Nrho,Ntheta,Method,Center,Shape);
  60.    if nargout==0,
  61.       imshow(r,g,b);
  62.       return;
  63.    elseif nargout==1,
  64.       if strcmp(ClassIn,'uint8');
  65.          rout = repmat(uint8(0),[size(r),3]);
  66.          rout(:,:,1) = uint8(round(r*255));
  67.          rout(:,:,2) = uint8(round(g*255));
  68.          rout(:,:,3) = uint8(round(b*255));
  69.       else
  70.          rout = zeros([size(r),3]);
  71.          rout(:,:,1) = r;
  72.          rout(:,:,2) = g;
  73.          rout(:,:,3) = b;
  74.       end
  75.    else % nargout==3
  76.       if strcmp(ClassIn,'uint8')
  77.          rout = uint8(round(r*255));
  78.          g = uint8(round(g*255));
  79.          b = uint8(round(b*255));
  80.       else
  81.          rout = r;        % g,b are already defined correctly above
  82.       end
  83.    end
  84. else
  85.    r = transformImage(Image,rows,cols,Nrho,Ntheta,Method,Center,Shape);
  86.    if nargout==0,
  87.       imshow(r);
  88.       return;
  89.    end
  90.    if strcmp(ClassIn,'uint8')
  91.       if islogical(Image)
  92.          r = im2uint8(logical(round(r)));   
  93.       else
  94.          r = im2uint8(r);
  95.       end
  96.    end
  97.    rout = r;
  98. 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. error(nargchk(3,6,nargin));

  10. A = varargin{1};

  11. Ar = size(A,1);     % Ar = number of rows of the input image
  12. Ac = size(A,2);     % Ac = number of columns of the input image

  13. Nrho = varargin{2};
  14. Ntheta = varargin{3};
  15. Class = class(A);

  16. if nargin < 4
  17.     Method = '';
  18. else
  19.     Method = varargin{4};
  20. end
  21. if isempty(Method)
  22.     Method = 'nearest';
  23. end
  24. Method = lower(Method);
  25. if ~any(strcmp(Method,{'nearest','bilinear','bicubic'}))
  26.     error('Method must be one of ''nearest'', ''bilinear'', or ''bicubic''.');
  27. end

  28. if nargin < 5
  29.     Center = [];
  30. else
  31.     Center = varargin{5};
  32. end
  33. if isempty(Center)
  34.     Center = [(Ac+1)/2 (Ar+1)/2];
  35. end
  36. if length(Center(:))~=2
  37.     error('Center should be 1x2 array.');
  38. end
  39. if
  40. any(Center(:)>[Ac;Ar] | Center(:)<1) % THIS LINE USED TO READ 'if
  41. any(Center(:)>[Ar;Ac] | Center(:)<1)' but Ar and Ac should be
  42. swapped round -- look at line 40 for whty this should be. A.I.Wilmer,
  43. 12th Oct 2002
  44. num2str(['Center is
  45. ',num2str(Center(1)),',',num2str(Center(2)) ' with size of image =
  46. ',num2str(Ar),'x',num2str(Ac),' (rows,columns)'])
  47.     warning('Center supplied is not within image boundaries.');
  48. end

  49. if nargin < 6
  50.     Shape = '';
  51. else
  52.     Shape = varargin{6};
  53. end
  54. if isempty(Shape)
  55.     Shape = 'valid';
  56. end
  57. Shape = lower(Shape);
  58. if ~any(strcmp(Shape,{'full','valid'}))
  59.     error('Shape must be one of ''full'' or ''valid''.');
  60. end

  61. if isa(A, 'uint8'),     % Convert A to Double grayscale for interpolation
  62.    if islogical(A)
  63.       A = double(A);
  64.    else
  65.       A = double(A)/255;
  66.    end
  67. 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. global rho;

  10. theta = linspace(0,2*pi,Ntheta+1); theta(end) = [];

  11. switch Shape
  12. case 'full'
  13.     corners = [1 1;Ar 1;Ar Ac;1 Ac];
  14.     d = max(sqrt(sum((repmat(Center(:)',4,1)-corners).^2,2)));
  15. case 'valid'
  16.     d = min([Ac-Center(1) Center(1)-1 Ar-Center(2) Center(2)-1]);
  17. end
  18. minScale = 1;
  19. rho = logspace(log10(minScale),log10(d),Nrho)';  % default 'base 10' logspace - play with d to change the scale of the log axis

  20. % convert polar coordinates to cartesian coordinates and center
  21. xx = rho*cos(theta+pi) + Center(1);
  22. yy = rho*sin(theta+pi) + Center(2);

  23. if nargout==3
  24.   if strcmp(Method,'nearest'), % Nearest neighbor interpolation
  25.     r=interp2(A(:,:,1),xx,yy,'nearest');
  26.     g=interp2(A(:,:,2),xx,yy,'nearest');
  27.     b=interp2(A(:,:,3),xx,yy,'nearest');
  28.   elseif strcmp(Method,'bilinear'), % Linear interpolation
  29.     r=interp2(A(:,:,1),xx,yy,'linear');
  30.     g=interp2(A(:,:,2),xx,yy,'linear');
  31.     b=interp2(A(:,:,3),xx,yy,'linear');
  32.   elseif strcmp(Method,'bicubic'), % Cubic interpolation
  33.     r=interp2(A(:,:,1),xx,yy,'cubic');
  34.     g=interp2(A(:,:,2),xx,yy,'cubic');
  35.     b=interp2(A(:,:,3),xx,yy,'cubic');
  36.   else
  37.     error(['Unknown interpolation method: ',method]);
  38.   end
  39.   % any pixels outside , pad with black
  40.   mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1);
  41.   r(mask)=NaN;
  42.   g(mask)=NaN;
  43.   b(mask)=NaN;
  44. else
  45.   if strcmp(Method,'nearest'), % Nearest neighbor interpolation
  46.     r=interp2(A,xx,yy,'nearest');
  47.   elseif strcmp(Method,'bilinear'), % Linear interpolation
  48.     r=interp2(A,xx,yy,'linear');
  49.   elseif strcmp(Method,'bicubic'), % Cubic interpolation
  50.     r=interp2(A,xx,yy,'cubic');
  51.   else
  52.     error(['Unknown interpolation method: ',method]);
  53.   end
  54.   % any pixels outside warp, pad with black
  55.   mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1);
  56.   r(mask)=NaN;
  57. end
复制代码

评分

1

查看全部评分

 楼主| 发表于 2007-6-24 09:59 | 显示全部楼层

回复

呵呵,谢谢楼上两位,这个我看过,可是有点不明白,极坐标下求得一个dleta函数,而dleta位置理论值为非整数值,用find函数得到的都是整数,在matlab中怎么达到这个要求啊?
谢谢了!
发表于 2007-6-24 20:16 | 显示全部楼层
原帖由 whatman 于 2007-6-24 09:59 发表
呵呵,谢谢楼上两位,这个我看过,可是有点不明白,极坐标下求得一个dleta函数,而dleta位置理论值为非整数值,用find函数得到的都是整数,在matlab中怎么达到这个要求啊?
谢谢了!


搞清楚find函数就不成问题了
 楼主| 发表于 2007-6-26 12:08 | 显示全部楼层
好的,谢谢楼上
发表于 2007-10-18 17:37 | 显示全部楼层

回复 #1 whatman 的帖子

我也遇到同样的问题,能不能指教一下啊
 楼主| 发表于 2007-10-27 09:42 | 显示全部楼层
我也只是用了一下上面的程序计算了一下,深入的还没有研究明白,可以共同切磋啊!

[ 本帖最后由 无水1324 于 2007-11-7 13:12 编辑 ]
发表于 2007-10-31 21:26 | 显示全部楼层
高手解释一下:d = max(sqrt(sum((repmat(Center(:)',4,1)-corners).^2,2)));这实现的是什么功能啊?
发表于 2007-10-31 21:34 | 显示全部楼层

回复 #9 zhanghongling 的帖子

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

[ 本帖最后由 eight 于 2007-10-31 22:00 编辑 ]
发表于 2007-11-7 11:03 | 显示全部楼层

回复 #8 whatman 的帖子

我引用上面的程序计算出来的旋转度数不是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) '度' ] );
是我的调用程序的问题吗?
 楼主| 发表于 2007-11-8 11:27 | 显示全部楼层
这个程序不是多大的角度都可以的,只有小角度的才可以,并且跟你所选图像的大小有关系,调用函数中间有个分辨率的关系
发表于 2007-11-22 11:43 | 显示全部楼层

回复 #12 whatman 的帖子

谢谢楼上,大角度也是可以准确求出来的。可是缩放参数该怎么求啊,经过photoshop处理过的图像求出的缩放参数总是不对,到底问题在哪儿呢?:@o
发表于 2007-12-20 22:35 | 显示全部楼层

有结果了吗?

楼主做的怎么样了,有结果了吗?
发表于 2009-9-21 19:54 | 显示全部楼层
我怎么调用有问题呢??现在我要求的是对图像的频谱图映射到对数极坐标,并将谱投影到角度轴上,然后通过频谱峰值计算角度,请问该怎么做??
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-5 19:35 , Processed in 0.078841 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表