From 0b00e11ffb855253501d3b1ed1e0e2fcc79635f3 Mon Sep 17 00:00:00 2001
From: Yoann Guilhem <yoann.guilhem@esrf.fr>
Date: Fri, 7 Dec 2012 11:06:05 +0000
Subject: [PATCH] Update and add some VTK utility functions
git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@947 4c865b51-4357-4376-afb4-474e03ccb993
---
zUtil_VTK/gtVTKConvertTypeVTK2Matlab.m | 50 ++++++++++
zUtil_VTK/gtVTKMeshReader.m | 52 ++++++----
zUtil_VTK/gtVTKMeshWriter.m | 35 ++++---
zUtil_VTK/gtVTKStructGridWriter.m | 131 +++++++++++++++++++++++++
4 files changed, 234 insertions(+), 34 deletions(-)
create mode 100644 zUtil_VTK/gtVTKConvertTypeVTK2Matlab.m
create mode 100644 zUtil_VTK/gtVTKStructGridWriter.m
diff --git a/zUtil_VTK/gtVTKConvertTypeVTK2Matlab.m b/zUtil_VTK/gtVTKConvertTypeVTK2Matlab.m
new file mode 100644
index 00000000..cf484518
--- /dev/null
+++ b/zUtil_VTK/gtVTKConvertTypeVTK2Matlab.m
@@ -0,0 +1,50 @@
+function [mType, inFMT, inType] = gtVTKConvertTypeVTK2Matlab(vtkType)
+% GTVTKCONVERTTYPEVTK2MATLAB Convert VTK datatype to Matlab datatype.
+
+switch vtkType
+case 'unsigned_char'
+ mType = 'uint8';
+ inFMT = '%d';
+ inType = 'uchar';
+case 'unsigned_short'
+ mType = 'uint16';
+ inFMT = '%d';
+ inType = 'ushort';
+case 'unsigned_int'
+ mType = 'uint32';
+ inFMT = '%d';
+ inType = 'uint';
+case 'unsigned_long'
+ mType = 'uint64';
+ inFMT = '%d';
+ inType = 'uint64';
+case 'char'
+ mType = 'int8';
+ inFMT = '%d';
+ inType = 'schar';
+case 'short'
+ mType = 'int16';
+ inFMT = '%d';
+ inType = 'short';
+case 'int'
+ mType = 'int32';
+ inFMT = '%d';
+ inType = 'int';
+case 'long'
+ mType = 'int64';
+ inFMT = '%d';
+ inType = 'int64';
+case 'float'
+ mType = 'single';
+ inFMT = '%f';
+ inType = 'float32';
+case 'double'
+ mType = 'double';
+ inFMT = '%f';
+ inType = 'double';
+ otherwise
+ gtError('gtVTKConvertTypeVTK2Matlab:wrong_datatype', ...
+ ['Unimplemented datatype "', vtkType, '", sorry...']);
+end
+
+end % end of function
diff --git a/zUtil_VTK/gtVTKMeshReader.m b/zUtil_VTK/gtVTKMeshReader.m
index 48cf1bd3..ae3b77eb 100644
--- a/zUtil_VTK/gtVTKMeshReader.m
+++ b/zUtil_VTK/gtVTKMeshReader.m
@@ -14,6 +14,18 @@ function mesh = gtVTKMeshReader(filename, varargin)
%
% Version 001 20-10-2012 by YGuilhem yoann.guilhem@esrf.fr
+%% Guessing output directory, filename and extension
+[fpath, fname, fext] = fileparts(filename);
+if isempty(fext)
+ fext = '.vtk';
+end
+% Trick to handle ~ path in unix
+if ~isempty(fpath) && isunix && fpath(1)=='~'
+ fpath(1) = '';
+ fpath = fullfile(getenv('HOME'), fpath);
+end
+filename = fullfile(fpath, [fname fext]);
+
fid = fopen(filename, 'r');
if fid ==-1
gtError('gtVTKMeshReader:bad_input_file', 'Can''t open the file.');
@@ -24,7 +36,7 @@ if ~strcmp(file.desc(3:5), 'vtk')
gtError('gtVTKMeshReader:bad_input_file', 'Not a valid VTK file.');
end
-% Read header, mode and dataset
+%% Read header, mode and dataset
file.header = fgets(fid);
str = fgets(fid);
file.mode = sscanf(str, '%s', 1);
@@ -39,7 +51,7 @@ if ~strcmp(file.dataset, 'POLYDATA')
'This vtk file does not contain PolyData dataset!');
end
-% Read number of vertices and coordinate datatype
+%% Read number of vertices and coordinate datatype
str = fgets(fid);
nVert = sscanf(str, 'POINTS %d', 1);
file.nVertices = nVert;
@@ -50,7 +62,8 @@ if binary
case {'float', 'double'}
coordFMT = coordType;
otherwise
- true;
+ gtError('gtVTKMeshReader:wrong_coord_type', ...
+ ['cannot handle coordinate with type: ' coordType]);
end
else
switch coordType
@@ -63,41 +76,42 @@ else
end
end
-% Read vertices
+%% Read vertices
nCoord = 3*nVert;
if binary
coord = fread(fid, nCoord, ['3*' coordFMT '=>' coordFMT], 0 , 'ieee-be');
else
[coord, sz] = fscanf(fid, '%f %f %f', nCoord);
if sz ~= nCoord
- warning('Problem in reading vertices.');
+ warning('gtVTKMeshReader:vertice_warning', ...
+ 'Problem in reading vertices.');
end
end
mesh.vertices = reshape(coord, 3, nVert)';
-% Read polygons
+%% Read polygons
str = fgets(fid);
while str(1) ~= 'P' %&& str(1) ~= 'V'
str = fgets(fid);
end
if str(1) == 'P'
- numbers = sscanf(str,'%*s %d %d', 2);
+ numbers = sscanf(str, '%*s %d %d', 2);
nFace = numbers(1);
nValue = numbers(2);
nNode = 3*nFace;
if binary
nbytes = 4;
- type = 'uint32';
fseek(fid, nbytes, 0);
connectivity = fread(fid, nValue, '3*int32=>uint32', nbytes, 'ieee-be');
else
- connectivity = uint32(fscanf(fid, '%*d %d %d %d\n', nNode));
- %[connectivity, cnt] = fscanf(fid, '%*d %d %d %d\n', nNode);
- %if cnt ~= nNode
- % warning('Problem in reading faces.');
- %end
- %connectivity = uint32(connectivity);
+ [connectivity, cnt] = fscanf(fid, '%*d %d %d %d\n', nNode);
+ if cnt ~= nNode
+ warning('gtVTKMeshReader:face_warning', ...
+ 'Problem in reading faces.');
+ end
+ connectivity = uint32(connectivity);
+ class(connectivity)
end
mesh.faces = reshape(connectivity+1, 3, nFace)';
else
@@ -106,21 +120,23 @@ end
% read vertex ids
%if info == 'V'
-% nv = sscanf(str,'%*s %d %*s', 1);
+% nv = sscanf(str, '%*s %d %*s', 1);
-% [A,cnt] = fscanf(fid,'%d %d \n', 2*nv);
+% [A, cnt] = fscanf(fid, '%d %d \n', 2*nv);
% if cnt~=2*nv
-% warning('Problem in reading faces.');
+% warning('gtVTKMeshReader:vertex_warning', ...
+% 'Problem in reading vertex.');
% end
% A = reshape(A, 2, cnt/2);
-% faces = repmat(A(2,:)+1, 3, 1);
+% faces = repmat(A(2, :)+1, 3, 1);
%end
%if str(1) ~= 'P' %&& str(1) ~= 'V'
% faces = 0;
%end
+%% Close the file end terminate
fclose(fid);
end % end of function
diff --git a/zUtil_VTK/gtVTKMeshWriter.m b/zUtil_VTK/gtVTKMeshWriter.m
index 1fb68a62..da60e742 100644
--- a/zUtil_VTK/gtVTKMeshWriter.m
+++ b/zUtil_VTK/gtVTKMeshWriter.m
@@ -34,7 +34,7 @@ par.facefield = [];
par.facecolor = [];
par = parse_pv_pairs(par, varargin);
-% Check coordinate data type
+%% Check coordinate data type
switch par.coordType
case {'float', 'double'}
coordFMT = '%f %f %f\n';
@@ -43,23 +43,23 @@ otherwise
['Unsupported coordinates datatype: ' par.coordType]);
end
-% Get file parts and check extension
+%% Get file parts and check extension
[fpath, fname, fext] = fileparts(filename);
if isempty(fext)
fext = '.vtk';
+ filename = fullfile(fpath, [fname fext]);
elseif ~strcmp(fext, '.vtk')
gtError('gtVTKMeshWriter:wrong_file_extension', ...
'Output file extension should be ''.vtk''!');
- fext = '.vtk';
end
-% Open file for writing
+%% Open file for writing
fid = fopen(filename, 'w');
if fid ==-1
gtError('gtVTKMeshWriter:bad_output_file', 'Can''t open the file.');
end
-% Header
+%% Header
fprintf(fid, '# vtk DataFile Version 3.0\n');
fprintf(fid, '%s\n', par.header);
if par.binary
@@ -69,20 +69,20 @@ else
end
fprintf(fid, 'DATASET POLYDATA\n');
-% Points / Vertices
+%% Points / Vertices
fprintf(fid, 'POINTS %d %s\n', size(mesh.vertices, 1), par.coordType);
if par.binary
count = fwrite(fid, mesh.vertices', par.coordType, 0, 'ieee-be');
- %if count ~= size(mesh.vertices, 1)
if count ~= numel(mesh.vertices)
- warning('Problem in writing vertices.');
+ warning('gtVTKMeshWriter:vertice_warning', ...
+ 'Problem in writing vertices.');
end
fprintf(fid, '\n');
else
fprintf(fid, coordFMT, mesh.vertices');
end
-% Polygons / Faces / Cells
+%% Polygons / Faces / Cells
nFaces = size(mesh.faces, 1);
nValues = 4*nFaces;
fprintf(fid, 'POLYGONS %d %d\n', nFaces, nValues);
@@ -90,14 +90,14 @@ output = [3*ones(1, nFaces) ; mesh.faces'-1];
if par.binary
count = fwrite(fid, output, 'uint32', 0, 'ieee-be');
if count ~= nValues
- warning('Problem in writing faces.');
+ warning('gtVTKMeshWriter:face_warning', 'Problem in writing faces.');
end
fprintf(fid, '\n');
else
fprintf(fid, '%d %d %d %d\n', output);
end
-% Face field
+%% Face field
if ~isempty(par.facefield)
for ifield=1:length(par.facefield)
fField = par.facefield(ifield);
@@ -127,7 +127,8 @@ if ~isempty(par.facefield)
if par.binary
count = fwrite(fid, fField.data, outType, 0, 'ieee-be');
if count ~= numel(fField.data)
- warning(['Problem in writing field of face field ''' fField.name '']);
+ warning('gtVTKMeshWriter:face_warning', ...
+ ['Problem in writing field of face field ''' fField.name '']);
end
fprintf(fid, '\n');
else
@@ -140,7 +141,8 @@ if ~isempty(par.facefield)
if par.binary
count = fwrite(fid, uint8(255*output), 'uchar', 0, 'ieee-be');
if count ~= numel(output)
- warning(['Problem in writing lookup table of face field ''' fField.name '']);
+ warning('gtVTKMeshWriter:lookup_table_warning', ...
+ ['Problem in writing lookup table of face field ''' fField.name '']);
end
fprintf(fid, '\n');
else
@@ -150,7 +152,7 @@ if ~isempty(par.facefield)
end
end
-% Face color
+%% Face color
if ~isempty(par.facecolor)
for ifield=1:length(par.facecolor)
fField = par.facecolor(ifield);
@@ -171,7 +173,8 @@ if ~isempty(par.facecolor)
if par.binary
count = fwrite(fid, uint8(255*output), 'uchar', 0, 'ieee-be');
if count ~= numel(output)
- warning(['Problem in writing field of face field ''' fField.name '']);
+ warning('gtVTKMeshWriter:face_field_warning', ...
+ ['Problem in writing field of face field ''' fField.name '']);
end
fprintf(fid, '\n');
else
@@ -180,7 +183,7 @@ if ~isempty(par.facecolor)
end
end
-% Close the file end terminate
+%% Close the file end terminate
fclose(fid);
end % end of function
diff --git a/zUtil_VTK/gtVTKStructGridWriter.m b/zUtil_VTK/gtVTKStructGridWriter.m
new file mode 100644
index 00000000..6a949395
--- /dev/null
+++ b/zUtil_VTK/gtVTKStructGridWriter.m
@@ -0,0 +1,131 @@
+function gtVTKStructGridWriter(vol, filename, varargin)
+% GTVTKSTRUCTGRIDWRITER Write a volume in VTK file.
+%
+% gtVTKStructGridWriter(vol, filename, varargin)
+% --------------------------------------------------------------------------
+%
+% INPUT:
+% vol =
+% filename = <string> Path to the output VTK file
+%
+% OPTIONAL INPUT:
+% 'binary' = <bool> Write output in binary format (true by default)
+% 'spacing' = <double> Vector defining spacing in each space dimension
+% ([1 1 1] by default)
+% 'origin' = <double> Vector defining volume origin
+% ([0 0 0] by default)
+% 'header' = <string> VTK file header description
+% 'datatype' = <string> Data type to consider for input volume
+% 'cmap' = <colormap>
+%
+%
+% Version 001 20-10-2012 by YGuilhem yoann.guilhem@esrf.fr
+
+par.spacing = ones(1, 3);
+par.origin = zeros(1, 3);
+par.dataname = 'Values';
+par.header = 'VTK structured volume created by matlabDCT';
+par.datatype = '';
+par.cmap = [];
+par.opacity = 1;
+par.verbose = false;
+par.binary = true;
+par = parse_pv_pairs(par, varargin);
+
+if iscell(vol)
+ error('Not implemented: only one volume can be handled by now!');
+else
+ [nx, ny, nz] = size(vol);
+end
+ntot = nx*ny*nz;
+
+if ~isempty(par.cmap) && size(par.cmap, 1) < max(max(max(vol)))
+ warning('gtVTKStructGridWriter:cmap_warning', ...
+ 'You may have not enough values in your cmap...');
+end
+
+%% Guess input & output datatype
+if isempty(par.datatype)
+ par.datatype = class(vol);
+end
+
+[outputTypeVTK, outputType, outputFMT] = gtVTKConvertTypeMatlab2VTK(par.datatype);
+
+if par.verbose
+ disp(['Spacing : ' num2str(par.spacing)]);
+ disp(['Origin : ' num2str(par.origin)]);
+ disp(['Color map : ' num2str(par.cmap)]);
+ disp(['Dataname : ' par.dataname]);
+ disp(['Header : ' par.header]);
+ disp(['OutputType : ' outputType]);
+ disp(['OutputTypeVTK : ' outputTypeVTK]);
+end
+
+%% Open the file
+fid = fopen(filename, 'w');
+if fid ==-1
+ gtError('gtVTKStructGridWriter:bad_output_file', 'Can''t open the file.');
+end
+
+%% Write ASCII file header
+fprintf(fid, '# vtk DataFile Version 3.0\n');
+fprintf(fid, '%s\n', par.header);
+if par.binary
+ fprintf(fid, 'BINARY\n');
+else
+ fprintf(fid, 'ASCII\n');
+end
+fprintf(fid, 'DATASET STRUCTURED_POINTS\n');
+fprintf(fid, 'DIMENSIONS %d %d %d\n', nx, ny, nz);
+fprintf(fid, 'ORIGIN %4.3f %4.3f %4.3f\n', par.origin);
+fprintf(fid, 'SPACING %4.3f %4.3f %4.3f\n', par.spacing);
+
+fprintf(fid, 'POINT_DATA %d\n', ntot);
+fprintf(fid, 'SCALARS %s %s\n', par.dataname, outputTypeVTK);
+
+if ~isempty(par.cmap)
+ fprintf(fid, 'LOOKUP_TABLE map\n');
+else
+ fprintf(fid, 'LOOKUP_TABLE default\n');
+end
+
+%% Write volume data
+if par.binary
+ count = fwrite(fid, reshape(vol, 1, ntot), par.datatype, 0, 'ieee-be');
+ if count ~= numel(vol)
+ warning('gtVTKStructGridWriter:structgrid_volume_warning', ...
+ 'Problem in writing structured grid volume.');
+ end
+ fprintf(fid, '\n');
+else
+ fprintf(fid, [outputFMT '\n'], reshape(vol, 1, ntot));
+end
+
+%% Write lookup table if needed
+if ~isempty(par.cmap)
+ mapSize = size(par.cmap);
+ if mapSize(2) == 3
+ output = [par.cmap'; par.opacity*ones(1, mapSize(1))];
+ elseif mapSize(2) == 4
+ output = par.cmap';
+ else
+ gtError('gtVTKStructGridWriter:wrong_cmap_size', ...
+ 'Color map sould be N row vector of size 3 or 4!');
+ end
+ fprintf(fid, 'LOOKUP_TABLE map %d\n', mapSize(1));
+ if par.binary
+ count = fwrite(fid, uint8(225*output), 'uchar', 0, 'ieee-be');
+ if count ~= numel(output)
+ warning('gtVTKMeshWriter:cmap_warning', ...
+ 'Problem in writing color map!');
+ end
+ fprintf(fid, '\n');
+ else
+ fprintf(fid, '%f %f %f %f\n', output);
+ end
+end
+
+%% Close file and terminate
+fclose(fid);
+
+end % end of function
--
GitLab