On Sat, Oct 5, 2019 at 10:02 PM Nicholas Jankowski <

[hidden email]> wrote:

>

> > With current functions implemented in Octave, is my only option to create such a grid a la meshgrid and griddata?

> > Matlab 2013 introduced a scatteredInterpolant function

>

well, potentially solving my own problem, I came across one answer to

this stackexchange:

"How to integrate over a discrete 2D surface in MATLAB?"

https://stackoverflow.com/a/27987120/4598449In particular, the last part of the answer:

"If your integration area is rather complicated or has holes, you

should consider triangulating your data.

Let's say your integration area is given by the triangulation trep,

which for example could be obtained by trep =

delaunayTriangulation(x(:), y(:)). If you have your values z

corresponding to z(i) = f(trep.Points(i,1), trep.Points(i,2)), you can

use the following integration routine. It computes the exact integral

of the linear interpolant. This is done by evaluating the areas of all

the triangles and then using these areas as weights for the

midpoint(mean)-value on each triangle.

modifying the script a bit (octave doesnt' have delaunayTriangulation

implemented as a class, but it does have delaunay)

given a set of data z(x,y) as three vectors x, y, z:

function int = tri_integ (x,y,z)

P = [x(:),y(:)]; T = delaunay(x(:), y(:)); z = z(:);

d21 = P(T(:,2),:)-P(T(:,1),:);

d31 = P(T(:,3),:)-P(T(:,1),:);

areas = abs(1/2*(d21(:,1).*d31(:,2)-d21(:,2).*d31(:,1)));

int = areas'*mean(z(T),2);

endfunction

testing it out on an irregular array:

r = logspace(-3,0)(2:end-1); theta = linspace(0,2*pi,21)(1:end-1);

x = r.*cos(theta'); y = r.*sin(theta');

[x2,y2]=meshgrid(linspace(-1,1,21);

x =[x(:);x2(:)];y = [y(:);y2(:)];

f = @(x,y) 1-sin(x.^2+y.^2);

z = f(x,y);

tri_integ(x,y,z)

ans = 1.754037103916947

(running 100x = 1.7seconds)

integral2(f,-1,1,-1,1)

ans = 1.754838406712362

(running 100x = 0.6 seconds)

trying the griddata approach:

[x3,y3]=meshgrid(sort(x),sort(y));

z3 = griddata(x,y,z,x3,y3);

trapz(sort(y),trapz(sort(x),z3,2))

ans = 1.753869939077482

(running 100x = 241 seconds)

not much less accurate, but goes from handling >7500 element arrays to

~1.9M element arrays. the delay comes from the griddata step. That

would explain hitting memory limits with large meshes.

Anyway, sorry for the message clutter. At least this will be in the

message archive when I need to look it up again. Discussion of

Matlab's scatteredInterpolant mentioned triangulation a few times, so

I'm wondering if it might have done something similar (rather than

just wrapping griddata) for efficiency. Might revisit that patch once

I'm done this current project.

nickj