Wolfgang Schwanghart

While this function has a reasonable help, it lacks various features that I'd wish to see in the future to make this function valuable for others. First, provide a h1 line. It's the first line in your help block and can be found by lookfor. Currently, your h1 line is empty and thus the function can never be found.

You allocate a 3d-array of zeros which you subsequently fill with ones in 3 nested while loops. If there are only ones and zeros, why don't you use a logical array. It requires much less memory than using double precision.

In line 9 you than call the function zeros in following manner:
XYZ = zeros(size(XY,1), size(XY,2), max(max(XY)));

This may produce a warning (actually I was wondering that it not even returns an error) for DEMs with a maximum elevation that is not an integer value.
Warning: Size vector should be a row vector with integer elements.

As such, you implicitly assume that everyone uses integer representations of DEMs, that may often be not the case. Moreover, wouln't it be better to give the user more control on how your function handles pixel to voxel conversion? I think, you should allow the user to define a base level (now it will always be the sea level or another fixed reference value) and vertical resolution of the generated 3d-array. Moreover, I'd like to see that the function returns a coordinate system for the 3d array (for example by calling meshgrid(x,y,z) or so.

Currently the function is heavily looped and I am sure there are ways to make it computationally more efficient.
 

Sorry, this function is absolutely useless. It can be easily replaced by a one-liner
x2 = x(~isnan(x));

Your function uses a loop and is incredibly slow since it does not preallocate the result vector.

Thanks for this great contribution! Would it, however, be possible that ipdm (in case only data1 is provided) only returns the row and column indices and distances so that they are arranged in the order (distance(rowix,colix)): (2,1), (3,1), ..., (n,1), (3,2), ..., (n,2), ..., (n,n–1)). That is the way pdist (part of the statistics toolbox) returns the distance matrix in a vector to be more memory efficient. pdist, however, does not support distance subsets.

So far, I use this workaround in my code.

% remove distances were d.columnindex = d.rowindex
I = d.columnindex==d.rowindex;
d.columnindex(I) = [];
d.rowindex(I) = [];
d.distance(I) = [];

% remove double entry of pairs
[dtemp, m] = unique(sort([d.rowindex d.columnindex],2),'rows');
d.rowindex = dtemp(:,1);
d.columnindex = dtemp(:,2);
d.distance = d.distance(m);

This reduces the size of d by more than a half. But it doesn't avoid the generation of the huge struct returned by ipdm, which is a real memory bottleneck for me. Would it be possible to calculate this more efficient output within ipdm? So far I don't see a straightforward way, but perhaps you do?

Best regards,
Wolfgang

To all: TopoToolbox is now available here: http://physiogeo.unibas.ch/topotoolbox/

TopoToolbox can handle larger grids than wflowacc can and has many more utilities related to digital elevation model analysis.

that was what I meant, thank you :)
one more question, how to pass the data that has been plot? i mean here:
plot(out.distance,out.val,marker)
I need to do some curve fitting on those values...