Re: help with 3D grid data

In article <dkfo5e$1sm$1@xxxxxxxxxxxxxxxx>, Kamran <kamran@xxxxxx>
wrote:

> Thanks John for answering
>
> they are regular but not uniform if one may say that.
> What I mean is that the inter-spacing between the points
> is different.
> I am not sure what you mean by Z=0. Every point in

By Z = 0, I meant that Z moves between two
planes. I arbitrarily assigned levels for Z
on those two planes: Z = 0 on the lower plane,
and Z = 1 on the upper plane.


> the top or bottom surface has 3 coordinates, x, y and z.
> So for number_of_xgrid = 40, number_of_ygrid =100 I will
> have 41*101*3 coordinates (3 being x,y,z for each point).
> The same number for the bottom surface
>
> > Then a second set of independent x-y pairs lie in the
> > plane Z = 1.
> >
> > Now you wish to triangulate the domain between the
> > two planes, at 25 levels of Z?
> >
>
> Um, right ! I am not well versed in 3D grafics. What I
> am to do is to make a 3D grid constructed of cubes each
> having 8 corner points and each corner point having 3 coordinates
> (x,y,z). Now I have only the z-coordinates for the cubes corner-points
> and not the x and y coordinates. Since it is not a uniform grid
> I can not use the x and y coordinates of the top and layer surfaces
> but have to find these by interpolating between
> the boundary layers(surfaces) for all my z points.

You do still want to tessellate that grid? Or just
leave it as cubes?



> Though the two boundary layers are parallel they may not be
> completely aligned in z-direction. One may be slightly shifted in
> respect to the other making a somehow slant 3D shape if we were to
> connect all the corresponding points on the two. (the
> number of grid points in the two surfaces are equal).
> Maybe cube is not the right word since they are different
> in size, better call them cells.
> So having 25 z-layers in between the two boundaries I will then
> have to construct 40*100*25 cubes each represented by
> the 8 corner points. After that I thought of using matlab's
> 'patch' with a fourth dimension to represent some charactristic
> of the geometry like porosity or ...


I think I understand now what you want. There should
be several ways to accomplish it, none via delaunay.
(The easiest way would use code I cannot give away.)

I've included a function below called tess_lat which
will tessellate an arbitrary domain on a lattice.
The lattice needs not be uniform. But it does need to
be a rectangular lattice. Simplest is just to create
a rectangular lattice and then distort it. I'll give
an example. I've made a small example below.

Choose the nodes for our object on the lower Z plane
as:

XL = [0 1 3 5 8 11]';
YL = [1 2 4 8 16]';

% Note that X and Y are not uniform, nor are they the
% same lengths. I'll set Z at 11 levels. Initially I'll
% use floating point numbers between 0 and 1, and then
% transform them too. The initial scaling from 0 to 1
% is useful in the transformation process. The final
% vector of Z levels will be 10:2:30.

nz = 11;
zmin = 10;
zmax = 30;

% At the upper plane in Z, I'll modify the grid to be:

XU = 2*XL;
YU = (1:2:9)';

% First, build a tessellation of a regular lattice,
% based on XL, YL, and Z.

nx = length(XL);
ny = length(YL);
[v0,tess] = tess_lat(1:nx,1:ny,linspace(0,1,nz));

% Note that v0 is a 330*3 array. Each row is one vertex
% of a uniform lattice.

plot3(v0(:,1),v0(:,2),v0(:,3),'o')

% The array tess is what you want delaunay to return.
% A list of tetrahedra.

% We need to transform the vertices array using XL, XU,
% YL, and YU.

vertices = zeros(size(v0));

vertices(:,1) = XL(v0(:,1)).*(1-v0(:,3))+...
XU(v0(:,1)).*v0(:,3);

vertices(:,2) = YL(v0(:,2)).*(1-v0(:,3))+...
YU(v0(:,2)).*v0(:,3);

vertices(:,3) = zmin*(1-v0(:,3)) + zmax*v0(:,3);

% Plot out the transformed array

figure
plot3(vertices(:,1),vertices(:,2),vertices(:,3),'o')

% rotate this grid around to verify that it is as
% I've claimed it to be.


The simplicial complex describes by the array tess is
still valid on that transformed array of vertices, at
least so under most transformations.

Its easy enough to add the simplex edges to this last
plot. But I left something to keep you awake. ;-)

HTH,
John



% =========================================================
% Begin tess_lat code
% =========================================================
function [vertices,tess]=tess_lat(varargin)
% tess_lat: simplicial tessellation of a rectangular lattice
% usage: [tessellation,vertices]=tess_lat(p1,p2,p3,...)
%
% arguments: input
% p1,p2,p3,... - numeric vectors defining the lattice in
% each dimension.
% Each vector must be of length >= 1
%
% arguments: (output)
% vertices - factorial lattice created from (p1,p2,p3,...)
% Each row of this array is one node in the lattice
% tess - integer array defining simplexes as references to
% rows of "vertices".

% dimension of the lattice
n = length(varargin);

% create a single n-d hypercube
% list of vertices of the cube itself
vhc=('1'==dec2bin(0:(2^n-1)));
% permutations of the integers 1:n
p=perms(1:n);
nt=factorial(n);
thc=zeros(nt,n+1);
for i=1:nt
thc(i,:)=find(all(diff(vhc(:,p(i,:)),[],2)>=0,2))';
end

% build the complete lattice
nodecount = cellfun('length',varargin);
if any(nodecount<2)
error 'Each dimension must be of size 2 or more.'
end
vertices = lattice(varargin{:});

% unrolled index into each hyper-rectangle in the lattice
ind = cell(1,n);
for i=1:n
ind{i} = 0:(nodecount(i)-2);
end
ind = lattice(ind{:});
k = cumprod([1,nodecount(1:(end-1))]);
ind = 1+ind*k';
nind = length(ind);

offset=vhc*k';
tess=zeros(nt*nind,n+1);
L=(1:nind)';
for i=1:nt
tess(L,:)=repmat(ind,1,n+1)+repmat(offset(thc(i,:))',nind,1);
L=L+nind;
end

% ======== subfunction ========
function g = lattice(varargin)
% generate a factorial lattice in n variables
n=nargin;
sizes = cellfun('length',varargin);
c=cell(1,n);
[c{1:n}]=ndgrid(varargin{:});
g=zeros(prod(sizes),n);
for i=1:n
g(:,i)=c{i}(:);
end
.

Relevant Pages