Boxcount in a DEM ALOS image

Hello. I am trying to make the boxcount in a satellite image (DEM ALOS) and this was the result:

Error using logical

NaN values cannot be converted to logicals.

Error in boxcount (line 64)

c = logical(squeeze(c));

Error in Alos_clip01 (line 6)

boxcount(c)

c = imread('Alos_clip01.tif');

imagesc(c)

colormap gray

axis square

%%

boxcount(c)

%%

[n, r] = boxcount(c);

loglog(r, n,'bo-', r, (r/r(end)).^(-2), 'r--')

xlabel('r')

ylabel('n(r)')

legend('actual box-count','space-filling box-count');

%%

boxcount(c, 'slope')

Boxcount code:

function [n,r] = boxcount(c,varargin)

%BOXCOUNT Box-Counting of a D-dimensional array (with D=1,2,3).

% [N, R] = BOXCOUNT(C), where C is a D-dimensional array (with D=1,2,3),

% counts the number N of D-dimensional boxes of size R needed to cover

% the nonzero elements of C. The box sizes are powers of two, i.e.,

% R = 1, 2, 4 ... 2^P, where P is the smallest integer such that

% MAX(SIZE(C)) <= 2^P. If the sizes of C over each dimension are smaller

% than 2^P, C is padded with zeros to size 2^P over each dimension (e.g.,

% a 320-by-200 image is padded to 512-by-512). The output vectors N and R

% are of size P+1. For a RGB color image (m-by-n-by-3 array), a summation

% over the 3 RGB planes is done first.

%

% The Box-counting method is useful to determine fractal properties of a

% 1D segment, a 2D image or a 3D array. If C is a fractal set, with

% fractal dimension DF < D, then N scales as R^(-DF). DF is known as the

% Minkowski-Bouligand dimension, or Kolmogorov capacity, or Kolmogorov

% dimension, or simply box-counting dimension.

%

% BOXCOUNT(C,'plot') also shows the log-log plot of N as a function of R

% (if no output argument, this option is selected by default).

%

% BOXCOUNT(C,'slope') also shows the semi-log plot of the local slope

% DF = - dlnN/dlnR as a function of R. If DF is contant in a certain

% range of R, then DF is the fractal dimension of the set C. The

% derivative is computed as a 2nd order finite difference (see GRADIENT).

%

% The execution time depends on the sizes of C. It is fastest for powers

% of two over each dimension.

%

% Examples:

%

% % Plots the box-count of a vector containing randomly-distributed

% % 0 and 1. This set is not fractal: one has N = R^-2 at large R,

% % and N = cste at small R.

% c = (rand(1,2048)<0.2);

% boxcount(c);

%

% % Plots the box-count and the fractal dimension of a 2D fractal set

% % of size 512^2 (obtained by RANDCANTOR), with fractal dimension

% % DF = 2 + log(P) / log(2) = 1.68 (with P=0.8).

% c = randcantor(0.8, 512, 2);

% boxcount(c);

% figure, boxcount(c, 'slope');

%

% F. Moisy

% Revision: 2.10, Date: 2008/07/09

% History:

% 2006/11/22: v2.00, joined into a single file boxcountn (n=1,2,3).

% 2008/07/09: v2.10, minor improvements

% control input argument

narginchk(1,2);

% check for true color image (m-by-n-by-3 array)

if ndims(c)==3

if size(c,3)==3 && size(c,1)>=8 && size(c,2)>=8

c = sum(c,3);

end

end

warning off

c = logical(squeeze(c));

warning on

dim = ndims(c); % dim is 2 for a vector or a matrix, 3 for a cube

if dim>3

error('Maximum dimension is 3.');

end

% transpose the vector to a 1-by-n vector

if length(c)==numel(c)

dim=1;

if size(c,1)~=1

c = c';

end

end

width = max(size(c)); % largest size of the box

p = log(width)/log(2); % nbre of generations

% remap the array if the sizes are not all equal,

% or if they are not power of two

% (this slows down the computation!)

if p~=round(p) || any(size(c)~=width)

p = ceil(p);

width = 2^p;

switch dim

case 1

mz = zeros(1,width);

mz(1:length(c)) = c;

c = mz;

case 2

mz = zeros(width, width);

mz(1:size(c,1), 1:size(c,2)) = c;

c = mz;

case 3

mz = zeros(width, width, width);

mz(1:size(c,1), 1:size(c,2), 1:size(c,3)) = c;

c = mz;

end

end

n=zeros(1,p+1); % pre-allocate the number of box of size r

switch dim

case 1 %------------------- 1D boxcount ---------------------%

n(p+1) = sum(c);

for g=(p-1):-1:0

siz = 2^(p-g);

siz2 = round(siz/2);

for i=1:siz:(width-siz+1)

c(i) = ( c(i) || c(i+siz2));

end

n(g+1) = sum(c(1:siz:(width-siz+1)));

end

case 2 %------------------- 2D boxcount ---------------------%

n(p+1) = sum(c(:));

for g=(p-1):-1:0

siz = 2^(p-g);

siz2 = round(siz/2);

for i=1:siz:(width-siz+1)

for j=1:siz:(width-siz+1)

c(i,j) = ( c(i,j) || c(i+siz2,j) || c(i,j+siz2) || c(i+siz2,j+siz2) );

end

end

n(g+1) = sum(sum(c(1:siz:(width-siz+1),1:siz:(width-siz+1))));

end

case 3 %------------------- 3D boxcount ---------------------%

n(p+1) = sum(c(:));

for g=(p-1):-1:0

siz = 2^(p-g);

siz2 = round(siz/2);

for i=1:siz:(width-siz+1)

for j=1:siz:(width-siz+1)

for k=1:siz:(width-siz+1)

c(i,j,k)=( c(i,j,k) || c(i+siz2,j,k) || c(i,j+siz2,k) ...

|| c(i+siz2,j+siz2,k) || c(i,j,k+siz2) || c(i+siz2,j,k+siz2) ...

|| c(i,j+siz2,k+siz2) || c(i+siz2,j+siz2,k+siz2));

end

end

end

n(g+1) = sum(sum(sum(c(1:siz:(width-siz+1),1:siz:(width-siz+1),1:siz:(width-siz+1)))));

end

end

n = n(end:-1:1);

r = 2.^(0:p); % box size (1, 2, 4, 8...)

if any(strncmpi(varargin,'slope',1))

s=-gradient(log(n))./gradient(log(r));

semilogx(r, s, 's-');

ylim([0 dim]);

xlabel('r, box size'); ylabel('- d ln n / d ln r, local dimension');

title([num2str(dim) 'D box-count']);

elseif nargout==0 || any(strncmpi(varargin,'plot',1))

loglog(r,n,'s-');

xlabel('r, box size'); ylabel('n(r), number of boxes');

title([num2str(dim) 'D box-count']);

end

if nargout==0

clear r n

end