diff options
Diffstat (limited to 'pkg/plot/getdata.x')
| -rw-r--r-- | pkg/plot/getdata.x | 212 |
1 files changed, 212 insertions, 0 deletions
diff --git a/pkg/plot/getdata.x b/pkg/plot/getdata.x new file mode 100644 index 00000000..2e64b184 --- /dev/null +++ b/pkg/plot/getdata.x @@ -0,0 +1,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 |