-
Notifications
You must be signed in to change notification settings - Fork 65
/
plotimage2.m
135 lines (130 loc) · 3.17 KB
/
plotimage2.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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
%-*-octave-*--
%
% PLOTIMAGE2
%
% plotimage2( fil, machineformat, cvals )
%
% Plots the image on file 'fil' with contour, using the contour levels cvals.
% If cvals is ommitted, 21 countour levels will be obtained through the imageinfo fcn.
%
% The boundary of each grid patch is outlined in black. A vector cvals can be
% obtained from function imageinfo.
%
% Input:
% fil: Name of image file
% machineformat (optional): Passed to fopen to read big endian, little endian, etc
% cvals (optional): Vector of countour levels to plot
%
function plotimage2( fil, machineformat, cvals )
if nargin < 2
machineformat='native';
end;
if nargin < 3
nc=21;
[cv0, np, plane, mx, mn] = imageinfo(fil,nc,0,machineformat);
% printf("data min=%e, max=%e\n", mn, mx);
mn = max(0.006,mx/20);
if mx<5*mn
mx = 5*mn;
end
cvals = linspace(mn,mx,nc);
else
mn = cvals(1);
mx = cvals(end);
end;
fd=fopen(fil,'r',machineformat);
pr=fread(fd,1,'int');
nb=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');
fclose(fd);
npatches = nb;
firstCurviPatch = npatches - gridinfo + 1;
x1min=1e9;
x2min=1e9;
x1max=-1e9;
x2max=-1e9;
for b=1:nb
[im,x,y,z] = readimage(fil,b,0,machineformat);
if plane==0
contour(y,z,im,cvals);
elseif plane==1
contour(x,z,im,cvals);
elseif plane==2
contour(x,y,im,cvals);
end
if b==1
hold on;
end;
[n m]=size(im);
% make a frame
if plane==0
if (b >= firstCurviPatch)
plot(y(1,:), z(1,:),'k')
plot(y(n,:), z(n,:),'k')
plot(y(:,1), z(:,1),'k')
plot(y(:,m), z(:,m),'k')
x1mi=min(min(y));
x1ma=max(max(y));
x2mi=min(min(z));
x2ma=max(max(z));
else
plot(y, z(1)*ones(size(y)),'k')
plot(y, z(n)*ones(size(y)),'k')
plot(y(1)*ones(size(z)), z,'k')
plot(y(m)*ones(size(z)), z,'k')
x1mi=min(y);
x1ma=max(y);
x2mi=min(z);
x2ma=max(z);
end;
elseif plane==1
% if (b==nb) && (gridinfo == 1) % update
if (b >= firstCurviPatch)
plot(x(1,:), z(1,:),'k')
plot(x(n,:), z(n,:),'k')
plot(x(:,1), z(:,1),'k')
plot(x(:,m), z(:,m),'k')
x1mi=min(min(x));
x1ma=max(max(x));
x2mi=min(min(z));
x2ma=max(max(z));
else
plot(x, z(1)*ones(size(x)),'k')
plot(x, z(n)*ones(size(x)),'k')
plot(x(1)*ones(size(z)), z,'k')
plot(x(m)*ones(size(z)), z,'k')
x1mi=min(x);
x1ma=max(x);
x2mi=min(z);
x2ma=max(z);
end;
elseif plane==2
plot(x, y(1)*ones(size(x)),'k')
plot(x, y(n)*ones(size(x)),'k')
plot(x(1)*ones(size(y)), y,'k')
plot(x(m)*ones(size(y)), y,'k')
x1mi=min(x);
x1ma=max(x);
x2mi=min(y);
x2ma=max(y);
end;
% update global min/max
x1min = min(x1min,x1mi);
x1max = max(x1max,x1ma);
x2min = min(x2min,x2mi);
x2max = max(x2max,x2ma);
end;
axis([x1min x1max x2min x2max]);
hold off;
axis ij; % flip z-axis to point downwards
axis equal;
%caxis([0.0, mx]);
caxis([mn, mx]);
colorbar;
set(gca,"fontsize",16);
#title_str=sprintf("Time=%g", t);
#title(title_str);