diff options
| author | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
|---|---|---|
| committer | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
| commit | fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch) | |
| tree | bdda434976bc09c864f2e4fa6f16ba1952b1e555 /pkg/images/imgeom | |
| download | iraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz | |
Initial commit
Diffstat (limited to 'pkg/images/imgeom')
46 files changed, 8115 insertions, 0 deletions
diff --git a/pkg/images/imgeom/Revisions b/pkg/images/imgeom/Revisions new file mode 100644 index 00000000..5e4ad630 --- /dev/null +++ b/pkg/images/imgeom/Revisions @@ -0,0 +1,2026 @@ +.help revisions Jan97 images.imgeom +.nf + +=============================== +Package Reorganization +=============================== + +pkg/images/imarith/t_imsum.x +pkg/images/imarith/t_imcombine.x +pkg/images/doc/imsum.hlp +pkg/images/doc/imcombine.hlp + Provided options for USHORT data. (12/10/96, Valdes) + +pkg/images/imarith/icsetout.x +pkg/images/doc/imcombine.hlp + A new option for computing offsets from the image WCS has been added. + (11/30/96, Valdes) + +pkg/images/imarith/t_imcombine.x +pkg/images/imarith/icombine.gx + Changed the error checking to catch additional errors relating to too + many files. (11/12/96, Valdes) + +pkg/images/imarith/icsort.gx + There was an error in the ic_2sort routine when there are exactly + three images that one of the explicit cases did not properly keep + the image identifications. See buglog 344. (8/1/96, Valdes) + +pkg/images/filters/median.x + The routine mde_yefilter was being called with the wrong number of + arguments. + (7/18/96, Davis) + +pkg/images/imarith/t_imcombine.x +pkg/images/imarith/icombine.gx +pkg/images/imarith/icimstack.x + +pkg/images/imarith/iclog.x +pkg/images/imarith/mkpkg +pkg/images/doc/imcombine.hlp + The limit on the maximum number of images that can be combined, set by + the maximum number of logical file descriptors, has been removed. If + the condition of too many files is detected the task now automatically + stacks all the images in a temporary image and then combines them with + the project option. + (5/14/96, Valdes) + +pkg/images/geometry/xregister/rgxfit.x + Changed several Memr[] references to Memi[] in the rg_fit routine. + This bug was causing a floating point error in the xregister task + on the Dec Alpha if the coords file was defined, and could potentially + cause problems on other machines. + (Davis, April 3, 1996) + +pkg/images/geometry/t_geotran.x +pkg/images/geometry/geograph.x +pkg/images/doc/geomap.hlp + Corrected the definition of skew in the routines which compute a geometric + interpretation of the 6-coefficient fit, which compute the coefficients + from the geometric parameters, and in the relevant help pages. + (2/19/96, Davis) + +pkg/images/median.par +pkg/images/rmedian.par +pkg/images/mode.par +pkg/images/rmode.par +pkg/images/fmedian.par +pkg/images/frmedian.par +pkg/images/fmode.par +pkg/images/frmode.par +pkg/images/doc/median.hlp +pkg/images/doc/rmedian.hlp +pkg/images/doc/mode.hlp +pkg/images/doc/rmode.hlp +pkg/images/doc/fmedian.hlp +pkg/images/doc/frmedian.hlp +pkg/images/doc/fmode.hlp +pkg/images/doc/frmode.hlp +pkg/images/filters/t_median.x +pkg/images/filters/t_rmedian.x +pkg/images/filters/t_mode.x +pkg/images/filters/t_rmode.x +pkg/images/filters/t_fmedian.x +pkg/images/filters/t_frmedian.x +pkg/images/filters/t_fmode.x +pkg/images/filters/t_frmode.x + Added a verbose parameter to the median, rmedian, mode, rmode, fmedian, + frmedian, fmode, and frmode tasks. (11/27/95, Davis) + +pkg/images/geometry/doc/geotran.hlp + Fixed an error in the help page for geotran. The default values for + the xscale and yscale parameters were incorrectly listed as INDEF, + INDEF instead of 1.0, 1.0. (11/14/95, Davis) + +pkg/images/imarith/icpclip.gx + Fixed a bug where a variable was improperly used for two different + purposes causing the algorithm to fail (bug 316). (10/19/95, Valdes) + +pkg/images/doc/imcombine.hlp + Clarified a point about how the sigma is calculated with the SIGCLIP + option. (10/11/95, Valdes) + +pkg/images/imarith/icombine.gx + To deal with the case of readnoise=0. and image data which has points with + negative mean or median and very small minimum readnoise is set + internally to avoid computing a zero sigma and dividing by it. This + applies to the noise model rejection options. (8/11/95, Valdes) + +pkg/images/frmedian.hlp +pkg/images/frmode.hlp +pkg/images/rmedian.hlp +pkg/images/rmode.hlp +pkg/images/frmedian.par +pkg/images/frmode.par +pkg/images/rmedian.par +pkg/images/rmode.par +pkg/images/filters/frmedian.h +pkg/images/filters/frmode.h +pkg/images/filters/rmedian.h +pkg/images/filters/rmode.h +pkg/images/filters/t_frmedian.x +pkg/images/filters/t_frmode.x +pkg/images/filters/t_rmedian.x +pkg/images/filters/t_rmode.x +pkg/images/filters/frmedian.x +pkg/images/filters/frmode.x +pkg/images/filters/rmedian.x +pkg/images/filters/rmode.x +pkg/images/filters/med_utils.x + Added new ring median and modal filtering tasks frmedian, rmedian, + frmode, and rmode to the images package. + (6/20/95, Davis) + +pkg/images/fmedian.hlp +pkg/images/fmode.hlp +pkg/images/median.hlp +pkg/images/mode.hlp +pkg/images/fmedian.par +pkg/images/fmode.par +pkg/images/median.par +pkg/images/mode.par +pkg/images/filters/fmedian.h +pkg/images/filters/fmode.h +pkg/images/filters/median.h +pkg/images/filters/mode.h +pkg/images/filters/t_fmedian.x +pkg/images/filters/t_fmode.x +pkg/images/filters/t_median.x +pkg/images/filters/t_mode.x +pkg/images/filters/fmedian.x +pkg/images/filters/fmode.x +pkg/images/filters/median.x +pkg/images/filters/mode.x +pkg/images/filters/fmd_buf.x +pkg/images/filters/fmd_hist.x +pkg/images/filters/fmd_maxmin.x +pkg/images/filters/med_buf.x +pkg/images/filters/med_sort.x + Added minimum and maximum good data parameters to the fmedian, fmode, + median, and mode filtering tasks. Removed the 64X64 kernel size limit + in the median and mode tasks. Replaced the common blocks with structures + and .h files. + (6/20/95, Davis) + +pkg/images/geometry/t_geotran.x +pkg/images/geometry/geotran.x +pkg/images/geometry/geotimtran.x + Fixed a bug in the buffering of the x and y coordinate surface interpolants + which can cause a memory corruption error if, nthe nxsample or nysample + parameters are > 1, and the nxblock or nyblock parameters are less + than the x and y dimensions of the input image. Took the opportunity + to clean up the code. + (6/13/95, Davis) + +======= +V2.10.4 +======= + +pkg/images/geometry/t_geomap.x + Corrected a harmless typo in the code which determines the minimum + and maximum x values and improved the precision of the test when the + input is double precision. + (4/18/95, Davis) + +pkg/images/doc/fit1d.hlp + Added a description of the interactive parameter to the fit1d help page. + (4/17/95, Davis) + +pkg/images/imarith/t_imcombine.x + If an error occurs while opening an input image header the error + recovery will close all open images and then propagate the error. + For the case of running out of file descriptors with STF format + images this will allow the error message to be printed rather + than the error code. (4/3/95, Valdes) + +pkg/images/geometry/xregister/t_xregister.x + Added a test on the status code returned from the fitting routine so + the xregister tasks does not go ahead and write an output image when + the user quits the task in in interactive mode. + (3/31/95, Davis) + +pkg/images/imarith/icscale.x +pkg/images/doc/imcombine.hlp + The behavior of the weights when using both multiplicative and zero + point scaling was incorrect; the zero levels have to account for + the scaling. (3/27/95, Valdes) + +pkg/images/geometry/xregister/rgxtools.x + Changed some amovr and amovi calls to amovkr and amovki calls. + (3/15/95, Davis) + +pkg/images/geometry/t_imshift.x +pkg/images/geometry/t_magnify.x +pkg/images/geometry/geotran.x +pkg/images/geometry/xregister/rgximshift.x + The buffering margins set for the bicubic spline interpolants were + increased to improve the flux conservation properties of the interpolant + in cases where the data is undersampled. (12/6/94, Davis) + +pkg/images/xregister/rgxbckgrd.x + In several places the construct array[1++nx-wborder] was being used + instead of array[1+nx-wborder]. Apparently caused by a typo which + propagated through the code, the Sun compilers did not catch this, but + the IBM/RISC6000 compilers did. (11/16/94, Davis) + + +pkg/images/xregister.par +pkg/images/doc/xregister.hlp +pkg/images/geometry/xregister/t_xregister.x +pkg/images/geometry/xregister/rgxcorr.x +pkg/images/geometry/xregister/rgxicorr.x +pkg/images/geometry/xregister/rgxcolon.x +pkg/images/geometry/xregister/rgxdbio.x + The xregister task was modified to to write the output shifts file + in either text database format (the current default) or in simple text + format. The change was made so that the output of xregister could + both be edited more easily by the user and be used directly with the + imshift task. (11/11/94, Davis) + +pkg/images/imfit/fit1d.x + A Memc in the ratio output option was incorrectly used instead of Memr + when the bug fix of 11/16/93 was made. (10/14/94, Valdes) + +pkg/images/geometry/xregister/rgxcorr.x + The procedure rg_xlaplace was being incorrectly declared as an integer + procedure. + (8/1/94, Davis) + +pkg/images/geometry/xregister/rgxregions.x + The routine strncmp was being called (with a missing number of characters + argument) instead of strcmp. This was causing a bus error under solaris + but not sun os whenever the user set regions to "grid ...". (7/27/94 LED) + +pkg/images/tv/imexaine/ierimexam.x + The Gaussian fitting can return a negative sigma**2 which would cause + an FPE when the square root is taken. This will only occur when + there is no reasonable signal. The results of the gaussian fitting + are now set to INDEF if this unphysical result occurs. (7/7/94, Valdes) + +pkg/images/geometry/geofit.x + A routine expecting two char arrays was being passed two real arrays + instead resulting in a segmentation violation if calctype=real + and reject > 0. + (6/21/94, Davis) + +pkg/images/imarith/t_imarith.x + IMARITH now deletes the CCDMEAN keyword if present. (6/21/94, Valdes) + +pkg/images/imarith/icaclip.gx +pkg/images/imarith/iccclip.gx +pkg/images/imarith/icpclip.gx +pkg/images/imarith/icsclip.gx + 1. The restoration of deleted pixels to satisfy the nkeep parameter + was being done inside the iteration loop causing the possiblity + of a non-terminating loop; i.e. pixels are rejected, they are + restored, and the number left then does not statisfy the termination + condition. The restoration step was moved following the iterative + rejection. + 2. The restoration was also incorrectly when mclip=no and could + lead to a segmentation violation. + (6/13/94, Valdes) + +pkg/images/geometry/xregister/rgxicorr.x + The path names to the xregister task interactive help files was incorrect. + (6/13/94, Davis) + +pkg/images/imarith/iccclip.gx +pkg/images/imarith/icsclip.gx + Found and fixed another typo bug. (6/7/94, Valdes/Zhang) + +pkg/images/imarith/icscale.x + The sigma scaling flag, doscale1, would not be set in the case of + a mean offset of zero though the scale factors could be different. + (5/25/94, Valdes/Zhang) + +pkg/images/imarith/icsclip.gx + There was a missing line: l = Memi[mp1]. (5/25/94, Valdes/Zhang) + +pkg/images/imarith/icaclip.gx +pkg/images/imarith/iccclip.gx +pkg/images/imarith/icpclip.gx +pkg/images/imarith/icsclip.gx + The reordering step when a central median is used during rejection + but the final combining is average was incorrect if the number + of rejected low pixels was greater than the number of pixel + number of pixels not rejected. (5/25/94, Valdes) + +pkg/images/geometry/t_geotran.x + In cases where there was no input geomap database, geotran was + unnecessarily overiding the size of the input image requested by the + user if the size of the image was bigger than the default output size + (the size of the output image which would include all the input image + pixels is no user shifts were applied). + (5/10/94, Davis) + +pkg/images/imarith/icscale.x +pkg/images/imarith/t_imcombine.x + 1. There is now a warning error if the scale, zero, or weight type + is unknown. + 2. An sfree was being called before the allocated memory was finished + being used. + (5/2/94, Valdes) + +pkg/images/tv/imexaine/ierimexam.x + For some objects the moment analysis could fail producing a floating + overflow error in imexamine, because the code was trying to use + INDEF as the initial value of the object fwhm. Changed the gaussian + fitting code to use a fraction of the fitting radius as the initial value + for the fitted full-width half-maximum in cases where the moment analysis + cannot compute an initial value. + (4/15/94 LED) + +pkg/images/imarith/iclog.x + Changed the mean, median, mode, and zero formats from 6g to 7.5g to + insure 5 significant digits regardless of signs and decimal points. + (4/13/94, Valdes) + +pkg/images/doc/imcombine.hlp + Tried again to clarify the scaling as multiplicative and the offseting + as additive for file input and for log output. (3/22/94, Valdes) + +pkg/images/imarith/iacclip.gx +pkg/images/imarith/iccclip.gx +pkg/images/imarith/iscclip.gx + The image sigma was incorrectly computed when an offset scaling is used. + (3/8/94, Valdes) + +pkg/images/doc/imcombine.hlp + The MINMAX example confused low and high. (3/7/94, Valdes) + +pkg/images/geometry/t_geomap.x +pkg/images/geometry/geofit.x +pkg/images/geometry/geograph.x + Fixed a bug in the geomap code which caused the linear portion of the transformation + to be computed incorrectly if the x and y fits had a different functional form. + (12/29/93, Davis) + +pkg/images/imarith/t_imcombine.x +pkg/images/imcombine.par +pkg/images/do/imcombine.hlp + The output pixel datatypes now include unsigned short integer. + (12/4/93, Valdes) + +pkg/images/doc/imcombine.hlp + Fixed an error in the example of offseting. (11/23/93, Valdes) + +pkg/images/imfit/fit1d.x + When doing operations in place the input and output buffers are the + same and the difference and ratio operations assumed they were not + causing the final results to be wrong. (11/16/93, Valdes) + +pkg/images/imarith/t_imarith.x +pkg/images/doc/imarith.hlp + If no calculation type is specified then it will be at least real + for a division. Since the output pixel type defaults to the + calculation type if not specified this will also result in a + real output if dividing two integer images. (11/12/93, Valdes) + +pkg/images/imarith/icgrow.gx +pkg/images/imarith/icpclip.gx +pkg/images/imarith/icsclip.gx +pkg/images/imarith/icaclip.gx +pkg/images/imarith/iccclip.gx +pkg/images/imarith/t_imcombine.x +pkg/images/doc/imcombine.hlp + If there were fewer initial pixels than specified by nkeep then the + task would attempt to add garbage data to achieve nkeep pixels. This + could occur when using offsets, bad pixel masks, or thresholds. The + code was changed to check against the initial number of pixels rather + than the number of images. Also a negative nkeep is no longer + converted to a positive value based on the number of images. Instead + it specifies the maximum number of pixels to reject from the initial + set of pixels. (11/8/93, Valdes) + +======= +V2.10.2 +======= + +pkg/images/imarith/icsetout.x + Added MWCS calls to update the axis mapping when using the project + option in IMCOMBINE. (10/8/93, Valdes) + +pkg$images/imarith/icscale.x +pkg$images/doc/imcombine.hlp + The help indicated that user input scale or zero level factors + by an @file or keyword are multiplicative and additive while the + task was using then as divisive and subtractive. This was + corrected to agree with the intend of the documentation. + Also the factors are no longer normalized. (9/24/93, Valdes) + +pkg$images/imarith/icsetout.x + The case in which absolute offsets are specified but the offsets are + all the same did not work correctly. (9/24/93, Valdes) + +pkg$images/imfit/imsurfit.h +pkg$images/imfit/t_imsurfit.x +pkg$images/imfit/imsurfit.x +pkg$images/lib/ranges.x + Fixed two bugs in the imsurfit task bad pixel rejection code. For low + k-sigma rejections factors the bad pixel list could overflow resulting + in a segmentation violation or a hung task. Overlapping ranges were + not being decoded into a bad pixel list properly resulting in + oscillating bad pixel rejection behavior where certain groups of + bad pixels were alternately being included and excluded from the fit. + Both bugs are fixed in iraf 2.10.3 + (9/21/93, Davis) + +pkg$images/doc/imcombine.hlp + Clarified how bad pixel masks work with the "project" option. + (9/13/93, Valdes) + +pkg$images/imfit/fit1d.x + When the input and output images are the same there was an typo error + such that the output was opened separately but then never unmapped + resulting in the end of the image not being updated. (8/6/93, Valdes) + +pkg$images/imarith/t_imcombine.x + The algorithm for making sure there are enough file descriptors failed + to account for the need to reopen the output image header for an + update. Thus when the number of input images + output images + logfile + was exactly 60 the task would fail. The update occurs when the output + image is unmapped so the solution was to close the input images first + except for the first image whose pointer is used in the new copy of the + output image. (8/4/93, Valdes) + +pkg$images/filters/t_mode.x +pkg$images/filters/t_median.x + Fixed a bug in the error trapping code in the median and mode tasks. + The call to eprintf contained an extra invalid error code agument. + (7/28/93, Davis) + +pkg$images/geometry/geomap.par +pkg$images/geometry/t_geomap.x +pkg$images/geometry/geogmap.x +pkg$images/geometry/geofit.x + Fixed a bug in the error handling code in geomap which was producing + a segmentation violation on exit if the user's coordinate list + had fewer than 3 data points. Also improved the error messages + presented to the user in both interactive and non-interactive mode. + (7/7/93, Davis) + +pkg$images/imarith/icgdata.gx + There was an indexing error in setting up the ID array when using + the grow option. This caused the CRREJECT/CCDCLIP algorithm to + fail with a floating divide by zero error when there were non-zero + shifts. (5/26/93, Valdes) + +pkg$images/imarith/icmedian.gx + The median calculation is now done so that the original input data + is not lost. This slightly greater inefficiency is required so + that an output sigma image may be computed if desired. (5/10/93, Valdes) + +pkg$images/geometry/t_imshift.x + Added support for type ushort to the imshift task in cases where the + pixel shifts are integral. + (5/8/93, Davis) + +pkg$images/doc/rotate.hlp + Fixed a bug in the rotate task help page which implied that automatic + image size computation would occur if ncols or nlines were set no 0 + instead of ncols and nlines. + (4/17/93, Davis) + +pkg$images/imarith/imcombine.gx + There was no error checking when writing to the output image. If + an error occurred (the example being when an imaccessible imdir was + set) obscure messages would result. Errchks were added. + (4/16/93, Valdes) + +pkg$images/doc/gauss.hlp + Fixed 2 sign errors in the equations in the documentation describing + the elliptical gaussian fucntion. + (4/13/92, Davis) + +pkg/images/imutil/t_imslice.x + Removed an error check in the imslice task, which was preventing it from + being used to reduce the dimensionality of images where the length of + the slice dimension is 1.0. + (2/16/83, Davis) + +pkg/images/filters/fmedian.x + The fmedian task was printing debugging information under iraf 2.10.2. + (1/25/93, Davis) + +pkg/images/imarith/icaclip.gx +pkg/images/imarith/iccclip.gx +pkg/images/imarith/icpclip.gx +pkg/images/imarith/icsclip.gx + When using mclip=yes and when more pixels are rejected than allowed by + the nkeep parameter there was a subtle bug in how the pixels are added + back which can result in a segmentation violation. + if (nh == n2) ==> if (nh == n[i]) + (1/20/93, Valdes) + + +======= +V2.10.1 +======= + +pkg/images/imarith/t_imcombine.x +pkg/images/imarith/icaclip.gx +pkg/images/imarith/iccclip.gx +pkg/images/imarith/icgrow.gx +pkg/images/imarith/iclog.x +pkg/images/imarith/icombine.com +pkg/images/imarith/icombine.gx +pkg/images/imarith/icombine.h +pkg/images/imarith/icpclip.gx +pkg/images/imarith/icscale.x +pkg/images/imarith/icsclip.gx +pkg/images/imarith/icsetout.x +pkg/images/imcombine.par +pkg/images/doc/combine.hlp + The weighting was changed from using the square root of the exposure time + or image statistics to using the values directly. This corresponds + to variance weighting. Other options for specifying the scaling and + weighting factors were added; namely from a file or from a different + image header keyword. The \fInkeep\fR parameter was added to allow + controlling the maximum number of pixels to be rejected by the clipping + algorithms. The \fIsnoise\fR parameter was added to include a sensitivity + or scale noise component to the noise model. Errors will now delete + the output image. + (9/30/92, Valdes) + +pkg/images/imutil/imcopy.x + Added a call to flush after the status line printout so that the output + will appear immediately. (8/19/92, Davis) + +pkg/images/filters/mkpkg +pkg/images/filters/t_fmedian.x +pkg/images/filters/fmedian.x +pkg/images/filters/fmd_buf.x +pkg/images/filters/fmd_maxmin.x + The fmedian task could crash with a segmentation violation if mapping + was turned off (hmin = zmin and hmax = zmax) and the input image + contained data outside the range defined by zmin and zmax. (8/18/92, Davis) + +pkg/images/imarith/icaclip.gx +pkg/images/imarith/iccclip.gx +pkg/images/imarith/icpclip.gx +pkg/images/imarith/icsclip.gx + There was a very unlikely possibility that if all the input pixels had + exactly the same number of rejected pixels the weighted average would + be done incorrectly because the dflag would not be set. (8/11/92, Valdes) + +pkg/images/imarith/icmm.gx + This procedure failed to set the dflag resulting in the weighted average + being computed in correctly. (8/11/92, Valdes) + +pkg/images/imfit/fit1d.x + At some point changes were made but not documented dealing with image + sections on the input/output. The changes seem to have left off the + final step of opening the output image using the appropriate image + sections. Because of this it is an error to use an image section + on an input image when the output image is different; i.e. + + cl> fit1d dev$pix[200:400,*] junk + + This has now been fixed. (8/10/92, Valdes) + +pkg/images/imarith/icscales.x + The zero levels were incorrectly scaled twice. (8/10/92, Valdes) + +pkg/images/imarith/icstat.gx + Contained the statement + nv = max (1., (Memi[v2+i] - Memi[v1+i]) / Memi[dv+i] + 1) + which is max(real,int). Changed the 1. to a 1. (8/10/92, Valdes) + +pkg$images/imarith/icaclip.gx +pkg$images/imarith/iccclip.gx +pkg$images/imarith/icsclip.gx + These files contained multiple cases (ten or so) of constructs such as + "max (1., ...)" or "max (0., ...)" where the ... could be either real + or double. In the double cases the DEC compiler complained about a + type mismatch since 1. is real. (8/10/92, Valdes) + +pkg$images/imfit/t_imsurfit.x + Fixed a bug in the section reading code. Imsurfit is supposed to switch + the order of the section delimiters in x and y if x2 < x1 or y2 < 1. + Unfortunately the y test was actually "if (y2 < x1)" instead of + "if (y2 < y1)". Whether or not the code actually works correctly + depends on the value of x1 relative to y2. This bug was not present + in 2.9.1 but is present in subsequent releases. (7/30/92 LED) + +======= +V2.10.1 +======= + +pkg$images/filters/t_gauss.x + The case theta=90 and ratio > 0.0 but < 1.0 was producing an incorrect + convolution if bilinear=yes, because the major axis sigmas being + input along the x and y axes were sigma and ratio * sigma respectively + instead of ratio * sigma and sigma in this case. + +pkg$images/imutil/imcopy.x + Modified imcopy to write its verbose output to STDOUT instead of + STDERR. (6/24/92, Davis) + +pkg$images/imarith/imcombine.gx + The step where impl1$t is called to check if there is enough memory + did not set the return buffer because the values are irrelevant for + this check. However, depending on history, this buffer could have + arbitrary values and later when IMIO attempts to flush this buffer, + at least in the case of image type coersion, cause arithmetic errors. + The fix was to clear the returned buffers. (4/27/92, Valdes) + +pkg$images/imutil/t_imstack.x + Modified the imslice task to read the old and write a new axis map. + (4/23/92, Davis) + +pkg$images/geometry/t_imslice.x + Modified the imslice task to read the old and write a new axis map. + (4/23/92, Davis) + +pkg$images/geometry/t_blkavg.x +pkg$images/geometry/t_blkrep.x + Modified the calls to mw_shift and mw_scale to explicitly set the + number of logical axes instead of using the default of 0. + (4/23/92, Davis) + +pkg$images/geometry/t_imtrans.x + Modified imtranspose so that it correctly picks up the axis map + and writes it to the output image wcs. (4/23/92, Davis) + +pkg$images/register.par +pkg$images/geotran.par +pkg$images/doc/register.hlp +pkg$images/doc/geotran.hlp + Changed the default values of the parameters xscale and yscale in + the register and geotran tasks from INDEF to 1.0 (4/23/92, Davis) + +pkg$images/geometry/t_imtrans.x +pkg$images/doc/imtranspose.hlp + Modified the imtranspose task so it does a true transpose of the + axes instead of simply modifying the lterm. (4/8/92, Davis) + +pkg$images/iminfo/listpixels.x + Added the formats parameter for formatting the output pixel coordinates + to the listpixels task. These formats take precedence over the formats + stored in the WCS in the image header and the previous default format. + (4/7/92, Davis) + +pkg$images/imutil/t_imstack.x + Added wcs support to the imstack task. (4/2/92, Davis) + +pkg$images/iminfo/listpixels.x + Modified listpixels so that it will work correctly if the dimension + of the wcs is less than the dimension of the image. (3/16/92, Davis) + +pkg$images/geometry/t_geotran.x + Modified the rotate, imlintran, register and geotran tasks wcs updating + code to deal correclty with dimensionally reduced data. (3/16/92, Davis) + +pkg$images/imarith/icalip.gx +pkg$images/imarith/icclip.gx +pkg$images/imarith/ipslip.gx +pkg$images/imarith/icslip.gx +pkg$images/imarith/icmedian.gx + The median calculation with an even number of points for short data + could overflow (addition of two short values) and be incorrect. + (3/16/92, Valdes) + +pkg$images/geometry/t_blkavg.x +pkg$images/geometry/t_blkrep.x + 1. Improved the precision of the blkavg task wcs updating code. + 2. Changed the blkrep task wcs updating code so that it is consistent + with blkavg. This means that a blkrep command followed by a blkavg + command or vice versa will return the original coordinate system + to within machine precision. (3/16/92, Davis) + +pkg$images/iminfo/listpixels.x + Modified listpixels to print out an error if it could not open the + wcs in the image. (3/15/92, Davis) + +pkg$images/geometry/t_magnify.x + Fixed a bug in the magnify task wcs updating code which was not + working correctly for dimensionally reduced images. (3/15/92, Davis) + +pkg$images/geometry/t_imtrans.x + Fixed a bug in the imtranspose task wcs updating code which was not + working correctly for dimensionally reduced images. (3/14/92, Davis) + +pkg$images/imarith/icalip.gx +pkg$images/imarith/icclip.gx +pkg$images/imarith/icslip.gx + There was a bug allowing the number of valid pixels counter to become + negative. Also there was a step which should not be done if the + number of valid pixels is less than 1; i.e. all pixels rejected. + A test was put in to skip this step. (3/13/92, Valdes) + +pkg$images/iminfo/t_imslice.x +pkg$images/doc/imslice.hlp + Added wcs support to the imslice task. + (3/12/92, Davis) + +pkg$images/iminfo/t_imstat.x + Fixed a bug in the code for computing the standard deviation, kurtosis, + and skew, wherein precision was being lost because two of the intermediate + variables in the computation were real instead of double precision. + (3/10/92, Davis) + +pkg$images/iminfo/listpixels.x + 1. Modified listpixels task to use the MWCS axis "format" attributes + if they are present in the image header. + 2. Added support for dimensionally reduced images, i.e. + images which are sections of larger images and whose coordinate + transformations depend on the reduced axes, to the listpixels task. + (3/6/92, Davis) + +pkg$images/imarith/t_imcombine.x +pkg$images/imarith/icsetout.x + Changed error messages to say IMCOMBINE instead of ICOMBINE. + (3/2/92, Valdes) + +pkg$images/imarith/iclog.x + Added listing of read noise and gain. (2/10/92, Valdes) + +pkg$images/imarith/icscale.x +pkg$images/imarith/icpclip.gx + 1. Datatype declaration for asumi was incorrect. + 2. Reduced the minimum number of images allowed for PCLIP to 3. + (1/7/92, Valdes) + +pkg$images/imarith/icgrow.gx + The first pixel to be checked was incorrectly set to 0 instead of 1 + resulting in a segvio when using the grow option. (12/6/91, Valdes) + +pkg$images/imarith/icgdata.gx +pkg$images/imarith/icscale.x + Fixed datatype declaration errors found by SPPLINT. (11/22/91, Valdes) + +pkg$images/iminfo/t_imstat.x + Fixed a bug in the kurtosis computation found by ST. + (Davis 10/11/91) + +pkg$images/iminfo/t_imstat.x +pkg$images/doc/imstat.hlp + Corrected a bug in the mode computation in imstatistics. The parabolic + interpolation correction for computing the histogram peak was being + applied in the wrong direction. Note that for dev$pix the wrong answer + is actually closer to the expected answer than the correct answer + due to binning effects. + (Davis 9/24/91) + +pkg$images/filters/t_gauss.x + The code which computes the gaussian kernel was producing a divide by + zero error if ratio=0.0 and bilinear=yes (2.10 version only). + (Davis 9/18/91) + +pkg$images/doc/magnify.hlp + Corrected a bug in the magnify help page. + (Davis 9/18/91) + +pkg$images/imarith/icsclip.gx +pkg$images/imarith/icaclip.gx +pkg$images/imarith/iccclip.gx + There was a typo, Memr[d[k]+k] --> Memr[d[j]+k]. (9/17/91, Valdes) + +pkg$images/imarith/icstat.gx +pkg$images/imarith/icmask.x + The offsets were used improperly in computing image statistics. + (Valdes, 9/17/91) + +pkg$images/geometry/t_imshift.x + The shifts file pointer was not being correctly initialized to NULL + in the case where no shifts file was declared. When the task + was invoked repeatedly from a script, this could result in an array being + referenced, for which space had not been previously allocated. + (Davis 7/29/91) + +pkg$images/imarith/imc* - +pkg$images/imarith/ic* + +pkg$images/imarith/t_imcombine.x +pkg$images/imarith/mkpkg +pkg$images/imarith/generic/mkpkg +pkg$images/imcombine.par +pkg$images/doc/imcombine.hlp + Replaced old version of IMCOMBINE with new version supporting masks, + offsets, and new algorithms. (Valdes 7/19/91) + +pkg$images/iminfo/imhistogram.x + Imhistogram has been modified to print the value of the middle of + histogram bin instead of the left edge if the output type is list + instead of plot. (Davis 6/11/91) + +pkg$images/t_imsurfit.x + Modified the sections file reading code to check the order of the + x1 x2 y1 y2 parameters and switch (x1,x2) or (y1,y2) if x2 < x1 or + y2 < y1 respectively. (Davis 5/28/91) + +pkg$images/listpixels.par +pkg$images/iminfo/listpixels.x +pkg$images/doc/listpixels.hlp + Modified the listpixels task to be able to print the pixel coordinates + in logical, physical or world coordinates. The default coordinate + system is still logical as before. (Davis 5/17/91) + +pkg$images/images.par +pkg$images/doc/minmax.hlp +pkg$images/imutil/t_minmax.x +pkg$images/imutil/minmax.x + Minmax was modified to do the minimum and maximum values computations + in double precision or complex instead of real if the input image + pixel type is double precision or complex. Note that the minimum and + maximum header values are still stored as real however. + (Davis 5/16/91) + +imarith/t_imarith.x + There was a missing statement to set the error flag if the image + dimensions did not match. (5/14/91, Valdes) + +doc/imarith.hlp + Fixed some formatting problems in the imarith help page. (5/2/91 Davis) + +imarith$imcombine.x + Changed the order in which images are unmapped to have the output images + closed last. This is to allow file descriptors for the temporary image + used when updating STF headers. (4/22/91, Valdes) + +pkg$images/geometry/t_blkavg.x +pkg$images/geometry/blkavg.gx +pkg$images/geometry/blkavg.x + The blkavg task was partially modified to support complex image data. + The full modifications cannot be made because of an error in abavx.x + and the missing routine absux.x. + (4/18/91 Davis) + +pkg$images/geometry/geofit.x + The x and y fits cross-terms switch was not being set correctly to "yes" + in the case where xxorder=2 and xyorder=2 or in the case where yxorder=2 + and yyorder=2. + (4/9/91 Davis) + +pkg$images/geometry/geogmap.x + Modified the line which prints the geometric parameters to use the + variable name xshift and yshift instead of delx and dely. + (4/9/91 Davis) + +pkg$images/imfit/imsurfit.x + Fixed a bug in the pixel rejection code which occurred when upper was > + 0.0 and lower = 0.0 or lower > 0 and upper = 0.0. The problem was that + the code was simply setting the rejection limits to the computed sigma + times the upper and lower parameters without checking for the 0.0 + condition first. In the first case this results in all points with + negative residuals being rejected and in the latter all points with + positive residuals are rejected. + (2/25/91 Davis) + +pkg$images/doc/hedit.hlp +pkg$images/doc/hselect.hlp +pkg$images/doc/imheader.hlp +pkg$images/doc/imgets.hlp + Added a reference to imgets in the SEE ALSO sections of the hedit and + hselect tasks. + Added a reference to hselect and hedit in the SEE ALSO sections of the + imheader and imgets tasks. + (2/22/91 Davis) + +pkg$images/gradient.hlp +pkg$images/laplace.hlp +pkg$images/gauss.hlp +pkg$images/convolve.hlp +pkg$images/gradient.par +pkg$images/laplace.par +pkg$images/gauss.par +pkg$images/convolve.par +pkg$images/t_gradient.x +pkg$images/t_laplace.x +pkg$images/t_gauss.x +pkg$images/t_convolve.x +pkg$images/convolve.x +pkg$images/xyconvolve.x +pkg$images/radcnv.x + The convolution operators were modified to run more efficiently in + certain cases. The LAPLACE task was modified to make use of the + radial symmetry of the convolution kernel in the y direction as well + as the x direction resulting in a modest speedup in execution time. + A new parameter bilinear was added to the GAUSS and CONVOLVE tasks. + By default and if appropriate mathematically, GAUSS now makes use of + the bilinearity or separability of the Gaussian function, + to separate the 2D convolution in x and y into two equivalent + 1D convolutions in x and y, resulting in a considerable speedup + in execution time. Similarly the user can know program CONVOLVE to + compute a bilinear convolution instead of a full 2D 1 if appropriate. + (1/29/91 Davis) + +pkg$images/filters/t_convolve.x + CONVOLVE was not decoding the legal 1D kernel "1.0 2.0 1.0" correctly + although the alternate form "1.0 2.0 1.0;" worked. Leading + blanks in string kernels as in for example " 1.0 2.0 1.0" also generated + and error. Fixed these bugs and added some additional error checking code. + (11/28/90 Davis) + +pkg$images/doc/gauss.hlp + Added a detailed mathematical description of the gaussian kernel used + by the GAUSS task to the help page. + +pkg$images/images.hd +pkg$images/rotate.cl +pkg$images/imlintran.cl +pkg$images/register.cl +pkg$images/register.par + Added src="script file name" entries to the IMAGES help database + for the tasks ROTATE, IMLINTRAN, and REGISTER. Changed the CL + script for REGISTER to a procedure script to remove the ugly + local variable declarations. Added a few comments to the scripts. + (12/11/90, Davis) + +pkg$images/iminfo/imhistogram.x + Added a new parameter binwidth to imhistogram. If binwidth is defined + it determines the histogram resolution in intensity units, otherwise + nbins determines the resolution as before. (10/26/90, Davis) + +pkg$images/doc/sections.hlp + Clarified what is meant by an image template and that the task itself + does not check whether the specified names are actually images. + The examples were improved. (10/3/90, Valdes) + +pkg$images/doc/fit1d.hlp + Changed lines to columns in example 2. (10/3/90, Valdes) + +pkg$images/imarith/imcscales.x + When an error occured while parsing the mode section the untrapped error + caused further problems downstream. Because it would require adding + lots of errchks to cause the program to gracefully abort I instead made + it a warning. (10/2/90, Valdes) + +pkg$images/imutil/hedit.x + Hedit was computing but not using min_lenarea. If the user specified + a min_lenuserarea greater than the default of 28800 then the default + was being used instead of the larger number. + +pkg$imarith/imasub.gx + The case of subtracting an image from the constant zero had a bug + which is now fixed. (8/14/90, Valdes) + +pkg$images/t_imtrans.x + Modified the imtranspose task so it will work on type ushort images. + (6/6/90 Davis) + +pkg$images + Added world coordinate system support to the following tasks: imshift, + shiftlines, magnify, imtranspose, blkrep, blkavg, rotate, imlintran, + register and geotran. The only limitation is that register and geotran + will only support simple linear transformations. + (2/24/90 Davis) + +pkg$images/geometry/geotimtran.x + Fixed a problem in the boundary extension "reflect" option code for small + images which was causing odd values to be inserted at the edges of the + image. + (2/14/90 Davis) + +pkg$images/iminfo/imhistogram.x + A new parameter "hist_type" was added to the imhistogram task giving + the user the option of plotting the integral, first derivative and + second derivative of the histogram as well as the normal histogram. + Code was contributed by Rob Seaman. + (2/2/90 Davis) + +pkg$images/geometry/geogmap.x + The path name of the help file was being erroneously renamed with + the result that when users ran the double precision version of the + code they could not find the help file. + (26/1/90 Davis) + +pkg$images/filters/t_boxcar.x,t_convolve.x + Added some checks for 1-D images. + (1/20/90 Davis) + +pkg$images/iminfo/t_imstat.x,imstat.h + Made several minor bug fixes and alterations in the imstatistics task + in response to user complaints and suggestions. + + 1. Changed the verbose parameter to the format parameter. If format is + "yes" (the default) then the selected fields are printed in fixed format + with column labels. Other wise the fields are printed in free format + separated by 2 blanks. This fixes the problem of fields running together. + + 2. Fixed a bug in the code which estimates the median from the image + histogram by linearly interpolating around the midpt of the integrated + histogram. The bug occurred when more than half the pixels were in the + first bin. + + 3. Added a check to ensure that the number of fields did not overflow + the fields array. + + 4. Removed the extraneous blank line printed after the title. + + 5. The pound sign is now printed at the beginning of the column header + string regardless of which field is printed first. In the previous + versions it was only being printed if the image name field was + printed first. + + 6. Changed the name of the median field to midpt in response to user + confusions about how the median is computed. + + (1/20/90, Davis) + +pkg$images/imutil/t_imslice.hlp + The imslice was not correctly computing the number of lines in the + output image in the case where the slice dimension was 1. + (12/4/89, Davis) + +pkg$images/doc/imcombine.hlp + Clarified and documented definitions of the scale, offset, and weights. + (11/30/89, Valdes) + +pkg$images/geometry/geotran.x + High order surfaces of a certain functional form could occasionally + produce out of bounds pixel errors. The bug was caused by not properly + computing the distortion of the image boundary for higher order + surfaces. + (11/21/89, Davis) + +pkg$images/geometry/imshift.x + The circulating buffer space was not being freed after each execution + of IMSHIFT. This did not cause an error in execution but for a long + list of frames could result in alot of memory being tied up. + (10/25/89, Davis) + +pkg$images/imarith/t_imarith.x + IMARITH is not prepared to deal with images sections in the output. + It used to look for '[' to decide if the output specification included + and image section. This has been changed to call the IMIO procedure + imgsection and check if a non-null section string is returned. + Thus it is up to IMIO to decide what part of the image name is + an image section. (9/5/89, Valdes) + +pkg$images/imarith/imcmode.gx + Fixed bug causing infinite loop when computing mode of constant value + section. (8/14/89, Valdes) + +==== +V2.8 +==== + +pkg$images/iminfo/t_imstat.x + Davis, Jun 15, 1989 + Added a couple of switches to that skew and kurtosis are not computed + if they are not to be printed. + +pkg$images/iminfo/t_imstat.x + Davis, Jun 14, 1989 + A simple mod was made to the skew and kurtosis computation to avoid + divide by zero errors in case of underflow. + +pkg$images/imutil/chpixtype.par + Davis, Jun 13, 1989 + The parameter file has been modified to accept an output pixel + type of ushort. + +pkg$images/imarith/imcombine.gx + Valdes, Jun 2, 1989 + A new scheme to detect file errors is now used. + +pkg$images/imfit/t_imsurfit.x + Davis, Jun 1, 1989 + 1. If the user set regions = "sections" but the sections file + did not exist the task would go into an infinite loop. The problem + was a missing error check on the open statement. + +pkg$images/iminfo/imhistogram.x,imhistogram.par + Davis, May 31, 1989 + A new version of imhistogram has been installed. These mods have + been made over a period of a month by Doug Tody and Rob Seaman. + The mods include + 1. An option to turn off log scaling of the y axis of the histogram plot. + 2. A new autoscale parameter which avoids aliasing problems for integer + data. + 3. A new parameter top_close which resolves the ambiguity in the top + bin of the histogram. + +pkg$images/imarith/imcombine.gx + Valdes, May 9, 1989 + Because a file descriptor was not reserved for string buffer operations + and a call to stropen in cnvdate was not error checked the task would + hang when more than 115 images were combined. Better error checking + was added and now an error message is printed when the maximum number + of images that can be combined is exceeded. + +pkg$images/imarith/t_imarith.x + Valdes, May 6, 1989 + Operations in which the output image has an image section are now + skipped with a warning message. + +pkg$images/imarith/sigma.gx +pkg$images/imarith/imcmode.gx + Valdes, May 6, 1989 + 1. The weighted sigma was being computed incorrectly. + 2. The argument declarations were wrong for integer input images. + Namely the mean vector is always real. + 3. Minor change to imcmode.gx to return correct datatype. + +pkg$images/imstack,imslice + Davis, April 1, 1989 + The proto images tasks imstack and imslice have been moved from the + proto package to the images package. Imstack is unchanged except that + it now supports the image data types USHORT and COMPLEX. Imslice has + been modified to allow slicing along any dimension of the image instead + of just the highest dimension. + +pkg$images/imstatistics. + Davis, Mar 31, 1989 + 1. A totally new version of the imstatistics task has been written + and replaces the old version. The new task allows the user to select + which statistical parameters to compute and print. These include + the mean, median, mode, standard deviation, skew, kurtosis and the + minimum and maximum pixel values. + +pkg$images/imhistogram.par +pkg$images/iminfo/imhistogram.x +pkg$images/doc/imhistogram.hlp + Davis, Mar 31, 1989 + 1. The imhistogram task has been modified to plot "box" style histograms + as well as "line" type histograms. Type "line" remains the default. + +pkg$images/geometry/geotran.par,register.par,geomap.par +pkg$images/doc/geomap.hlp,register.hlp,geotran.hlp + Davis, Mar 6, 1989 + 1. Improved the parameter prompting in GEOMAP, REGISTER and GEOTRAN + and improved the help pages. + 2. Changed GEOMAP database quantities "xscale" and "yscale" to "xmag" + and "ymag" for consistency . Geotran was changed appropriately. + +pkg$images/imarith/imcmode.gx + For short data a short variable was wraping around when there were + a significant number of saturated pixels leading to an infinite loop. + The variables were made real regardless of the image datatype. + (3/1/89, Valdes) + +pkg$images/imutil/imcopy.x + Davis, Feb 28, 1989 + 1. Added support for type USHORT to the imcopy task. This is a merged + ST modification. + +pkg$images/imarith/imcthreshold.gx +pkg$images/imcombine.par +pkg$images/doc/imcombine.hlp +pkg$images/imarith/imcscales.x + Valdes, Feb 16, 1989 + 1. Added provision for blank value when all pixels are rejected by the + threshold. + 2. Fixed a bug that was improperly scaling images in the threshold option. + 3. The offset printed in the log now has the opposite sign so that it + is the value "added" to bring images to a common level. + +pkg$images/imfit/imsurfit.x + Davis, Feb 23, 1989 + Fixed a bug in the median fitting code which could cause the porgram + to go into an infinite loop if the region to be fitted was less than + the size of the whole image. + +pkg$images/geometry/t_magnify.x + Davis, Feb 16, 1989 + Modified magnify to work on 1D images as well as 2D images. The + documentation has been updated. + +pkg$images/geometry/t_geotran.x + Davis, Feb 15, 1989 + Modified the GEOTRAN and REGISTER tasks so that they can handle a list + of transform records one for each input image. + +pkg$images/imarith/imcmode.gx + Valdes, Feb 8, 1989 + Added test for nx=1. + +pkg$images/imarith/t_imcombine.x + Valdes, Feb 3, 1989 + The test for the datatype of the output sigma image was wrong. + +pkg$images/iminfo/listpixels.x,listpixels.par + Davis, Feb 6, 1989 + The listpixels task has been modified to print out the pixels for a + list of images instead of a single image only. A title line for each + image listed can optionally be printed on the standard output if + the new parameter verbose is set to yes. + +pkg$images/geometry/t_imshift.x + Davis, Feb 2, 1989 + Added a new parameter shifts_file to the imshift task. Shifts_file + is the name of a text file containing the the x and yshifts for + each input image to be shifted. The number of input shifts must + equal the number of input images. + +pkg$images/geometry/t_geomap.x + Davis, Jan 17, 1989 + Added an error message for the case where the coordinates is empty + of there are no points in the specified data range. Previously the + task would proceed to the next coordinate file without any message. + +pkg$images/geometry/t_magnify.x + Davis, Jan 14, 1989 + Added the parameter flux conserve to the magnify task to bring it into + line with all the other geometric transformation tasks. + +pgk$images/geometry/geotran.x,geotimtran.x + Davis, Jan 2, 1989 + A bug was fixed in the flux conserve code. If the x and y reference + coordinates are not in pixel units and are not 1 then + the computed flux per pixel was too small by xscale * yscale. + +pkg$images/filters/acnvrr.x,convolve.x,boxcar.x,aboxcar.x + Davis, Dec 27, 1988 + I changed the name of the acnvrr procedure to cnv_radcnvr to avoid + a name conflict with a vops library procedure. This only showed + up when shared libraries were implemented. I also changed the name + of the aboxcarr procedure to cnv_aboxr to avoid conflict with the + vops naming conventions. + +pkg$images/imarith/imcaverage.gx + Davis, Dec 22, 1988 + Added an errchk statement for imc_scales and imgnl$t to stop the + program bombing with segmentation violations when mode <= 0. + +pkg$images/imarith/imcscales.x + Valdes, Dec 8, 1988 + 1. IMCOMBINE now prints the scale as a multiplicative quantity. + 2. The combined exposure time was not being scaled by the scaling + factors resulting in a final exposure time inconsistent with the + data. + +pkg$images/iminfo/imhistogram.x + Davis, Nov 30, 1988 + Changed the list+ mode so that bin value and count are printed out instead + of bin count and value. This makes the plot and list modes compatable. + +pkg$images/iminfo/t_imstat.x + Davis, Nov 17, 1988 + Added the n=n+1 back into the inner loop of imstat. + +pkg$images/geotran.par,register.par + Davis, Nov 11 , 1988 + Fixed to glaring errors in the parameter files for register and geotran. + Xscale and yscale were described as pixels per reference unit when + they should be reference units per pixel. The appropriate bug fix has been + made. + +pkg$images/geometry/t_geotran.x + Davis, November 7, 1988 + The routine gsrestore was not being error checked. If either of the + input x or y coordinate surface was linear and the other was not, + the message came back GSRESTORE: Illegal x coordinate. This bug has + been fixed. + +pkg$images/imarith/imcombine.gx + Valdes, October 19, 1988 + A vops clear routine was not called generically causing a crash with + double images. + +pkg$images/filters/t_fmedian.x,t_median.x,t_fmode.x,t_mode.x,t_gradient.x + t_gauss.x,t_boxcar.x,t_convolve.x,t_laplace.x + Davis, October 4, 1988 + I fixed a bug in the error handling code for the filters tasks. If + and error occurred during task execution and the input image name was + the same as the output image name then the input image was trashed. + +pkg$images/imarith/imcscales.gx + Valdes, September 28, 1988 + It is now an error for the mode to be nonpositive when scaling or weighting. + +pkg$images/imarith/imcmedian.gx + Valdes, August 16, 1988 + The median option was selecting the n/2 value instead of (n+1)/2. Thus, + for an odd number of images the wrong value was being determined for the + median. + +pkg$images/geometry/t_imshift.x + Davis, August 11, 1988 + 1. Imshift has been modified to uses the optimized code if nearest + neighbour interpolation is requested. A nint is done on the shifts + before calling the quick shift routine. + 2. If the requested pixel shift is too large imshift will now + clean up any pixelless header files before continuing execution. + +pkg$images/geometry/blkavg.gx + Davis, July 13, 1988 + Blkavg has been fixed so that it will work on 1D images. + +pkg$images/geometry/t_imtrans.x,imtrans.x + Davis, July 12, 1988 + Imtranspose has been modified to work on complex images. + +pkg$images/imutil/t_chpix.x + Davis, June 29, 1988 + A new task chpixtype has been added to the images package. Chpixtype + changes the pixel types of a list of images to a specified output pixel + type. Seven data types are supported "short", "ushort", "int", "long" + "real" and "double". + +pkg$images/geometry/rotate.cl,imlintran.cl,t_geotran.x + Davis, June 10, 1988 + The rotate and imlintran scripts have been rewritten to use procedure + scripts. This removes all the annoying temporary cl variables which + appear when the user does an lpar. In previous versions of these + two tasks the output was restricted to being the same size as the input + image. This is still the default case, but the user can now set the + ncols and nrows parameters to the desired output size. I ncols or nlines + < 0 then then the task compute the output image size required to contain + the whole input image. + +pkg$images/filters/t_convolve.x,t_laplace.x,t_gradient.x,t_gauss.x,convolve.x + Davis, June 1, 1988 + The convolution operators laplace, gauss and convolve have been modified + to make use of radial symmetry in the convolution kernel. In gauss and + laplace the change is transparent to the user. For the convolve operator + the user must indicate that the kernel is radially symmetric by setting + the parameter radsym. For kernels of 7 by 7 or greater the speedup + in timings is on the order of 30% on the Vax 750 with the fpa. + +pkg$images/imarith/imcmode.gx + Valdes, Apr 11, 1988 + 1. The use of a mode sections was handled incorrectly. + +pkg$images/imfit/fit1d.x + Valdes, Jan 4, 1988 + 1. Added an error check for a failure in IMMAP. The missing error check + caused FIT1D to hang when a bad input image was specified. + +pkg$images/magnify.par +pkg$images/imcombine.par +pkg$images/imarith/imcmode.gx +pkg$images/doc/imarith.hlp + Valdes, Dec 7, 1987 + 1. Added option list to parameter prompts. + 2. Fixed minor typo in help page + 3. The mode calculation in IMCOMBINE would go into an infinite loop + if all the pixel values were the same. If all the pixels are the + same them it skips searching for the mode and returns the constant + number. + +pkg$images/geometry/geotimtran.x + Davis, Nov 25, 1987 + 1. A bug in the boundary extension = wrap option was found in the + IMLINTRAN task. The problem occured in computing values for out of + bounds pixels in the range 0.0 < x < 1.0, ncols < x < ncols + 1.0, + 0.0 < y < 1.0 and nlines < y < nlines + 1. The computed coordinates + were falling outside the boundaries of the interpolation array. + +pkg$images/geometry/t_geomap.x,geograph.x + Davis, Nov 19, 1987 + 1. The geomap task now writes the name of the output file into the database. + 2. Rotation angles of 360. degrees have been altered to 0 degrees. + +pkg$images/imfit/t_imsurfit.x,imsurfit.x +pkg$images/lib/ranges.x + Davis, Nov 2, 1987 + A bug in the regions fitting option of the IMSURFIT task has been found + and fixed. This bug would occur when the user set the regions parameter + to sections and then listed section which overlapped each other. The + modified ranges package was not handling the overlap correctly and + computing a number of points which was incorrect. + +pkg$images/imarith/* + + Valdes, Sep 30, 1987 + The directory was reorganized to put generic code in the subdirectory + generic. + + A new task called IMCOMBINE has been added. It provides for combining + images by a number of algorithms, statistically weighting the images + when averaging, scaling or offsetting the images by the exposure time + or image mode before combining, and rejecting deviant pixels. It is + almost fully generic including complex images and works on images of + any dimension. + +pkg$images/geometry/geotran.x + Davis, Sept 3, 1987 + A bug in the flux conserving algorithm was found in the geotran code. + The symptom was that the flux of the output image occasionally was + negative. This would happen when two conditions were met, the transformation + was of higher order than a simple rotation, magnification, translation + and an axis flip was involved. The mathematical interpretation of this + bug is that the coordinate surface had turned upside down. The solution + for people running systems with this bug is to multiply there images + by -1. + +pkg$images/imfit/imsurfit.h,t_imsurfit.x + Davis, Aug 6, 1987 + A new option was added to the parameter regions in the imsurfit task. + Imsurfit will now fit a surface to a single circular region defined + by an x and y center and a radius. + +pkg$images/geometry/geotimtran.x + Davis, Jun 15, 1987 + Geotran and register were failing when the output image number of rows + and columns was different from the input number of rows and columns. + Geotran was mistakenly using the input images sizes to determine the + number of output lines that should be produced. The same problem occurred + when the values of the boundary pixels were being computed. The program + was using the output image dimensions to compute the boundary pixels + instead of the input image dimensions. + +pkg$images/geometry/geofit.x,geogmap.x + Davis, Jun 11, 1987 + A bug in the error checking code in the geomap task was fixed. The + condition of too few points for a reasonable was not being trapped + correctly. The appropriate errchk statements were added. + +pkg$images/geomap.par + Davis, Jun 10, 1987 + The default fitting function was changed to polynomial. This will satisfy + most users who wish to do shifts, rotations, and magnifications and + avoid the neccessity of correctly setting the xmin, xmax, ymin, and ymax + parameters. For the chebyshev and legendre polynomial functions these + parameters must be explicitly set. For reference coordinates in pixel + units the normal settings are 1, ncols, 1 and nlines respectively. + +pkg$images/iminfo/hselect.x,imheader.x,images$/imutil/hselect.x + Davis, Jun 8, 1987 + Imheader has been modified to open an image with the default min_lenuserarea + Hselect and hedit will now open the image setting the user area to the + maximum of 28800 chars or the min_lenuser environment variable. + +pkg$images/iminfo/t_imstat.x + Davis, May 22, 1987 + An error in the image minimum computation was corrected. This error + would show up most noiticeably if imstat was run on a 1 pixel image. + The min value would be left set to MAX_REAL. + +pkg$images/filters/mkpkg + Davis, May 22, 1987 + I added mach.h to the dependency file list of t_fmedian.x and + recompiled. The segmentation violations I had been getting in the + program disappeared. + +pkg$images/t_shiftlines.x,shiftlines.x + Davis, April 15, 1987 + 1. I changed the names of the procedures shiftlines and shiftlinesi + to sh_lines and sh_linesi. When the original names were contracted + to 6 letter fortran names they became shifti and shifts which just + so happens to collide with shifti and shifts in the subdirectory + osb. On VMS this was causing problems with the shareable libraries. + If images was linked with -z there was no problem. + +pkg$images/imarith/t_imsum.x + Valdes, March 24, 1987 + 1. IMSUM was failing to unmap images opened to check image dimensions + in a quick first pass through the image list. This is probably + the source of the out of files problem with STF images. It may + be the source of the out of memory problem reported from AOS/IRAF. + +pkg$images/imfit/fit1d.x +pkg$images/imfit/mkpkg + Valdes, March 17, 1987 + 1. Added error checking for the illegal operation in which both input + and output image had an image section. This was causing the task + to crash. The task now behaves properly in this circumstance and + even allows the fitted output to be placed in an image section of + an existing output image (even different than the input image + section) provided the input and output images have the same sizes. + +pkg$images/t_convolve.x + Davis, March 3, 1987 + 1. Fixed the kernel decoding routine in the convolve task so that + it now recognizes the row delimter character in string entry mode. + +pkg$images/geometry,filters + Davis, February 27, 1987 + 1. Changed all the imseti (im, TY_BNDRYPIXVAL, value) calls to imsetr. + +pkg$images/t_minmax.x,minmax.x + Davis, February 24, 1987 + 1. Minmax has been changed to compute the minimum and maximum pixel + as well as the minimum and maximum pixel values. The pixels are output + in section notation and stored in the minmax parameter file. + +pkg$images/t_magnify.x + Davis, February 19, 1987 + 1. Magnify was aborting with the error MSIFIT: Too few datapoints + when trying to reduce an image using the higher order interpolants + poly3, poly5 and spline3. I increased the NEDGE defined constant + from 2 to three and modified the code to use the out of bounds + imio. + +pkg$images/geograph.x,geogmap.x + Davis, February 17, 1987 + 1. Geomap now uses the gpagefile routine to page the .keys file. + The :show command deactivates the workstation before printing a + block of text and reactivates it when it is finished. + +pkg$images/geometry/geomap,geotran + Davis, January 26, 1987 + 1. There have been substantial changes to the geomap, and geotrans + tasks and those tasks rotate, imlintran and register which depend + on them. + 2. Geomap has been changed to be able to compute a transformation + in both single and double precision. + 3. The geotran code has been speeded up considerably. A simple rotate + now takes 70 seconds instead of 155 seconds using bilinear interpolation. + 4. Two new cl parameters nxblock and nyblock have been added to the + rotate, imlintran, register and geotran tasks. If the output image + is smaller than these parameters then the entire output image + is computed at once. Otherwise the output image is computed in blocks + nxblock by nyblock in size. + 5. The 3 geotran parameters rotation, scangle and flip have been replaced + with two parameters xrotation and yrotation which serve the same purpose. + +pkg$images/geometry/t_shiftlines.x,shiftlines.x + Davis, January 19, 1987 + 1. The shiftlines task has been completely rewritten. The following + are the major changes. + 2. Shiftlines now makes use of the imio boundary extension operations. + Therefore the four options: nearest pixel, reflect, wrap and constant + boundary extension are available. + 3. The interpolation code has been vectorised. The previous version + was using the function call asieval for every output pixel evaluated. + The asieval call were replaced with asivector calls. + 4. An extra CL parameter constant to support constant boundary + exension was added. + 5. The shiftlines help page was modified and the date changed to + January 1987. + +pkg$images/imfit/imsurfit.x + Davis, January 12, 1987 + 1. I changed the amedr call to asokr calls. For my application it did + not matter whether the input array is left partially sorted and the asokr + routine is more efficient. + +pkg$images/lib/pixlist.x + Davis, December 12, 1986 + 1. A bug in the pl_get_ranges routine caused the routine to fail when the + number of ranges got too large. The program could not detect the end of + the ranges and would go into an infinite loop. + +pkg$images/iminfo/t_imstat.x + Davis, December 3, 1986 + 1. Imstat was failing on constant images because finite machine precision + could result in a negative sigma squared. Added a check for this condition. + +pkg$images/filters/fmode.x + Davis, October 27, 1986 + 1. Added a check for 0 data range before calling amapr. + +pkg$images/imarith/imsum.gx + Valdes, October 20, 1986 + 1. Found and fixed bug in this routine which caused pixel rejection + to fail some fraction of the time. + +pkg$images/geometry/blkrp.gx + Valdes, October 13, 1986 + 1. There was a bug when the replication factor for axis 1 was 1. + +pkg$images/iminfo/imhistogram.x + Hammond, October 8, 1986 + 1. Running imhistogram on a constant valued image would result in + a "floating divide by zero fault" in ahgm. This condition is + now trapped and a warning printed if there is no range in the data. + +pkg$images/tv/doc/cvl.hlp + Valdes, October 7, 1986 + 1. Typo in V2.3 documentation fixed: "zcale" -> "zscale". + +pkg$images/fit1d.par + Valdes, October 7, 1986 + 1. When querying for the output type the query was: + +Type of output (fit, difference, ratio) (fit|difference|ratio) (): + + The enumerated values were removed since they are given in the + prompt string. + +pkg$images/imarith/t_imsum.x +pkg$images/imarith/imsum.gx +pkg$images/do/imsum.hlp + Valdes, October 7, 1986 + 1. Medians or pixel rejection with more than 15 images is now + correct. There was an error in buffering. + 2. Averages of integer datatype images are now correct. The error + was caused by summing the pixel values divided by the number + of images instead of summing the pixel values and then dividing + by the number of images. + 3. Option keywords may now be abbreviated. + 4. The output pixel datatype now defaults to the calculation datatype + as is done in IMARITH. The help page was modified to indicate this. + 5. Dynamic memory is now used throughout to reduce the size of the + executable. + 6. The bugs 1-2 are present in V2.3 and not in V2.2. + +pkg$images/imarith/t_imarith.x +pkg$images/imarith.par +pkg$images/doc/imarith.hlp + Valdes, October 6, 1986 + 1. The parameter "debug" was changed to "noact". "debug" is reserved + for debugging information. + 2. The output pixel type now defaults to the calculation datatype. + 3. The datatype of constant operands is determined with LEXNUM. This + fixes a bug in which a constant such as "1." was classified as an + integer. + 4. Trailing whitespace in the string for a constant operand is allowed. + This fixes a bug with using "@" files created with the task FIELDS + from a table of numbers. Trailing whitespace in image names is + not checked for since this should be taken care of by lower level + system services. + 5. The reported bug with the "max" operation not creating a pixel file + was the result of the previous round of changes. This has been + corrected. This problem does not exist in the released version. + 6. All strings are now dynamically allocated. Also IMTOPENP is used + to open a CL list directly. + 7. The help page was revised for points (1) and (2). + +pkg$images/fmode.par +pkg$images/fmd_buf.x +pkg$images/med_sort.x + Davis, September 29, 1986 + 1. Changed the default value of the unmap parameter in fmode to yes. The + documentation was changed and the date modified. + 2. Added a test to make sure that the input image was not a constant + image in fmode and fmedian. + 3. Fixed the recently added swap macro in the sort routines which + was giving erroneous results for small boxes in tasks median and mode. + +pkg$images/imfit/fit1d.x + Valdes, September 24, 1986 + 1. Changed subroutine name with a VOPS prefix to one with a FIT1D + prefix. + +pkg$images/imarith/t_imdivide.x +pkg$images/doc/imdivide.hlp +pkg$images/imdivide.par + Valdes, September 24, 1986 + 1. Modified this ancient and obsolete task to remove redundant + subroutines now available in the VOPS library. + 2. The option to select action on zero divide was removed since + there was only one option. Parameter file changed. + 3. Help page revised. + +pkg$images/geometry/t_blkrep.x + +pkg$images/geometry/blkrp.gx + +pkg$images/geometry/blkrep.x + +pkg$images/doc/blkrep.hlp + +pkg$images/doc/mkpkg +pkg$images/images.cl +pkg$images/images.men +pkg$images/images.hd +pkg$images/x_images.x + Valdes, September 24, 1986 + 1. A new task called BLKREP for block replicating images has been added. + This task is a complement to BLKAVG and performs a function not + available in any other way. + 2. Help for BLKREP has been added. + +pkg$images/imarith/t_imarith.x +pkg$images/imarith/imadiv.gx +pkg$images/doc/imarith.hlp +pkg$images/imarith.par + Valdes, September 24, 1986 + 1. IMARITH has been modified to provide replacement of divisions + by zero with a constant parameter value. + 2. The documentation has been revised to include this change and to + clarify and emphasize areas of possible confusion. + +pkg$images/doc/magnify.hlp +pkg$images/doc/blkavg.hlp + Valdes, September 18, 1986 + 1. The MAGNIFY help document was expanded to clarify that images with axis + lengths of 1 cannot be magnified. Also a discussion of the output + size of a magnified image. This has been misunderstood often. + 2. Minor typo fix for BLKAVG. + +images$geometry/blkav.gx: Davis, September 7, 1986 + 1. The routine blkav$t was declared a function but called everywhere as + a procedure. Removed the function declaration. + +images$filters/med_sort.x: Davis, August 14, 1986 + 1. A bug in the sorting routine for MEDIAN and MODE in which the doop + loop increment was being set to zero has been fixed. This bug was + causing MEDIAN and MODE to fail on class 6 for certain sized windows. + +images$imfit/fit1d.x: Davis, July 24, 1986 + 1. A bug in the type=ratio option of fit1d was fixed. The iferr call + on the vector operator adivr was not trapping a divide by zero + condition. Changed adivr to adivzr. + +images$iminfo/listpixels.x: Davis, July 21, 1986 + 1. I changed a pargl to pargi for writing out the column number of the + pixels. + +images$iminfo/t_imstat.x: Davis, July 21, 1986 + 1. I changed a pargr to a pargd for the double precision quantitiies + sum(MIN) and sum(MAX). + +images$imfit/t_lineclean.x: Davis, July 14, 1986 + 1. Bug in the calling sequence for ic_clean fixed. The ic pointer + was not being passed to ic_clean causing access violation and/or + segmentation violation errors. + +images$imfit/fit1d.x, lineclean.x: Valdes, July 3, 1986 + 1. FIT1D and LINECLEAN modified to use new ICFIT package. + +From Valdes June 19, 1986 + +1. The help page for IMSUM was modified to explicitly state what the +median of an even number of images does. + +----------------------------------------------------------------------------- + +From Davis June 13, 1986 + +1. A bug in CONVOLVE in which insufficient space was being allocated for +long (> 161 elements) 1D kernels has been fixed. CONVOLVE was not +allocating sufficent extra space. + +----------------------------------------------------------------------------- + +From Davis June 12, 1986 + +1. I have changed the default value of parameter unmap in task FMEDIAN to +yes to preserve the original data range. + +2. I have changed the value of parameter row_delimiter from \n to ;. + +----------------------------------------------------------------------------- + +From Davis May 12, 1986 + +1. Changed the angle convention in GAUSS so that theta is the angle of the +major axis with respect to the x axis measured counter-clockwise as specified +in the help page instead of the negative of that angle. + +----------------------------------------------------------------------------- + +From Davis Apr 28, 1986 + +1. Moved geomap.key to lib$scr and made redefined HELPFILE in geogmap.x +appropriately. + +------------------------------------------------------------------------------ + +images$imarith/imsum.gx: Valdes Apr 25, 1986 + 1. Fixed bug in generic code which called the real VOPS operator + regardless of the datatype. This caused IMSUM to fail on short + images. + +From Davis Apr 17, 1986 + +1. Changed constructs of the form boolean == false in the file imdelete.x +to ! boolean. + +------------------------------------------------------------------------------ + +images$imarith: Valdes, April 8, 1986 + 1. IMARITH has been modified to also operate on a list of specified + header parameters. This is primarily used when adding images to + also added the exposure times. A new parameter was added and the + help page modified. + 2. IMSUM has been modified to also operate on a list of specified + header parameters. This is primarily used when summing images to + also sum the exposure times. A new parameter was added and the + help page modified. + +------------------------------------------------------------------------------ + +From Valdes Mar 24, 1986: + +1. When modifying IMARITH to handle mixed dimensions the output image header +was made a copy of the image with the higher dimension. However, the default +when the images were of the same dimension changed to be a copy of the second +operand. This has been changed back to being a copy of the first operand +image. + +------------------------------------------------------------------------------ + +From Davis Mar 21, 1986: + +1. A NULL pointer bug in the subroutine plfree inside IMSURFIT was causing +segmentation violation errors. A null pointer test was added to plfree. + +------------------------------------------------------------------------------ + +From Davis Mar 20, 1986: + +1. A bug involving in place operations in several image tasks has been fixed. + +------------------------------------------------------------------------------ + +From Davis Mar 19, 1986: + +1. IMSURFIT no longer permits the input image to be replaced by the output +image. + +2. The tasks IMSHIFT, IMTRANSPOSE, SHIFTLINES, and GEOTRAN have been modified +to use the images tools xt_mkimtemp and xt_delimtemp for in place +calculations. + +------------------------------------------------------------------------------- + +From Valdes Mar 13, 1986: + +1. Bug dealing with type coercion in short datatype images in IMARITH and IMSUM +which occurs on the SUN has been fixed. +------ +From Valdes Mar 10, 1986: + +1. IMSUM has been modified to work on any number of images. + +2. Modified the help page +------ +From Valdes Feb 25, 1986: + +There have been two changes to IMARITH: + +1. A bug preventing use of image sections has been removed. + +2. An improvement allowing use of images of different dimension. +The algorithm is as follow: + + a. Check if both operands are images. If not the output + image is a copy of the operand image. + + b. Check that the axes lengths are the same for the dimensions + in common. For example a 3D and 2D image must have the same + number of columns and lines. + + c. Set the output image to be a copy of the image with the + higher dimension. + + d. Repeat the operation over the lower dimensions for each of + the higher dimensions. + +For example, consider subtracting a 2D image from a 3D image. The output +image will be 3D and the 2D image is subtracted from each band of the +3D image. This will work for any combination of dimensions. Another +example is dividing a 3D image by a 1D image. Then each line of each +plane and each band will be divided by the 1D image. Likely applications +will be subtracting biases and darks and dividing by response calibrations +in stacked observations. + +3. Modified the help page +=========== +Release 2.2 +=========== +From Davis Mar 6, 1986: + +1. A serious bug had crept into GAUSS after I made some changes. For 2D +images the sense of the sigma was reversed, i.e sigma = 2.0 was actually +sigma = 0.5. This bug has now been fixed. + +--------------------------------------------------------------------------- + +From Davis Jan 13, 1986: + +1. Listpixels will now print out complex pixel values correctly. + +--------------------------------------------------------------------------- + +From Davis Dec 12, 1985: + +1. The directional gradient operator has been added to the images package. + +--------------------------------------------------------------------------- + +From Valdes Dec 11, 1985: + +1. IMARITH has been modified to first check if an operand is an existing +file. This allows purely numeric image names to be used. + +--------------------------------------------------------------------------- + +From Davis Dec 11, 1985: + +1. A Laplacian (second derivatives) operator has been added to the images +package. + +--------------------------------------------------------------------------- + +From Davis Dec 10, 1985: + +1. The new convolution tasks boxcar, gauss and convolve have been added +to the images package. Convolve convolves an image with an arbitrary +user supplied rectangular kernel. Gauss convolves an image with a 2D +Gaussian of arbitrary size. Boxcar will smooth an image using a smoothing +window of arbitrary size. + +2. The images package source code has been reorganized into the following +subdirectories: 1) filters 2) geometry 3) imfit 4) imarith 4) iminfo and +5) imutil 6) lib. Lib contains routines which may be of use to several IRAF +tasks such as ranges. The imutil subdirectory contains tasks which modify +images in some way such as hedit. The iminfo subdirectory contains code +for displaying header and pixel values and other image characteristics +such as the histogram. Image arithmetic and fitting routines are found +in imarith and imfit respectively. Filters contains the convolution and +median filtering routines and geometry contains the geometric distortion +corrections routines. + +3. The documentation of the main images package has been brought into +conformity with the new IRAF standards. + +4. Documentation for imdelete, imheader, imhistogram, listpixels and +sections has been added to the help database. + +5. The parameter structure for imhistogram has been simplified. The +redundant parameters sections and setranges have been removed. + +--------------------------------------------------------------------------- + + +From Valdes Nov 4, 1985: + +1. IMCOPY modified so that the output image may be a directory. Previously +logical directories were not correctly identified. +------ + +From Davis Oct 21, 1985: + +1. A bug in the pixel rejection cycle of IMSURFIT was corrected. The routine +make_ranges in ranges.x was not successfully converting a sorted list of +rejected pixels into a list of ranges in all cases. + +2. Automatic zero divide error checking has been added to IMSURFIT. +------ +From Valdes Oct 17, 1985: + +1. Fit1d now allows averaging of image lines or columns when interactively +setting the fitting parameters. The syntax is "Fit line = 10 30"; i.e. +blank separated line or column numbers. A single number selects just one +line or column. Be aware however, that the actual fitting of the image +is still done on each column or line individually. + +2. The zero line in the interactive curve fitting graphs has been removed. +This zero line interfered with fitting data near zero. +------ +From Rooke Oct 10, 1985: + +1. Blkaverage was changed to "blkavg" and modified to support any allowed +number of dimensions. It was also made faster in most cases, depending on +the blocking factors in each dimension. +------ +From Valdes Oct 4, 1985: + +1. Fit1d and lineclean modified to allow separate low and high rejection +limits and rejection iterations. +------ +From Davis Oct 3, 1985: + +1. Minmax was not calculating the minimum correctly for integer images. +because the initial values were not being set correctly. +------ +From Valdes Oct 1, 1985: + +1. Imheader was modified to print the image history. Though the history +mechanism is little used at the moment it should become an important part +of any image. + +2. Task revisions renamed to revs. +------ +From Davis Sept 30, 1985: + +1. Two new tasks median and fmedian have been added to the images package. +Fmedian is a fast median filtering algorithm for integer data which uses +the histogram of the image to calculate the median at each window. Median +is a slower but more general algorithm which performs the same task. +------ +From Valdes August 26, 1985: + +1. Blkaverage has been modified to include an new parameter called option. +The current options are to average the blocks or sum the blocks. +------ +From Valdes August 7, 1985 + +1. Fit1d and lineclean wer recompiled with the modified icfit package. +The new package contains better labeling and graph documentation. + +2. The two tasks now have parameters for setting the graphics device +and reading cursor input from a file. +______ +From: /u2/davis/ Tue 08:27:09 06-Aug-85 +Package: images +Title: imshift bug + +Imshift was shifting incorrectly when an integral pixel shift in x and +a fractional pixel shift in y was requested. The actual x shift was +xshift + 1. The bug has been fixed and imshift will now work correctly for +any combination of fractional and integral pixel shifts +------ +From: /u2/davis/ Fri 18:14:12 02-Aug-85 +Package: images +Title: new images task + +A new task GEOMAP has been added to the images package. GEOMAP calculates +the spatial transformation required to map one image onto another. +------ +From: /u2/davis/ Thu 16:47:49 01-Aug-85 +Package: images +Title: new images tasks + +The tasks ROTATE, IMLINTRAN and GEODISTRAN have been added to the images +package. ROTATE rotates and shifts an image. IMLINTRAN will rotate, rescale +and shift an an image. GEODISTRAN corrects an image for geometric distortion. +------ +From Valdes July 26, 1985: + +1. The task revisions has been added to page revisions to the images +package. The intent is that each package will have a revisions task. +Note that this means there may be multiple tasks named revisions loaded +at one time. Typing revisions alone will give the revisions for the +current package. To get the system revisions type system.revisions. + +2. A new task called fit1d replaces linefit. It is essentially the same +as linefit except for an extra parameter "axis" which selects the axis along +which the functions are to be fit. Axis 1 is lines and axis 2 is columns. +The advantages of this change are: + + a. Column fitting can now be done without transposing the image. + This allows linefit to be used with image sections along + both axes. + b. For 1D images there is no prompt for the line number. +.endhelp diff --git a/pkg/images/imgeom/blkavg.par b/pkg/images/imgeom/blkavg.par new file mode 100644 index 00000000..eb5f1d82 --- /dev/null +++ b/pkg/images/imgeom/blkavg.par @@ -0,0 +1,12 @@ +# BLKAVG -- Block average or sum on n-dimensional images. + +input,s,a,,,,Input images +output,s,a,,,,Output block averaged images +option,s,h,average,"average|sum",,Type of operation +b1,i,a,1,1,,blocking factor in dimension 1 (x or column) +b2,i,a,1,1,,blocking factor in dimension 2 (y or line) +b3,i,a,1,1,,blocking factor in dimension 3 (z or band) +b4,i,a,1,1,,blocking factor in dimension 4 +b5,i,a,1,1,,blocking factor in dimension 5 +b6,i,a,1,1,,blocking factor in dimension 6 +b7,i,a,1,1,,blocking factor in dimension 7 diff --git a/pkg/images/imgeom/blkrep.par b/pkg/images/imgeom/blkrep.par new file mode 100644 index 00000000..befcc002 --- /dev/null +++ b/pkg/images/imgeom/blkrep.par @@ -0,0 +1,11 @@ +# BLKREP -- Block replicate an n-dimensional images. + +input,s,a,,,,Input images +output,s,a,,,,Output block replicated images +b1,i,a,1,1,,block replication factor in dimension 1 (x or column) +b2,i,a,1,1,,block replication factor in dimension 2 (y or line) +b3,i,a,1,1,,block replication factor in dimension 3 (z or band) +b4,i,a,1,1,,block replication factor in dimension 4 +b5,i,a,1,1,,block replication factor in dimension 5 +b6,i,a,1,1,,block replication factor in dimension 6 +b7,i,a,1,1,,block replication factor in dimension 7 diff --git a/pkg/images/imgeom/doc/blkavg.hlp b/pkg/images/imgeom/doc/blkavg.hlp new file mode 100644 index 00000000..c4491788 --- /dev/null +++ b/pkg/images/imgeom/doc/blkavg.hlp @@ -0,0 +1,65 @@ +.help blkavg Oct85 images.imgeom +.ih +NAME +blkavg -- block average or sum n-dimensional images +.ih +USAGE +.nf +blkavg input output b1 b2 b3 b4 b5 b6 b7 +.fi +.ih +PARAMETERS +.ls input +List of images to be block averaged. Image sections are allowed. +.le +.ls output +List of output image names. If the output image name is the same as the input +image name then the block averaged image replaces the input image. +.le +.ls b1 +The number of columns to be block averaged (dimension 1, or x). +.le +.ls b2 +The number of lines to be block averaged (dimension 2, or y). +.le +.ls b3 +The number of bands to be block averaged (dimension 3, or z). +.le +.ls b4 +The number of pixels to be block averaged in dimension 4 (... etc. for b5-b7). +.le +.ls option = "average" +Type of block average. The choices are "average" and "sum". +.le +.ih +DESCRIPTION +The list of input images are block averaged or summed to form the output images. +The output image names are specified by the output list. The number of +output image names must be the same as the number of input images. +An output image name may be the same +as the corresponding input image in which case the block averaged image replaces +the input image. The last column, line, etc. of the output image may be +a partial block. The option parameter selects whether to block average +or block sum. +.ih +TIMINGS +It requires approximately 10 cpu seconds to block average a 512 by 512 +short image by a factor of 8 in each direction (Vax 11/750 with fpa). +.ih +EXAMPLES +1. To block average a 2-d image in blocks of 2 by 3: + + cl> blkavg imagein imageout 2 3 + +2. To block sum two 2-d images in blocks of 5 by 5: + + cl> blkavg image1,image2 out1,out2 5 5 op=sum + +3. To block average a 3-d image by 4 in x and y and 2 in z: + + cl> blkavg imagein imageout 4 4 2 + + or + + cl> blkavg imagein imageout b1=4 b2=4 b3=2 +.endhelp diff --git a/pkg/images/imgeom/doc/blkrep.hlp b/pkg/images/imgeom/doc/blkrep.hlp new file mode 100644 index 00000000..7f72616b --- /dev/null +++ b/pkg/images/imgeom/doc/blkrep.hlp @@ -0,0 +1,103 @@ +.help blkrep Sep86 images.imgeom +.ih +NAME +blkrep -- block replicate n-dimensional images +.ih +USAGE +.nf +blkrep input output b1 [b2 b3 b4 b5 b6 b7] +.fi +.ih +PARAMETERS +.ls input +List of images to be block replicated. Image sections are allowed. +.le +.ls output +List of output image names. If the output image name is the same as the input +image name then the block replicated image replaces the input image. +.le +.ls b1, b2, b3, b4, b5, b6, b7 +Block replication factor for dimensions 1 - 7. Only the block factors for +the dimensions of the input image are required. Dimension 1 is the column +or x axis, dimension 2 is the line or y axis. +.le +.ih +DESCRIPTION +The list of input images are block replicated by the specified factors +to form the output images. The output image names are specified by the +output list. The number of output image names must be the same as the +number of input images. An output image name may be the same as the +corresponding input image in which case the block averaged image +replaces the input image. Only the block factors for the dimensions +of the input images are used. + +This task is a complement to \fBblkavg\fR which block averages or sums +images. Another related task is \fBmagnify\fR which interpolates +images to arbitrary sizes. This task, however, is only applicable to +two dimensional images with at least two pixels in each dimension. +Finally, in conjunction with \fBimstack\fR a lower dimensional image +can be replicated to higher dimensions. +.ih +TIMINGS +VAX 11/750 with FPA running UNIX 4.3BSD and IRAF V2.4: + +.nf + SIZE DATATYPE REPLICATION CPU CLOCK + 100 short 5 0.5s 1s + 100 real 5 0.5s 1s + 100x100 short 5x5 1.7s 5s + 100x100 real 5x5 2.1s 6s + 100x100x1 real 5x5x5 9.7s 33s +.fi +.ih +EXAMPLES +.ls 4 1. +To block replicate a 1D image in blocks of 3: + +cl> blkrep imagein imageout 3 +.le +.ls 4 2. +To block replicate a 2D image in blocks of 2 by 3: + +cl> blkrep imagein imageout 2 3 +.le +.ls 4 3. +To block replicate two 2D images in blocks of 5 by 5: + +cl> blkrep image1,image2 out1,out2 5 5 +.le +.ls 4 4. +To block replicate a 3D image in place by factors of 2: + +cl> blkrep image1 image1 2 2 2 +.le +.ls 4 5. +To smooth an image by block averaging and expanding by a factor of 2: + +.nf +cl> blkavg imagein imageout 2 2 +cl> blkrep imageout imageout 2 2 +.fi +.le +.ls 4 6. +To take a 1D image and create a 2D image in which each line is the same: + +.nf +cl> imstack image1d image2d +cl> blkrep image2d image2d 1 100 +.fi +.le +.ls 4 7. +To take a 1D image and create a 2D image in which each column is the same: + +.nf +cl> imstack image1d image2d +cl> imtranspose image2d image2d +cl> blkrep image2d image2d 100 1 +.fi +.le + +.ih +SEE ALSO +blkavg, imstack, magnify +.endhelp diff --git a/pkg/images/imgeom/doc/im3dtran.hlp b/pkg/images/imgeom/doc/im3dtran.hlp new file mode 100644 index 00000000..68a2b0cd --- /dev/null +++ b/pkg/images/imgeom/doc/im3dtran.hlp @@ -0,0 +1,94 @@ +.help im3dtran Jan97 images.imgeom +.ih +NAME +im3dtran -- Transpose a 3D image +.ih +USAGE +im3dtran input output +.ih +PARAMETERS +.ls input +The input 3d image. +.le +.ls output +The output transposed 3D image. If the output image name is the same as +the input image name then the original image will be overwritten. The number +of output images must equal the number of input images. +.le +.ls new_x = 3 +The new x axis. The default (3) replaces new x with old z. +.le +.ls new_y = 2 +The new y axis = old axis. The default (2) does not change the y axis. +.le +.ls new_z = 1 +The new z axis. The default (1) replaces new z with old x. +.le +.ls len_blk = 128 +The size in pixels of the linear internal subraster. Im3dtran will try to +transpose a subraster up to len_blk cubed at one time. Runtime is much +faster with larger \fBlen_blk\fR, but the task may run out of memory. +.le +.ls verbose = yes +Print messages about actions taken by the task. +.le + +.ih +DESCRIPTION + +IM3DTRAN transposes the input images \fIinput\fR in 3 dimensions and +writes the transposed images to \fIoutput\fR. The 6 possible axis +mappings are specified by setting the parameters \fInew_x\fR, \fInew_y\fR, +and \fInew_z\fR. + +IM3DTRAN can be used to rotate a datacube 90 degrees in any direction by +combining the transpose operation with an axis flip. For +example, Consider a datacube is visualized with its origin at the lower +left front +when seen by the viewer, with its abscissa being the x axis, its ordinate +the y axis, and its depth the z axis, with z increasing away from the viewer +or into the datacube [this +is a left-handed coordinate system]. To rotate the datacube +by 90 degrees clockwise about the y axis when viewed from the +y direction; +the old z axis must become the new x axis, and the old x axis must become +the new z axis, while the new y axis remains old y axis. In this case the +axis that must be flipped prior to transposition is the x axis as shown +in example 2. + +The parameter \fBlen_blk\fR controls how much memory is used during the +transpose operation. \fBlen_blk\fR elements are used in each axis at a +time, or a cube len_blk elements on a side. If \fBlen_blk\fR is too large, +the task will abort with an "out of memory" error. If it is too small, +the task can take a very long time to run. The maximum size of len_blk +depends on how much memory is available at the time IM3DTRAN is run, +and the size and datatype of the image to be transposed. + +.ih +EXAMPLES + +1. Transpose axes 1 2 and 3 of a list of input images to axes 2 1 and 3 of +a list of output images. + +.nf + cl> im3dtran image1,image2,image3 tr1,tr2,tr3 2 1 3 +.fi + +2. For an input datacube with columns = x = abscissa, lines = y = ordinate, +and bands = z = depth increasing away from viewer, and with the image +origin at the lower left front, rotate datacube 90 degrees clockwise +around the y axis when viewed from +y (top): + +.nf + cl> im3dtran input[-*,*,*] output 3 2 1 +.fi + +.ih +TIMINGS + +.ih +BUGS + +.ih +SEE ALSO +imtranspose, imjoin, imstack, imslice +.endhelp diff --git a/pkg/images/imgeom/doc/imlintran.hlp b/pkg/images/imgeom/doc/imlintran.hlp new file mode 100644 index 00000000..f595ffda --- /dev/null +++ b/pkg/images/imgeom/doc/imlintran.hlp @@ -0,0 +1,184 @@ +.help imlintran Dec98 images.imgeom +.ih +NAME +imlintran -- shift, scale, rotate a list of images +.ih +USAGE +imlintran input output xrotation yrotation xmag ymag +.ih +PARAMETERS +.ls input +List of images to be transformed. +.le +.ls output +List of output images. +.le +.ls xrotation, yrotation +Angle of rotation of points on the image axes in degrees. +Positive angles rotate in a counter-clockwise sense around the x axis. +For a normal coordinate axes rotation xrotation and yrotation should +be the same. A simple y axis flip can be introduced by make yrotation +equal to xrotation plus 180 degrees. An axis skew can be introduced by +making the angle between xrotation and yrotation other than a +multiple of 90 degrees. +.le +.ls xmag, ymag +The number of input pixels per output pixel in x and y. The magnifications +must always be positive numbers. Numbers less than 1 magnify the image; +numbers greater than one reduce the image. +.le +.ls xin = INDEF, yin = INDEF +The origin of the input picture in pixels. Xin and yin default to the center of +the input image. +.le +.ls xout = INDEF, yout = INDEF +The origin of the output image. Xout and yout default to the center of the +output image. +.le +.ls ncols = INDEF, nlines = INDEF +The number of columns and rows in the output image. The default is to +keep the dimensions the same as the input image. If ncols and nrows are +less than or equal to zero then the task computes the number of rows and +columns required to include the whole input image, excluding the effects +of any origin shift. +.le +.ls interpolant = "linear" +The choices are the following. +.ls nearest +Nearest neighbor. +.le +.ls linear +Bilinear interpolation in x and y. +.le +.ls poly3 +Third order interior polynomial in x and y. +.le +.ls poly5 +Fifth order interior polynomial in x and y. +.le +.ls spline3 +Bicubic spline. +.le +.ls sinc +2D sinc interpolation. Users can specify the sinc interpolant width by +appending a width value to the interpolant string, e.g. sinc51 specifies +a 51 by 51 pixel wide sinc interpolant. The sinc width will be rounded up to +the nearest odd number. The default sinc width is 31 by 31. +.le +.ls lsinc +Look-up table sinc interpolation. Users can specify the look-up table sinc +interpolant width by appending a width value to the interpolant string, e.g. +lsinc51 specifies a 51 by 51 pixel wide look-up table sinc interpolant. The user +supplied sinc width will be rounded up to the nearest odd number. The default +sinc width is 31 by 31 pixels. Users can specify the resolution of the lookup +table sinc by appending the look-up table size in square brackets to the +interpolant string, e.g. lsinc51[20] specifies a 20 by 20 element sinc +look-up table interpolant with a pixel resolution of 0.05 pixels in x and y. +The default look-up table size and resolution are 20 by 20 and 0.05 pixels +in x and y respectively. +.le +.ls drizzle +2D drizzle resampling. Users can specify the drizzle pixel fraction in x and y +by appending a value between 0.0 and 1.0 in square brackets to the +interpolant string, e.g. drizzle[0.5]. The default value is 1.0. +The value 0.0 is increased internally to 0.001. Drizzle resampling +with a pixel fraction of 1.0 in x and y is equivalent to fractional pixel +rotated block summing (fluxconserve = yes) or averaging (flux_conserve = no) if +xmag and ymag are > 1.0. +.le +.le +.ls boundary = "nearest" +The choices are: +.ls nearest +Use the value of the nearest boundary pixel. +.le +.ls constant +Use a constant value. +.le +.ls reflect +Generate value by reflecting about the boundary. +.le +.ls wrap +Generate a value by wrapping around to the opposite side of the image. +.le +.le +.ls constant = 0. +The value of the constant for boundary extension. +.le +.ls fluxconserve = yes +Preserve the total image flux? +.le +.ls nxblock = 512, nyblock = 512 +If the size of the output image is less than nxblock by nyblock then +the entire image is transformed at once. Otherwise the output image +is computed in blocks of nxblock by nxblock pixels. +.le +.ih +DESCRIPTION + +IMLINTRAN linearly transforms a the list of images in input using rotation +angles and magnification factors supplied by the user and writes the output +images into output. The coordinate transformation from input to output +image is described below. + +.nf + 1. subtract the origin + + xt = x(input) - xin + yt = y(input) - yin + + 2. scale the image + + xt = xt / xmag + yt = xt / xmag + + 3. rotate the image + + xt = xt * cos (xrotation) - yt * sin (yrotation) + yt = xt * sin (yrotation) + yt * cos (yrotation) + + 4. new orgin + + x(output) = xt + xout + y(output) = yt + yout + +.fi + +The output image gray levels are determined by interpolating in the input +image at the positions of the transformed output pixels using the inverse +of the above transformation. +IMLINTRAN uses the routines in the 2-D interpolation package. + +.ih +TIMINGS +It requires approximately 70 and 290 cpu seconds respectively to linearly +transform a 512 by 512 real image using bilinear and biquintic +interpolation respectively (Vax 11/750 fpa). + +.ih +EXAMPLES + +.nf +1. Rotate an image 45 degrees around its center and magnify + the image by a factor of 2. in each direction. + + cl> imlintran n4151 n4151rm 45.0 45.0 0.50 0.50 + +2. Rotate the axes of an image by 45 degrees around 100. and 100., + shift the orgin to 150. and 150. and flip the y axis. + + cl> imlintran n1068 n1068r 45.0 225.0 1.0 1.0 xin=100. yin=100. \ + >>> xout=150. yout=150. + +3. Rotate an image by 45 degrees and reduce the scale in x and y + by a factor of 1.5 + + cl> imlintran n7026 n7026rm 45.0 45.0 1.5 1.5 +.fi + +.ih +BUGS +.ih +SEE ALSO +imshift, magnify, rotate, lintran, register, geotran, geomap +.endhelp diff --git a/pkg/images/imgeom/doc/imshift.hlp b/pkg/images/imgeom/doc/imshift.hlp new file mode 100644 index 00000000..5d3b3dd5 --- /dev/null +++ b/pkg/images/imgeom/doc/imshift.hlp @@ -0,0 +1,125 @@ +.help imshift Dec98 images.imgeom +.ih +NAME +imshift -- shift a set of images in x and y +.ih +USAGE +imshift input output xshift yshift +.ih +PARAMETERS +.ls input +List of images to be transformed. +.le +.ls output +List of output images. +.le +.ls xshift, yshift +Fractional pixel shift in x and y such that xout = xin + xshift and +yout = yin + yshift. +.le +.ls shifts_file = "" +The name of the text file containing the shifts for each input image. If no +file name is supplied each input image is shifted by \fIxshift\fR and +\fIyshift\fR. Shifts are listed in the text file, 1 set of shifts per image, +with the x and y shift in columns 1 and 2 respectively. The number of +shifts in the file must equal the number of input images. +.le +.ls interp_type = "linear" +The interpolant type use to computed the output shifted image. +The choices are the following: +.ls nearest +nearest neighbor. +.le +.ls linear +bilinear interpolation in x and y. +.le +.ls poly3 +third order interior polynomial in x and y. +.le +.ls poly5 +fifth order interior polynomial in x and y. +.le +.ls spline3 +bicubic spline. +.le +.ls sinc +2D sinc interpolation. Users can specify the sinc interpolant width by +appending a width value to the interpolant string, e.g. sinc51 specifies +a 51 by 51 pixel wide sinc interpolant. The sinc width input by the +user will be rounded up to the nearest odd number. The default sinc width +is 31 by 31. +.le +.ls drizzle +2D drizzle resampling. Users can specify the drizzle pixel fractions in x and y +by appending values between 0.0 and 1.0 in square brackets to the +interpolant string, e.g. drizzle[0.5]. The default value is 1.0. The +value 0.0 is increased to 0.001. Drizzle resampling with a pixel fraction +of 1.0 in x and y is identical to bilinear interpolation. +.le +.le +.ls boundary_type = "nearest" +The choices are: +.ls nearest +Use the value of the nearest boundary pixel. +.le +.ls constant +Use a constant value. +.le +.ls reflect +Generate value by reflecting about the boundary. +.le +.ls wrap +Generate a value by wrapping around to the opposite side of the image. +.le +.le +.ih +DESCRIPTION + +IMSHIFT will shift an image in x and y such that: + +.nf + xout = xin + xshift + yout = yin + yshift + +.fi + +The output image gray levels are determined by interpolating in the input +image at the positions of the shifted output pixels. +IMSHIFT uses the routines in the 2-D interpolator package. + +.ih +EXAMPLES + +1. Shift an image by (+3.2, -4.5) using a biquintic interior polynomial + interpolant and boundary extension. + + cl> imshift vys70 vys70shift 3.2 -4.5 inter=poly5 bound=neare + +2. Shift an image by (-6., 1.2) using bilinear interpolation and + boundary extension. + + cl> imshift ugc1040 ugc1040shift -6.0 1.2 bound=neare + +3. Shift a set of images using shifts listed in the textfile "shifts". + + cl> page shifts + + 3.5 4.86 + -2. 8.9 + 10.1 7.8 + + cl> imshift im1,im2,im3 im1.s,im2.s,im3.s shifts_file=shifts + +.ih +TIMINGS +The time required to shift a 512 by 512 real image by fractional pixel +amounts in x and y is approximately 10, 20, 70, 120, and 120 cpu seconds for the +nearest neighbor, bilinear, bicubic, biquintic and bicubic spline +interpolants respectively (Vax 11/750 fpa). + +.ih +BUGS +.ih +SEE ALSO +shiftlines, magnify, rotate, geomap, geotran, imlintran +.endhelp diff --git a/pkg/images/imgeom/doc/imtrans.hlp b/pkg/images/imgeom/doc/imtrans.hlp new file mode 100644 index 00000000..332cf040 --- /dev/null +++ b/pkg/images/imgeom/doc/imtrans.hlp @@ -0,0 +1,69 @@ +.help imtranspose Aug84 images.imgeom +.ih +NAME +imtranspose -- transpose two dimensional images +.ih +USAGE +.nf +imtranspose input output +.fi +.ih +PARAMETERS +.ls input +List of images to be transposed. +.le +.ls output +List of output transposed images. If the output image name is the same as +the input image name then the output image will replace the input image. +The number of output images must be the same as the number of input images. +.le +.ls len_blk = 512 +The one dimensional length of the transpose blocks. +.le +.ih +DESCRIPTION +Imtranspose transposes the list of images in input by interchanging +their rows and columns and writing the results to images specified in +output. The number of input and output images must be the same. + +The transpose is done in square blocks whose dimensions are equal \fIlen_blk\fR. + +The imtranspose tasks can be used to perform counter-clockwise or +clockwise ninety degree rotations by flipping the y or x axis respectively +in the input image section specification. + +.ih +EXAMPLES +1. To transpose an image: + + cl> imtranspose image1 image2 + +2. To transpose an image in place: + + cl> imtranspose image1 image1 + +3. To rotate an image 90 degrees counter-clockwise and clockwise: + + cl> imtranspose image1[*,-*] image2 + + cl> imtranspose image1[-*,*] image2 + +3. To transpose a set of 3 images listed 1 per line in the file imlist to +the new images trans001, trans002, and trans003: + + cl> imtranspose @imlist trans001,trans002,trans003 + +4. To transpose a set of images in place: + + cl> imtranspose frame* frame* + +5. To rotate an image 90 degrees counter-clockwise in place: + + cl> imtranspose image[*,-*] image +.ih +BUGS + +It is currently not legal to transpose images with a wcs type of MULTISPEC. +.ih +SEE ALSO +.endhelp diff --git a/pkg/images/imgeom/doc/magnify.hlp b/pkg/images/imgeom/doc/magnify.hlp new file mode 100644 index 00000000..8b31817e --- /dev/null +++ b/pkg/images/imgeom/doc/magnify.hlp @@ -0,0 +1,202 @@ +.help magnify Dec98 images.imgeom +.ih +NAME +magnify -- interpolate two dimensional images +.ih +USAGE +.nf +magnify input output xmag ymag +.fi +.ih +PARAMETERS +.ls input +List of one or two dimensional images to be magnified. Image sections are +allowed. Images with an axis containing only one pixel are not magnified. +.le +.ls output +List of output image names. If the output image name is the same as the input +image name then the magnified image replaces the input image. +.le +.ls xmag, ymag +The magnification factors for the first and second image dimensions +respectively. The magnifications need not be integers. Magnifications +greater than 1 increase the size of the image while negative magnifications +less than -1 decrease the size by the specified factor. Magnifications +between -1 and 1 are interpreted as reciprocal magnifications. +.le +.ls x1 = INDEF, x2 = INDEF +The starting and ending coordinates in x in the input image which become +the first and last pixel in x in the magnified image. The values need not +be integers. If indefinite the values default to the first and last pixel +in x of the input image; i.e. a value of 1 and nx. +.le +.ls y1 = INDEF, y2 = INDEF +The starting and ending coordinates in y in the input image which become +the first and last pixel in y in the magnified image. The values need not +be integers. If indefinite the values default to the first and last pixel +in y of the input image; i.e. a value of 1 and ny. +.le +.ls dx = INDEF, dy = INDEF +The intervals between the output pixels in terms of the input image. +The values need not be integers. If these values are specified they +override the magnification factors. +.le +.ls interpolant = "linear" +The interpolant used for rebinning the image. +The choices are the following. +.ls nearest +Nearest neighbor. +.le +.ls linear +Bilinear interpolation in x and y. +.le +.ls poly3 +Third order polynomial in x and y. +.le +.ls poly5 +Fifth order polynomial in x and y. +.le +.ls spline3 +Bicubic spline. +.le +.ls sinc +2D sinc interpolation. Users can specify the sinc interpolant width by +appending a width value to the interpolant string, e.g. sinc51 specifies +a 51 by 51 pixel wide sinc interpolant. The sinc width will be rounded up to +the nearest odd number. The default sinc width is 31 by 31. +.le +.ls lsinc +Look-up table sinc interpolation. Users can specify the look-up table sinc +interpolant width by appending a width value to the interpolant string, e.g. +lsinc51 specifies a 51 by 51 pixel wide look-up table sinc interpolant. The user +supplied sinc width will be rounded up to the nearest odd number. The default +sinc width is 31 by 31 pixels. Users can specify the resolution of the lookup +table sinc by appending the look-up table size in square brackets to the +interpolant string, e.g. lsinc51[20] specifies a 20 by 20 element sinc +look-up table interpolant with a pixel resolution of 0.05 pixels in x and y. +The default look-up table size and resolution are 20 by 20 and 0.05 pixels +in x and y respectively. +.le +.ls drizzle +2D drizzle resampling. Users can specify the drizzle pixel fraction in x and y +by appending a value between 0.0 and 1.0 in square brackets to the +interpolant string, e.g. drizzle[0.5]. The default value is 1.0. +The value 0.0 is increased internally to 0.001. Drizzle resampling +with a pixel fraction of 1.0 in x and y is equivalent to fractional pixel +block summing (fluxconserve = yes) or averaging (flux_conserve = no) if +xmag and ymag are < 1.0. +.le +.le +.ls boundary = "nearest" +Boundary extension type for references to pixels outside the bounds of the +input image. The choices are: +.ls nearest +Use the value of the nearest boundary pixel. +.le +.ls constant +Use a constant value. +.le +.ls reflect +Generate value by reflecting about the boundary. +.le +.ls wrap +Generate a value by wrapping around to the opposite side of the image. +.le +.le +.ls constant = 0. +Constant value for constant boundary extension. +.le +.ls fluxconserve = yes +Preserve the total image flux. +.le +.ls logfile = STDOUT +Log file for recording information about the magnification. A null +logfile may be used to turn off log information. +.le +.ih +DESCRIPTION +The list of input images are expanded or contracted by interpolation +to form the output images. The output image names are specified by the +output list. The number of output image names must be the +same as the number of input images. An output image name may be the same +as the corresponding input image in which case the magnified image replaces +the input image. The input images must be one or two dimensional and each +axis must be of at least length 2 (i.e. there have to be distinct +endpoints between which to interpolate). + +The magnification factor determines the pixel step size or interval. +Positive magnifications are related to the step size as the reciprocal; +for example a magnification of 2.5 implies a step size of .4 and a +magnification of .2 implies a step size of 5. Negative magnifications +are related to the step size as the absolute value; for example a +magnification of -2.2 implies a step size of 2.2. This definition +frees the user from dealing with reciprocals and irrational numbers. +Note that the step size may be specified directly with the parameters +\fIdx\fR and \fIdy\fR, in which case the magnification factor is +not required. + +If fluxconserve = yes, the magnification is approximately flux conserving +in that the image values are scaled by the ratio of the output to the input +pixel areas; i.e dx * dy. + +In the default case with only the magnifications specified the full +image is expanded or contracted. By specifying additional parameters +the size and origin of the output image may be changed. Only those +parameters to be fixed need to be specified and the values of the +remaining parameters are either determined from these values or +default as indicated in the PARAMETERS section. + +The user may select the type of two dimensional interpolation and boundary +extension to be used. Note that the image interpolation is performed on +the boundary extended input image. Thus, boundary extensions which are +discontinuous (constant and wrap) may introduce interpolation errors. +.ih +EXAMPLES +1. To expand an image by a factor of 2.5: + + cl> magnify imagein imageout 2.5 2.5 + +2. To subsample the lines of an image in steps of 3.5: + + cl> magnify imagein imageout dx=3.5 dy=1 + +3. To magnify the central part of an image by 2 into a 11 by 31 image: + +.nf + cl> magnify imagein imageout 2 2 x1=25.3 x2=30.3 \ + >>> y1=20 y2=35 +.fi + +4. To use a higher order interpolator with wrap around boundary extension: + +.nf + cl> magnify imagein imageout 2 2 x1=-10 y1=-10 \ + >>> interpolation=spline3 boundary=wrap +.fi + +It is important to remember that the magnification affects the pixel intervals! +This means that the number of pixels in an expanded image is not simply +a multiple of the original number. The following example illustrates this +point. Begin with an image which is 100 by 10. This means the +x coordinates run between 1 and 100 and the y coordinates run between 1 and +10 with a pixel interval of 1. + +Let's magnify the x axis by 0.5 and the y axis by 2. +The output pixel intervals, in terms of the input pixel intervals, +are then 2 and 0.5. This means the output x pixels are at +1, 3, 5, etc. and output y pixels are at 1, 1.5, 2, 2.5, etc., again in +terms of the input pixel coordinates. The last output x pixel is then +at 99 in the input coordinates and the number of pixels is 50. For the +y axis the last output pixel is at 10 in the input coordinates and the +number of pixels between 1 and 10 in intervals of 0.5 is 19! Thus, the +final image is 50 by 19 and not 50 by 20 which you would get if you +multiplied the axis lengths by the magnification factors. + +A more complex example is given above in which x1=25.3, +x2=30.3, y1=20, and y2=35 with magnification factors of 2. +It is important to understand why the output image is 11 by 31 and +what the pixel coordinates are in terms of the input pixel coordinates. +.ih +SEE ALSO +imshift, blkavg, rotate, imlintran, register, geotran, geomap +.endhelp diff --git a/pkg/images/imgeom/doc/rotate.hlp b/pkg/images/imgeom/doc/rotate.hlp new file mode 100644 index 00000000..2b534738 --- /dev/null +++ b/pkg/images/imgeom/doc/rotate.hlp @@ -0,0 +1,164 @@ +.help rotate Dec98 images.imgeom +.ih +NAME +rotate -- rotate and shift a list of images +.ih +USAGE +rotate input output rotation +.ih +PARAMETERS +.ls input +List of images to be rotated. +.le +.ls output +List of output images. +.le +.ls rotation +Angle of rotation of the image in degrees. Positive angles will rotate +the image counter-clockwise from the x axis. +.le +.ls xin = INDEF, yin = INDEF +The origin of the rotation in pixels. Xin and yin default to the center of +the input image. +.le +.ls xout = INDEF, yout = INDEF +The origin of the output image. Xout and yout default to the center of the +output image. +.le +.ls ncols = INDEF, nlines = INDEF +The number of columns and rows in the output image. The default is to +keep the dimensions the same as the input image. If ncols and nrows is +less then or equal to zero the program will compute the number of columns +and rows needed to include the whole image, excluding the effects of +any origin shifts. +.le +.ls interpolant = "linear" +The interpolant. The options are the following: +.ls nearest +Nearest neighbor. +.le +.ls linear +Bilinear interpolation in x and y. +.le +.ls poly3 +Third order polynomial in x and y. +.le +.ls poly5 +Fifth order polynomial in x and y. +.le +.ls spline3 +Bicubic spline. +.le +.ls sinc +2D sinc interpolation. Users can specify the sinc interpolant width by +appending a width value to the interpolant string, e.g. sinc51 specifies +a 51 by 51 pixel wide sinc interpolant. The sinc width will be rounded up to +the nearest odd number. The default sinc width is 31 by 31. +.le +.ls lsinc +Look-up table sinc interpolation. Users can specify the look-up table sinc +interpolant width by appending a width value to the interpolant string, e.g. +lsinc51 specifies a 51 by 51 pixel wide look-up table sinc interpolant. The user +supplied sinc width will be rounded up to the nearest odd number. The default +sinc width is 31 by 31 pixels. Users can specify the resolution of the lookup +table sinc by appending the look-up table size in square brackets to the +interpolant string, e.g. lsinc51[20] specifies a 20 by 20 element sinc +look-up table interpolant with a pixel resolution of 0.05 pixels in x and y. +The default look-up table size and resolution are 20 by 20 and 0.05 pixels +in x and y respectively. +.le +.ls drizzle +2D drizzle resampling. Users can specify the drizzle pixel fraction in x and y +by appending a value between 0.0 and 1.0 in square brackets to the +interpolant string, e.g. drizzle[0.5]. The default value is 1.0. +The value 0.0 is increased internally to 0.001. Drizzle resampling +with a pixel fraction of 1.0 in x and y is equivalent to fractional pixel +rotated block summing (fluxconserve = yes) or averaging (flux_conserve = no) if +xmag and ymag are > 1.0. +.le +.le +.ls boundary = "nearest" +The choices are: +.ls nearest +Use the value of the nearest boundary pixel. +.le +.ls constant +Use a constant value. +.le +.ls reflect +Generate a value by reflecting around the boundary. +.le +.ls wrap +Generate a value by wrapping around to the opposite side of the image. +.le +.le +.ls constant = 0. +The value of the constant for constant boundary extension. +.le +.ls nxblock = 512, nyblock = 512 +If the dimensions of the output image are less than nxblock and nyblock +then the entire image is rotated at once. Otherwise nxblock by nyblock +segments of the image are rotated. +.le +.ih +DESCRIPTION + +ROTATE rotates the list of images in input by rotation degrees and writes +the output to the images specified by output. The origins of the input and +output images may be specified by setting xin, yin, xout and yout. The +transformation is described below. + +.nf + xt = (x - xin) * cos (rotation) - (y - yin) * sin (rotation) + xout + yt = (x - xin) * sin (rotation) + (y - yin) * cos (rotation) + yout + +.fi + +The output image gray levels are determined by interpolating in the input +image at the positions of the transformed output pixels. ROTATE uses the +routines in the 2-D interpolation package. + +.ih +EXAMPLES + +.nf +1. Rotate an image 45 degrees around its center. + + cl> rotate m51 m51r45 45.0 + +2. Rotate an image by 45 degrees around (100., 100.) and + shift the origin to (150., 150.0) using bicubic interpolation. + + cl> rotate m92 m92r45 45.0 xin=100. yin=100. xout=150. yout=150.\ + >>> interp=poly3 + +3. Rotate an image 90 degrees counter-clockwise and clockwise around its + center. Note the use of imtranspose and image section notation. + + cl> imtranspose m92[*,-*] m92d90 + + cl> imtranspose m92[-*,*] m92d270 + +4. Rotate an image 180 degrees counter-clockwise. Note the use of imcopy + and image section notation. + + cl> imcopy m92[-*,-*] m92d180 +.fi + +.ih +TIMINGS +It requires approximately 70 and 290 cpu seconds to rotate a 512 by 512 +real image using bilinear and biquintic interpolation respectively +(Vax 11/750 fpa). +.ih +BUGS +The interpolation operation is done in real arithmetic. However the output +type of the pixels is set equal to the input type. This can lead to truncation +problems for integer images. + +Simple 90, 180, 270 etc degree rotations are best performed using the +imtranspose task and/or image section notation. +.ih +SEE ALSO +imtranspose, imshift, magnify, lintran, geotran, geomap +.endhelp diff --git a/pkg/images/imgeom/doc/shiftlines.hlp b/pkg/images/imgeom/doc/shiftlines.hlp new file mode 100644 index 00000000..856ab3a8 --- /dev/null +++ b/pkg/images/imgeom/doc/shiftlines.hlp @@ -0,0 +1,119 @@ +.help shiftlines Dec98 images.imgeom +.ih +NAME +shiftlines -- shift lines in a list of images +.ih +USAGE +.nf +shiftlines input output shift +.fi +.ih +PARAMETERS +.ls input +List of images to be shifted. Image sections are allowed. +.le +.ls output +List of output image names. If the output image name is the same as the input +image name then the shifted image replaces the input image. +.le +.ls shift +Shift in pixels. +.le +.ls interp_type = "linear" +The interpolant type use to computed the output shifted image. +The choices are the following: +.ls nearest +nearest neighbor interpolation. +.le +.ls linear +linear interpolation in x. +.le +.ls poly3 +third order interior polynomial in x. +.le +.ls poly5 +fifth order interior polynomial in x. +.le +.ls spline3 +cubic spline in x. +.le +.ls sinc +sinc interpolation in x. Users can specify the sinc interpolant width by +appending a width value to the interpolant string, e.g. sinc51 specifies +a 51 pixel wide sinc interpolant. The sinc width input by the user will +be rounded up to the nearest odd number. The default sinc width +is 31 pixels. +.le +.ls drizzle +1D drizzle resampling. Users can specify the drizzle pixel fraction +by appending a value between 0.0 and 1.0 in square brackets to the +interpolant string, e.g. drizzle[0.5]. The default value is 1.0. The +value 0.0 is increased to 0.001. Drizzle resampling with a pixel fraction +of 1.0 is identical to linear interpolation. +.le +.le +.ls boundary_type = "nearest" +Boundary condition for shifts outside the input image. +The minimum match abbreviated choices are: +.ls "nearest" +Use the values of the nearest boundary pixel. +.le +.ls "wrap" +Generate a value by wrapping around to the opposite boundary. +.le +.ls "reflect" +Generate a value by reflecting around the boundary +.le +.ls "constant" +Use a user supplied constant pixel value. +.le +.le +.ls constant = "0.0" +The constant for constant boundary extension. +.le +.ih +DESCRIPTION +The list of images in \fIinput\fR is shifted by the amount \fIshift\fR +and copied to the list of output images \fIoutput\fR. +The number of output image names must be the same as the number of input +images. An output image name may be the same as the corresponding +input image in which case the shifted image replaces the input image. + +The shift is defined by the following relation. + + xout = xint + shift + +Features in the input image are moved to higher columns when the shift +is positive and to lower columns when the shift is negative. For example, +to shift a feature at column 10 to column 12 the shift is 2.0. The task +has been optimized for integral pixel shifts. + +There are five choices for the one dimensional image interpolation +which is selected with the parameter \fIinterp_type\fR. +The value of the output pixels corresponding to input pixel positions +outside the boundaries of the image is determined by the parameter +\fIboundary_type\fR. + +.ih +EXAMPLES + +1. Shift the lines of an image by 0.25 pixels to the right. + + cl> shiftlines imagein imageout 0.25 + +2. Shift the lines of an image by -.3 pixels using cubic spline interpolation +and replace the input image by the output image. + + cl> shiftlines image image -.3 interp=spline3 + +.ih +TIMINGS +It requires approximately 28 and 59 seconds to shift a 512 square image +using linear and cubic spline interpolation respectively +(Vax 11/750 with fpa). +.ih +BUGS +.ih +SEE ALSO +imshift, magnify, rotate, imlintran, blkrep, blkav, geotran +.endhelp diff --git a/pkg/images/imgeom/im3dtran.par b/pkg/images/imgeom/im3dtran.par new file mode 100644 index 00000000..3c993815 --- /dev/null +++ b/pkg/images/imgeom/im3dtran.par @@ -0,0 +1,9 @@ +# Parameter list for the IM3DTRAN task + +input,s,a,,,,Input 3d image +output,s,a,,,,Output 3d image +new_x,i,a,3,1,3,"New x axis" +new_y,i,a,2,1,3,"New y axis" +new_z,i,a,1,1,3,"New z axis" +len_blk,i,h,128,,,Working block size in pixels +verbose,b,h,yes,,,Print messages about actions taken by the task ? diff --git a/pkg/images/imgeom/imgeom.cl b/pkg/images/imgeom/imgeom.cl new file mode 100644 index 00000000..8642c8f7 --- /dev/null +++ b/pkg/images/imgeom/imgeom.cl @@ -0,0 +1,30 @@ +#{ IMGEOM -- The Image Geometric Transformation Package. + +set imgeom = "images$imgeom/" + +package imgeom + +# Tasks. + +task blkavg, + blkrep, + imshift, + imtranspose, + im3dtran, + magnify, + shiftlines = "imgeom$x_images.e" + +# Tasks in other packages + +# Geotran is used by the imlintran and rotate tasks. + +task geotran = "immatch$x_images.e" +hidetask geotran + +# Scripts + +task imlintran = "imgeom$imlintran.cl" +task rotate = "imgeom$rotate.cl" + + +clbye() diff --git a/pkg/images/imgeom/imgeom.hd b/pkg/images/imgeom/imgeom.hd new file mode 100644 index 00000000..c9ed726a --- /dev/null +++ b/pkg/images/imgeom/imgeom.hd @@ -0,0 +1,16 @@ +# Help directory for the IMGEOM package + +$doc = "images$imgeom/doc/" +$src = "images$imgeom/src/" + +blkavg hlp=doc$blkavg.hlp, src=src$t_blkavg.x +blkrep hlp=doc$blkrep.hlp, src=src$t_blkrep.x +imlintran hlp=doc$imlintran.hlp, src=imgeom$imlintran.cl +imshift hlp=doc$imshift.hlp, src=src$t_imshift.x +imtranspose hlp=doc$imtrans.hlp, src=src$t_imtrans.x +im3dtran hlp=doc$im3dtran.hlp, src=src$t_im3dtran.x +magnify hlp=doc$magnify.hlp, src=src$t_magnify.x +rotate hlp=doc$rotate.hlp, src=imgeom$rotate.cl +shiftlines hlp=doc$shiftlines.hlp, src=src$t_shiftlines.x +revisions sys=Revisions + diff --git a/pkg/images/imgeom/imgeom.men b/pkg/images/imgeom/imgeom.men new file mode 100644 index 00000000..565f5b6d --- /dev/null +++ b/pkg/images/imgeom/imgeom.men @@ -0,0 +1,9 @@ + blkavg - Block average or sum a list of N-D images + blkrep - Block replicate a list of N-D images + imlintran - Linearly transform a list of 2-D images + imshift - Shift a list of 1-D or 2-D images + imtranspose - Transpose a list of 2-D images + im3dtran - Transpose a list of 3-D images + magnify - Magnify a list of 1-D or 2-D images + rotate - Rotate and shift a list of 2-D images + shiftlines - Shift the lines of a list of N-D images diff --git a/pkg/images/imgeom/imgeom.par b/pkg/images/imgeom/imgeom.par new file mode 100644 index 00000000..cef3f3ff --- /dev/null +++ b/pkg/images/imgeom/imgeom.par @@ -0,0 +1 @@ +version,s,h,"Jan97" diff --git a/pkg/images/imgeom/imlintran.cl b/pkg/images/imgeom/imlintran.cl new file mode 100644 index 00000000..59db414f --- /dev/null +++ b/pkg/images/imgeom/imlintran.cl @@ -0,0 +1,50 @@ +# IMLINTRAN -- Linearly transform and image by calling the GEOTRAN task +# with the appropriate parameters. + +procedure imlintran (input, output, xrotation, yrotation, xmag, ymag, xin, yin, + xout, yout, ncols, nlines, interpolant, boundary, constant, + fluxconserve, nxblock, nyblock, verbose) + +string input +string output +real xrotation +real yrotation +real xmag +real ymag +real xin +real yin +real xout +real yout +real ncols +real nlines +string interpolant +string boundary +real constant +bool fluxconserve +int nxblock +int nyblock +bool verbose + + +begin + # Declare local variables. + string tinput, toutput + real txrotation, tyrotation + + # Get the parameters. + tinput = input + toutput = output + txrotation = xrotation + tyrotation = yrotation + + # Call GEOTRAN. + geotran (input=tinput, output=toutput, database="", + xrotation=txrotation, yrotation=tyrotation, xin=xin, yin=yin, + xout=xout, yout=yout, xshift=INDEF, yshift=INDEF, xmin=1.0, + xmax=ncols, ymin=1.0, ymax=nlines, xscale=1.0, yscale=1.0, + ncols=INDEF, nlines=INDEF, xmag=xmag, ymag=ymag, + interpolant=interpolant, boundary=boundary, constant=constant, + xsample=1., ysample=1., fluxconserve=fluxconserve, nxblock=nxblock, + nyblock=nyblock, verbose=verbose) +end + diff --git a/pkg/images/imgeom/imlintran.par b/pkg/images/imgeom/imlintran.par new file mode 100644 index 00000000..a45ed454 --- /dev/null +++ b/pkg/images/imgeom/imlintran.par @@ -0,0 +1,30 @@ +# IMLINTRAN Parameters + +# Required parameters. +input,f,a,,,,Input data +output,f,a,,,,Output data +xrotation,r,a,,,,X rotation angle in degrees +yrotation,r,a,,,,Y rotation angle in degrees +xmag,r,a,1.0,0.0,,X output pixels per input pixel +ymag,r,a,1.0,0.0,,Y output pixels per input pixel + +# Change transformation parameters. +xin,r,h,INDEF,1.,,X origin of input image in pixels +yin,r,h,INDEF,1.,,Y origin of input image in pixels +xout,r,h,INDEF,1.,,X origin of output image in pixels +yout,r,h,INDEF,1.,,Y origin of output image in pixels +ncols,r,h,INDEF,,,Number of columns in the output image +nlines,r,h,INDEF,,,Number of lines in the output image + +# Coordinate surface and image interpolation parameters. +interpolant,s,h,'linear',,,'Interpolant (nearest,linear,poly3,poly5,spline3,sinc,lsinc,drizzle)' +boundary,s,h,'nearest',|nearest|constant|reflect|wrap|,,'Boundary extension (nearest,constant,reflect,wrap)' +constant,r,h,0.,,,Constant for constant boundary extension +fluxconserve,b,h,yes,,,Preserve image flux ? + +# Transformation blocking factors. +nxblock,i,h,512,,,X dimension of working block size in pixels +nyblock,i,h,512,,,Y dimension of working block size in pixels +verbose,b,h,yes,,,Print messages about the progress of the task ? + +mode,s,h,'ql' diff --git a/pkg/images/imgeom/imshift.par b/pkg/images/imgeom/imshift.par new file mode 100644 index 00000000..bcd630c1 --- /dev/null +++ b/pkg/images/imgeom/imshift.par @@ -0,0 +1,11 @@ +# IMSHIFT + +input,f,a,,,,Input images to be fit +output,f,a,,,,Output images +xshift,r,a,,,,Fractional pixel shift in x +yshift,r,a,,,,Fractional pixel shift in y +shifts_file,f,h,"",,,Text file containing shifts for each image +interp_type,s,h,'linear',,,'Interpolant (nearest,linear,poly3,poly5,spline3,sinc,drizzle)' +boundary_type,s,h,'nearest',"|nearest|constant|reflect|wrap|",,'Boundary (constant,nearest,reflect,wrap)' +constant,r,h,0,,,Constant for boundary extension +mode,s,h,'ql' diff --git a/pkg/images/imgeom/imtranspose.par b/pkg/images/imgeom/imtranspose.par new file mode 100644 index 00000000..d77fe522 --- /dev/null +++ b/pkg/images/imgeom/imtranspose.par @@ -0,0 +1,3 @@ +input,s,a,,,,Images to be transposed +output,s,a,,,,Output image names +len_blk,i,h,512,,,Size in pixels of internal subraster diff --git a/pkg/images/imgeom/junk.cl b/pkg/images/imgeom/junk.cl new file mode 100644 index 00000000..59db414f --- /dev/null +++ b/pkg/images/imgeom/junk.cl @@ -0,0 +1,50 @@ +# IMLINTRAN -- Linearly transform and image by calling the GEOTRAN task +# with the appropriate parameters. + +procedure imlintran (input, output, xrotation, yrotation, xmag, ymag, xin, yin, + xout, yout, ncols, nlines, interpolant, boundary, constant, + fluxconserve, nxblock, nyblock, verbose) + +string input +string output +real xrotation +real yrotation +real xmag +real ymag +real xin +real yin +real xout +real yout +real ncols +real nlines +string interpolant +string boundary +real constant +bool fluxconserve +int nxblock +int nyblock +bool verbose + + +begin + # Declare local variables. + string tinput, toutput + real txrotation, tyrotation + + # Get the parameters. + tinput = input + toutput = output + txrotation = xrotation + tyrotation = yrotation + + # Call GEOTRAN. + geotran (input=tinput, output=toutput, database="", + xrotation=txrotation, yrotation=tyrotation, xin=xin, yin=yin, + xout=xout, yout=yout, xshift=INDEF, yshift=INDEF, xmin=1.0, + xmax=ncols, ymin=1.0, ymax=nlines, xscale=1.0, yscale=1.0, + ncols=INDEF, nlines=INDEF, xmag=xmag, ymag=ymag, + interpolant=interpolant, boundary=boundary, constant=constant, + xsample=1., ysample=1., fluxconserve=fluxconserve, nxblock=nxblock, + nyblock=nyblock, verbose=verbose) +end + diff --git a/pkg/images/imgeom/magnify.par b/pkg/images/imgeom/magnify.par new file mode 100644 index 00000000..8f1fbf5f --- /dev/null +++ b/pkg/images/imgeom/magnify.par @@ -0,0 +1,17 @@ +# Magnify parameters + +input,s,a,,,,Input two dimensional images +output,s,a,,,,Output magnified images +xmag,r,a,,,,X magnification factor +ymag,r,a,,,,Y magnification factor +x1,r,h,INDEF,,,X window origin relative to input image +x2,r,h,INDEF,,,X window end point relative to input image +dx,r,h,INDEF,0.,,X Pixel interval relative to input image +y1,r,h,INDEF,,,Y window origin relative to input image +y2,r,h,INDEF,,,Y window end point relative to input image +dy,r,h,INDEF,0.,,Y Pixel interval relative to input image +interpolation,s,h,linear,,,"Interpolation type(nearest,linear,poly3,poly5,spline3,sinc,lsinc,drizzle" +boundary,s,h,nearest,"|nearest|constant|reflect|wrap|",,"Boundary extension type (constant,nearest,reflect,wrap)" +constant,r,h,0.,,,Boundary extension constant +fluxconserve,b,h,yes,,,Preserve total image flux? +logfile,f,h,STDOUT,,,Log file diff --git a/pkg/images/imgeom/mkpkg b/pkg/images/imgeom/mkpkg new file mode 100644 index 00000000..b196aebc --- /dev/null +++ b/pkg/images/imgeom/mkpkg @@ -0,0 +1,5 @@ +# MKPKG for the IMGEOTRAN Package + +libpkg.a: + @src + ; diff --git a/pkg/images/imgeom/rotate.cl b/pkg/images/imgeom/rotate.cl new file mode 100644 index 00000000..249fb50e --- /dev/null +++ b/pkg/images/imgeom/rotate.cl @@ -0,0 +1,43 @@ +# ROTATE -- Rotate an image by calling the GEOTRAN task with the appropriate +# parameters. + +procedure rotate (input, output, rotation, xin, yin, xout, yout, ncols, nlines, + interpolant, boundary, constant, nxblock, nyblock, verbose) + +string input +string output +real rotation +real xin +real yin +real xout +real yout +real ncols +real nlines +string interpolant +string boundary +real constant +int nxblock +int nyblock +bool verbose + + +begin + # Declare local variables. + string tinput, toutput + real trotation + + # Get the parameters. + tinput = input + toutput = output + trotation = rotation + + # Call GEOTRAN. + geotran (input=tinput, output=toutput, database="", xrotation=trotation, + yrotation=trotation, xin=xin, yin=yin, xout=xout, yout=yout, + xshift=INDEF, yshift=INDEF, xmin=1.0, xmax=ncols, ymin=1.0, + ymax=nlines, xscale=1.0, yscale=1.0, ncols=INDEF, nlines=INDEF, + xmag=INDEF, ymag=INDEF, interpolant=interpolant, + boundary=boundary, constant=constant, xsample=1., ysample=1., + fluxconserve=no, nxblock=nxblock, nyblock= nyblock, verbose=verbose) +end + diff --git a/pkg/images/imgeom/rotate.par b/pkg/images/imgeom/rotate.par new file mode 100644 index 00000000..9e44b99b --- /dev/null +++ b/pkg/images/imgeom/rotate.par @@ -0,0 +1,24 @@ +# ROTATE Parameters + +# Required parameters. +input,f,a,,,,Input data +output,f,a,,,,Output data +rotation,r,a,,,,Rotation angle in degrees + +# Transformation parameters. +xin,r,h,INDEF,,,X origin of input image in pixels +yin,r,h,INDEF,,,Y origin of input image in pixels +xout,r,h,INDEF,,,X origin of output image in pixels +yout,r,h,INDEF,,,Y origin of output image in pixels +ncols,r,h,INDEF,,,Number of columns in the output image +nlines,r,h,INDEF,,,Number of lines in the output image + +# Image interpolation parameters. +interpolant,s,h,'linear',,,'Interpolant (nearest,linear,poly3,poly5,spline3,sinc,lsinc,drizzle)' +boundary,s,h,'nearest',|nearest|constant|reflect|wrap|,,'Boundary extension (nearest,constant,reflect,wrap)' +constant,r,h,0.,,,Constant for constant boundary extension +nxblock,i,h,512,,,X dimension of working block size in pixels +nyblock,i,h,512,,,Y dimension of working block size in pixels +verbose,b,h,yes,,,Print messages about the progress of the task ? + +mode,s,h,'ql' diff --git a/pkg/images/imgeom/shiftlines.par b/pkg/images/imgeom/shiftlines.par new file mode 100644 index 00000000..957f8b98 --- /dev/null +++ b/pkg/images/imgeom/shiftlines.par @@ -0,0 +1,9 @@ +# Task parameters for SHIFTLINES. + +input,s,a,,,,Input images +output,s,a,,,,Output images +shift,r,a,,,,Line shift in pixels +interp_type,s,h,"linear",,,'Interpolant (nearest,linear,poly3,poly5,spline3,sinc,drizzle)' +boundary_type,s,h,"nearest","|nearest|constant|reflect|wrap|",,'Boundary (constant,nearest,reflect,wrap)' +constant,r,h,0.0,,,Constant for boundary extension +mode,s,h,ql diff --git a/pkg/images/imgeom/src/blkav.gx b/pkg/images/imgeom/src/blkav.gx new file mode 100644 index 00000000..9c83ebab --- /dev/null +++ b/pkg/images/imgeom/src/blkav.gx @@ -0,0 +1,131 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <error.h> + +define AVG 1 # operation = arithmetic average +define SUM 2 # operation = arithmetic sum + +# change to (lrdx) in future +$for (lrd) + +# BLKAVG -- Block average or sum on n-dimensional images. +# For example, block averaging by a factor of 2 means that pixels 1 and 2 +# are averaged to produce the first output pixel, 3 and 4 are averaged to +# produce the second output pixel, and so on. Remainder pixels at the end +# of a line, col, etc. are averaged correctly. + +procedure blkav$t (im1, im2, blkfac, option) + +pointer im1 # input image +pointer im2 # output image +int blkfac[IM_MAXDIM] # blocking factors +int option # block operation (average, sum, ...) + +int num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx +int junk, num_ilines, num_olines, oline +long blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM] +long vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM] +$if (datatype != l) +PIXEL sum +$else +real sum +$endif +pointer sp, accum_ptr, iline_ptr, oline_ptr + +int blkcomp(), imggs$t(), impnl$t() +errchk imggs$t(), impnl$t() + +begin + call smark (sp) + call salloc (accum_ptr, IM_LEN(im1, 1), TY_PIXEL) + + # Initialize; input and output vectors, block counters. + ndim = IM_NDIM(im1) + nfull_blkx = IM_LEN(im1, 1) / blkfac[1] + blkin_s[1] = long(1) + blkin_e[1] = long(IM_LEN(im1, 1)) + vin_s[1] = blkin_s[1] + vin_e[1] = blkin_e[1] + + do i = 1, ndim { + num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i] + IM_LEN(im2, i) = num_oblks[i] + } + num_olines = 1 + do i = 2, ndim + num_olines = num_olines * num_oblks[i] + call amovkl (long(1), vout, IM_MAXDIM) + + # For each sequential output-image line, ... + do oline = 1, num_olines { + + call aclr$t (Mem$t[accum_ptr], IM_LEN(im1, 1)) + nlines_in_sum = 0 + + # Compute block vector; initialize dim>1 input vector. + num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e, + vin_s, vin_e) + + # Output pointer; note impnl$t returns vout for NEXT output line. + junk = impnl$t (im2, oline_ptr, vout) + + # For all input lines mapping to current output line, ... + do i = 1, num_ilines { + # Get line from input image. + iline_ptr = imggs$t (im1, vin_s, vin_e, ndim) + + # Increment general section input vector between block bounds. + do dim = 2, ndim { + if (vin_s[dim] < blkin_e[dim]) { + vin_s[dim] = vin_s[dim] + 1 + vin_e[dim] = vin_s[dim] + break + } else { + vin_s[dim] = blkin_s[dim] + vin_e[dim] = vin_s[dim] + } + } + + # Accumulate line into block sum. Keep track of no. of + # lines in sum so that we can compute block average later. + + call aadd$t (Mem$t[iline_ptr], Mem$t[accum_ptr], + Mem$t[accum_ptr], IM_LEN(im1,1)) + nlines_in_sum = nlines_in_sum + 1 + } + + # We now have a templine of sums; average/sum into output buffer + # first the full blocks using a vop. + if (option == AVG) + call abav$t (Mem$t[accum_ptr], Mem$t[oline_ptr], nfull_blkx, + blkfac[1]) + else + call absu$t (Mem$t[accum_ptr], Mem$t[oline_ptr], nfull_blkx, + blkfac[1]) + + # Now average/sum the final partial block in x, if any. + if (nfull_blkx < num_oblks[1]) { + sum = 0$f + count = 0 + do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) { + sum = sum + Mem$t[accum_ptr+i-1] + count = count + 1 + } + if (option == AVG) + Mem$t[oline_ptr+num_oblks[1]-1] = sum / count + else + Mem$t[oline_ptr+num_oblks[1]-1] = sum + } + + # Block average into output line from the sum of all lines block + # averaged in X. + if (option == AVG) + call adivk$t (Mem$t[oline_ptr], PIXEL(nlines_in_sum), + Mem$t[oline_ptr], num_oblks[1]) + } + + call sfree (sp) +end + +$endfor diff --git a/pkg/images/imgeom/src/blkcomp.x b/pkg/images/imgeom/src/blkcomp.x new file mode 100644 index 00000000..814be4ee --- /dev/null +++ b/pkg/images/imgeom/src/blkcomp.x @@ -0,0 +1,38 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# BLKCOMP -- compute starting and ending input vectors for blocks in each +# dimension. Initialize input-line vectors to block vectors. Return total +# number of input lines mapping to current output line. + +int procedure blkcomp (im1, blkfac, vout, blkin_s, blkin_e, + vin_s, vin_e) + +pointer im1 # pointer to input image descriptor +int blkfac[IM_MAXDIM] # blocking factors for each dimension +long vout[IM_MAXDIM] # output image vectors for each dimension +long blkin_s[IM_MAXDIM] # index of starting block for each dimension +long blkin_e[IM_MAXDIM] # index of ending block for each dimension +long vin_s[IM_MAXDIM] # initial starting input vector +long vin_e[IM_MAXDIM] # initial ending input vector + +int num_ilines, dim + +begin + num_ilines = 1 + + # Compute starting and ending indices of input image pixels in each + # dimension mapping to current output line. + + do dim = 2, IM_NDIM(im1) { + blkin_s[dim] = long(1 + (vout[dim] - 1) * blkfac[dim]) + blkin_e[dim] = long(min (IM_LEN(im1,dim), vout[dim] * blkfac[dim])) + vin_s[dim] = blkin_s[dim] + vin_e[dim] = blkin_s[dim] + num_ilines = num_ilines * (blkin_e[dim] - blkin_s[dim] + 1) + } + + return (num_ilines) +end + diff --git a/pkg/images/imgeom/src/blkrp.gx b/pkg/images/imgeom/src/blkrp.gx new file mode 100644 index 00000000..fb297633 --- /dev/null +++ b/pkg/images/imgeom/src/blkrp.gx @@ -0,0 +1,103 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +$for (dlrs) + +# BLKRP -- Block replicate an image. + +procedure blkrp$t (in, out, blkfac) + +pointer in # Input IMIO pointer +pointer out # Output IMIO pointer +int blkfac[ARB] # Block replication factors + +int i, j, ndim, nin, nout +pointer sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout +pointer imgl1$t(), impl1$t(), imgnl$t(), impnl$t() + +begin + call smark (sp) + + ndim = IM_NDIM(in) + nin = IM_LEN(in, 1) + nout = nin * blkfac[1] + IM_LEN(out,1) = nout + + if (ndim == 1) { + # For one dimensional images do the replication directly. + + buf1 = imgl1$t (in) + buf2 = impl1$t (out) + ptrin = buf1 + ptrout = buf2 + do i = 1, nin { + do j = 1, blkfac[1] { + Mem$t[ptrout] = Mem$t[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + + } else { + # For higher dimensional images use line access routines. + + do i = 2, ndim + IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i] + + # Allocate memory. + call salloc (buf, nout, TY_PIXEL) + call salloc (v1, IM_MAXDIM, TY_LONG) + call salloc (v2, IM_MAXDIM, TY_LONG) + call salloc (v3, IM_MAXDIM, TY_LONG) + + # Initialize the input line vector and the output section vectors. + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovkl (long(1), Meml[v2], IM_MAXDIM) + call amovkl (long(1), Meml[v3], IM_MAXDIM) + + # For each output line compute a block replicated line from the + # input image. If the line replication factor is greater than + # 1 then simply repeat the input line. This algorithm is + # sequential in both the input and output image though the + # input image will be recycled for each repetition of higher + # dimensions. + + while (impnl$t (out, buf2, Meml[v2]) != EOF) { + # Get the input vector corresponding to the output line. + do i = 2, ndim + Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1 + i = imgnl$t (in, buf1, Meml[v1]) + + # Block replicate the columns. + if (blkfac[1] == 1) + buf3 = buf1 + else { + ptrin = buf1 + ptrout = buf + do i = 1, nin { + do j = 1, blkfac[1] { + Mem$t[ptrout] = Mem$t[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + buf3 = buf + } + + # Copy the input line to the output line. + call amov$t (Mem$t[buf3], Mem$t[buf2], nout) + + # Repeat for each repetition of the input line. + for (i=2; i <= blkfac[2]; i=i+1) { + j = impnl$t (out, buf2, Meml[v2]) + call amov$t (Mem$t[buf3], Mem$t[buf2], nout) + } + + call amovl (Meml[v2], Meml[v3], IM_MAXDIM) + } + } + + call sfree (sp) +end +$endfor diff --git a/pkg/images/imgeom/src/generic/blkav.x b/pkg/images/imgeom/src/generic/blkav.x new file mode 100644 index 00000000..5a7df840 --- /dev/null +++ b/pkg/images/imgeom/src/generic/blkav.x @@ -0,0 +1,361 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <error.h> + +define AVG 1 # operation = arithmetic average +define SUM 2 # operation = arithmetic sum + +# change to (lrdx) in future + + +# BLKAVG -- Block average or sum on n-dimensional images. +# For example, block averaging by a factor of 2 means that pixels 1 and 2 +# are averaged to produce the first output pixel, 3 and 4 are averaged to +# produce the second output pixel, and so on. Remainder pixels at the end +# of a line, col, etc. are averaged correctly. + +procedure blkavl (im1, im2, blkfac, option) + +pointer im1 # input image +pointer im2 # output image +int blkfac[IM_MAXDIM] # blocking factors +int option # block operation (average, sum, ...) + +int num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx +int junk, num_ilines, num_olines, oline +long blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM] +long vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM] +real sum +pointer sp, accum_ptr, iline_ptr, oline_ptr + +int blkcomp(), imggsl(), impnll() +errchk imggsl(), impnll() + +begin + call smark (sp) + call salloc (accum_ptr, IM_LEN(im1, 1), TY_LONG) + + # Initialize; input and output vectors, block counters. + ndim = IM_NDIM(im1) + nfull_blkx = IM_LEN(im1, 1) / blkfac[1] + blkin_s[1] = long(1) + blkin_e[1] = long(IM_LEN(im1, 1)) + vin_s[1] = blkin_s[1] + vin_e[1] = blkin_e[1] + + do i = 1, ndim { + num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i] + IM_LEN(im2, i) = num_oblks[i] + } + num_olines = 1 + do i = 2, ndim + num_olines = num_olines * num_oblks[i] + call amovkl (long(1), vout, IM_MAXDIM) + + # For each sequential output-image line, ... + do oline = 1, num_olines { + + call aclrl (Meml[accum_ptr], IM_LEN(im1, 1)) + nlines_in_sum = 0 + + # Compute block vector; initialize dim>1 input vector. + num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e, + vin_s, vin_e) + + # Output pointer; note impnl$t returns vout for NEXT output line. + junk = impnll (im2, oline_ptr, vout) + + # For all input lines mapping to current output line, ... + do i = 1, num_ilines { + # Get line from input image. + iline_ptr = imggsl (im1, vin_s, vin_e, ndim) + + # Increment general section input vector between block bounds. + do dim = 2, ndim { + if (vin_s[dim] < blkin_e[dim]) { + vin_s[dim] = vin_s[dim] + 1 + vin_e[dim] = vin_s[dim] + break + } else { + vin_s[dim] = blkin_s[dim] + vin_e[dim] = vin_s[dim] + } + } + + # Accumulate line into block sum. Keep track of no. of + # lines in sum so that we can compute block average later. + + call aaddl (Meml[iline_ptr], Meml[accum_ptr], + Meml[accum_ptr], IM_LEN(im1,1)) + nlines_in_sum = nlines_in_sum + 1 + } + + # We now have a templine of sums; average/sum into output buffer + # first the full blocks using a vop. + if (option == AVG) + call abavl (Meml[accum_ptr], Meml[oline_ptr], nfull_blkx, + blkfac[1]) + else + call absul (Meml[accum_ptr], Meml[oline_ptr], nfull_blkx, + blkfac[1]) + + # Now average/sum the final partial block in x, if any. + if (nfull_blkx < num_oblks[1]) { + sum = 0 + count = 0 + do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) { + sum = sum + Meml[accum_ptr+i-1] + count = count + 1 + } + if (option == AVG) + Meml[oline_ptr+num_oblks[1]-1] = sum / count + else + Meml[oline_ptr+num_oblks[1]-1] = sum + } + + # Block average into output line from the sum of all lines block + # averaged in X. + if (option == AVG) + call adivkl (Meml[oline_ptr], long(nlines_in_sum), + Meml[oline_ptr], num_oblks[1]) + } + + call sfree (sp) +end + + + +# BLKAVG -- Block average or sum on n-dimensional images. +# For example, block averaging by a factor of 2 means that pixels 1 and 2 +# are averaged to produce the first output pixel, 3 and 4 are averaged to +# produce the second output pixel, and so on. Remainder pixels at the end +# of a line, col, etc. are averaged correctly. + +procedure blkavr (im1, im2, blkfac, option) + +pointer im1 # input image +pointer im2 # output image +int blkfac[IM_MAXDIM] # blocking factors +int option # block operation (average, sum, ...) + +int num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx +int junk, num_ilines, num_olines, oline +long blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM] +long vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM] +real sum +pointer sp, accum_ptr, iline_ptr, oline_ptr + +int blkcomp(), imggsr(), impnlr() +errchk imggsr(), impnlr() + +begin + call smark (sp) + call salloc (accum_ptr, IM_LEN(im1, 1), TY_REAL) + + # Initialize; input and output vectors, block counters. + ndim = IM_NDIM(im1) + nfull_blkx = IM_LEN(im1, 1) / blkfac[1] + blkin_s[1] = long(1) + blkin_e[1] = long(IM_LEN(im1, 1)) + vin_s[1] = blkin_s[1] + vin_e[1] = blkin_e[1] + + do i = 1, ndim { + num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i] + IM_LEN(im2, i) = num_oblks[i] + } + num_olines = 1 + do i = 2, ndim + num_olines = num_olines * num_oblks[i] + call amovkl (long(1), vout, IM_MAXDIM) + + # For each sequential output-image line, ... + do oline = 1, num_olines { + + call aclrr (Memr[accum_ptr], IM_LEN(im1, 1)) + nlines_in_sum = 0 + + # Compute block vector; initialize dim>1 input vector. + num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e, + vin_s, vin_e) + + # Output pointer; note impnl$t returns vout for NEXT output line. + junk = impnlr (im2, oline_ptr, vout) + + # For all input lines mapping to current output line, ... + do i = 1, num_ilines { + # Get line from input image. + iline_ptr = imggsr (im1, vin_s, vin_e, ndim) + + # Increment general section input vector between block bounds. + do dim = 2, ndim { + if (vin_s[dim] < blkin_e[dim]) { + vin_s[dim] = vin_s[dim] + 1 + vin_e[dim] = vin_s[dim] + break + } else { + vin_s[dim] = blkin_s[dim] + vin_e[dim] = vin_s[dim] + } + } + + # Accumulate line into block sum. Keep track of no. of + # lines in sum so that we can compute block average later. + + call aaddr (Memr[iline_ptr], Memr[accum_ptr], + Memr[accum_ptr], IM_LEN(im1,1)) + nlines_in_sum = nlines_in_sum + 1 + } + + # We now have a templine of sums; average/sum into output buffer + # first the full blocks using a vop. + if (option == AVG) + call abavr (Memr[accum_ptr], Memr[oline_ptr], nfull_blkx, + blkfac[1]) + else + call absur (Memr[accum_ptr], Memr[oline_ptr], nfull_blkx, + blkfac[1]) + + # Now average/sum the final partial block in x, if any. + if (nfull_blkx < num_oblks[1]) { + sum = 0.0 + count = 0 + do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) { + sum = sum + Memr[accum_ptr+i-1] + count = count + 1 + } + if (option == AVG) + Memr[oline_ptr+num_oblks[1]-1] = sum / count + else + Memr[oline_ptr+num_oblks[1]-1] = sum + } + + # Block average into output line from the sum of all lines block + # averaged in X. + if (option == AVG) + call adivkr (Memr[oline_ptr], real(nlines_in_sum), + Memr[oline_ptr], num_oblks[1]) + } + + call sfree (sp) +end + + + +# BLKAVG -- Block average or sum on n-dimensional images. +# For example, block averaging by a factor of 2 means that pixels 1 and 2 +# are averaged to produce the first output pixel, 3 and 4 are averaged to +# produce the second output pixel, and so on. Remainder pixels at the end +# of a line, col, etc. are averaged correctly. + +procedure blkavd (im1, im2, blkfac, option) + +pointer im1 # input image +pointer im2 # output image +int blkfac[IM_MAXDIM] # blocking factors +int option # block operation (average, sum, ...) + +int num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx +int junk, num_ilines, num_olines, oline +long blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM] +long vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM] +double sum +pointer sp, accum_ptr, iline_ptr, oline_ptr + +int blkcomp(), imggsd(), impnld() +errchk imggsd(), impnld() + +begin + call smark (sp) + call salloc (accum_ptr, IM_LEN(im1, 1), TY_DOUBLE) + + # Initialize; input and output vectors, block counters. + ndim = IM_NDIM(im1) + nfull_blkx = IM_LEN(im1, 1) / blkfac[1] + blkin_s[1] = long(1) + blkin_e[1] = long(IM_LEN(im1, 1)) + vin_s[1] = blkin_s[1] + vin_e[1] = blkin_e[1] + + do i = 1, ndim { + num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i] + IM_LEN(im2, i) = num_oblks[i] + } + num_olines = 1 + do i = 2, ndim + num_olines = num_olines * num_oblks[i] + call amovkl (long(1), vout, IM_MAXDIM) + + # For each sequential output-image line, ... + do oline = 1, num_olines { + + call aclrd (Memd[accum_ptr], IM_LEN(im1, 1)) + nlines_in_sum = 0 + + # Compute block vector; initialize dim>1 input vector. + num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e, + vin_s, vin_e) + + # Output pointer; note impnl$t returns vout for NEXT output line. + junk = impnld (im2, oline_ptr, vout) + + # For all input lines mapping to current output line, ... + do i = 1, num_ilines { + # Get line from input image. + iline_ptr = imggsd (im1, vin_s, vin_e, ndim) + + # Increment general section input vector between block bounds. + do dim = 2, ndim { + if (vin_s[dim] < blkin_e[dim]) { + vin_s[dim] = vin_s[dim] + 1 + vin_e[dim] = vin_s[dim] + break + } else { + vin_s[dim] = blkin_s[dim] + vin_e[dim] = vin_s[dim] + } + } + + # Accumulate line into block sum. Keep track of no. of + # lines in sum so that we can compute block average later. + + call aaddd (Memd[iline_ptr], Memd[accum_ptr], + Memd[accum_ptr], IM_LEN(im1,1)) + nlines_in_sum = nlines_in_sum + 1 + } + + # We now have a templine of sums; average/sum into output buffer + # first the full blocks using a vop. + if (option == AVG) + call abavd (Memd[accum_ptr], Memd[oline_ptr], nfull_blkx, + blkfac[1]) + else + call absud (Memd[accum_ptr], Memd[oline_ptr], nfull_blkx, + blkfac[1]) + + # Now average/sum the final partial block in x, if any. + if (nfull_blkx < num_oblks[1]) { + sum = 0.0D0 + count = 0 + do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) { + sum = sum + Memd[accum_ptr+i-1] + count = count + 1 + } + if (option == AVG) + Memd[oline_ptr+num_oblks[1]-1] = sum / count + else + Memd[oline_ptr+num_oblks[1]-1] = sum + } + + # Block average into output line from the sum of all lines block + # averaged in X. + if (option == AVG) + call adivkd (Memd[oline_ptr], double(nlines_in_sum), + Memd[oline_ptr], num_oblks[1]) + } + + call sfree (sp) +end + + diff --git a/pkg/images/imgeom/src/generic/blkrp.x b/pkg/images/imgeom/src/generic/blkrp.x new file mode 100644 index 00000000..bc43a3e5 --- /dev/null +++ b/pkg/images/imgeom/src/generic/blkrp.x @@ -0,0 +1,397 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + + + +# BLKRP -- Block replicate an image. + +procedure blkrpd (in, out, blkfac) + +pointer in # Input IMIO pointer +pointer out # Output IMIO pointer +int blkfac[ARB] # Block replication factors + +int i, j, ndim, nin, nout +pointer sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout +pointer imgl1d(), impl1d(), imgnld(), impnld() + +begin + call smark (sp) + + ndim = IM_NDIM(in) + nin = IM_LEN(in, 1) + nout = nin * blkfac[1] + IM_LEN(out,1) = nout + + if (ndim == 1) { + # For one dimensional images do the replication directly. + + buf1 = imgl1d (in) + buf2 = impl1d (out) + ptrin = buf1 + ptrout = buf2 + do i = 1, nin { + do j = 1, blkfac[1] { + Memd[ptrout] = Memd[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + + } else { + # For higher dimensional images use line access routines. + + do i = 2, ndim + IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i] + + # Allocate memory. + call salloc (buf, nout, TY_DOUBLE) + call salloc (v1, IM_MAXDIM, TY_LONG) + call salloc (v2, IM_MAXDIM, TY_LONG) + call salloc (v3, IM_MAXDIM, TY_LONG) + + # Initialize the input line vector and the output section vectors. + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovkl (long(1), Meml[v2], IM_MAXDIM) + call amovkl (long(1), Meml[v3], IM_MAXDIM) + + # For each output line compute a block replicated line from the + # input image. If the line replication factor is greater than + # 1 then simply repeat the input line. This algorithm is + # sequential in both the input and output image though the + # input image will be recycled for each repetition of higher + # dimensions. + + while (impnld (out, buf2, Meml[v2]) != EOF) { + # Get the input vector corresponding to the output line. + do i = 2, ndim + Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1 + i = imgnld (in, buf1, Meml[v1]) + + # Block replicate the columns. + if (blkfac[1] == 1) + buf3 = buf1 + else { + ptrin = buf1 + ptrout = buf + do i = 1, nin { + do j = 1, blkfac[1] { + Memd[ptrout] = Memd[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + buf3 = buf + } + + # Copy the input line to the output line. + call amovd (Memd[buf3], Memd[buf2], nout) + + # Repeat for each repetition of the input line. + for (i=2; i <= blkfac[2]; i=i+1) { + j = impnld (out, buf2, Meml[v2]) + call amovd (Memd[buf3], Memd[buf2], nout) + } + + call amovl (Meml[v2], Meml[v3], IM_MAXDIM) + } + } + + call sfree (sp) +end + + +# BLKRP -- Block replicate an image. + +procedure blkrpl (in, out, blkfac) + +pointer in # Input IMIO pointer +pointer out # Output IMIO pointer +int blkfac[ARB] # Block replication factors + +int i, j, ndim, nin, nout +pointer sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout +pointer imgl1l(), impl1l(), imgnll(), impnll() + +begin + call smark (sp) + + ndim = IM_NDIM(in) + nin = IM_LEN(in, 1) + nout = nin * blkfac[1] + IM_LEN(out,1) = nout + + if (ndim == 1) { + # For one dimensional images do the replication directly. + + buf1 = imgl1l (in) + buf2 = impl1l (out) + ptrin = buf1 + ptrout = buf2 + do i = 1, nin { + do j = 1, blkfac[1] { + Meml[ptrout] = Meml[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + + } else { + # For higher dimensional images use line access routines. + + do i = 2, ndim + IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i] + + # Allocate memory. + call salloc (buf, nout, TY_LONG) + call salloc (v1, IM_MAXDIM, TY_LONG) + call salloc (v2, IM_MAXDIM, TY_LONG) + call salloc (v3, IM_MAXDIM, TY_LONG) + + # Initialize the input line vector and the output section vectors. + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovkl (long(1), Meml[v2], IM_MAXDIM) + call amovkl (long(1), Meml[v3], IM_MAXDIM) + + # For each output line compute a block replicated line from the + # input image. If the line replication factor is greater than + # 1 then simply repeat the input line. This algorithm is + # sequential in both the input and output image though the + # input image will be recycled for each repetition of higher + # dimensions. + + while (impnll (out, buf2, Meml[v2]) != EOF) { + # Get the input vector corresponding to the output line. + do i = 2, ndim + Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1 + i = imgnll (in, buf1, Meml[v1]) + + # Block replicate the columns. + if (blkfac[1] == 1) + buf3 = buf1 + else { + ptrin = buf1 + ptrout = buf + do i = 1, nin { + do j = 1, blkfac[1] { + Meml[ptrout] = Meml[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + buf3 = buf + } + + # Copy the input line to the output line. + call amovl (Meml[buf3], Meml[buf2], nout) + + # Repeat for each repetition of the input line. + for (i=2; i <= blkfac[2]; i=i+1) { + j = impnll (out, buf2, Meml[v2]) + call amovl (Meml[buf3], Meml[buf2], nout) + } + + call amovl (Meml[v2], Meml[v3], IM_MAXDIM) + } + } + + call sfree (sp) +end + + +# BLKRP -- Block replicate an image. + +procedure blkrpr (in, out, blkfac) + +pointer in # Input IMIO pointer +pointer out # Output IMIO pointer +int blkfac[ARB] # Block replication factors + +int i, j, ndim, nin, nout +pointer sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout +pointer imgl1r(), impl1r(), imgnlr(), impnlr() + +begin + call smark (sp) + + ndim = IM_NDIM(in) + nin = IM_LEN(in, 1) + nout = nin * blkfac[1] + IM_LEN(out,1) = nout + + if (ndim == 1) { + # For one dimensional images do the replication directly. + + buf1 = imgl1r (in) + buf2 = impl1r (out) + ptrin = buf1 + ptrout = buf2 + do i = 1, nin { + do j = 1, blkfac[1] { + Memr[ptrout] = Memr[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + + } else { + # For higher dimensional images use line access routines. + + do i = 2, ndim + IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i] + + # Allocate memory. + call salloc (buf, nout, TY_REAL) + call salloc (v1, IM_MAXDIM, TY_LONG) + call salloc (v2, IM_MAXDIM, TY_LONG) + call salloc (v3, IM_MAXDIM, TY_LONG) + + # Initialize the input line vector and the output section vectors. + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovkl (long(1), Meml[v2], IM_MAXDIM) + call amovkl (long(1), Meml[v3], IM_MAXDIM) + + # For each output line compute a block replicated line from the + # input image. If the line replication factor is greater than + # 1 then simply repeat the input line. This algorithm is + # sequential in both the input and output image though the + # input image will be recycled for each repetition of higher + # dimensions. + + while (impnlr (out, buf2, Meml[v2]) != EOF) { + # Get the input vector corresponding to the output line. + do i = 2, ndim + Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1 + i = imgnlr (in, buf1, Meml[v1]) + + # Block replicate the columns. + if (blkfac[1] == 1) + buf3 = buf1 + else { + ptrin = buf1 + ptrout = buf + do i = 1, nin { + do j = 1, blkfac[1] { + Memr[ptrout] = Memr[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + buf3 = buf + } + + # Copy the input line to the output line. + call amovr (Memr[buf3], Memr[buf2], nout) + + # Repeat for each repetition of the input line. + for (i=2; i <= blkfac[2]; i=i+1) { + j = impnlr (out, buf2, Meml[v2]) + call amovr (Memr[buf3], Memr[buf2], nout) + } + + call amovl (Meml[v2], Meml[v3], IM_MAXDIM) + } + } + + call sfree (sp) +end + + +# BLKRP -- Block replicate an image. + +procedure blkrps (in, out, blkfac) + +pointer in # Input IMIO pointer +pointer out # Output IMIO pointer +int blkfac[ARB] # Block replication factors + +int i, j, ndim, nin, nout +pointer sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout +pointer imgl1s(), impl1s(), imgnls(), impnls() + +begin + call smark (sp) + + ndim = IM_NDIM(in) + nin = IM_LEN(in, 1) + nout = nin * blkfac[1] + IM_LEN(out,1) = nout + + if (ndim == 1) { + # For one dimensional images do the replication directly. + + buf1 = imgl1s (in) + buf2 = impl1s (out) + ptrin = buf1 + ptrout = buf2 + do i = 1, nin { + do j = 1, blkfac[1] { + Mems[ptrout] = Mems[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + + } else { + # For higher dimensional images use line access routines. + + do i = 2, ndim + IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i] + + # Allocate memory. + call salloc (buf, nout, TY_SHORT) + call salloc (v1, IM_MAXDIM, TY_LONG) + call salloc (v2, IM_MAXDIM, TY_LONG) + call salloc (v3, IM_MAXDIM, TY_LONG) + + # Initialize the input line vector and the output section vectors. + call amovkl (long(1), Meml[v1], IM_MAXDIM) + call amovkl (long(1), Meml[v2], IM_MAXDIM) + call amovkl (long(1), Meml[v3], IM_MAXDIM) + + # For each output line compute a block replicated line from the + # input image. If the line replication factor is greater than + # 1 then simply repeat the input line. This algorithm is + # sequential in both the input and output image though the + # input image will be recycled for each repetition of higher + # dimensions. + + while (impnls (out, buf2, Meml[v2]) != EOF) { + # Get the input vector corresponding to the output line. + do i = 2, ndim + Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1 + i = imgnls (in, buf1, Meml[v1]) + + # Block replicate the columns. + if (blkfac[1] == 1) + buf3 = buf1 + else { + ptrin = buf1 + ptrout = buf + do i = 1, nin { + do j = 1, blkfac[1] { + Mems[ptrout] = Mems[ptrin] + ptrout = ptrout + 1 + } + ptrin = ptrin + 1 + } + buf3 = buf + } + + # Copy the input line to the output line. + call amovs (Mems[buf3], Mems[buf2], nout) + + # Repeat for each repetition of the input line. + for (i=2; i <= blkfac[2]; i=i+1) { + j = impnls (out, buf2, Meml[v2]) + call amovs (Mems[buf3], Mems[buf2], nout) + } + + call amovl (Meml[v2], Meml[v3], IM_MAXDIM) + } + } + + call sfree (sp) +end + diff --git a/pkg/images/imgeom/src/generic/im3dtran.x b/pkg/images/imgeom/src/generic/im3dtran.x new file mode 100644 index 00000000..ae3153fe --- /dev/null +++ b/pkg/images/imgeom/src/generic/im3dtran.x @@ -0,0 +1,583 @@ + + +# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z. The arrays need not be +# identical. + +procedure txyz3s (a, b, nx, ny, nz) + +short a[nx, ny, nz], b[nx, ny, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, y, z] = a[x, y, z] +end + + +# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y. The arrays need not be +# identical. + +procedure txzy3s (a, b, nx, ny, nz) + +short a[nx, ny, nz], b[nx, nz, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, z, y] = a[x, y, z] +end + + +# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z. The arrays need not be +# identical. + +procedure tyxz3s (a, b, nx, ny, nz) + +short a[nx, ny, nz], b[ny, nx, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, x, z] = a[x, y, z] +end + + +# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x. The arrays need not be +# identical. + +procedure tyzx3s (a, b, nx, ny, nz) + +short a[nx, ny, nz], b[ny, nz, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, z, x] = a[x, y, z] +end + + +# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y. The arrays need not be +# identical. + +procedure tzxy3s (a, b, nx, ny, nz) + +short a[nx, ny, nz], b[nz, nx, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, x, y] = a[x, y, z] +end + + +# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x. The arrays need not be +# identical. + +procedure tzyx3s (a, b, nx, ny, nz) + +short a[nx, ny, nz], b[nz, ny, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, y, x] = a[x, y, z] +end + + + +# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z. The arrays need not be +# identical. + +procedure txyz3i (a, b, nx, ny, nz) + +int a[nx, ny, nz], b[nx, ny, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, y, z] = a[x, y, z] +end + + +# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y. The arrays need not be +# identical. + +procedure txzy3i (a, b, nx, ny, nz) + +int a[nx, ny, nz], b[nx, nz, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, z, y] = a[x, y, z] +end + + +# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z. The arrays need not be +# identical. + +procedure tyxz3i (a, b, nx, ny, nz) + +int a[nx, ny, nz], b[ny, nx, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, x, z] = a[x, y, z] +end + + +# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x. The arrays need not be +# identical. + +procedure tyzx3i (a, b, nx, ny, nz) + +int a[nx, ny, nz], b[ny, nz, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, z, x] = a[x, y, z] +end + + +# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y. The arrays need not be +# identical. + +procedure tzxy3i (a, b, nx, ny, nz) + +int a[nx, ny, nz], b[nz, nx, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, x, y] = a[x, y, z] +end + + +# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x. The arrays need not be +# identical. + +procedure tzyx3i (a, b, nx, ny, nz) + +int a[nx, ny, nz], b[nz, ny, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, y, x] = a[x, y, z] +end + + + +# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z. The arrays need not be +# identical. + +procedure txyz3l (a, b, nx, ny, nz) + +long a[nx, ny, nz], b[nx, ny, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, y, z] = a[x, y, z] +end + + +# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y. The arrays need not be +# identical. + +procedure txzy3l (a, b, nx, ny, nz) + +long a[nx, ny, nz], b[nx, nz, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, z, y] = a[x, y, z] +end + + +# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z. The arrays need not be +# identical. + +procedure tyxz3l (a, b, nx, ny, nz) + +long a[nx, ny, nz], b[ny, nx, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, x, z] = a[x, y, z] +end + + +# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x. The arrays need not be +# identical. + +procedure tyzx3l (a, b, nx, ny, nz) + +long a[nx, ny, nz], b[ny, nz, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, z, x] = a[x, y, z] +end + + +# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y. The arrays need not be +# identical. + +procedure tzxy3l (a, b, nx, ny, nz) + +long a[nx, ny, nz], b[nz, nx, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, x, y] = a[x, y, z] +end + + +# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x. The arrays need not be +# identical. + +procedure tzyx3l (a, b, nx, ny, nz) + +long a[nx, ny, nz], b[nz, ny, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, y, x] = a[x, y, z] +end + + + +# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z. The arrays need not be +# identical. + +procedure txyz3r (a, b, nx, ny, nz) + +real a[nx, ny, nz], b[nx, ny, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, y, z] = a[x, y, z] +end + + +# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y. The arrays need not be +# identical. + +procedure txzy3r (a, b, nx, ny, nz) + +real a[nx, ny, nz], b[nx, nz, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, z, y] = a[x, y, z] +end + + +# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z. The arrays need not be +# identical. + +procedure tyxz3r (a, b, nx, ny, nz) + +real a[nx, ny, nz], b[ny, nx, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, x, z] = a[x, y, z] +end + + +# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x. The arrays need not be +# identical. + +procedure tyzx3r (a, b, nx, ny, nz) + +real a[nx, ny, nz], b[ny, nz, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, z, x] = a[x, y, z] +end + + +# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y. The arrays need not be +# identical. + +procedure tzxy3r (a, b, nx, ny, nz) + +real a[nx, ny, nz], b[nz, nx, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, x, y] = a[x, y, z] +end + + +# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x. The arrays need not be +# identical. + +procedure tzyx3r (a, b, nx, ny, nz) + +real a[nx, ny, nz], b[nz, ny, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, y, x] = a[x, y, z] +end + + + +# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z. The arrays need not be +# identical. + +procedure txyz3d (a, b, nx, ny, nz) + +double a[nx, ny, nz], b[nx, ny, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, y, z] = a[x, y, z] +end + + +# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y. The arrays need not be +# identical. + +procedure txzy3d (a, b, nx, ny, nz) + +double a[nx, ny, nz], b[nx, nz, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, z, y] = a[x, y, z] +end + + +# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z. The arrays need not be +# identical. + +procedure tyxz3d (a, b, nx, ny, nz) + +double a[nx, ny, nz], b[ny, nx, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, x, z] = a[x, y, z] +end + + +# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x. The arrays need not be +# identical. + +procedure tyzx3d (a, b, nx, ny, nz) + +double a[nx, ny, nz], b[ny, nz, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, z, x] = a[x, y, z] +end + + +# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y. The arrays need not be +# identical. + +procedure tzxy3d (a, b, nx, ny, nz) + +double a[nx, ny, nz], b[nz, nx, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, x, y] = a[x, y, z] +end + + +# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x. The arrays need not be +# identical. + +procedure tzyx3d (a, b, nx, ny, nz) + +double a[nx, ny, nz], b[nz, ny, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, y, x] = a[x, y, z] +end + + + +# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z. The arrays need not be +# identical. + +procedure txyz3x (a, b, nx, ny, nz) + +complex a[nx, ny, nz], b[nx, ny, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, y, z] = a[x, y, z] +end + + +# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y. The arrays need not be +# identical. + +procedure txzy3x (a, b, nx, ny, nz) + +complex a[nx, ny, nz], b[nx, nz, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, z, y] = a[x, y, z] +end + + +# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z. The arrays need not be +# identical. + +procedure tyxz3x (a, b, nx, ny, nz) + +complex a[nx, ny, nz], b[ny, nx, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, x, z] = a[x, y, z] +end + + +# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x. The arrays need not be +# identical. + +procedure tyzx3x (a, b, nx, ny, nz) + +complex a[nx, ny, nz], b[ny, nz, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, z, x] = a[x, y, z] +end + + +# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y. The arrays need not be +# identical. + +procedure tzxy3x (a, b, nx, ny, nz) + +complex a[nx, ny, nz], b[nz, nx, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, x, y] = a[x, y, z] +end + + +# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x. The arrays need not be +# identical. + +procedure tzyx3x (a, b, nx, ny, nz) + +complex a[nx, ny, nz], b[nz, ny, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, y, x] = a[x, y, z] +end + + diff --git a/pkg/images/imgeom/src/generic/imtrans.x b/pkg/images/imgeom/src/generic/imtrans.x new file mode 100644 index 00000000..26754fe6 --- /dev/null +++ b/pkg/images/imgeom/src/generic/imtrans.x @@ -0,0 +1,93 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + + + +# IMTR2 -- Generic transpose. The arrays need not be identical. + +procedure imtr2s (a, b, nx, ny) + +short a[nx, ny], b[ny, nx] +int nx, ny, x, y + +begin + do x = 1, nx + do y = 1, ny + b[y, x] = a[x, y] +end + + + +# IMTR2 -- Generic transpose. The arrays need not be identical. + +procedure imtr2i (a, b, nx, ny) + +int a[nx, ny], b[ny, nx] +int nx, ny, x, y + +begin + do x = 1, nx + do y = 1, ny + b[y, x] = a[x, y] +end + + + +# IMTR2 -- Generic transpose. The arrays need not be identical. + +procedure imtr2l (a, b, nx, ny) + +long a[nx, ny], b[ny, nx] +int nx, ny, x, y + +begin + do x = 1, nx + do y = 1, ny + b[y, x] = a[x, y] +end + + + +# IMTR2 -- Generic transpose. The arrays need not be identical. + +procedure imtr2r (a, b, nx, ny) + +real a[nx, ny], b[ny, nx] +int nx, ny, x, y + +begin + do x = 1, nx + do y = 1, ny + b[y, x] = a[x, y] +end + + + +# IMTR2 -- Generic transpose. The arrays need not be identical. + +procedure imtr2d (a, b, nx, ny) + +double a[nx, ny], b[ny, nx] +int nx, ny, x, y + +begin + do x = 1, nx + do y = 1, ny + b[y, x] = a[x, y] +end + + + +# IMTR2 -- Generic transpose. The arrays need not be identical. + +procedure imtr2x (a, b, nx, ny) + +complex a[nx, ny], b[ny, nx] +int nx, ny, x, y + +begin + do x = 1, nx + do y = 1, ny + b[y, x] = a[x, y] +end + + diff --git a/pkg/images/imgeom/src/generic/mkpkg b/pkg/images/imgeom/src/generic/mkpkg new file mode 100644 index 00000000..9fe8f222 --- /dev/null +++ b/pkg/images/imgeom/src/generic/mkpkg @@ -0,0 +1,13 @@ +# Library for the GEOMETRIC TRANSFORMATION tasks + +$checkout libpkg.a pkg$images/ +$update libpkg.a +$checkin libpkg.a pkg$images/ +$exit + +libpkg.a: + blkav.x <imhdr.h> <error.h> + blkrp.x <imhdr.h> + imtrans.x + im3dtran.x + ; diff --git a/pkg/images/imgeom/src/im3dtran.gx b/pkg/images/imgeom/src/im3dtran.gx new file mode 100644 index 00000000..1b683124 --- /dev/null +++ b/pkg/images/imgeom/src/im3dtran.gx @@ -0,0 +1,98 @@ +$for (silrdx) + +# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z. The arrays need not be +# identical. + +procedure txyz3$t (a, b, nx, ny, nz) + +PIXEL a[nx, ny, nz], b[nx, ny, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, y, z] = a[x, y, z] +end + + +# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y. The arrays need not be +# identical. + +procedure txzy3$t (a, b, nx, ny, nz) + +PIXEL a[nx, ny, nz], b[nx, nz, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[x, z, y] = a[x, y, z] +end + + +# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z. The arrays need not be +# identical. + +procedure tyxz3$t (a, b, nx, ny, nz) + +PIXEL a[nx, ny, nz], b[ny, nx, nz] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, x, z] = a[x, y, z] +end + + +# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x. The arrays need not be +# identical. + +procedure tyzx3$t (a, b, nx, ny, nz) + +PIXEL a[nx, ny, nz], b[ny, nz, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[y, z, x] = a[x, y, z] +end + + +# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y. The arrays need not be +# identical. + +procedure tzxy3$t (a, b, nx, ny, nz) + +PIXEL a[nx, ny, nz], b[nz, nx, ny] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, x, y] = a[x, y, z] +end + + +# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x. The arrays need not be +# identical. + +procedure tzyx3$t (a, b, nx, ny, nz) + +PIXEL a[nx, ny, nz], b[nz, ny, nx] +int nx, ny, nz, x, y, z + +begin + do x = 1, nx + do y = 1, ny + do z = 1, nz + b[z, y, x] = a[x, y, z] +end + +$endfor diff --git a/pkg/images/imgeom/src/imtrans.gx b/pkg/images/imgeom/src/imtrans.gx new file mode 100644 index 00000000..3749e314 --- /dev/null +++ b/pkg/images/imgeom/src/imtrans.gx @@ -0,0 +1,18 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +$for (silrdx) + +# IMTR2 -- Generic transpose. The arrays need not be identical. + +procedure imtr2$t (a, b, nx, ny) + +PIXEL a[nx, ny], b[ny, nx] +int nx, ny, x, y + +begin + do x = 1, nx + do y = 1, ny + b[y, x] = a[x, y] +end + +$endfor diff --git a/pkg/images/imgeom/src/mkpkg b/pkg/images/imgeom/src/mkpkg new file mode 100644 index 00000000..2e784d54 --- /dev/null +++ b/pkg/images/imgeom/src/mkpkg @@ -0,0 +1,35 @@ +# Library for the GEOMETRIC TRANSFORMATION tasks + +$checkout libpkg.a ../../ +$update libpkg.a +$checkin libpkg.a ../../ +$exit + +generic: + $set GEN = "$$generic -k" + + $ifolder (generic/blkav.x, blkav.gx) + $(GEN) blkav.gx -o generic/blkav.x $endif + $ifolder (generic/blkrp.x, blkrp.gx) + $(GEN) blkrp.gx -o generic/blkrp.x $endif + $ifolder (generic/imtrans.x, imtrans.gx) + $(GEN) imtrans.gx -o generic/imtrans.x $endif + $ifolder (generic/im3dtran.x, im3dtran.gx) + $(GEN) im3dtran.gx -o generic/im3dtran.x $endif + ; + +libpkg.a: + $ifeq (USE_GENERIC, yes) $call generic $endif + + @generic + + blkcomp.x <imhdr.h> + shiftlines.x <imhdr.h> <imset.h> <math/iminterp.h> + t_blkavg.x <imhdr.h> + t_blkrep.x <imhdr.h> + t_imshift.x <error.h> <imhdr.h> <imset.h> <math/iminterp.h> + t_imtrans.x <imhdr.h> <error.h> <mwset.h> + t_im3dtran.x <imhdr.h> <error.h> <mwset.h> + t_magnify.x <imhdr.h> <imset.h> <error.h> <math/iminterp.h> + t_shiftlines.x <error.h> + ; diff --git a/pkg/images/imgeom/src/shiftlines.x b/pkg/images/imgeom/src/shiftlines.x new file mode 100644 index 00000000..e0bd6d9a --- /dev/null +++ b/pkg/images/imgeom/src/shiftlines.x @@ -0,0 +1,279 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <imset.h> +include <math/iminterp.h> + +define NMARGIN 3 + +# SH_LINES -- Shift image lines. +# +# If an integer shift is given then do an efficient integer shift +# without datatype I/O conversions and using array move operators. +# If the shift is non-integer then use image interpolation. + +procedure sh_lines (im1, im2, shift, boundary, constant, interpstr) + +pointer im1 # Input image descriptor +pointer im2 # Output image descriptor +real shift # Shift in pixels +int boundary # Boundary extension type +real constant # Constant boundary extension +char interpstr[ARB] # Interpolation type + +int i, nsinc, nincr, ncols, nimcols, nlines, nbpix, nmargin, interpolation +long v1[IM_MAXDIM], v2[IM_MAXDIM], vout[IM_MAXDIM] +real dx, deltax, cx +pointer sp, x, asi, junk, buf1, buf2 + +bool fp_equalr() +int imggsr(), impnlr(), asigeti() + +begin + # Check for out of bounds shifts. + ncols = IM_LEN(im1, 1) + if (shift < -ncols || shift > ncols) + call error (0, "SHIFTLINES: Shift out of bounds") + + # Compute the shift. + dx = abs (shift - int (shift)) + if (fp_equalr (dx, 0.0)) + deltax = 0.0 + else if (shift > 0.0) + deltax = 1.0 - dx + else + deltax = dx + + # Initialize the interpolation. + call asitype (interpstr, interpolation, nsinc, nincr, cx) + if (interpolation == II_LSINC || interpolation == II_SINC) + call asisinit (asi, II_LSINC, nsinc, 1, deltax - nint (deltax), + 0.0) + else + call asisinit (asi, interpolation, nsinc, 1, cx, 0.0) + + # Set up the image boundary conditions. + if (interpolation == II_SINC || interpolation == II_LSINC) + nmargin = asigeti (asi, II_ASINSINC) + else + nmargin = NMARGIN + nbpix = int (abs (shift) + 1.0) + nmargin + call imseti (im1, IM_TYBNDRY, boundary) + call imseti (im1, IM_NBNDRYPIX, nbpix) + call imsetr (im1, IM_BNDRYPIXVAL, constant) + + # Allocate space for and set up the interpolation coordinates. + call smark (sp) + call salloc (x, 2 * ncols, TY_REAL) + deltax = deltax + nmargin + if (interpolation == II_DRIZZLE) { + do i = 1, ncols { + Memr[x+2*i-2] = i + deltax - 0.5 + Memr[x+2*i-1] = i + deltax + 0.5 + } + } else { + do i = 1, ncols + Memr[x+i-1] = i + deltax + } + + # Initialize the input v vectors. + cx = 1. - nmargin - shift + if ((cx <= 0.0) && (! fp_equalr (dx, 0.0))) + v1[1] = long (cx) - 1 + else + v1[1] = long (cx) + v2[1] = ncols - shift + nmargin + 1 + nimcols = v2[1] - v1[1] + 1 + do i = 2, IM_NDIM(im1) { + v1[i] = long (1) + v2[i] = long (1) + } + + # Compute the number of output lines. + nlines = 1 + do i = 2, IM_NDIM(im1) + nlines = nlines * IM_LEN(im1, i) + + # Initialize the output v vector. + call amovkl (long(1), vout, IM_MAXDIM) + + # Shift the images. + do i = 1, nlines { + + # Get the input image data buffer. + buf1 = imggsr (im1, v1, v2, IM_NDIM(im1)) + if (buf1 == EOF) + call error (0, "SHIFTLINES: Error reading input image\n") + + # Get the output image data buffer. + junk = impnlr (im2, buf2, vout) + if (junk == EOF) + call error (0, "SHIFTLINES: Error writing output image\n") + + # Evaluate the interpolation at the shifted points. + call asifit (asi, Memr[buf1], nimcols) + call asivector (asi, Memr[x], Memr[buf2], ncols) + + # Increment the v vectors. + call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1)) + } + + call asifree (asi) + call sfree (sp) +end + + +# SH_LINESI -- Integer pixel shift. +# +# Shift the pixels in an image by an integer amount. Perform I/O without +# out type conversion and use array move operators. + +procedure sh_linesi (im1, im2, shift, boundary, constant) + +pointer im1 # Input image descriptor +pointer im2 # Output image descriptor +int shift # Integer shift +int boundary # Boundary extension type +real constant # Constant for boundary extension + +int i, ncols, nlines, junk +long v1[IM_MAXDIM], v2[IM_MAXDIM], vout[IM_MAXDIM] +pointer buf1, buf2 + +int imggss(), imggsi(), imggsl(), imggsr(), imggsd(), imggsx() +int impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx() + +begin + # Set the boundary extension parameters. + call imseti (im1, IM_NBNDRYPIX, abs (shift)) + call imseti (im1, IM_TYBNDRY, boundary) + call imsetr (im1, IM_BNDRYPIXVAL, constant) + + # Return if shift off image. + ncols = IM_LEN(im1, 1) + if (shift < -ncols || shift > ncols) + call error (0, "Shiftlinesi: shift out of bounds") + + # Setup start vector for sequential reads and writes. + v1[1] = max (-ncols + 1, -shift + 1) + v2[1] = min (2 * ncols, ncols - shift) + do i = 2, IM_NDIM(im1) { + v1[i] = long (1) + v2[i] = long (1) + } + call amovkl (long(1), vout, IM_MAXDIM) + + # Setup line counter. + nlines = 1 + do i = 2, IM_NDIM(im1) + nlines = nlines * IM_LEN(im1, i) + + + # Shift the image using appropriate datatype operators. + switch (IM_PIXTYPE(im1)) { + + case TY_SHORT: + do i = 1, nlines { + buf1 = imggss (im1, v1, v2, IM_NDIM(im1)) + if (buf1 == EOF) + call error (0, "Shiftlines: error reading input image") + junk = impnls (im2, buf2, vout) + if (junk == EOF) + call error (0, "Shiftlinesi: error writing out image") + call amovs (Mems[buf1], Mems[buf2], ncols) + call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1)) + } + + case TY_INT: + do i = 1, nlines { + buf1 = imggsi (im1, v1, v2, IM_NDIM(im1)) + if (buf1 == EOF) + call error (0, "Shiftlines: error reading input image") + junk = impnli (im2, buf2, vout) + if (junk == EOF) + call error (0, "Shiftlinesi: error writing out image") + call amovi (Memi[buf1], Memi[buf2], ncols) + call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1)) + } + + case TY_LONG: + do i = 1, nlines { + buf1 = imggsl (im1, v1, v2, IM_NDIM(im1)) + if (buf1 == EOF) + call error (0, "Shiftlines: error reading input image") + junk = impnll (im2, buf2, vout) + if (junk == EOF) + call error (0, "Shiftlinesi: error writing out image") + call amovl (Meml[buf1], Meml[buf2], ncols) + call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1)) + } + + case TY_REAL: + do i = 1, nlines { + buf1 = imggsr (im1, v1, v2, IM_NDIM(im1)) + if (buf1 == EOF) + call error (0, "Shiftlines: error reading input image") + junk = impnlr (im2, buf2, vout) + if (junk == EOF) + call error (0, "Shiftlinesi: error writing out image") + call amovr (Memr[buf1], Memr[buf2], ncols) + call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1)) + } + + case TY_DOUBLE: + do i = 1, nlines { + buf1 = imggsd (im1, v1, v2, IM_NDIM(im1)) + if (buf1 == EOF) + call error (0, "Shiftlines: error reading input image") + junk = impnld (im2, buf2, vout) + if (junk == EOF) + call error (0, "Shiftlinesi: error writing out image") + call amovd (Memd[buf1], Memd[buf2], ncols) + call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1)) + } + + case TY_COMPLEX: + do i = 1, nlines { + buf1 = imggsx (im1, v1, v2, IM_NDIM(im1)) + if (buf1 == EOF) + call error (0, "Shiftlines: error reading input image") + junk = impnlx (im2, buf2, vout) + if (junk == EOF) + call error (0, "Shiftlinesi: error writing out image") + call amovx (Memx[buf1], Memx[buf2], ncols) + call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1)) + } + + default: + call error (0, "Unknown pixel datatype") + } +end + + +# SH_LOOP -- Increment the vector V from VS to VE (nested do loops cannot +# be used because of the variable number of dimensions). Return LOOP_DONE +# when V exceeds VE. + +procedure sh_loop (vs, ve, szdims, ndim) + +long vs[ndim] # the start vector +long ve[ndim] # the end vector +long szdims[ndim] # the dimensions vector +int ndim # the number of dimensions + +int dim + +begin + for (dim=2; dim <= ndim; dim=dim+1) { + vs[dim] = vs[dim] + 1 + ve[dim] = vs[dim] + if (vs[dim] - szdims[dim] == 1) { + if (dim < ndim) { + vs[dim] = 1 + ve[dim] = 1 + } else + break + } else + break + } +end diff --git a/pkg/images/imgeom/src/t_blkavg.x b/pkg/images/imgeom/src/t_blkavg.x new file mode 100644 index 00000000..e621bc64 --- /dev/null +++ b/pkg/images/imgeom/src/t_blkavg.x @@ -0,0 +1,115 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +define OPTIONS "|average|sum|" # Values of task param options +define AVG 1 # Arithmetic average in block +define SUM 2 # Sum of pixels within block +define DEF_BLKFAC 1 # Default blocking factor + +# T_BLKAVG -- Block average or sum on n-dimensional images. +# +# The input and output images are given by image template lists. The +# number of output images must match the number of input images. Image +# sections are allowed in the input images and are ignored in the output +# images. If the input and output image names are the same then the +# blocking operation is performed to a temporary file which then replaces +# the input image. + +procedure t_blkavg() + +char imtlist1[SZ_LINE] # Input image list +char imtlist2[SZ_LINE] # Output image list +int option # Type of operation +int blkfac[IM_MAXDIM] # Block sizes + +char image1[SZ_FNAME] # Input image name +char image2[SZ_FNAME] # Output image name +char imtemp[SZ_FNAME] # Temporary file + +int list1, list2, i +pointer im1, im2, mw +real shifts[IM_MAXDIM], mags[IM_MAXDIM] + +bool envgetb() +int imtopen(), imtgetim(), imtlen(), clgeti(), clgwrd() +pointer immap(), mw_openim() + +string blk_param "bX" + +begin + # Get input and output image template lists and the block sizes. + + call clgstr ("input", imtlist1, SZ_LINE) + call clgstr ("output", imtlist2, SZ_LINE) + option = clgwrd ("option", image1, SZ_FNAME, OPTIONS) + call amovki (INDEFI, blkfac, IM_MAXDIM) + + # Expand the input and output image lists. + + list1 = imtopen (imtlist1) + list2 = imtopen (imtlist2) + + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + call error (0, "Number of input and output images not the same") + } + + # Do each set of input/output images. + + while ((imtgetim (list1, image1, SZ_FNAME) != EOF) && + (imtgetim (list2, image2, SZ_FNAME) != EOF)) { + + call xt_mkimtemp (image1, image2, imtemp, SZ_FNAME) + + im1 = immap (image1, READ_ONLY, 0) + im2 = immap (image2, NEW_COPY, im1) + + do i = 1, IM_NDIM(im1) { + if (IS_INDEFI(blkfac[i])) { + call sprintf (blk_param[2], SZ_CHAR, "%1d") + call pargi (i) + blkfac[i] = max (1, min (clgeti (blk_param), + IM_LEN(im1, i))) + } + } + + # Perform the block operation. + switch (IM_PIXTYPE (im1)) { + case TY_SHORT, TY_USHORT, TY_INT, TY_LONG: + call blkavl (im1, im2, blkfac, option) + case TY_REAL: + call blkavr (im1, im2, blkfac, option) + case TY_DOUBLE: + call blkavd (im1, im2, blkfac, option) + case TY_COMPLEX: + #call blkavx (im1, im2, blkfac, option) + call error (0, + "Blkavg does not currently support pixel data type complex.") + default: + call error (0, "Unknown pixel data type") + } + + # Update the world coordinate system. + if (!envgetb ("nomwcs")) { + mw = mw_openim (im1) + do i = 1, IM_NDIM(im1) + mags[i] = 1.0d0 / double (blkfac[i]) + call mw_scale (mw, mags, (2 ** IM_NDIM(im1) - 1)) + do i = 1, IM_NDIM(im1) + shifts[i] = 0.5d0 - 1.0d0 / double (blkfac[i]) / 2.0d0 + call mw_shift (mw, shifts, (2 ** IM_NDIM(im1) - 1)) + call mw_saveim (mw, im2) + call mw_close (mw) + } + + call imunmap (im2) + call imunmap (im1) + + call xt_delimtemp (image2, imtemp) + } + + call imtclose (list1) + call imtclose (list2) +end diff --git a/pkg/images/imgeom/src/t_blkrep.x b/pkg/images/imgeom/src/t_blkrep.x new file mode 100644 index 00000000..2f92a567 --- /dev/null +++ b/pkg/images/imgeom/src/t_blkrep.x @@ -0,0 +1,96 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# T_BLKREP -- Block replicate n-dimensional images. +# +# The input and output images are given by image template lists. The +# number of output images must match the number of input images. Image +# sections are allowed in the input images and are ignored in the output +# images. If the input and output image names are the same then the +# replication operation is performed to a temporary file which then replaces +# the input image. + +procedure t_blkrep() + +int i, list1, list2 +pointer sp, image1, image2, image3, blkfac, im1, im2, mw +real shifts[IM_MAXDIM], mags[IM_MAXDIM] + +bool envgetb() +int imtopenp(), imtgetim(), imtlen(), clgeti() +pointer immap(), mw_openim() +string blk_param "bX" + +begin + # Allocate memory. + call smark (sp) + call salloc (image1, SZ_LINE, TY_CHAR) + call salloc (image2, SZ_LINE, TY_CHAR) + call salloc (image3, SZ_LINE, TY_CHAR) + call salloc (blkfac, IM_MAXDIM, TY_INT) + + # Expand the input and output image lists. + + list1 = imtopenp ("input") + list2 = imtopenp ("output") + + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + call error (1, "Number of input and output images not the same") + } + + # Do each set of input/output images. + + call amovki (INDEFI, Memi[blkfac], IM_MAXDIM) + while ((imtgetim (list1, Memc[image1], SZ_LINE) != EOF) && + (imtgetim (list2, Memc[image2], SZ_LINE) != EOF)) { + + call xt_mkimtemp (Memc[image1], Memc[image2], Memc[image3], SZ_LINE) + + im1 = immap (Memc[image1], READ_ONLY, 0) + im2 = immap (Memc[image2], NEW_COPY, im1) + + do i = 1, IM_NDIM(im1) { + if (IS_INDEFI(Memi[blkfac+i-1])) { + call sprintf (blk_param[2], SZ_CHAR, "%1d") + call pargi (i) + Memi[blkfac+i-1] = clgeti (blk_param) + } + } + + # Perform the block operation. + switch (IM_PIXTYPE (im1)) { + case TY_SHORT: + call blkrps (im1, im2, Memi[blkfac]) + case TY_INT, TY_LONG: + call blkrpl (im1, im2, Memi[blkfac]) + case TY_DOUBLE: + call blkrpd (im1, im2, Memi[blkfac]) + default: + call blkrpr (im1, im2, Memi[blkfac]) + } + + if (!envgetb ("nomwcs")) { + mw = mw_openim (im1) + do i = 1, IM_NDIM(im1) + mags[i] = double (Memi[blkfac+i-1]) + call mw_scale (mw, mags, (2 ** IM_NDIM(im1) - 1)) + do i = 1, IM_NDIM(im1) + shifts[i] = 0.5d0 - double (Memi[blkfac+i-1]) / 2.0d0 + call mw_shift (mw, shifts, (2 ** IM_NDIM(im1) - 1)) + call mw_saveim (mw, im2) + call mw_close (mw) + } + + call imunmap (im2) + call imunmap (im1) + + call xt_delimtemp (Memc[image2], Memc[image3]) + } + + call imtclose (list1) + call imtclose (list2) + call sfree (sp) +end diff --git a/pkg/images/imgeom/src/t_im3dtran.x b/pkg/images/imgeom/src/t_im3dtran.x new file mode 100644 index 00000000..46d4a536 --- /dev/null +++ b/pkg/images/imgeom/src/t_im3dtran.x @@ -0,0 +1,719 @@ +include <imhdr.h> +include <error.h> +include <mwset.h> + + +# Define all possible tranpose operations. +define XYZ 1 # xyz -> xyz +define XZY 2 # xyz -> xzy +define YXZ 3 # xyz -> yxz +define YZX 4 # xyz -> yzx +define ZXY 5 # xyz -> zxy +define ZYX 6 # xyz -> zyx + + +# T_IM3DTRAN -- Transpose 3d images. +# +# The input and output images are given by image template lists. The +# number of output images must match the number of input images. Image +# sections are allowed in the input images and are ignored in the output +# images. If the input and output image names are the same then the transpose +# is performed to a temporary file which then replaces the input image. + +procedure t_im3dtran () + +bool verbose +int list1, list2, len_blk, new_ax[3], which3d +pointer sp, imtlist1, imtlist2, image1, image2, imtemp, im1, im2, mw + +bool clgetb(), envgetb() +int clgeti(), imtopen(), imtgetim(), imtlen(), whichtran() +pointer immap(), mw_openim() +errchk im3dtranpose(), mw_openim(), mw_saveim(), mw_close(), im3dtrmw() + +begin + # Get some working space. + call smark (sp) + call salloc (imtlist1, SZ_LINE, TY_CHAR) + call salloc (imtlist2, SZ_LINE, TY_CHAR) + call salloc (image1, SZ_FNAME, TY_CHAR) + call salloc (image2, SZ_FNAME, TY_CHAR) + call salloc (imtemp, SZ_FNAME, TY_CHAR) + + # Get input and output image template lists, the size of the transpose + # block, and the transpose mapping. + call clgstr ("input", Memc[imtlist1], SZ_LINE) + call clgstr ("output", Memc[imtlist2], SZ_LINE) + new_ax[1] = clgeti ("new_x") + new_ax[2] = clgeti ("new_y") + new_ax[3] = clgeti ("new_z") + len_blk = clgeti ("len_blk") + verbose = clgetb ("verbose") + + # Determine the type of 3d transpose. + which3d = whichtran (new_ax) + if (which3d <= 0) + call error (0, "Invalid mapping of new_x, new_y, or new_z") + + # Expand the input and output image lists. + + list1 = imtopen (Memc[imtlist1]) + list2 = imtopen (Memc[imtlist2]) + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + call error (1, "Number of input and output images not the same") + } + + # Do each set of input/output images. + while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF) && + (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF)) { + + + if (verbose) { + call printf ( + "Image: %s axes: [123] -> Image: %s axes: [%d%d%d]\n") + call pargstr (Memc[image1]) + call pargstr (Memc[image2]) + call pargi (new_ax[1]) + call pargi (new_ax[2]) + call pargi (new_ax[3]) + call flush (STDOUT) + } + + call xt_mkimtemp (Memc[image1], Memc[image2], Memc[imtemp], + SZ_FNAME) + im1 = immap (Memc[image1], READ_ONLY, 0) + im2 = immap (Memc[image2], NEW_COPY, im1) + + iferr { + + # Do the transpose. + call im3dtranspose (im1, im2, len_blk, which3d, new_ax) + + # Update the image WCS to reflect the transpose. + if (!envgetb ("nomwcs")) { + mw = mw_openim (im1) + call im3dtrmw (mw, which3d) + call mw_saveim (mw, im2) + call mw_close (mw) + } + + } then { + + call eprintf ("Error transposing image: %s\n") + call pargstr (Memc[image1]) + call erract (EA_WARN) + call imunmap (im2) + call imunmap (im1) + call imdelete (Memc[image2]) + + } else { + + # Finish up + call imunmap (im2) + call imunmap (im1) + call xt_delimtemp (Memc[image2], Memc[imtemp]) + } + } + + call imtclose (list1) + call imtclose (list2) + + call sfree (sp) +end + + +# IM3DTRANSPOSE -- Transpose a 3D image. +# +# Divide the image into square blocks of size len_blk by len_blk. +# Transpose each block with a generic array transpose operator. + +procedure im3dtranspose (im_in, im_out, len_blk, which3d, new_ax) + +pointer im_in #I Input image descriptor +pointer im_out #I Output image descriptor +int len_blk #I 1D length of transpose block +int which3d #I Parameterized transpose order +int new_ax[3] #I Map old axis[index] to new value + +int x1, x2, nx, y1, y2, ny, z1, z2, nz +pointer buf_in, buf_out +pointer imgs3s(), imps3s(), imgs3i(), imps3i(), imgs3l(), imps3l() +pointer imgs3r(), imps3r(), imgs3d(), imps3d(), imgs3x(), imps3x() + +begin + # Check that the image is 3D. + if (IM_NDIM(im_in) != 3) + call error (1, "image is not 3D.") + + # Output image is a copy of input image with dimensions transposed. + + IM_LEN (im_out, 1) = IM_LEN (im_in, new_ax[1]) + IM_LEN (im_out, 2) = IM_LEN (im_in, new_ax[2]) + IM_LEN (im_out, 3) = IM_LEN (im_in, new_ax[3]) + + # Break the input image into blocks of at most (len_blk) ** 3. + + do x1 = 1, IM_LEN (im_in, 1), len_blk { + x2 = x1 + len_blk - 1 + if (x2 > IM_LEN(im_in, 1)) + x2 = IM_LEN(im_in, 1) + nx = x2 - x1 + 1 + + do y1 = 1, IM_LEN (im_in, 2), len_blk { + y2 = y1 + len_blk - 1 + if (y2 > IM_LEN(im_in, 2)) + y2 = IM_LEN(im_in, 2) + ny = y2 - y1 + 1 + + do z1 = 1, IM_LEN (im_in, 3), len_blk { + z2 = z1 + len_blk - 1 + if (z2 > IM_LEN(im_in, 3)) + z2 = IM_LEN(im_in, 3) + nz = z2 - z1 + 1 + + # Switch on the pixel type to optimize IMIO. + + switch (IM_PIXTYPE (im_in)) { + + case TY_SHORT: + buf_in = imgs3s (im_in, x1, x2, y1, y2, z1, z2) + switch (which3d) { + case XYZ: + buf_out = imps3s (im_out, x1, x2, y1, y2, z1, z2) + call txyz3s (Mems[buf_in], Mems[buf_out], nx,ny,nz) + case XZY: + buf_out = imps3s (im_out, x1, x2, z1, z2, y1, y2) + call txzy3s (Mems[buf_in], Mems[buf_out], nx,ny,nz) + case YXZ: + buf_out = imps3s (im_out, y1, y2, x1, x2, z1, z2) + call tyxz3s (Mems[buf_in], Mems[buf_out], nx,ny,nz) + case YZX: + buf_out = imps3s (im_out, y1, y2, z1, z2, x1, x2) + call tyzx3s (Mems[buf_in], Mems[buf_out], nx,ny,nz) + case ZXY: + buf_out = imps3s (im_out, z1, z2, x1, x2, y1, y2) + call tzxy3s (Mems[buf_in], Mems[buf_out], nx,ny,nz) + case ZYX: + buf_out = imps3s (im_out, z1, z2, y1, y2, x1, x2) + call tzyx3s (Mems[buf_in], Mems[buf_out], nx,ny,nz) + } + + case TY_INT: + buf_in = imgs3i (im_in, x1, x2, y1, y2, z1, z2) + switch (which3d) { + case XYZ: + buf_out = imps3i (im_out, x1, x2, y1, y2, z1, z2) + call txyz3i (Memi[buf_in], Memi[buf_out], nx,ny,nz) + case XZY: + buf_out = imps3i (im_out, x1, x2, z1, z2, y1, y2) + call txzy3i (Memi[buf_in], Memi[buf_out], nx,ny,nz) + case YXZ: + buf_out = imps3i (im_out, y1, y2, x1, x2, z1, z2) + call tyxz3i (Memi[buf_in], Memi[buf_out], nx,ny,nz) + case YZX: + buf_out = imps3i (im_out, y1, y2, z1, z2, x1, x2) + call tyzx3i (Memi[buf_in], Memi[buf_out], nx,ny,nz) + case ZXY: + buf_out = imps3i (im_out, z1, z2, x1, x2, y1, y2) + call tzxy3i (Memi[buf_in], Memi[buf_out], nx,ny,nz) + case ZYX: + buf_out = imps3i (im_out, z1, z2, y1, y2, x1, x2) + call tzyx3i (Memi[buf_in], Memi[buf_out], nx,ny,nz) + } + + case TY_LONG, TY_USHORT: + buf_in = imgs3l (im_in, x1, x2, y1, y2, z1, z2) + switch (which3d) { + case XYZ: + buf_out = imps3l (im_out, x1, x2, y1, y2, z1, z2) + call txyz3l (Meml[buf_in], Meml[buf_out], nx,ny,nz) + case XZY: + buf_out = imps3l (im_out, x1, x2, z1, z2, y1, y2) + call txzy3l (Meml[buf_in], Meml[buf_out], nx,ny,nz) + case YXZ: + buf_out = imps3l (im_out, y1, y2, x1, x2, z1, z2) + call tyxz3l (Meml[buf_in], Meml[buf_out], nx,ny,nz) + case YZX: + buf_out = imps3l (im_out, y1, y2, z1, z2, x1, x2) + call tyzx3l (Meml[buf_in], Meml[buf_out], nx,ny,nz) + case ZXY: + buf_out = imps3l (im_out, z1, z2, x1, x2, y1, y2) + call tzxy3l (Meml[buf_in], Meml[buf_out], nx,ny,nz) + case ZYX: + buf_out = imps3l (im_out, z1, z2, y1, y2, x1, x2) + call tzyx3l (Meml[buf_in], Meml[buf_out], nx,ny,nz) + } + + case TY_REAL: + buf_in = imgs3r (im_in, x1, x2, y1, y2, z1, z2) + switch (which3d) { + case XYZ: + buf_out = imps3r (im_out, x1, x2, y1, y2, z1, z2) + call txyz3r (Memr[buf_in], Memr[buf_out], nx,ny,nz) + case XZY: + buf_out = imps3r (im_out, x1, x2, z1, z2, y1, y2) + call txzy3r (Memr[buf_in], Memr[buf_out], nx,ny,nz) + case YXZ: + buf_out = imps3r (im_out, y1, y2, x1, x2, z1, z2) + call tyxz3r (Memr[buf_in], Memr[buf_out], nx,ny,nz) + case YZX: + buf_out = imps3r (im_out, y1, y2, z1, z2, x1, x2) + call tyzx3r (Memr[buf_in], Memr[buf_out], nx,ny,nz) + case ZXY: + buf_out = imps3r (im_out, z1, z2, x1, x2, y1, y2) + call tzxy3r (Memr[buf_in], Memr[buf_out], nx,ny,nz) + case ZYX: + buf_out = imps3r (im_out, z1, z2, y1, y2, x1, x2) + call tzyx3r (Memr[buf_in], Memr[buf_out], nx,ny,nz) + } + + case TY_DOUBLE: + buf_in = imgs3d (im_in, x1, x2, y1, y2, z1, z2) + switch (which3d) { + case XYZ: + buf_out = imps3d (im_out, x1, x2, y1, y2, z1, z2) + call txyz3d (Memd[buf_in], Memd[buf_out], nx,ny,nz) + case XZY: + buf_out = imps3d (im_out, x1, x2, z1, z2, y1, y2) + call txzy3d (Memd[buf_in], Memd[buf_out], nx,ny,nz) + case YXZ: + buf_out = imps3d (im_out, y1, y2, x1, x2, z1, z2) + call tyxz3d (Memd[buf_in], Memd[buf_out], nx,ny,nz) + case YZX: + buf_out = imps3d (im_out, y1, y2, z1, z2, x1, x2) + call tyzx3d (Memd[buf_in], Memd[buf_out], nx,ny,nz) + case ZXY: + buf_out = imps3d (im_out, z1, z2, x1, x2, y1, y2) + call tzxy3d (Memd[buf_in], Memd[buf_out], nx,ny,nz) + case ZYX: + buf_out = imps3d (im_out, z1, z2, y1, y2, x1, x2) + call tzyx3d (Memd[buf_in], Memd[buf_out], nx,ny,nz) + } + + case TY_COMPLEX: + buf_in = imgs3x (im_in, x1, x2, y1, y2, z1, z2) + switch (which3d) { + case XYZ: + buf_out = imps3x (im_out, x1, x2, y1, y2, z1, z2) + call txyz3x (Memx[buf_in], Memx[buf_out], nx,ny,nz) + case XZY: + buf_out = imps3x (im_out, x1, x2, z1, z2, y1, y2) + call txzy3x (Memx[buf_in], Memx[buf_out], nx,ny,nz) + case YXZ: + buf_out = imps3x (im_out, y1, y2, x1, x2, z1, z2) + call tyxz3x (Memx[buf_in], Memx[buf_out], nx,ny,nz) + case YZX: + buf_out = imps3x (im_out, y1, y2, z1, z2, x1, x2) + call tyzx3x (Memx[buf_in], Memx[buf_out], nx,ny,nz) + case ZXY: + buf_out = imps3x (im_out, z1, z2, x1, x2, y1, y2) + call tzxy3x (Memx[buf_in], Memx[buf_out], nx,ny,nz) + case ZYX: + buf_out = imps3x (im_out, z1, z2, y1, y2, x1, x2) + call tzyx3x (Memx[buf_in], Memx[buf_out], nx,ny,nz) + } + + default: + call error (3, "unknown pixel type") + } + } + } + } +end + + +# WHICHTRAN -- Return the transpose type given the axes transpose list. + +int procedure whichtran (new_ax) + +int new_ax[3] #I the input axes transpose list. + +int which + +begin + which = 0 + + if (new_ax[1] == 1) { + if (new_ax[2] == 2) + which = XYZ + else if (new_ax[2] == 3) + which = XZY + } else if (new_ax[1] == 2) { + if (new_ax[2] == 1) + which = YXZ + else if (new_ax[2] == 3) + which = YZX + } else if (new_ax[1] == 3) { + if (new_ax[2] == 1) + which = ZXY + else if (new_ax[2] == 2) + which = ZYX + } + + return (which) +end + + +define LTM Memd[ltr+(($2)-1)*pdim+($1)-1] +define NCD Memd[ncd+(($2)-1)*pdim+($1)-1] +define swap {temp=$1;$1=$2;$2=temp} + +# IM3DTRMW -- Perform a transpose operation on the image WCS. + +procedure im3dtrmw (mw, which3d) + +pointer mw #I pointer to the mwcs structure +int which3d #I type of 3D transpose + +int i, axes[IM_MAXDIM], axval[IM_MAXDIM] +int naxes, pdim, nelem, axmap, ax1, ax2, ax3, szatstr +pointer sp, ltr, ltm, ltv, cd, r, w, ncd, nr +pointer attribute1, attribute2, attribute3, atstr1, atstr2, atstr3, mwtmp +double temp +int mw_stati(), itoc(), strlen() +pointer mw_open() +errchk mw_gwattrs(), mw_newsystem() + +begin + # Convert axis bitflags to the axis lists. + call mw_gaxlist (mw, 07B, axes, naxes) + if (naxes < 2) + return + + # Get the dimensions of the wcs and turn off axis mapping. + pdim = mw_stati (mw, MW_NPHYSDIM) + nelem = pdim * pdim + axmap = mw_stati (mw, MW_USEAXMAP) + call mw_seti (mw, MW_USEAXMAP, NO) + szatstr = SZ_LINE + + # Allocate working space. + call smark (sp) + call salloc (ltr, nelem, TY_DOUBLE) + call salloc (cd, nelem, TY_DOUBLE) + call salloc (r, pdim, TY_DOUBLE) + call salloc (w, pdim, TY_DOUBLE) + call salloc (ltm, nelem, TY_DOUBLE) + call salloc (ltv, pdim, TY_DOUBLE) + call salloc (ncd, nelem, TY_DOUBLE) + call salloc (nr, pdim, TY_DOUBLE) + call salloc (attribute1, SZ_FNAME, TY_CHAR) + call salloc (attribute2, SZ_FNAME, TY_CHAR) + call salloc (attribute3, SZ_FNAME, TY_CHAR) + + # Get the wterm which corresponds to the original logical to + # world transformation. + call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], pdim) + call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim) + call mwvmuld (Memd[ltm], Memd[r], Memd[nr], pdim) + call aaddd (Memd[nr], Memd[ltv], Memd[nr], pdim) + call mwinvertd (Memd[ltm], Memd[ltr], pdim) + call mwmmuld (Memd[cd], Memd[ltr], Memd[ncd], pdim) + + # Define which physical axes the logical axes correspond to. + # and recompute the above wterm to take into account the transpose. + ax1 = axes[1] + ax2 = axes[2] + ax3 = axes[3] + + switch (which3d) { + case XYZ: + # do nothing + + case XZY: + # switch axes 3 and 2 + call amovd (Memd[ncd], Memd[ltr], nelem) + NCD(ax1,ax1) = LTM(ax1,ax1) + NCD(ax2,ax1) = LTM(ax3,ax1) + NCD(ax3,ax1) = LTM(ax2,ax1) + NCD(ax1,ax2) = LTM(ax1,ax3) + NCD(ax2,ax2) = LTM(ax3,ax3) + NCD(ax3,ax2) = LTM(ax2,ax3) + NCD(ax1,ax3) = LTM(ax1,ax2) + NCD(ax2,ax3) = LTM(ax3,ax2) + NCD(ax3,ax3) = LTM(ax2,ax2) + swap (Memd[w+ax3-1], Memd[w+ax2-1]) + swap (Memd[nr+ax3-1], Memd[nr+ax2-1]) + + case YXZ: + # switch axes 1 and 2 + call amovd (Memd[ncd], Memd[ltr], nelem) + NCD(ax1,ax1) = LTM(ax2,ax2) + NCD(ax2,ax1) = LTM(ax1,ax2) + NCD(ax3,ax1) = LTM(ax3,ax2) + NCD(ax1,ax2) = LTM(ax2,ax1) + NCD(ax2,ax2) = LTM(ax1,ax1) + NCD(ax3,ax2) = LTM(ax3,ax1) + NCD(ax1,ax3) = LTM(ax2,ax3) + NCD(ax2,ax3) = LTM(ax1,ax3) + NCD(ax3,ax3) = LTM(ax3,ax3) + swap (Memd[w+ax1-1], Memd[w+ax2-1]) + swap (Memd[nr+ax1-1], Memd[nr+ax2-1]) + + case YZX: + # map axes 123 to 231 + call amovd (Memd[ncd], Memd[ltr], nelem) + NCD(ax1,ax1) = LTM(ax2,ax2) + NCD(ax2,ax1) = LTM(ax3,ax2) + NCD(ax3,ax1) = LTM(ax1,ax2) + NCD(ax1,ax2) = LTM(ax2,ax3) + NCD(ax2,ax2) = LTM(ax3,ax3) + NCD(ax3,ax2) = LTM(ax1,ax3) + NCD(ax1,ax3) = LTM(ax2,ax1) + NCD(ax2,ax3) = LTM(ax3,ax1) + NCD(ax3,ax3) = LTM(ax1,ax1) + call amovd (Memd[w], Memd[ltv], pdim) + Memd[w+ax1-1] = Memd[ltv+ax2-1] + Memd[w+ax2-1] = Memd[ltv+ax3-1] + Memd[w+ax3-1] = Memd[ltv+ax1-1] + call amovd (Memd[nr], Memd[ltv], pdim) + Memd[nr+ax1-1] = Memd[ltv+ax2-1] + Memd[nr+ax2-1] = Memd[ltv+ax3-1] + Memd[nr+ax3-1] = Memd[ltv+ax1-1] + + case ZXY: + # map axes 123 to 312 + call amovd (Memd[ncd], Memd[ltr], nelem) + NCD(ax1,ax1) = LTM(ax3,ax3) + NCD(ax2,ax1) = LTM(ax1,ax3) + NCD(ax3,ax1) = LTM(ax2,ax3) + NCD(ax1,ax2) = LTM(ax3,ax1) + NCD(ax2,ax2) = LTM(ax1,ax1) + NCD(ax3,ax2) = LTM(ax2,ax1) + NCD(ax1,ax3) = LTM(ax3,ax2) + NCD(ax2,ax3) = LTM(ax1,ax2) + NCD(ax3,ax3) = LTM(ax2,ax2) + call amovd (Memd[w], Memd[ltv], pdim) + Memd[w+ax1-1] = Memd[ltv+ax3-1] + Memd[w+ax2-1] = Memd[ltv+ax1-1] + Memd[w+ax3-1] = Memd[ltv+ax2-1] + call amovd (Memd[nr], Memd[ltv], pdim) + Memd[nr+ax1-1] = Memd[ltv+ax3-1] + Memd[nr+ax2-1] = Memd[ltv+ax1-1] + Memd[nr+ax3-1] = Memd[ltv+ax2-1] + + case ZYX: + # switch axes 3 and 1 + call amovd (Memd[ncd], Memd[ltr], nelem) + NCD(ax1,ax1) = LTM(ax3,ax3) + NCD(ax2,ax1) = LTM(ax2,ax3) + NCD(ax3,ax1) = LTM(ax1,ax3) + NCD(ax1,ax2) = LTM(ax3,ax2) + NCD(ax2,ax2) = LTM(ax2,ax2) + NCD(ax3,ax2) = LTM(ax1,ax2) + NCD(ax1,ax3) = LTM(ax3,ax1) + NCD(ax2,ax3) = LTM(ax2,ax1) + NCD(ax3,ax3) = LTM(ax1,ax1) + swap (Memd[w+ax1-1], Memd[w+ax3-1]) + swap (Memd[nr+ax1-1], Memd[nr+ax3-1]) + } + + # Perform the transpose of the lterm. + call mw_mkidmd (Memd[ltr], pdim) + switch (which3d) { + + case XYZ: + # do nothing + + case XZY: + # switch axes 3 and 2 + LTM(ax2,ax2) = 0.0d0 + LTM(ax3,ax2) = 1.0d0 + LTM(ax2,ax3) = 1.0d0 + LTM(ax3,ax3) = 0.0d0 + + case YXZ: + # switch axes 1 and 2 + LTM(ax1,ax1) = 0.0d0 + LTM(ax1,ax2) = 1.0d0 + LTM(ax2,ax1) = 1.0d0 + LTM(ax2,ax2) = 0.0d0 + + case YZX: + # map axes 123 to 231 + LTM(ax1,ax1) = 0.0d0 + LTM(ax1,ax2) = 1.0d0 + LTM(ax1,ax3) = 0.0d0 + LTM(ax2,ax1) = 0.0d0 + LTM(ax2,ax2) = 0.0d0 + LTM(ax2,ax3) = 1.0d0 + LTM(ax3,ax1) = 1.0d0 + LTM(ax3,ax2) = 0.0d0 + LTM(ax3,ax3) = 0.0d0 + + case ZXY: + # map axes 123 to 312 + LTM(ax1,ax1) = 0.0d0 + LTM(ax1,ax2) = 0.0d0 + LTM(ax1,ax3) = 1.0d0 + LTM(ax2,ax1) = 1.0d0 + LTM(ax2,ax2) = 0.0d0 + LTM(ax2,ax3) = 0.0d0 + LTM(ax3,ax1) = 0.0d0 + LTM(ax3,ax2) = 1.0d0 + LTM(ax3,ax3) = 0.0d0 + + case ZYX: + # switch axes 3 and 1 + LTM(ax3,ax3) = 0.0d0 + LTM(ax3,ax1) = 1.0d0 + LTM(ax1,ax3) = 1.0d0 + LTM(ax1,ax1) = 0.0d0 + + } + call aclrd (Memd[ltv], pdim) + call aclrd (Memd[r], pdim) + call mw_translated (mw, Memd[ltv], Memd[ltr], Memd[r], pdim) + + # Get the new lterm, recompute the wterm, and store it. + call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim) + call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], pdim) + call mwinvertd (Memd[ltm], Memd[ltr], pdim) + call asubd (Memd[nr], Memd[ltv], Memd[r], pdim) + call mwvmuld (Memd[ltr], Memd[r], Memd[nr], pdim) + call mw_swtermd (mw, Memd[nr], Memd[w], Memd[cd], pdim) + + # Make a new temporary wcs and set the system name. + mwtmp = mw_open (NULL, pdim) + call mw_gsystem (mw, Memc[attribute1], SZ_FNAME) + iferr (call mw_newsystem (mwtmp, Memc[attribute1], pdim)) + call mw_ssystem (mwtmp, Memc[attribute1]) + + # Copy the wterm and the lterm to it. + call mw_gwtermd (mw, Memd[r], Memd[w], Memd[ltr], pdim) + call mw_swtermd (mwtmp, Memd[r], Memd[w], Memd[ltr], pdim) + call mw_gltermd (mw, Memd[ltr], Memd[r], pdim) + call mw_sltermd (mwtmp, Memd[ltr], Memd[r], pdim) + + # Set the axis map and the axis types. + call mw_gaxmap (mw, axes, axval, pdim) + call mw_saxmap (mwtmp, axes, axval, pdim) + iferr (call mw_gwattrs (mw, ax1, "wtype", Memc[attribute1], SZ_FNAME)) + call strcpy ("linear", Memc[attribute1], SZ_FNAME) + iferr (call mw_gwattrs (mw, ax2, "wtype", Memc[attribute2], SZ_FNAME)) + call strcpy ("linear", Memc[attribute2], SZ_FNAME) + iferr (call mw_gwattrs (mw, ax3, "wtype", Memc[attribute3], SZ_FNAME)) + call strcpy ("linear", Memc[attribute3], SZ_FNAME) + + switch (which3d) { + case XYZ: + call mw_swtype (mwtmp, ax1, 1, Memc[attribute1], "") + call mw_swtype (mwtmp, ax2, 1, Memc[attribute2], "") + call mw_swtype (mwtmp, ax3, 1, Memc[attribute3], "") + case XZY: + call mw_swtype (mwtmp, ax1, 1, Memc[attribute1], "") + call mw_swtype (mwtmp, ax2, 1, Memc[attribute3], "") + call mw_swtype (mwtmp, ax3, 1, Memc[attribute2], "") + case YXZ: + call mw_swtype (mwtmp, ax1, 1, Memc[attribute2], "") + call mw_swtype (mwtmp, ax2, 1, Memc[attribute1], "") + call mw_swtype (mwtmp, ax3, 1, Memc[attribute3], "") + case YZX: + call mw_swtype (mwtmp, ax1, 1, Memc[attribute2], "") + call mw_swtype (mwtmp, ax2, 1, Memc[attribute3], "") + call mw_swtype (mwtmp, ax3, 1, Memc[attribute1], "") + case ZXY: + call mw_swtype (mwtmp, ax1, 1, Memc[attribute3], "") + call mw_swtype (mwtmp, ax2, 1, Memc[attribute1], "") + call mw_swtype (mwtmp, ax3, 1, Memc[attribute2], "") + case ZYX: + call mw_swtype (mwtmp, ax1, 1, Memc[attribute3], "") + call mw_swtype (mwtmp, ax2, 1, Memc[attribute2], "") + call mw_swtype (mwtmp, ax3, 1, Memc[attribute1], "") + } + + # Copy the axis attributes. + call malloc (atstr1, szatstr, TY_CHAR) + call malloc (atstr2, szatstr, TY_CHAR) + call malloc (atstr3, szatstr, TY_CHAR) + + for (i = 1; ; i = i + 1) { + + if (itoc (i, Memc[attribute1], SZ_FNAME) <= 0) + Memc[attribute1] = EOS + if (itoc (i, Memc[attribute2], SZ_FNAME) <= 0) + Memc[attribute2] = EOS + if (itoc (i, Memc[attribute3], SZ_FNAME) <= 0) + Memc[attribute3] = EOS + + repeat { + iferr (call mw_gwattrs (mw, ax1, Memc[attribute1], + Memc[atstr1], szatstr)) + Memc[atstr1] = EOS + iferr (call mw_gwattrs (mw, ax2, Memc[attribute2], + Memc[atstr2], szatstr)) + Memc[atstr2] = EOS + iferr (call mw_gwattrs (mw, ax3, Memc[attribute3], + Memc[atstr3], szatstr)) + Memc[atstr3] = EOS + if ((strlen (Memc[atstr1]) < szatstr) && + (strlen (Memc[atstr2]) < szatstr) && + (strlen (Memc[atstr3]) < szatstr)) + break + szatstr = szatstr + SZ_LINE + call realloc (atstr1, szatstr, TY_CHAR) + call realloc (atstr2, szatstr, TY_CHAR) + call realloc (atstr3, szatstr, TY_CHAR) + } + if ((Memc[atstr1] == EOS) && (Memc[atstr2] == EOS) && + (Memc[atstr3] == EOS)) + break + + switch (which3d) { + case XYZ: + if (Memc[atstr1] != EOS) + call mw_swattrs (mwtmp, ax1, Memc[attribute1], Memc[atstr1]) + if (Memc[atstr2] != EOS) + call mw_swattrs (mwtmp, ax2, Memc[attribute2], Memc[atstr2]) + if (Memc[atstr3] != EOS) + call mw_swattrs (mwtmp, ax3, Memc[attribute3], Memc[atstr3]) + case XZY: + if (Memc[atstr1] != EOS) + call mw_swattrs (mwtmp, ax1, Memc[attribute1], Memc[atstr1]) + if (Memc[atstr3] != EOS) + call mw_swattrs (mwtmp, ax2, Memc[attribute3], Memc[atstr3]) + if (Memc[atstr2] != EOS) + call mw_swattrs (mwtmp, ax3, Memc[attribute2], Memc[atstr2]) + case YXZ: + if (Memc[atstr2] != EOS) + call mw_swattrs (mwtmp, ax1, Memc[attribute2], Memc[atstr2]) + if (Memc[atstr1] != EOS) + call mw_swattrs (mwtmp, ax2, Memc[attribute1], Memc[atstr1]) + if (Memc[atstr3] != EOS) + call mw_swattrs (mwtmp, ax3, Memc[attribute3], Memc[atstr3]) + case YZX: + if (Memc[atstr2] != EOS) + call mw_swattrs (mwtmp, ax1, Memc[attribute2], Memc[atstr2]) + if (Memc[atstr3] != EOS) + call mw_swattrs (mwtmp, ax2, Memc[attribute3], Memc[atstr3]) + if (Memc[atstr1] != EOS) + call mw_swattrs (mwtmp, ax3, Memc[attribute1], Memc[atstr1]) + case ZXY: + if (Memc[atstr3] != EOS) + call mw_swattrs (mwtmp, ax1, Memc[attribute3], Memc[atstr3]) + if (Memc[atstr1] != EOS) + call mw_swattrs (mwtmp, ax2, Memc[attribute1], Memc[atstr1]) + if (Memc[atstr2] != EOS) + call mw_swattrs (mwtmp, ax3, Memc[attribute2], Memc[atstr2]) + case ZYX: + if (Memc[atstr3] != EOS) + call mw_swattrs (mwtmp, ax1, Memc[attribute3], Memc[atstr3]) + if (Memc[atstr2] != EOS) + call mw_swattrs (mwtmp, ax2, Memc[attribute2], Memc[atstr2]) + if (Memc[atstr1] != EOS) + call mw_swattrs (mwtmp, ax3, Memc[attribute1], Memc[atstr1]) + } + + } + call mfree (atstr1, TY_CHAR) + call mfree (atstr2, TY_CHAR) + call mfree (atstr3, TY_CHAR) + call mw_close (mw) + + # Delete the old wcs and set equal to the new one. + call sfree (sp) + mw = mwtmp + call mw_seti (mw, MW_USEAXMAP, axmap) +end diff --git a/pkg/images/imgeom/src/t_imshift.x b/pkg/images/imgeom/src/t_imshift.x new file mode 100644 index 00000000..4a940b99 --- /dev/null +++ b/pkg/images/imgeom/src/t_imshift.x @@ -0,0 +1,530 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include <imset.h> +include <math/iminterp.h> + +define NYOUT 16 # number of lines output at once +define NMARGIN 3 # number of boundary pixels required +define NMARGIN_SPLINE3 16 # number of spline boundary pixels required + + +# T_IMSHIFT -- Shift a 2-D image by an arbitrary amount in X and Y, using +# boundary extension to preserve the image size. + +procedure t_imshift() + +pointer imtlist1 # Input image list +pointer imtlist2 # Output image list + +pointer image1 # Input image +pointer image2 # Output image +pointer imtemp # Temporary file +pointer sfile # Text file containing list of shifts +pointer interpstr # Interpolant string + +int list1, list2, boundary_type, ixshift, iyshift, nshifts, interp_type +pointer sp, str, xs, ys, im1, im2, sf, mw +real constant, shifts[2] +double txshift, tyshift, xshift, yshift + +bool fp_equald(), envgetb() +int imtgetim(), imtlen(), clgwrd(), strdic(), open(), ish_rshifts() +pointer immap(), imtopen(), mw_openim() +real clgetr() +double clgetd() +errchk ish_ishiftxy, ish_gshiftxy, mw_openim, mw_saveim, mw_shift + +begin + call smark (sp) + call salloc (imtlist1, SZ_LINE, TY_CHAR) + call salloc (imtlist2, SZ_LINE, TY_CHAR) + call salloc (image1, SZ_LINE, TY_CHAR) + call salloc (image2, SZ_LINE, TY_CHAR) + call salloc (imtemp, SZ_LINE, TY_CHAR) + call salloc (sfile, SZ_FNAME, TY_CHAR) + call salloc (interpstr, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Get task parameters. + call clgstr ("input", Memc[imtlist1], SZ_FNAME) + call clgstr ("output", Memc[imtlist2], SZ_FNAME) + call clgstr ("shifts_file", Memc[sfile], SZ_FNAME) + + # Get the 2-D interpolation parameters. + call clgstr ("interp_type", Memc[interpstr], SZ_FNAME) + interp_type = strdic (Memc[interpstr], Memc[str], SZ_LINE, + II_BFUNCTIONS) + boundary_type = clgwrd ("boundary_type", Memc[str], SZ_LINE, + ",constant,nearest,reflect,wrap,") + if (boundary_type == BT_CONSTANT) + constant = clgetr ("constant") + + # Open the input and output image lists. + list1 = imtopen (Memc[imtlist1]) + list2 = imtopen (Memc[imtlist2]) + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + call error (1, "Number of input and output images not the same.") + } + + # Determine the source of the shifts. + if (Memc[sfile] != EOS) { + sf = open (Memc[sfile], READ_ONLY, TEXT_FILE) + call salloc (xs, imtlen (list1), TY_DOUBLE) + call salloc (ys, imtlen (list1), TY_DOUBLE) + nshifts = ish_rshifts (sf, Memd[xs], Memd[ys], imtlen (list1)) + if (nshifts != imtlen (list1)) + call error (2, + "The number of input images and shifts are not the same.") + } else { + sf = NULL + txshift = clgetd ("xshift") + tyshift = clgetd ("yshift") + } + + + # Do each set of input and output images. + nshifts = 0 + while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF) && + (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF)) { + + call xt_mkimtemp (Memc[image1], Memc[image2], Memc[imtemp], + SZ_FNAME) + + im1 = immap (Memc[image1], READ_ONLY, 0) + im2 = immap (Memc[image2], NEW_COPY, im1) + + if (sf != NULL) { + xshift = Memd[xs+nshifts] + yshift = Memd[ys+nshifts] + } else { + xshift = txshift + yshift = tyshift + } + + ixshift = int (xshift) + iyshift = int (yshift) + + iferr { + # Perform the shift. + if (interp_type == II_BINEAREST) { + call ish_ishiftxy (im1, im2, nint(xshift), nint(yshift), + boundary_type, constant) + } else if (fp_equald (xshift, double(ixshift)) && + fp_equald (yshift, double(iyshift))) { + call ish_ishiftxy (im1, im2, ixshift, iyshift, + boundary_type, constant) + } else { + call ish_gshiftxy (im1, im2, xshift, yshift, + Memc[interpstr], boundary_type, constant) + } + + # Update the image WCS to reflect the shift. + if (!envgetb ("nomwcs")) { + mw = mw_openim (im1) + shifts[1] = xshift + shifts[2] = yshift + call mw_shift (mw, shifts, 03B) + call mw_saveim (mw, im2) + call mw_close (mw) + } + + } then { + call eprintf ("Error shifting image: %s\n") + call pargstr (Memc[image1]) + call erract (EA_WARN) + call imunmap (im2) + call imunmap (im1) + call imdelete (Memc[image2]) + + } else { + # Finish up. + call imunmap (im2) + call imunmap (im1) + call xt_delimtemp (Memc[image2], Memc[imtemp]) + } + + nshifts = nshifts + 1 + } + + if (sf != NULL) + call close (sf) + call imtclose (list1) + call imtclose (list2) + call sfree (sp) +end + + +# ISH_ISHIFTXY -- Shift a 2-D image by integral pixels in x and y. + +procedure ish_ishiftxy (im1, im2, ixshift, iyshift, boundary_type, + constant) + +pointer im1 #I pointer to the input image +pointer im2 #I pointer to the output image +int ixshift #I shift in x and y +int iyshift #I +int boundary_type #I type of boundary extension +real constant #I constant for boundary extension + +pointer buf1, buf2 +long v[IM_MAXDIM] +int ncols, nlines, nbpix +int i, x1col, x2col, yline + +int impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx() +pointer imgs2s(), imgs2i(), imgs2l(), imgs2r(), imgs2d(), imgs2x() +errchk impnls, impnli, impnll, impnlr, impnld, impnlx +errchk imgs2s, imgs2i, imgs2l, imgs2r, imgs2d, imgs2x +string wrerr "ISHIFTXY: Error writing in image." + +begin + ncols = IM_LEN(im1,1) + nlines = IM_LEN(im1,2) + + # Cannot shift off image. + if (ixshift < -ncols || ixshift > ncols) + call error (3, "ISHIFTXY: X shift out of bounds.") + if (iyshift < -nlines || iyshift > nlines) + call error (4, "ISHIFTXY: Y shift out of bounds.") + + # Calculate the shift. + switch (boundary_type) { + case BT_CONSTANT,BT_REFLECT,BT_NEAREST: + ixshift = min (ncols, max (-ncols, ixshift)) + iyshift = min (nlines, max (-nlines, iyshift)) + case BT_WRAP: + ixshift = mod (ixshift, ncols) + iyshift = mod (iyshift, nlines) + } + + # Set the boundary extension values. + nbpix = max (abs (ixshift), abs (iyshift)) + call imseti (im1, IM_NBNDRYPIX, nbpix) + call imseti (im1, IM_TYBNDRY, boundary_type) + if (boundary_type == BT_CONSTANT) + call imsetr (im1, IM_BNDRYPIXVAL, constant) + + # Get column boundaries in the input image. + x1col = max (-ncols + 1, - ixshift + 1) + x2col = min (2 * ncols, ncols - ixshift) + + call amovkl (long (1), v, IM_MAXDIM) + + # Shift the image using the appropriate data type operators. + switch (IM_PIXTYPE(im1)) { + case TY_SHORT: + do i = 1, nlines { + if (impnls (im2, buf2, v) == EOF) + call error (5, wrerr) + yline = i - iyshift + buf1 = imgs2s (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (5, wrerr) + call amovs (Mems[buf1], Mems[buf2], ncols) + } + case TY_INT: + do i = 1, nlines { + if (impnli (im2, buf2, v) == EOF) + call error (5, wrerr) + yline = i - iyshift + buf1 = imgs2i (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (5, wrerr) + call amovi (Memi[buf1], Memi[buf2], ncols) + } + case TY_USHORT, TY_LONG: + do i = 1, nlines { + if (impnll (im2, buf2, v) == EOF) + call error (5, wrerr) + yline = i - iyshift + buf1 = imgs2l (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (5, wrerr) + call amovl (Meml[buf1], Meml[buf2], ncols) + } + case TY_REAL: + do i = 1, nlines { + if (impnlr (im2, buf2, v) == EOF) + call error (5, wrerr) + yline = i - iyshift + buf1 = imgs2r (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (5, wrerr) + call amovr (Memr[buf1], Memr[buf2], ncols) + } + case TY_DOUBLE: + do i = 1, nlines { + if (impnld (im2, buf2, v) == EOF) + call error (0, wrerr) + yline = i - iyshift + buf1 = imgs2d (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (0, wrerr) + call amovd (Memd[buf1], Memd[buf2], ncols) + } + case TY_COMPLEX: + do i = 1, nlines { + if (impnlx (im2, buf2, v) == EOF) + call error (0, wrerr) + yline = i - iyshift + buf1 = imgs2x (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (0, wrerr) + call amovx (Memx[buf1], Memx[buf2], ncols) + } + default: + call error (6, "ISHIFTXY: Unknown IRAF type.") + } +end + + +# ISH_GSHIFTXY -- Shift an image by fractional pixels in x and y. +# Unfortunately, this code currently performs the shift only on single +# precision real, so precision is lost if the data is of type double, +# and the imaginary component is lost if the data is of type complex. + +procedure ish_gshiftxy (im1, im2, xshift, yshift, interpstr, boundary_type, + constant) + +pointer im1 #I pointer to input image +pointer im2 #I pointer to output image +double xshift #I shift in x direction +double yshift #I shift in y direction +char interpstr[ARB] #I type of interpolant +int boundary_type #I type of boundary extension +real constant #I value of constant for boundary extension + +int lout1, lout2, nyout, nxymargin, interp_type, nsinc, nincr +int cin1, cin2, nxin, lin1, lin2, nyin, i +int ncols, nlines, nbpix, fstline, lstline +double xshft, yshft, deltax, deltay, dx, dy, cx, ly +pointer sp, x, y, msi, sinbuf, soutbuf + +pointer imps2r() +int msigeti() +bool fp_equald() +errchk msisinit(), msifree(), msifit(), msigrid() +errchk imgs2r(), imps2r() + +begin + ncols = IM_LEN(im1,1) + nlines = IM_LEN(im1,2) + + # Check for out of bounds shift. + if (xshift < -ncols || xshift > ncols) + call error (7, "GSHIFTXY: X shift out of bounds.") + if (yshift < -nlines || yshift > nlines) + call error (8, "GSHIFTXY: Y shift out of bounds.") + + # Get the real shift. + if (boundary_type == BT_WRAP) { + xshft = mod (xshift, real (ncols)) + yshft = mod (yshift, real (nlines)) + } else { + xshft = xshift + yshft = yshift + } + + # Allocate temporary space. + call smark (sp) + call salloc (x, 2 * ncols, TY_REAL) + call salloc (y, 2 * nlines, TY_REAL) + sinbuf = NULL + + # Define the x and y shifts for the interpolation. + dx = abs (xshft - int (xshft)) + if (fp_equald (dx, 0D0)) + deltax = 0.0 + else if (xshft > 0.) + deltax = 1. - dx + else + deltax = dx + dy = abs (yshft - int (yshft)) + if (fp_equald (dy, 0D0)) + deltay = 0.0 + else if (yshft > 0.) + deltay = 1. - dy + else + deltay = dy + + # Initialize the 2-D interpolation routines. + call msitype (interpstr, interp_type, nsinc, nincr, cx) + if (interp_type == II_BILSINC || interp_type == II_BISINC ) + call msisinit (msi, II_BILSINC, nsinc, 1, 1, + deltax - nint (deltax), deltay - nint (deltay), 0.0) + else + call msisinit (msi, interp_type, nsinc, nincr, nincr, cx, cx, 0.0) + + # Set boundary extension parameters. + if (interp_type == II_BISPLINE3) + nxymargin = NMARGIN_SPLINE3 + else if (interp_type == II_BISINC || interp_type == II_BILSINC) + nxymargin = msigeti (msi, II_MSINSINC) + else + nxymargin = NMARGIN + nbpix = max (int (abs(xshft)+1.0), int (abs(yshft)+1.0)) + nxymargin + call imseti (im1, IM_NBNDRYPIX, nbpix) + call imseti (im1, IM_TYBNDRY, boundary_type) + if (boundary_type == BT_CONSTANT) + call imsetr (im1, IM_BNDRYPIXVAL, constant) + + # Define the x interpolation coordinates + deltax = deltax + nxymargin + if (interp_type == II_BIDRIZZLE) { + do i = 1, ncols { + Memr[x+2*i-2] = i + deltax - 0.5 + Memr[x+2*i-1] = i + deltax + 0.5 + } + } else { + do i = 1, ncols + Memr[x+i-1] = i + deltax + } + + # Define the y interpolation coordinates. + deltay = deltay + nxymargin + if (interp_type == II_BIDRIZZLE) { + do i = 1, NYOUT { + Memr[y+2*i-2] = i + deltay - 0.5 + Memr[y+2*i-1] = i + deltay + 0.5 + } + } else { + do i = 1, NYOUT + Memr[y+i-1] = i + deltay + } + + # Define column ranges in the input image. + cx = 1. - nxymargin - xshft + if ((cx <= 0.) && (! fp_equald (dx, 0D0))) + cin1 = int (cx) - 1 + else + cin1 = int (cx) + cin2 = ncols - xshft + nxymargin + 1 + nxin = cin2 - cin1 + 1 + + # Loop over output sections. + for (lout1 = 1; lout1 <= nlines; lout1 = lout1 + NYOUT) { + + # Define range of output lines. + lout2 = min (lout1 + NYOUT - 1, nlines) + nyout = lout2 - lout1 + 1 + + # Define correspoding range of input lines. + ly = lout1 - nxymargin - yshft + if ((ly <= 0.0) && (! fp_equald (dy, 0D0))) + lin1 = int (ly) - 1 + else + lin1 = int (ly) + lin2 = lout2 - yshft + nxymargin + 1 + nyin = lin2 - lin1 + 1 + + # Get appropriate input section and calculate the coefficients. + if ((sinbuf == NULL) || (lin1 < fstline) || (lin2 > lstline)) { + fstline = lin1 + lstline = lin2 + call ish_buf (im1, cin1, cin2, lin1, lin2, sinbuf) + call msifit (msi, Memr[sinbuf], nxin, nyin, nxin) + } + + # Output the section. + soutbuf = imps2r (im2, 1, ncols, lout1, lout2) + if (soutbuf == EOF) + call error (9, "GSHIFTXY: Error writing output image.") + + # Evaluate the interpolant. + call msigrid (msi, Memr[x], Memr[y], Memr[soutbuf], ncols, nyout, + ncols) + } + + if (sinbuf != NULL) + call mfree (sinbuf, TY_REAL) + + call msifree (msi) + call sfree (sp) +end + + +# ISH_BUF -- Provide a buffer of image lines with minimum reads. + +procedure ish_buf (im, col1, col2, line1, line2, buf) + +pointer im #I pointer to input image +int col1, col2 #I column range of input buffer +int line1, line2 #I line range of input buffer +pointer buf #U buffer + +pointer buf1, buf2 +int i, ncols, nlines, nclast, llast1, llast2, nllast +errchk malloc, realloc +pointer imgs2r() + +begin + ncols = col2 - col1 + 1 + nlines = line2 - line1 + 1 + + # Make sure the buffer is large enough. + if (buf == NULL) { + call malloc (buf, ncols * nlines, TY_REAL) + llast1 = line1 - nlines + llast2 = line2 - nlines + } else if ((nlines != nllast) || (ncols != nclast)) { + call realloc (buf, ncols * nlines, TY_REAL) + llast1 = line1 - nlines + llast2 = line2 - nlines + } + + # The buffers must be contiguous. + if (line1 < llast1) { + do i = line2, line1, -1 { + if (i > llast1) + buf1 = buf + (i - llast1) * ncols + else + buf1 = imgs2r (im, col1, col2, i, i) + buf2 = buf + (i - line1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } else if (line2 > llast2) { + do i = line1, line2 { + if (i < llast2) + buf1 = buf + (i - llast1) * ncols + else + buf1 = imgs2r (im, col1, col2, i, i) + buf2 = buf + (i - line1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } + + llast1 = line1 + llast2 = line2 + nclast = ncols + nllast = nlines +end + + +# ISH_RSHIFTS -- Read shifts from a file. + +int procedure ish_rshifts (fd, x, y, max_nshifts) + +int fd #I shifts file +double x[ARB] #O x array +double y[ARB] #O y array +int max_nshifts #I the maximum number of shifts + +int nshifts +int fscan(), nscan() + +begin + nshifts = 0 + while (fscan (fd) != EOF && nshifts < max_nshifts) { + call gargd (x[nshifts+1]) + call gargd (y[nshifts+1]) + if (nscan () != 2) + next + nshifts = nshifts + 1 + } + + return (nshifts) +end diff --git a/pkg/images/imgeom/src/t_imtrans.x b/pkg/images/imgeom/src/t_imtrans.x new file mode 100644 index 00000000..04fa1d61 --- /dev/null +++ b/pkg/images/imgeom/src/t_imtrans.x @@ -0,0 +1,299 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <error.h> +include <mwset.h> + +# T_IMTRANSPOSE -- Transpose images. +# +# The input and output images are given by image template lists. The +# number of output images must match the number of input images. Image +# sections are allowed in the input images and are ignored in the output +# images. If the input and output image names are the same then the transpose +# is performed to a temporary file which then replaces the input image. + +procedure t_imtranspose () + +char imtlist1[SZ_LINE] # Input image list +char imtlist2[SZ_LINE] # Output image list +int len_blk # 1D length of transpose block + +char image1[SZ_FNAME] # Input image name +char image2[SZ_FNAME] # Output image name +char imtemp[SZ_FNAME] # Temporary file + +int list1, list2 +pointer im1, im2, mw + +bool envgetb() +int clgeti(), imtopen(), imtgetim(), imtlen() +pointer immap(), mw_openim() + +begin + # Get input and output image template lists and the size of + # the transpose block. + + call clgstr ("input", imtlist1, SZ_LINE) + call clgstr ("output", imtlist2, SZ_LINE) + len_blk = clgeti ("len_blk") + + # Expand the input and output image lists. + + list1 = imtopen (imtlist1) + list2 = imtopen (imtlist2) + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + call error (0, "Number of input and output images not the same") + } + + # Do each set of input/output images. + + while ((imtgetim (list1, image1, SZ_FNAME) != EOF) && + (imtgetim (list2, image2, SZ_FNAME) != EOF)) { + + call xt_mkimtemp (image1, image2, imtemp, SZ_FNAME) + + im1 = immap (image1, READ_ONLY, 0) + im2 = immap (image2, NEW_COPY, im1) + + # Do the transpose. + call imtranspose (im1, im2, len_blk) + + # Update the image WCS to reflect the shift. + if (!envgetb ("nomwcs")) { + mw = mw_openim (im1) + call imtrmw (mw) + call mw_saveim (mw, im2) + call mw_close (mw) + } + + # Unmap the input and output images. + call imunmap (im2) + call imunmap (im1) + + call xt_delimtemp (image2, imtemp) + } + + call imtclose (list1) + call imtclose (list2) +end + + +# IMTRANSPOSE -- Transpose an image. +# +# Divide the image into square blocks of size len_blk by len_blk. +# Transpose each block with a generic array transpose operator. + +procedure imtranspose (im_in, im_out, len_blk) + +pointer im_in # Input image descriptor +pointer im_out # Output image descriptor +int len_blk # 1D length of transpose block + +int x1, x2, nx +int y1, y2, ny +pointer buf_in, buf_out + +pointer imgs2s(), imps2s(), imgs2i(), imps2i(), imgs2l(), imps2l() +pointer imgs2r(), imps2r(), imgs2d(), imps2d(), imgs2x(), imps2x() + +begin + # Output image is a copy of input image with dims transposed. + + IM_LEN (im_out, 1) = IM_LEN (im_in, 2) + IM_LEN (im_out, 2) = IM_LEN (im_in, 1) + + # Break the input image into blocks of at most len_blk by len_blk. + + do x1 = 1, IM_LEN (im_in, 1), len_blk { + x2 = x1 + len_blk - 1 + if (x2 > IM_LEN(im_in, 1)) + x2 = IM_LEN(im_in, 1) + nx = x2 - x1 + 1 + + do y1 = 1, IM_LEN (im_in, 2), len_blk { + y2 = y1 + len_blk - 1 + if (y2 > IM_LEN(im_in, 2)) + y2 = IM_LEN(im_in, 2) + ny = y2 - y1 + 1 + + # Switch on the pixel type to optimize IMIO. + + switch (IM_PIXTYPE (im_in)) { + case TY_SHORT: + buf_in = imgs2s (im_in, x1, x2, y1, y2) + buf_out = imps2s (im_out, y1, y2, x1, x2) + call imtr2s (Mems[buf_in], Mems[buf_out], nx, ny) + case TY_INT: + buf_in = imgs2i (im_in, x1, x2, y1, y2) + buf_out = imps2i (im_out, y1, y2, x1, x2) + call imtr2i (Memi[buf_in], Memi[buf_out], nx, ny) + case TY_USHORT, TY_LONG: + buf_in = imgs2l (im_in, x1, x2, y1, y2) + buf_out = imps2l (im_out, y1, y2, x1, x2) + call imtr2l (Meml[buf_in], Meml[buf_out], nx, ny) + case TY_REAL: + buf_in = imgs2r (im_in, x1, x2, y1, y2) + buf_out = imps2r (im_out, y1, y2, x1, x2) + call imtr2r (Memr[buf_in], Memr[buf_out], nx, ny) + case TY_DOUBLE: + buf_in = imgs2d (im_in, x1, x2, y1, y2) + buf_out = imps2d (im_out, y1, y2, x1, x2) + call imtr2d (Memd[buf_in], Memd[buf_out], nx, ny) + case TY_COMPLEX: + buf_in = imgs2x (im_in, x1, x2, y1, y2) + buf_out = imps2x (im_out, y1, y2, x1, x2) + call imtr2x (Memx[buf_in], Memx[buf_out], nx, ny) + default: + call error (0, "unknown pixel type") + } + } + } +end + +define LTM Memd[ltr+(($2)-1)*pdim+($1)-1] +define NCD Memd[ncd+(($2)-1)*pdim+($1)-1] +define swap {temp=$1;$1=$2;$2=temp} + + +# IMTRMW -- Perform a transpose operation on the image WCS. + +procedure imtrmw (mw) + +pointer mw # pointer to the mwcs structure + +int i, axes[IM_MAXDIM], axval[IM_MAXDIM] +int naxes, pdim, nelem, axmap, ax1, ax2, szatstr +pointer sp, ltr, ltm, ltv, cd, r, w, ncd, nr +pointer attribute1, attribute2, atstr1, atstr2, mwtmp +double temp +int mw_stati(), itoc(), strlen() +pointer mw_open() +errchk mw_gwattrs(), mw_newsystem() + +begin + # Convert axis bitflags to the axis lists. + call mw_gaxlist (mw, 03B, axes, naxes) + if (naxes < 2) + return + + # Get the dimensions of the wcs and turn off axis mapping. + pdim = mw_stati (mw, MW_NPHYSDIM) + nelem = pdim * pdim + axmap = mw_stati (mw, MW_USEAXMAP) + call mw_seti (mw, MW_USEAXMAP, NO) + szatstr = SZ_LINE + + # Allocate working space. + call smark (sp) + call salloc (ltr, nelem, TY_DOUBLE) + call salloc (cd, nelem, TY_DOUBLE) + call salloc (r, pdim, TY_DOUBLE) + call salloc (w, pdim, TY_DOUBLE) + call salloc (ltm, nelem, TY_DOUBLE) + call salloc (ltv, pdim, TY_DOUBLE) + call salloc (ncd, nelem, TY_DOUBLE) + call salloc (nr, pdim, TY_DOUBLE) + call salloc (attribute1, SZ_FNAME, TY_CHAR) + call salloc (attribute2, SZ_FNAME, TY_CHAR) + + # Get the wterm which corresponds to the original logical to + # world transformation. + call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], pdim) + call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim) + call mwvmuld (Memd[ltm], Memd[r], Memd[nr], pdim) + call aaddd (Memd[nr], Memd[ltv], Memd[nr], pdim) + call mwinvertd (Memd[ltm], Memd[ltr], pdim) + call mwmmuld (Memd[cd], Memd[ltr], Memd[ncd], pdim) + + # Define which physical axes the logical axes correspond to. + # and recompute the above wterm to take into account the transpose. + ax1 = axes[1] + ax2 = axes[2] + swap (NCD(ax1,ax1), NCD(ax2,ax2)) + swap (NCD(ax1,ax2), NCD(ax2,ax1)) + swap (Memd[w+ax1-1], Memd[w+ax2-1]) + swap (Memd[nr+ax1-1], Memd[nr+ax2-1]) + + # Perform the transpose of the lterm. + call mw_mkidmd (Memd[ltr], pdim) + LTM(ax1,ax1) = 0.0d0 + LTM(ax1,ax2) = 1.0d0 + LTM(ax2,ax1) = 1.0d0 + LTM(ax2,ax2) = 0.0d0 + call aclrd (Memd[ltv], pdim) + call aclrd (Memd[r], pdim) + call mw_translated (mw, Memd[ltv], Memd[ltr], Memd[r], pdim) + + # Get the new lterm, recompute the wterm, and store it. + call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim) + call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], pdim) + call mwinvertd (Memd[ltm], Memd[ltr], pdim) + call asubd (Memd[nr], Memd[ltv], Memd[r], pdim) + call mwvmuld (Memd[ltr], Memd[r], Memd[nr], pdim) + call mw_swtermd (mw, Memd[nr], Memd[w], Memd[cd], pdim) + + # Make a new temporary wcs and set the system name. + mwtmp = mw_open (NULL, pdim) + call mw_gsystem (mw, Memc[attribute1], SZ_FNAME) + iferr (call mw_newsystem (mwtmp, Memc[attribute1], pdim)) + call mw_ssystem (mwtmp, Memc[attribute1]) + + # Copy the wterm and the lterm to it. + call mw_gwtermd (mw, Memd[r], Memd[w], Memd[ltr], pdim) + call mw_swtermd (mwtmp, Memd[r], Memd[w], Memd[ltr], pdim) + call mw_gltermd (mw, Memd[ltr], Memd[r], pdim) + call mw_sltermd (mwtmp, Memd[ltr], Memd[r], pdim) + + # Set the axis map and the axis types. + call mw_gaxmap (mw, axes, axval, pdim) + call mw_saxmap (mwtmp, axes, axval, pdim) + iferr (call mw_gwattrs (mw, ax1, "wtype", Memc[attribute1], SZ_FNAME)) + call strcpy ("linear", Memc[attribute1], SZ_FNAME) + iferr (call mw_gwattrs (mw, ax2, "wtype", Memc[attribute2], SZ_FNAME)) + call strcpy ("linear", Memc[attribute2], SZ_FNAME) + call mw_swtype (mwtmp, ax1, 1, Memc[attribute2], "") + call mw_swtype (mwtmp, ax2, 1, Memc[attribute1], "") + + # Copy the axis attributes. + call malloc (atstr1, szatstr, TY_CHAR) + call malloc (atstr2, szatstr, TY_CHAR) + for (i = 1; ; i = i + 1) { + + if (itoc (i, Memc[attribute1], SZ_FNAME) <= 0) + Memc[attribute1] = EOS + if (itoc (i, Memc[attribute2], SZ_FNAME) <= 0) + Memc[attribute2] = EOS + + repeat { + iferr (call mw_gwattrs (mw, ax1, Memc[attribute1], + Memc[atstr1], szatstr)) + Memc[atstr1] = EOS + iferr (call mw_gwattrs (mw, ax2, Memc[attribute2], + Memc[atstr2], szatstr)) + Memc[atstr2] = EOS + if ((strlen (Memc[atstr1]) < szatstr) && + (strlen (Memc[atstr2]) < szatstr)) + break + szatstr = szatstr + SZ_LINE + call realloc (atstr1, szatstr, TY_CHAR) + call realloc (atstr2, szatstr, TY_CHAR) + } + if ((Memc[atstr1] == EOS) && (Memc[atstr2] == EOS)) + break + + if (Memc[atstr2] != EOS) + call mw_swattrs (mwtmp, ax1, Memc[attribute2], Memc[atstr2]) + if (Memc[atstr1] != EOS) + call mw_swattrs (mwtmp, ax2, Memc[attribute1], Memc[atstr1]) + } + call mfree (atstr1, TY_CHAR) + call mfree (atstr2, TY_CHAR) + call mw_close (mw) + + # Delete the old wcs and set equal to the new one. + call sfree (sp) + mw = mwtmp + call mw_seti (mw, MW_USEAXMAP, axmap) +end diff --git a/pkg/images/imgeom/src/t_magnify.x b/pkg/images/imgeom/src/t_magnify.x new file mode 100644 index 00000000..ed797500 --- /dev/null +++ b/pkg/images/imgeom/src/t_magnify.x @@ -0,0 +1,624 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include <imset.h> +include <math/iminterp.h> + +# T_MAGNIFY -- Change coordinate origin and pixel interval in 2D images. +# +# The input and output images are given by image template lists. The +# number of output images must match the number of input images. Image +# sections are allowed in the input images and ignored in the output +# images. If the input and output image names are the same then the +# magnification is performed to a temporary file which then replaces +# the input image. + +# Interpolation types and boundary extension types. + +define BTYPES "|constant|nearest|reflect|wrap|project|" +define SZ_BTYPE 8 + +procedure t_magnify () + +pointer input # Pointer to input image list +pointer output # Pointer to output image list +pointer interp # Pointer to image interpolation type +pointer boundary # Pointer to boundary extension type +real bconst # Boundary extension pixel value +real xmag, ymag # Image magnifications +real dx, dy # Step size +real x1, y1 # Starting coordinates +real x2, y2 # Ending coordinates +int flux # Flux conserve + +int list1, list2, btype, logfd +pointer sp, in, out, image1, image2, image3, mw, errmsg +real a, b, c, d, shifts[2], scale[2] + +bool clgetb(), envgetb(), fp_equalr() +int clgwrd(), imtopen(), imtgetim(), imtlen(), open(), btoi(), errget() +pointer mw_openim(), immap() +real clgetr() +errchk open(), mg_magnify1(), mg_magnify2() + +begin + call smark (sp) + call salloc (input, SZ_LINE, TY_CHAR) + call salloc (output, SZ_LINE, TY_CHAR) + call salloc (interp, SZ_FNAME, TY_CHAR) + call salloc (boundary, SZ_BTYPE, TY_CHAR) + call salloc (image1, SZ_FNAME, TY_CHAR) + call salloc (image2, SZ_FNAME, TY_CHAR) + call salloc (image3, SZ_FNAME, TY_CHAR) + call salloc (errmsg, SZ_FNAME, TY_CHAR) + + # Get task parameters. + call clgstr ("input", Memc[input], SZ_LINE) + call clgstr ("output", Memc[output], SZ_LINE) + call clgstr ("interpolation", Memc[interp], SZ_FNAME) + btype = clgwrd ("boundary", Memc[boundary], SZ_BTYPE, BTYPES) + bconst = clgetr ("constant") + a = clgetr ("x1") + b = clgetr ("x2") + dx = clgetr ("dx") + c = clgetr ("y1") + d = clgetr ("y2") + dy = clgetr ("dy") + flux = btoi (clgetb ("fluxconserve")) + + # If the pixel interval INDEF then use the a magnification factor + # to determine the pixel interval. + + if (IS_INDEF (dx)) { + xmag = clgetr ("xmag") + if (xmag < 0.0) + dx = -xmag + else if (xmag > 0.0) + dx = 1.0 / xmag + else + dx = 0.0 + } + + if (IS_INDEF (dy)) { + ymag = clgetr ("ymag") + if (ymag < 0.0) + dy = -ymag + else if (ymag > 0.0) + dy = 1.0 / ymag + else + dy = 0.0 + } + + if (fp_equalr (dx, 0.0) || fp_equalr (dy, 0.0)) { + call error (0, "Illegal magnification") + } else { + xmag = 1.0 / dx + ymag = 1.0 / dy + } + + + # Open the log file. + call clgstr ("logfile", Memc[image1], SZ_FNAME) + iferr (logfd = open (Memc[image1], APPEND, TEXT_FILE)) + logfd = NULL + + # Expand the input and output image lists. + list1 = imtopen (Memc[input]) + list2 = imtopen (Memc[output]) + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + call error (0, "Number of input and output images not the same") + } + + # Magnify each set of input/output images with the 2D interpolation + # package. + + while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF) && + (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF)) { + + # Map the input and output images. + call xt_mkimtemp (Memc[image1], Memc[image2], Memc[image3], + SZ_FNAME) + in = immap (Memc[image1], READ_ONLY, 0) + out = immap (Memc[image2], NEW_COPY, in) + + # Set the limits of the output image. + x1 = a + x2 = b + y1 = c + y2 = d + + # Magnify the image making sure to update the wcs. + iferr { + if (IM_NDIM(in) == 1) { + call mg_magnify1 (in, out, Memc[interp], btype, bconst, + x1, x2, dx, flux) + if (!envgetb ("nomwcs")) { + mw = mw_openim (in) + scale[1] = xmag + shifts[1] = 1. - xmag * x1 + call mw_scale (mw, scale, 01B) + call mw_shift (mw, shifts, 01B) + call mw_saveim (mw, out) + call mw_close (mw) + } + } else if (IM_NDIM(in) == 2) { + call mg_magnify2 (in, out, Memc[interp], btype, bconst, + x1, x2, dx, y1, y2, dy, flux) + if (!envgetb ("nomwcs")) { + mw = mw_openim (in) + scale[1] = xmag + scale[2] = ymag + shifts[1] = 1. - xmag * x1 + shifts[2] = 1. - ymag * y1 + call mw_scale (mw, scale, 03B) + call mw_shift (mw, shifts, 03B) + call mw_saveim (mw, out) + call mw_close (mw) + } + } else { + call imunmap (out) + call imunmap (in) + call xt_delimtemp (Memc[image2], Memc[image3]) + if (logfd != NULL) { + call fprintf (logfd, "\n%s\n") + call pargstr (Memc[image3]) + call fprintf (logfd, + " Cannot magnify image %s to image %s.\n") + call pargstr (Memc[image1]) + call pargstr (Memc[image3]) + call fprintf (logfd, + " Dimensions are greater than 2.\n") + } + call imdelete (Memc[image3]) + next + } + + } then { + if (logfd != NULL) { + call fprintf (logfd, "\n%s\n") + call pargstr (Memc[image3]) + call fprintf (logfd, + " Cannot magnify image %s to image %s.\n") + call pargstr (Memc[image1]) + call pargstr (Memc[image3]) + if (errget (Memc[errmsg], SZ_FNAME) <= 0) + ; + call fprintf (logfd, "%s") + call pargstr (Memc[errmsg]) + } + call imdelete (Memc[image3]) + call imunmap (out) + call imunmap (in) + call xt_delimtemp (Memc[image2], Memc[image3]) + } else { + + if (logfd != NULL) { + call fprintf (logfd, "\n%s\n") + call pargstr (Memc[image3]) + call fprintf (logfd, " Magnify image %s to image %s.\n") + call pargstr (Memc[image1]) + call pargstr (Memc[image3]) + call fprintf (logfd, " Interpolation is %s.\n") + call pargstr (Memc[interp]) + call fprintf (logfd, " Boundary extension is %s.\n") + call pargstr (Memc[boundary]) + if (btype == 1) { + call fprintf (logfd, + " Boundary pixel constant is %g.\n") + call pargr (bconst) + } + call fprintf (logfd, + " Output coordinates in terms of input coordinates:\n") + call fprintf (logfd, + " x1 = %10.4g, x2 = %10.4g, dx = %10.6g\n") + call pargr (x1) + call pargr (x2) + call pargr (dx) + if (IM_NDIM(in) == 2) { + call fprintf (logfd, + " y1 = %10.4g, y2 = %10.4g, dy = %10.6g\n") + call pargr (y1) + call pargr (y2) + call pargr (dy) + } + } + + call imunmap (out) + call imunmap (in) + call xt_delimtemp (Memc[image2], Memc[image3]) + } + + } + + call imtclose (list1) + call imtclose (list2) + call close (logfd) + call sfree (sp) +end + + +define NYOUT2 16 # Number of input lines to use for interpolation +define NMARGIN 3 # Number of edge lines to add for interpolation +define NMARGIN_SPLINE3 16 # Number of edge lines to add for interpolation + + +# MG_MAGNIFY1 -- Magnify the input input image to create the output image. + +procedure mg_magnify1 (in, out, interp, btype, bconst, x1, x2, dx, flux) + +pointer in # pointer to the input image +pointer out # pointer to the output image +char interp[ARB] # Interpolation type +int btype # Boundary extension type +real bconst # Boundary extension constant +real x1, x2 # Starting and ending points of output image +real dx # Pixel interval +int flux # Conserve flux? + +int i, nxin, nxout, nxymargin, itype, nsinc, nincr, col1, col2 +pointer sp, x, z, buf, asi +real xshift +pointer imgs1r(), impl1r() +int asigeti() + +begin + # Set the default values for the output image limits if they are INDEF + # and calculate the number of output pixels. + + if (IS_INDEF (x1)) + x1 = 1. + if (IS_INDEF (x2)) + x2 = IM_LEN (in, 1) + if (x1 > x2) + call error (0, " X1 cannot be greater than X2\n") + + # Set the number of output pixels in the image header. + + nxout = (x2 - x1) / dx + 1 + IM_LEN(out, 1) = nxout + + # Initialize the interpolator. + + call asitype (interp, itype, nsinc, nincr, xshift) + call asisinit (asi, itype, nsinc, nincr, xshift, 0.0) + + # Round the coordinate limits to include the output image coordinate + # limits and the set boundary. + + col1 = x1 + col2 = nint (x2) + if (itype == II_SPLINE3) + nxymargin = NMARGIN_SPLINE3 + else if (itype == II_SINC || itype == II_LSINC) + nxymargin = asigeti (asi, II_ASINSINC) + else if (itype == II_DRIZZLE) + nxymargin = max (nint (dx), NMARGIN) + else + nxymargin = NMARGIN + call mg_setboundary1 (in, col1, col2, btype, bconst, nxymargin) + col1 = col1 - nxymargin + if (col1 <= 0) + col1 = col1 - 1 + col2 = col2 + nxymargin + 1 + + # Allocate memory for the interpolation coordinates. + # Also initialize the image data buffer. + + call smark (sp) + call salloc (x, 2 * nxout, TY_REAL) + + # Set the x interpolation coordinates. The coordinates are relative + # to the boundary extended input image. + + if (itype == II_DRIZZLE) { + do i = 1, nxout { + Memr[x+2*i-2] = x1 + (i - 1.5) * dx - col1 + 1 + Memr[x+2*i-1] = x1 + (i - 0.5) * dx - col1 + 1 + } + } else { + do i = 1, nxout + Memr[x+i-1] = x1 + (i - 1) * dx - col1 + 1 + } + + # Fit the output image. + nxin = col2 - col1 + 1 + buf = imgs1r (in, col1, col2) + call asifit (asi, Memr[buf], nxin) + + # Evaluate the output image pixel values. + z = impl1r (out) + call asivector (asi, Memr[x], Memr[z], nxout) + #if (itype != II_DRIZZLE && flux == YES) + if (flux == YES) + call amulkr (Memr[z], dx, Memr[z], nxout) + + # Free memory and unmap the images. + call asifree (asi) + call sfree (sp) +end + + +# MG_MAGNIFY2 -- Magnify the input input image to create the output image. + +procedure mg_magnify2 (in, out, interp, btype, bconst, x1, x2, dx, y1, y2, + dy, flux) + +pointer in # pointer to the input image +pointer out # pointer to the output image +char interp[ARB] # Interpolation type +int btype # Boundary extension type +real bconst # Boundary extension constant +real x1, y1 # Starting point of output image +real x2, y2 # Ending point of output image +real dx, dy # Pixel interval +int flux # Conserve flux? + +int i, nxin, nxout, nyout, nxymargin, itype, nsinc, nincr +int l1out, l2out, nlout, l1in, l2in, nlin, fstline, lstline +int col1, col2, line1, line2 +real shift +pointer msi +pointer sp, x, y, z, buf + +pointer imps2r() +int msigeti() + +begin + # Set the default values for the output image limits if they are INDEF + # and calculate the number of output pixels. + + if (IS_INDEF (x1)) + x1 = 1. + if (IS_INDEF (x2)) + x2 = IM_LEN (in, 1) + if (IS_INDEF (y1)) + y1 = 1. + if (IS_INDEF (y2)) + y2 = IM_LEN (in, 2) + if (x1 > x2) + call error (0, " X1 cannot be greater than X2\n") + if (y1 > y2) + call error (0, " Y1 cannot be greater than Y2\n") + nxout = (x2 - x1) / dx + 1 + nyout = (y2 - y1) / dy + 1 + + # Set the number of output pixels in the image header. + + IM_LEN(out, 1) = nxout + IM_LEN(out, 2) = nyout + + # Initialize the interpolator. + + call msitype (interp, itype, nsinc, nincr, shift) + call msisinit (msi, itype, nsinc, nincr, nincr, shift, shift, 0.0) + + # Compute the number of margin pixels required + + if (itype == II_BISPLINE3) + nxymargin = NMARGIN_SPLINE3 + else if (itype == II_BISINC || itype == II_BILSINC) + nxymargin = msigeti (msi, II_MSINSINC) + else if (itype == II_BIDRIZZLE) + nxymargin = max (nint (dx), nint(dy), NMARGIN) + else + nxymargin = NMARGIN + + # Round the coordinate limits to include the output image coordinate + # limits and the set boundary. + + col1 = x1 + col2 = nint (x2) + line1 = y1 + line2 = nint (y2) + call mg_setboundary2 (in, col1, col2, line1, line2, btype, bconst, + nxymargin) + + # Compute the input image column limits. + col1 = col1 - nxymargin + if (col1 <= 0) + col1 = col1 - 1 + col2 = col2 + nxymargin + 1 + nxin = col2 - col1 + 1 + + # Allocate memory for the interpolation coordinates. + # Also initialize the image data buffer. + + call smark (sp) + call salloc (x, 2 * nxout, TY_REAL) + call salloc (y, 2 * NYOUT2, TY_REAL) + buf = NULL + fstline = 0 + lstline = 0 + + # Set the x interpolation coordinates which do not change from + # line to line. The coordinates are relative to the boundary + # extended input image. + + if (itype == II_BIDRIZZLE) { + do i = 1, nxout { + Memr[x+2*i-2] = x1 + (i - 1.5) * dx - col1 + 1 + Memr[x+2*i-1] = x1 + (i - 0.5) * dx - col1 + 1 + } + } else { + do i = 1, nxout + Memr[x+i-1] = x1 + (i - 1) * dx - col1 + 1 + } + + # Loop over the image sections. + for (l1out = 1; l1out <= nyout; l1out = l1out + NYOUT2) { + + # Define the range of output lines. + l2out = min (l1out + NYOUT2 - 1, nyout) + nlout = l2out - l1out + 1 + + # Define the corresponding range of input lines. + l1in = y1 + (l1out - 1) * dy - nxymargin + if (l1in <= 0) + l1in = l1in - 1 + l2in = y1 + (l2out - 1) * dy + nxymargin + 1 + nlin = l2in - l1in + 1 + + # Get the apporiate image section and compute the coefficients. + if ((buf == NULL) || (l1in < fstline) || (l2in > lstline)) { + fstline = l1in + lstline = l2in + call mg_bufl2r (in, col1, col2, l1in, l2in, buf) + call msifit (msi, Memr[buf], nxin, nlin, nxin) + } + + # Output the section. + z = imps2r (out, 1, nxout, l1out, l2out) + + # Compute the y values. + if (itype == II_BIDRIZZLE) { + do i = l1out, l2out { + Memr[y+2*(i-l1out)] = y1 + (i - 1.5) * dy - fstline + 1 + Memr[y+2*(i-l1out)+1] = y1 + (i - 0.5) * dy - fstline + 1 + } + } else { + do i = l1out, l2out + Memr[y+i-l1out] = y1 + (i - 1) * dy - fstline + 1 + } + + # Evaluate the interpolant. + call msigrid (msi, Memr[x], Memr[y], Memr[z], nxout, nlout, nxout) + if (flux == YES) + call amulkr (Memr[z], dx * dy, Memr[z], nxout * nlout) + } + + # Free memory and buffers. + call msifree (msi) + call mfree (buf, TY_REAL) + call sfree (sp) +end + + +# MG_BUFL2R -- Maintain buffer of image lines. A new buffer is created when +# the buffer pointer is null or if the number of lines requested is changed. +# The minimum number of image reads is used. + +procedure mg_bufl2r (im, col1, col2, line1, line2, buf) + +pointer im # Image pointer +int col1 # First image column of buffer +int col2 # Last image column of buffer +int line1 # First image line of buffer +int line2 # Last image line of buffer +pointer buf # Buffer + +int i, ncols, nlines, nclast, llast1, llast2, nllast +pointer buf1, buf2 + +pointer imgs2r() + +begin + ncols = col2 - col1 + 1 + nlines = line2 - line1 + 1 + + # If the buffer pointer is undefined then allocate memory for the + # buffer. If the number of columns or lines requested changes + # reallocate the buffer. Initialize the last line values to force + # a full buffer image read. + + if (buf == NULL) { + call malloc (buf, ncols * nlines, TY_REAL) + llast1 = line1 - nlines + llast2 = line2 - nlines + } else if ((nlines != nllast) || (ncols != nclast)) { + call realloc (buf, ncols * nlines, TY_REAL) + llast1 = line1 - nlines + llast2 = line2 - nlines + } + + # Read only the image lines with are different from the last buffer. + + if (line1 < llast1) { + do i = line2, line1, -1 { + if (i > llast1) + buf1 = buf + (i - llast1) * ncols + else + buf1 = imgs2r (im, col1, col2, i, i) + + buf2 = buf + (i - line1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } else if (line2 > llast2) { + do i = line1, line2 { + if (i < llast2) + buf1 = buf + (i - llast1) * ncols + else + buf1 = imgs2r (im, col1, col2, i, i) + + buf2 = buf + (i - line1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } + + # Save the buffer parameters. + + llast1 = line1 + llast2 = line2 + nclast = ncols + nllast = nlines +end + + +# MG_SETBOUNDARY1 -- Set boundary extension for a 1D image. + +procedure mg_setboundary1 (im, col1, col2, btype, bconst, nxymargin) + +pointer im # IMIO pointer +int col1, col2 # Range of columns +int btype # Boundary extension type +real bconst # Constant for constant boundary extension +int nxymargin # Number of margin pixels + +int btypes[5] +int nbndrypix + +data btypes /BT_CONSTANT, BT_NEAREST, BT_REFLECT, BT_WRAP, BT_PROJECT/ + +begin + nbndrypix = 0 + nbndrypix = max (nbndrypix, 1 - col1) + nbndrypix = max (nbndrypix, col2 - IM_LEN(im, 1)) + + call imseti (im, IM_TYBNDRY, btypes[btype]) + call imseti (im, IM_NBNDRYPIX, nbndrypix + nxymargin + 1) + if (btypes[btype] == BT_CONSTANT) + call imsetr (im, IM_BNDRYPIXVAL, bconst) +end + + +# MG_SETBOUNDARY2 -- Set boundary extension for a 2D image. + +procedure mg_setboundary2 (im, col1, col2, line1, line2, btype, bconst, + nxymargin) + +pointer im # IMIO pointer +int col1, col2 # Range of columns +int line1, line2 # Range of lines +int btype # Boundary extension type +real bconst # Constant for constant boundary extension +int nxymargin # Number of margin pixels to allow + +int btypes[5] +int nbndrypix + +data btypes /BT_CONSTANT, BT_NEAREST, BT_REFLECT, BT_WRAP, BT_PROJECT/ + +begin + nbndrypix = 0 + nbndrypix = max (nbndrypix, 1 - col1) + nbndrypix = max (nbndrypix, col2 - IM_LEN(im, 1)) + nbndrypix = max (nbndrypix, 1 - line1) + nbndrypix = max (nbndrypix, line2 - IM_LEN(im, 2)) + + call imseti (im, IM_TYBNDRY, btypes[btype]) + call imseti (im, IM_NBNDRYPIX, nbndrypix + nxymargin + 1) + if (btypes[btype] == BT_CONSTANT) + call imsetr (im, IM_BNDRYPIXVAL, bconst) +end diff --git a/pkg/images/imgeom/src/t_shiftlines.x b/pkg/images/imgeom/src/t_shiftlines.x new file mode 100644 index 00000000..36ce5dec --- /dev/null +++ b/pkg/images/imgeom/src/t_shiftlines.x @@ -0,0 +1,102 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> + +# T_SHIFTLINES -- Shift image lines. +# +# The input and output images are given by image template lists. The +# number of output images must match the number of input images. Image +# sections are allowed in the input images and are ignored in the output +# images. If the input and output image names are the same then the shift +# is performed to a temporary file which then replaces the input image. + + +procedure t_shiftlines() + +char imtlist1[SZ_LINE] # Input image list +char imtlist2[SZ_LINE] # Output image list +real shift # Amount of pixel shift +int boundary # Type of boundary extension +real constant # Constant for boundary extension + +char image1[SZ_FNAME] # Input image name +char image2[SZ_FNAME] # Output image name +char imtemp[SZ_FNAME] # Temporary file + +char str[SZ_LINE], interpstr[SZ_FNAME] +int list1, list2, ishift +pointer im1, im2, mw + +bool fp_equalr(), envgetb() +int clgwrd(), imtopen(), imtgetim(), imtlen() +pointer immap(), mw_openim() +real clgetr() +errchk sh_lines, sh_linesi, mw_openim, mw_shift, mw_saveim, mw_close + +begin + # Get input and output image template lists. + call clgstr ("input", imtlist1, SZ_LINE) + call clgstr ("output", imtlist2, SZ_LINE) + + # Get the shift, interpolation type, and boundary condition. + shift = clgetr ("shift") + call clgstr ("interp_type", interpstr, SZ_LINE) + boundary = clgwrd ("boundary_type", str, SZ_LINE, + ",constant,nearest,reflect,wrap,") + constant = clgetr ("constant") + + # Expand the input and output image lists. + list1 = imtopen (imtlist1) + list2 = imtopen (imtlist2) + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + call error (0, "Number of input and output images not the same") + } + + ishift = shift + + # Do each set of input/output images. + while ((imtgetim (list1, image1, SZ_FNAME) != EOF) && + (imtgetim (list2, image2, SZ_FNAME) != EOF)) { + + call xt_mkimtemp (image1, image2, imtemp, SZ_FNAME) + + im1 = immap (image1, READ_ONLY, 0) + im2 = immap (image2, NEW_COPY, im1) + + # Shift the image. + iferr { + + if (fp_equalr (shift, real (ishift))) + call sh_linesi (im1, im2, ishift, boundary, constant) + else + call sh_lines (im1, im2, shift, boundary, constant, + interpstr) + + # Update the image WCS to reflect the shift. + if (!envgetb ("nomwcs")) { + mw = mw_openim (im1) + call mw_shift (mw, shift, 1B) + call mw_saveim (mw, im2) + call mw_close (mw) + } + + } then { + call eprintf ("Error shifting image: %s\n") + call pargstr (image1) + call erract (EA_WARN) + } + + # Unmap images. + call imunmap (im2) + call imunmap (im1) + + # If in place operation replace the input image with the temporary + # image. + call xt_delimtemp (image2, imtemp) + } + + call imtclose (list1) + call imtclose (list2) +end |