-
Notifications
You must be signed in to change notification settings - Fork 65
/
readimage3d.m
90 lines (89 loc) · 3.18 KB
/
readimage3d.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
%
% READIMAGE3D
%
% Read volumeimage produced by sw4 on format with the *.sw4img extension.
% This script will not read the older WPP image format.
%
% [im,x,y,z,t,timestring]=readimage3d( imfile, pnr, verbose )
%
% Input: imfile - Image file name
% pnr - Patch number, if more than one grid is used.
% verbose - Set to 1 to display file header information.
% Output: im - The image, as a 2D array.
% x, y, z - The spatial coordinates of the image, one of these is a scalar.
% t - Simulation time at which the image was output.
% timestring - String holding creation date
function [im,x,y,z,t,timestring]=readimage3d( imfile, pnr, verbose )
if nargin < 3
verbose= 0;
end;
if nargin < 2
pnr = 1;
end;
fd=fopen(imfile,'r');
if fd ~= -1
% Read header
prec =fread(fd,1,'int');
npatches=fread(fd,1,'int');
t =fread(fd,1,'double');
plane =fread(fd,1,'int');
coord =fread(fd,1,'double');
mode =fread(fd,1,'int');
gridinfo=fread(fd,1,'int');
timecreated=fread(fd,[1 25],'uchar');
timestring=num2str(timecreated,'%c');
mstr=getimagemodestr3d(mode);
% Display header
if verbose == 1
disp(['Found: prec = ' num2str(prec) ' t= ' num2str(t) ' plane= ' num2str(plane)]);
disp([' coord = ' num2str(coord) ' mode= ' mstr ' npatches= ' num2str(npatches) ]);
disp([' gridinfo = ' num2str(gridinfo) ]);
disp([' file created ' timecreated(1:24)]);
end;
for p=1:npatches
h(p) = fread(fd,1,'double');
zmin(p) = fread(fd,1,'double');
ib(p) = fread(fd,1,'int');
ni(p) = fread(fd,1,'int');
jb(p) = fread(fd,1,'int');
nj(p) = fread(fd,1,'int');
kb(p) = fread(fd,1,'int');
nk(p) = fread(fd,1,'int');
if verbose == 1
disp([' patch nr ' num2str(p) ' has h = ' num2str(h(p)) ' zmin = ' num2str(zmin(p))]);
end;
end;
% Read data
readz = 0;
if pnr <= npatches
for p=1:pnr-1
fseek(fd,(ni(p)-ib(p)+1)*(nj(p)-jb(p)+1)*(nk(p)-kb(p)+1)*prec,'cof');
end;
if prec == 4
im0 = fread(fd,(ni(pnr)-ib(pnr)+1)*(nj(pnr)-jb(pnr)+1)*(nk(pnr)-kb(pnr)+1),'float');
else
im0 = fread(fd,(ni(pnr)-ib(pnr)+1)*(nj(pnr)-jb(pnr)+1)*(nk(pnr)-kb(pnr)+1),'double');
end;
im = reshape(im0,ni(pnr)-ib(pnr)+1,nj(pnr)-jb(pnr)+1,nk(pnr)-kb(pnr)+1);
% read grid z-coordinates
if (pnr == npatches) && (gridinfo == 1 )
if prec == 4
z0 = fread(fd,(ni(pnr)-ib(pnr)+1)*(nj(pnr)-jb(pnr)+1)*(nk(pnr)-kb(pnr)+1),'float');
else
z0 = fread(fd,(ni(pnr)-ib(pnr)+1)*(nj(pnr)-jb(pnr)+1)*(nk(pnr)-kb(pnr)+1),'double');
end;
z = reshape(z0,ni(pnr)-ib(pnr)+1,nj(pnr)-jb(pnr)+1,nk(pnr)-kb(pnr)+1);
readz = 1;
end;
x = (ib(pnr)-1:ni(pnr)-1)*h(pnr);
y = (jb(pnr)-1:nj(pnr)-1)*h(pnr);
if readz == 0
z = (kb(pnr)-1:nk(pnr)-1)*h(pnr);
end;
fclose(fd);
else
disp(['Error: number of patches on file ' num2str(npatches) ' is smaller than input pnr = ' num2str(pnr)]);
end;
else
disp(['Error: could not open file ' imfile ]);
end;