Matlab讨论区 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

whatman 发表于 2007-6-22 09:29

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

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

xjzuo 发表于 2007-6-23 13:00

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

happy 发表于 2007-6-23 21:30

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

function = imlogpolar(varargin)
%IMLOGPOLAR Compute logarithmic polar transformation of image.
% B = IMLOGPOLAR(A,NRHO,NTHETA,METHOD) computes the logarithmic
% polar transformation of image A, generating a log polar image
% of size NRHO by NTHETA.METHOD describes the interpolation
% method.METHOD is a string that can have one of these values:
%
%    'nearest'(default) nearest neighbor interpolation
%
%    'bilinear' bilinear interpolation
%
%    'bicubic'bicubic interpolation
%
% If you omit the METHOD argument, IMLOGPOLAR uses the default
% method of 'nearest'.
%
% B = IMLOGPOLAR(A,NRHO,NTHETA,METHOD,CTR) assumes that the 2x1
% vector CTR contains the coordinates of the origin in image A.
% If CTR is not supplied, the default is CTR = [(m+1)/2,(n+1)/2],
% where A has n rows and m columns.
%
% B = IMLOGPOLAR(A,NRHO,NTHETA,METHOD,CTR,SHAPE) where SHAPE is a
% string that can have one of these values:
%
%    'full' - returns log polar transformation containing ALL
%             pixels from image A (the circumscribed circle
%             centered at CTR)
%
%    'valid' - returns log polar transformation containing only
%             pixels from the largest inscribed circle in image A
%             centered at CTR.
%
% If you omit the SHAPE argument, IMLOGPOLAR uses the default shape
% of 'valid'.If you specify the shape 'full', invalid values on the
% periphery of B are set to NaN.
%
% Class Support
% -------------
% The input image can be of class uint8 or double. The output
% image is of the same class as the input image.
%
% Example
% -------
%    I = imread('ic.tif');
%    J = imlogpolar(I,64,64,'bilinear');
%    imshow(I), figure, imshow(J)
%
% See also IMCROP, IMRESIZE, IMROTATE.

% Nathan D. Cahill 8-16-01, modified from:
% Clay M. Thompson 8-4-92
% Copyright 1993-1998 The MathWorks, Inc.All Rights Reserved.
% $Revision: 5.10 $$Date: 1997/11/24 15:35:33 $

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

= parse_inputs(varargin{:});

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

if threeD,
= transformImage(Image,rows,cols,Nrho,Ntheta,Method,Center,Shape);
if nargout==0,
   imshow(r,g,b);
   return;
elseif nargout==1,
   if strcmp(ClassIn,'uint8');
      rout = repmat(uint8(0),);
      rout(:,:,1) = uint8(round(r*255));
      rout(:,:,2) = uint8(round(g*255));
      rout(:,:,3) = uint8(round(b*255));
   else
      rout = zeros();
      rout(:,:,1) = r;
      rout(:,:,2) = g;
      rout(:,:,3) = b;
   end
else % nargout==3
   if strcmp(ClassIn,'uint8')
      rout = uint8(round(r*255));
      g = uint8(round(g*255));
      b = uint8(round(b*255));
   else
      rout = r;    % g,b are already defined correctly above
   end
end
else
r = transformImage(Image,rows,cols,Nrho,Ntheta,Method,Center,Shape);
if nargout==0,
   imshow(r);
   return;
end
if strcmp(ClassIn,'uint8')
   if islogical(Image)
      r = im2uint8(logical(round(r)));
   else
      r = im2uint8(r);
   end
end
rout = r;
end
function = parse_inputs(varargin)
% Outputs:A    the input image
%       Nrho the desired number of rows of transformed image
%       Nthetathe desired number of columns of transformed image
%       Methodinterpolation method (nearest,bilinear,bicubic)
%       Centerorigin of input image
%       Shape output size (full,valid)
%       Class storage class of A

error(nargchk(3,6,nargin));

A = varargin{1};

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

Nrho = varargin{2};
Ntheta = varargin{3};
Class = class(A);

if nargin < 4
Method = '';
else
Method = varargin{4};
end
if isempty(Method)
Method = 'nearest';
end
Method = lower(Method);
if ~any(strcmp(Method,{'nearest','bilinear','bicubic'}))
error('Method must be one of ''nearest'', ''bilinear'', or ''bicubic''.');
end

