如何用matlab计算体积

spano

知道了平面上一系列点的坐标和高程,用什么洋的算法可以算出由这些点形成的空间体的坐标(就是形成不规则棱形拼接成的空间棱形).

huright

既然是平面上的一系列的点，那么底面积就比较好求了。采用分割法，求出每一个三棱体的高，再相加，不就是总的体积了吗？

spano

但是是不规则的棱体啊，怎么找出其中最小的三棱体组合？
或者说是有一个点，怎么得到离它最近的2个点？
另外算的时候怎么控制运算的顺序？如果是方格的话可以按行优先或者列优先的顺序来，可是这个没有办法找到规律啊

twomao

我也不会。一点不成熟的想法，希望对你能有帮助。
给步长设置一个比较合理的范围，以该步长向外搜索，不在范围之内的点舍弃，如果范围之内有多个点，只取最近的两个。计算以起点和两个新点为顶点的三角形的面积。再以新点为起点，继续下一轮。用过的点不再用。

spano

恩，可以按照你的思路写下，谢谢

spano

换的是插值的算法，按照符合要求的间隔密化点。
然后按照点的位置对点的高程H值求权P，最后求H*P的总和就可以了

因为不是用来算标准几何体的体积，这种精度够用了，谢谢上面的指导。

kangtj

能把你写的程序帖出来供大家学习吗？！

spano

原帖由 kangtj 于 2007-11-11 19:32 发表

                               
登录/注册后可看大图

能把你写的程序帖出来供大家学习吗？！

clear;
[filename,pathname] = uigetfile('*.xyz','Select a file for read...');
InitialData = load(strcat(pathname,filename));
X=InitialData(:,1);
Y=InitialData(:,2);
Z=InitialData(:,3);
minX=min(X);
maxX=max(X);
minY=min(Y);
maxY=max(Y);
step=1;%设置单元格大小
Xi=minX:step:maxX;
Yi=minY:step:maxY;
[X2,Y2]=meshgrid(Xi,Yi);
Z2=griddata(X,Y,Z,X2,Y2,'cubic');
[m,n]=size(Z2);
top_bottom = ones(1,n)*NaN;
left_right = ones(1,m+2)'*NaN;
Z2 =[top_bottom;Z2;top_bottom];
Z2 =[left_right,Z2,left_right];
volume = 0;
for i =2:m+1
    for j=2:n+1
        if ~isnan(Z2(i,j))%为中间的每个点分配权
           counter = 0;
           if ~isnan(Z2(i-1,j-1))
               counter = counter+1;
           end
           if ~isnan(Z2(i-1,j))
               counter = counter+1;
           end
           if ~isnan(Z2(i-1,j+1))
               counter =counter+1;
           end
           if ~isnan(Z2(i,j-1))
               counter =counter+1;
           end
           if ~isnan(Z2(i,j+1))
               counter =counter+1;
           end
           if ~isnan(Z2(i+1,j-1))
               counter =counter+1;
           end
           if ~isnan(Z2(i+1,j))
               counter =counter+1;
           end
           if ~isnan(Z2(i+1,j+1))
               counter =counter+1;
           end
           cell_volume = counter/8*Z2(i,j)*step^2;
           volume = cell_volume+volume;
        end
    end
end



以上计算的都很基础，现在我想设定计算的范围又不知道怎么弄了。只能知道多边形定点的坐标，怎么判断点是否在多边形范围内呢？

[ 本帖最后由 spano 于 2007-11-12 14:38 编辑 ]

eight

只能知道多边形定点的坐标，怎么判断点是否在多边形范围内呢？

help inpolygon

这个问题我已经回答过多次，下次请先搜索一下

spano

原帖由 eight 于 2007-11-16 17:12 发表

                               
登录/注册后可看大图

help inpolygon

这个问题我已经回答过多次，下次请先搜索一下

...
没有找到,我后来自己写了个判断的函数.

spano

原帖由 spano 于 2007-11-16 18:56 发表

                               
登录/注册后可看大图

...
没有找到,我后来自己写了个判断的函数.

function mark=determine(polygon,point)
%令PI为a,a=8
a=8;
polypoint_num= size(polygon,1);
polygon =[polygon;polygon(1,:)];
t=zeros(size(polygon,1));
sum = 0;
x=polygon(:,1)-point(i,1);
y=polygon(:,2)-point(i,2);
for j=1:polypoint_num+1
    if x(j)>=0&y(j)>=0%判断象限
        t(j)=1;
    elseif x(j)<0&y(j)>=0
        t(j)=2;
    elseif x(j)<0&y(j)<0
        t(j)=3;
    elseif x(j)>=0&y(j)<0
        t(j)=4;
    end
end
for j =1:polypoint_num
    if x(j)==0&y(j)==0 %点为端点
        sum=16;
        break;
    end
    if x(j)*y(j+1)-x(j+1)*y(j)==0&x(j)*x(j+1)<=0&y(j)*y(j+1)<=0 %点在边上
        sum =16;
        break;
    end
    if mod(t(j+1),4)==mod(t(j)+1,4)
        sum=sum+4;
    elseif mod(t(j+1),4)==mod(t(j)+3,4)
        sum =sum -4;
    elseif mod(t(j+1),4)==mod(t(j)+2,4)
        if x(j)*y(j+1)-x(j+1)*y(j)>0
            sum =sum + 8;
        else
            sum = sum-8;
        end
    end
end
if sum==0
    mark=0;
else
    mark=1;
end

spano

当将土质分为多层，然后算每层的量的时候发现计算量很大。比如10*10的范围，分为6层，然后格网大小设置为0.1，计算时间大约要10分钟。当用1的时候计算很快。但是我想精度优先，尽量得到准确的结果。
有什么办法缩短计算时间没有?

[ 本帖最后由 spano 于 2007-11-20 10:37 编辑 ]

eight

原帖由 spano 于 2007-11-20 10:32 发表

                               
登录/注册后可看大图

当将土质分为多层，然后算每层的量的时候发现计算量很大。比如10*10的范围，分为6层，然后格网大小设置为0.1，计算时间大约要10分钟。当用1的时候计算很快。但是我想精度优先，尽量得到准确的结果。
有什么办法 ...

你的文件那么多循环，当然慢啦。用 inpolygon 命令不行吗？

spano

原帖由 eight 于 2007-11-20 15:58 发表

                               
登录/注册后可看大图

你的文件那么多循环，当然慢啦。用 inpolygon 命令不行吗？

用到了inpolygon函数，不过时间上并没有缩短多少。后来想可能是网格太密了，回去把程序发上来，各位帮忙分析下

eight

原帖由 spano 于 2007-11-20 16:21 发表

                               
登录/注册后可看大图

用到了inpolygon函数，不过时间上并没有缩短多少。后来想可能是网格太密了，回去把程序发上来，各位帮忙分析下

自己 profile 一下吧