Code covered by the BSD License  

Highlights from
Surface2VTK

Surface2VTK

by

 

Interpolates scattered data and saves the result to a vtk file.

surface2vtk(x,y,z,r,vtkname)
function surface2vtk(x,y,z,r,vtkname)
% surface2vtk. Given the scattered data r and its co-ordinates x,y,z,
% surface2vtk interpolates the data to a regular grid and saves the
% resulting surface as a legacy vtk file which can be visualized in vtk.
% e.g. surface2vtk(x,y,z,r,'surface.vtk')
tic;

%---read data in 'name'
%data=importdata(name);
%x=data(:,1);
%y=data(:,2);
%z=data(:,3);
%r=log10(data(:,4));


%---create xy grid for interpolation
res=input('Enter Resolution : ');

xmin=min(x);
xmax=max(x);
ymin=min(y);
ymax=max(y);
xx=xmin:res:xmax;
yy=ymin:res:ymax;
[xi, yi] = meshgrid(xx,yy);

%---interpolate elevation then data
zi = griddata(x,y,z,xi,yi,'cubic');
ri = griddata(x,y,r,xi,yi,'cubic');


%---get cells
fprintf('Arranging cells\n')
ncellX = size(zi,1)-1;
ncellY = size(zi,2)-1;

%---we save the data (ri) as point data, but need to know the cell data to
%correctly create the cells, so find average ri for each cell.
rri = meanRes(ri,ncellY,ncellX,1);

%---it may occur that a cell may have a NaN value in just 1-3 vertices and 
%not all. In this case the cell is used so those vertices need an
%elevation.
zi= fillElevation(zi,ncellX,ncellY,0,rri); 
index  = zeros(ncellX+1,ncellY+1);
k      = 0;

%---index matrix. 
for j = 1:ncellY+1
    for i = 1:ncellX+1
        if isfinite(zi(i,j))
        index(i,j) = k;
        k=k+1;
        end
    end
end


%---work through cells. If there is no value for a cell, discard it.
cells = cell(ncellX,ncellY);
for i = 1:ncellX
    for j = 1:ncellY
        if(isfinite(rri(i,j)))
	        cP1=index(i,j);
            cP2=index(i+1,j);
            cP3=index(i+1,j+1);
            cP4=index(i,j+1);
            cells(i,j) = {[4 cP1 cP2 cP3 cP4]};  
        end
    end
end

%---shuffle the data to vtk format.
cells=cells(:);
cells=cell2mat(cells);
ncells = size(cells,1);
xi(isnan(zi))=[];
yi(isnan(zi))=[];
zi(isnan(zi))=[];
points = [xi(:) yi(:) zi(:)];
npoints = length(zi(:));
ri(isnan(ri))=[];

%---write vtk file. remeber to change this to binary for faster reading
fprintf('Writing VTK file\n')
vtk_file=fopen(vtkname,'w');

fprintf(vtk_file,'# vtk DataFile Version 3.0\n3D export of map data\nASCII\nDATASET UNSTRUCTURED_GRID\nPOINTS %d float\n',npoints);
fprintf(vtk_file,'%i %i %.1f\n',points.');
fprintf(vtk_file,'CELLS %i %i\n',ncells,ncells*5);
fprintf(vtk_file,[repmat('%i ', 1, size(cells, 2)), '\n'],cells.');
fprintf(vtk_file,'CELL_TYPES %i\n',ncells);
fprintf(vtk_file,'%i\n',(ones(ncells,1)*9));
fprintf(vtk_file,'POINT_DATA %u \nSCALARS data float 1 \nLOOKUP_TABLE default \n ',npoints);
fprintf(vtk_file,[repmat('%.2f ', 1, size(ri, 2)), '\n'],ri.');
fclose all;
fprintf('VTK file written\n')
toc

Contact us