该算是一个利用MATLAB中的FFT程序开发的一个计算分形维算的程序。.zip

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)]);