if nargin < 5
Center = [];
else
Center = varargin{5};
end
if isempty(Center)
Center = [(Ac+1)/2 (Ar+1)/2];
end
if length(Center(:))~=2
error('Center should be 1x2 array.');
end
if
any(Center(:)> | Center(:)<1) % THIS LINE USED TO READ 'if
any(Center(:)> | Center(:)<1)' but Ar and Ac should be
swapped round -- look at line 40 for whty this should be. A.I.Wilmer,
12th Oct 2002
num2str(['Center is
',num2str(Center(1)),',',num2str(Center(2)) ' with size of image =
',num2str(Ar),'x',num2str(Ac),' (rows,columns)'])
warning('Center supplied is not within image boundaries.');
end

if nargin < 6
Shape = '';
else
Shape = varargin{6};
end
if isempty(Shape)
Shape = 'valid';
end
Shape = lower(Shape);
if ~any(strcmp(Shape,{'full','valid'}))
error('Shape must be one of ''full'' or ''valid''.');
end

if isa(A, 'uint8'), % Convert A to Double grayscale for interpolation
if islogical(A)
   A = double(A);
else
   A = double(A)/255;
end
end
function = transformImage(A,Ar,Ac,Nrho,Ntheta,Method,Center,Shape)
% Inputs: A    the input image
%       Nrho the desired number of rows of transformed image
%       Nthetathe desired number of columns of transformed image
%       Methodinterpolation method (nearest,bilinear,bicubic)
%       Centerorigin of input image
%       Shape output size (full,valid)
%       Class storage class of A

global rho;

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

switch Shape
case 'full'
corners = ;
d = max(sqrt(sum((repmat(Center(:)',4,1)-corners).^2,2)));
case 'valid'
d = min();
end
minScale = 1;
rho = logspace(log10(minScale),log10(d),Nrho)';% default 'base 10' logspace - play with d to change the scale of the log axis

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

if nargout==3
if strcmp(Method,'nearest'), % Nearest neighbor interpolation
r=interp2(A(:,:,1),xx,yy,'nearest');
g=interp2(A(:,:,2),xx,yy,'nearest');
b=interp2(A(:,:,3),xx,yy,'nearest');
elseif strcmp(Method,'bilinear'), % Linear interpolation
r=interp2(A(:,:,1),xx,yy,'linear');
g=interp2(A(:,:,2),xx,yy,'linear');
b=interp2(A(:,:,3),xx,yy,'linear');
elseif strcmp(Method,'bicubic'), % Cubic interpolation
r=interp2(A(:,:,1),xx,yy,'cubic');
g=interp2(A(:,:,2),xx,yy,'cubic');
b=interp2(A(:,:,3),xx,yy,'cubic');
else
error(['Unknown interpolation method: ',method]);
end
% any pixels outside , pad with black
mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1);
r(mask)=NaN;
g(mask)=NaN;
b(mask)=NaN;
else
if strcmp(Method,'nearest'), % Nearest neighbor interpolation
r=interp2(A,xx,yy,'nearest');
elseif strcmp(Method,'bilinear'), % Linear interpolation
r=interp2(A,xx,yy,'linear');
elseif strcmp(Method,'bicubic'), % Cubic interpolation
r=interp2(A,xx,yy,'cubic');
else
error(['Unknown interpolation method: ',method]);
end
% any pixels outside warp, pad with black
mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1);
r(mask)=NaN;
end

whatman 发表于 2007-6-24 09:59

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

happy 发表于 2007-6-24 20:16

原帖由 whatman 于 2007-6-24 09:59 发表 http://www.chinavib.com/forum/images/common/back.gif
呵呵，谢谢楼上两位，这个我看过，可是有点不明白，极坐标下求得一个dleta函数，而dleta位置理论值为非整数值，用find函数得到的都是整数，在matlab中怎么达到这个要求啊？
谢谢了！

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

whatman 发表于 2007-6-26 12:08

好的，谢谢楼上

zhanghongling 发表于 2007-10-18 17:37

回复 #1 whatman 的帖子

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

whatman 发表于 2007-10-27 09:42

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

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

zhanghongling 发表于 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 编辑 ]

zhanghongling 发表于 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,);
g1=imcrop(g,);
R=medfilt2(f1);
b=medfilt2(g1);
figure,
subplot(221);imshow(R);
subplot(222); imshow(b);
=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);
= 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 发表于 2007-11-8 11:27

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

zhanghongling 发表于 2007-11-22 11:43

回复 #12 whatman 的帖子

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

lmeagle01 发表于 2007-12-20 22:35

有结果了吗？

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

dorawan 发表于 2009-9-21 19:54

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

页: [1] 2