function [slope, intercept, averslope, averIC] = fdsurfft(im)
% FDSURFFT Compute fractal dimension (slope) of surface image im and draw rose plots
% of slope and intercept
% im: input array of surface image (grayvalue or range image)
% slope: an array of size 24 which stores the average slopes in 24
% directions
% intercept: an array of size 24 which stores the average intercepts in 24
% directions
% averslope: average slope for all directions
% averIC: average intercept for all directions
% This is a matlab version of John C. Russ 's program
% Written by Mr. Jianbo Zhang, J.Zhang@ewi.utwente.nl
% 2 Nov, 2004
tic
NUM_DIR = 24; % number of directions that the frequency space is uniformally divided
NUM_RAD = 30; % number of points that the radius is uniformally divided
if nargin < 1,
error('Missed input argument which must be an array!')
end
[M N] = size(im);
xctr = 1 + bitshift(N, -1); % x coordinate of center point
yctr = 1 + bitshift(M, -1); % y coordinate of center point
imMean = mean(im(:));
fim = fftshift(fft2(double(im) - imMean));
% power spectrum
mag = log(fim .* conj(fim)+ 10 ^ (-6));
sumBrite = zeros(NUM_DIR, NUM_RAD); %accumulation magnitude for each direction and radius
nCount = zeros(NUM_DIR, NUM_RAD); %number of magnitude
radius = zeros(2 * NUM_RAD,1);%accumulation magnitude for all directions
radCount = zeros(2 * NUM_RAD,1);% number of magintude for all directions
%Compute phase image and phase histogram
phaseim = zeros(M,N);
phase = zeros(180);
for j = 1:M
for i = 1:N
realv = real(fim(j,i));
imagv = imag(fim(j,i));
if realv == 0
value = pi/2;
else
value = atan((imagv / realv));
phaseim(j, i) = value;
ang = floor(180 * (pi / 2 + value) / pi);
end
if ang < 0
ang = 0;
end
if ang > 179
ang = 179;
end
phase(ang + 1) = phase(ang + 1) + 1;
end
end
maxphase = max(phase);
figure; imshow(phaseim, []);
title('Phase image');
figure; plot(phase / maxphase);
title('Phase histogram (0...2 \pi)');%
axis off
%accumulation of magnitude for each dirction and radius
rmax = log(min(M,N)/2);% maximum radius
for j = 1:M
if j ~= yctr
yval = yctr - j;
y2 = yval * yval;
for i = 1:N
if i ~= xctr
xval = i - xctr;
rho = log(sqrt(y2 + xval * xval));
if rho > 0 & rho <= rmax
mval = mag(j,i);
temp = yval /xval;
theta = atan(temp);
if xval < 0
theta = theta + pi;
end
if theta < 0
theta = theta + 2 * pi;
end
ang = floor(NUM_DIR * theta /(2* pi));
if ang > NUM_DIR - 1 | ang < 0
ang = NUM_DIR - 1;
end
k = floor(2*NUM_RAD * rho / rmax);
h = floor(k/2);
if k > 2 * NUM_RAD - 1
h = NUM_RAD - 1;
k = 2 * NUM_RAD - 1;
end
if h >= 5
sumBrite(ang + 1, h + 1) = sumBrite(ang + 1, h + 1) + mval;
nCount(ang + 1, h + 1) = nCount(ang + 1, h + 1) + 1;
end
if k >= 5
radius(k + 1) = radius(k + 1) + mval;
radCount(k + 1) = radCount(k + 1) + 1;
end
end
end
end
end
end
%linear regression
for ang = 1:NUM_DIR
sumx = 0;
sumy = 0;
sumx2 = 0;
sumxy = 0;
sumn = 0;
for range = 6:NUM_RAD
if nCount(ang, range) > 0
yval = sumBrite(ang, range)/nCount(ang, range);
xval = (range -1) * rmax / NUM_RAD;
sumx = sumx + xval;
sumy = sumy + yval;
sumx2 = sumx2 + xval * xval;
sumxy = sumxy + xval * yval;
sumn = sumn + 1;
end
end
slope(ang) = (sumn * sumxy - sumx * sumy) / (sumn * sumx2 - sumx * sumx);
intercept(ang) = (sumy - slope(ang) * sumx) / sumn;
end
%compute average slope over all directions and scales
sumn = 0;
for k = 6:(2 * NUM_RAD)
if radCount(k) > 0
sumn = sumn + 1;
yval(sumn) = radius(k) / radCount(k);
tempr(sumn) = (k -1) * rmax / (2 * NUM_RAD);
end
end
p = polyfit(tempr,yval,1);
averslope = p(1);
averIC = p(2);
fitln = polyval(p,tempr);
figure; plot(tempr,yval,tempr,fitln,'r-');
title('Log Log plot of Magn. vs Freq.');
ylabel('Log Magnitude');
xlabel('Log Frequency');
slope(NUM_DIR + 1) = slope(1);
intercept(NUM_DIR + 1) = intercept(1);
%draw rose plot of slope and intercept
ang = 1: (NUM_DIR + 1);
figure;
polar(pi/NUM_DIR + (ang -1) * 2* pi / NUM_DIR, intercept(ang), 'r-');
title('Rose plot of intercept');
figure;
polar(pi/NUM_DIR + (ang - 1)* 2 * pi / NUM_DIR, abs(slope(ang)), '-');
title('Rose plot of slope');
disp(['Elapsed time: ' num2str(toc)]);
没有合适的资源?快使用搜索试试~ 我知道了~
该算是一个利用MATLAB中的FFT程序开发的一个计算分形维算的程序 .zip
共1个文件
m:1个
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
0 下载量 198 浏览量
2023-03-10
16:24:55
上传
评论
收藏 2KB ZIP 举报
温馨提示
该算是一个利用MATLAB中的FFT程序开发的一个计算分形维算的程序。.zip
资源推荐
资源详情
资源评论
收起资源包目录
该算是一个利用MATLAB中的FFT程序开发的一个计算分形维算的程序。.zip (1个子文件)
该算是一个利用MATLAB中的FFT程序开发的一个计算分形维算的程序。
Calculationofractaldimensionoffractalsurfacesusingfft.m 5KB
共 1 条
- 1
资源评论
处处清欢
- 粉丝: 150
- 资源: 2471
下载权益
C知道特权
VIP文章
课程特权
开通VIP
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功