aboutsummaryrefslogtreecommitdiff
path: root/pkg/plot/getdata.x
blob: 2e64b184bf53170ad42919db25199a470e10dce5 (plain) (blame)
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
include <imhdr.h>

# PLT_GETDATA -- Decrease resolution of image by either subsampling
# or block averaging.  A pointer to the data values to be plotted is
# returned, as well as the number of columns and lines in the data matrix.

pointer procedure plt_getdata (im, sub, pre, xres, yres, nx, ny)

pointer	im			# Pointer to image header
bool	sub			# Subsample versus block average (yes/no)?
bool	pre			# Preserve aspect ratio (yes/no)?
int	xres, yres		# desired resolution
int	nx, ny			# dimensions of output array

int	nxin, nyin
pointer subras, data
pointer plt_blkaverage(), plt_subsample(), imgs2r()
errchk	plt_blkaverage, plt_subsample, calloc

begin
	# First, determine if the image resolution is to be decreased or not.
	nxin = IM_LEN(im,1)
	nyin = IM_LEN(im,2)

	if ((nxin > xres && xres != 0) || (nyin > yres && yres != 0)) {
	    # Need to decrease resolution; image can be either subsampled
	    # or block averaged.

	    if (sub) {
	        data = plt_subsample (im, xres, yres, pre, nx, ny)
		return (data)
	    } else { 
	        data = plt_blkaverage (im, xres, yres, pre, nx, ny)
		return (data)
	    }

	} else {
	    # Return entire image as data matrix.
	    nx = nxin
	    ny = nyin
	    call calloc (data, nx * ny, TY_REAL)
	    subras = imgs2r (im, 1, nxin, 1, nyin)
	    call amovr (Memr[subras], Memr[data], nx * ny)
	    return (data)
	}
end


# PLT_SUBSAMPLE -- The input data array is subsampled by the specified
# factors.  The reduced data array and its dimensions are returned.

pointer procedure plt_subsample (im, xres, yres, pre, nxout, nyout)

pointer	im			# input image
int	xres, yres		# desired output resolution
bool	pre			# preserve aspect ratio?
int	nxout, nyout		# dimensions of output array (returned)

pointer	sp, xvec, data
int	x_factor, y_factor, xrf, yrf
int	nxin, nyin, yin, ii, jj, index, nop
errchk	imgl2r, calloc
pointer	imgl2r()

begin
	call smark (sp)

	nxin = IM_LEN(im,1)
	nyin = IM_LEN(im,2)

	# User could have disabled subsampling in x or y, but not both.
	if (xres == 0)
	    xres = nxin
	if (yres == 0)
	    yres = nyin

	x_factor = (nxin + xres - 1) / xres
	y_factor = (nyin + yres - 1) / yres
	if (pre) {
	    xrf = max (x_factor, y_factor)
	    yrf = max (x_factor, y_factor)
	} else {
	    xrf = x_factor
	    yrf = y_factor
	}

	nxout = nxin / xrf
	nyout = nyin / yrf

	call eprintf ("Image will be subsampled by %d in x and %d in y\n")
	    call pargi (xrf)
	    call pargi (yrf)

	call salloc (xvec, nxin, TY_REAL)
	call calloc (data, nxout * nyout, TY_REAL)

	yin = 1
	do jj = 1, nyout {
	    call amovr (Memr[imgl2r (im, yin)], Memr[xvec],  nxin)
	    nop = 1
	    do ii = 1, nxin {
		index = data + ((jj-1) * nxout) + nop - 1
		if (mod (ii-1, xrf) == 0) {
		    Memr[index] = Memr[xvec+ii-1]
		    nop = nop + 1
		}
	    }
	    yin = yin + yrf
	}

	call sfree (sp)
	return (data)
end


# PLT_BLKAVERAGE -- Block average the input image by the specified reduction
# factors.  The reduced array and its dimensions are returned.

pointer procedure plt_blkaverage (im, xres, yres, pre, nx, ny)

pointer	im			# input image
int	xres, yres		# blocking factors
bool	pre			# preserve aspect ratio?
int	nx, ny			# dimensions of output array (returned)

real	sum
pointer	sp, xvec, data
int	nxin, nyin, nxout, nyout, nxout_full, nyout_full
int	jj, ii, index, nxcols, yin, nxrows, bfx, bfy, x_factor, y_factor
errchk 	abavr, aclrr, aaddr, salloc, calloc, imgl2r
pointer	imgl2r()

begin
	call smark (sp)

	nxin = IM_LEN(im,1)
	nyin = IM_LEN(im,2)

	# User could have disabled blockaveraging in x or y, but not both.
	if (xres == 0)
	    xres = nxin
	if (yres == 0)
	    yres = nyin

	bfx = nxin / xres
	bfy = nyin / yres
	if (pre) {
	    x_factor = max (bfx, bfy)
	    y_factor = max (bfx, bfy)
	} else{
	    x_factor = bfx
	    y_factor = bfy
	}

	call eprintf ("Image will be block averaged by %d in x and %d in y\n")
	    call pargi (x_factor)
	    call pargi (y_factor)

	nxout = (nxin + x_factor - 1) / x_factor
	nyout = (nyin + y_factor - 1) / y_factor
	nxout_full = nxin / x_factor
	nyout_full = nyin / y_factor
	nxcols = nxin - (nxout_full * x_factor)
	nxrows = nyin - (nyout_full * y_factor)

	call salloc (xvec, nxin, TY_REAL)
	call calloc (data, nxout * nyout, TY_REAL)

	yin = 1
	do jj = 1, nyout_full {
	    call aclrr (Memr[xvec], nxin)
	    do ii = 1, y_factor {
		call aaddr (Memr[imgl2r(im,yin)], Memr[xvec], Memr[xvec], nxin)
		yin = yin + 1
	    }
	    call adivkr (Memr[xvec], real(y_factor), Memr[xvec], nxin)

	    index = data + (jj-1) * nxout 
	    call abavr (Memr[xvec], Memr[index], nxout_full, x_factor) 

	    if (nxcols != 0) {
		sum = 0.0
		# deal with trailing column pixels
		do ii = 1, nxcols
		    sum = sum + Memr[xvec+(nxout_full*x_factor)+ii-1]
		Memr[index+nxout_full] = sum / real (nxcols)
	    }

	}

	if (nxrows != 0) {
	    call aclrr (Memr[xvec], nxin)
	    do ii = yin, nyin 
		call aaddr (Memr[imgl2r(im,ii)], Memr[xvec], Memr[xvec], nxin)
	    call adivkr (Memr[xvec], real (nxrows), Memr[xvec], nxin)

	    index = data + (nyout - 1) * nxout
	    call abavr (Memr[xvec], Memr[index], nxout_full, x_factor)

	    if (nxcols != 0) {
	        # Deal with trailing column pixels.
	        do ii = 1, nxcols
		    sum = sum + Memr[xvec+(nxout_full*x_factor)+ii-1]
		Memr[index+nxout_full] = sum / real (nxcols)
	    }
	}

	nx = nxout
	ny = nyout
	call sfree (sp)
	return (data)
end