aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/immatch/src/psfmatch
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/immatch/src/psfmatch')
-rw-r--r--pkg/images/immatch/src/psfmatch/mkpkg21
-rw-r--r--pkg/images/immatch/src/psfmatch/psfmatch.h274
-rw-r--r--pkg/images/immatch/src/psfmatch/psfmatch.key50
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpbckgrd.x70
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpcolon.x501
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpconvolve.x106
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpfft.x443
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpfilter.x502
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpisfm.x556
-rw-r--r--pkg/images/immatch/src/psfmatch/rgppars.x124
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpregions.x464
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpsfm.x815
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpshow.x116
-rw-r--r--pkg/images/immatch/src/psfmatch/rgptools.x641
-rw-r--r--pkg/images/immatch/src/psfmatch/t_psfmatch.x365
15 files changed, 5048 insertions, 0 deletions
diff --git a/pkg/images/immatch/src/psfmatch/mkpkg b/pkg/images/immatch/src/psfmatch/mkpkg
new file mode 100644
index 00000000..da3951dc
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/mkpkg
@@ -0,0 +1,21 @@
+# Make the PSFMATCH task
+
+$checkout libpkg.a ../../../
+$update libpkg.a
+$checkin libpkg.a ../../../
+$exit
+
+libpkg.a:
+ rgpbckgrd.x <math.h> <math/gsurfit.h> "psfmatch.h"
+ rgpcolon.x <imhdr.h> <imset.h> <error.h> "psfmatch.h"
+ rgpconvolve.x <error.h> <imhdr.h> <imset.h>
+ rgpisfm.x <imhdr.h> <gset.h> <ctype.h> "psfmatch.h"
+ rgpfft.x
+ rgpfilter.x <math.h>
+ rgppars.x "psfmatch.h"
+ rgpregions.x <imhdr.h> <fset.h> "psfmatch.h"
+ rgpsfm.x <imhdr.h> <math/gsurfit.h> "psfmatch.h"
+ rgpshow.x "psfmatch.h"
+ rgptools.x "psfmatch.h"
+ t_psfmatch.x <fset.h> <imhdr.h> "psfmatch.h"
+ ;
diff --git a/pkg/images/immatch/src/psfmatch/psfmatch.h b/pkg/images/immatch/src/psfmatch/psfmatch.h
new file mode 100644
index 00000000..c6b7d563
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/psfmatch.h
@@ -0,0 +1,274 @@
+# Header file for PSFMATCH
+
+define LEN_PSFSTRUCT (45 + 12 * SZ_FNAME + 12)
+
+# Define the psf fitting structure
+
+define PM_RC1 Memi[$1] # pointer to first column of region
+define PM_RC2 Memi[$1+1] # pointer to last column of region
+define PM_RL1 Memi[$1+2] # pointer to first line of region
+define PM_RL2 Memi[$1+3] # pointer to last line of region
+define PM_RZERO Memi[$1+4] # pointer to zero point of ref regions
+define PM_RXSLOPE Memi[$1+5] # pointer to x slopes of ref regions
+define PM_RYSLOPE Memi[$1+6] # pointer to y slopes of ref regions
+define PM_NREGIONS Memi[$1+7] # total number of regions
+define PM_CNREGION Memi[$1+8] # the current region
+
+define PM_CENTER Memi[$1+9] # the the psf objects
+define PM_BACKGRD Memi[$1+10] # type of background subtraction
+define PM_BVALUER Memr[P2R($1+11)] # reference background value
+define PM_BVALUE Memr[P2R($1+12)] # image background value
+define PM_LOREJECT Memr[P2R($1+13)] # low side rejection
+define PM_HIREJECT Memr[P2R($1+14)] # high side rejection
+define PM_APODIZE Memr[P2R($1+15)] # fraction of region to be apodized
+
+define PM_CONVOLUTION Memi[$1+16] # the convolution type
+define PM_DNX Memi[$1+17] # x dimension of kernel
+define PM_DNY Memi[$1+18] # y dimension of kernel
+define PM_PNX Memi[$1+19] # x dimension of user kernel
+define PM_PNY Memi[$1+20] # y dimension of user kernel
+define PM_KNX Memi[$1+21] # x size of kernel
+define PM_KNY Memi[$1+22] # x size of kernel
+
+define PM_POWER Memi[$1+23] # save power spectrum of kernel ?
+
+define PM_UFLUXRATIO Memr[P2R($1+24)] # the user ref / input flux ratio
+define PM_FLUXRATIO Memr[P2R($1+25)] # ref / input flux ratio
+define PM_FILTER Memi[$1+26] # background filtering
+define PM_SXINNER Memr[P2R($1+27)] # inner radius for cosine bell
+define PM_SXOUTER Memr[P2R($1+28)] # outer radius for cosine bell
+define PM_SYINNER Memr[P2R($1+29)] # inner radius for cosine bell
+define PM_SYOUTER Memr[P2R($1+30)] # outer radius for cosine bell
+define PM_RADSYM Memi[$1+31] # radial symmetry in convolution
+define PM_THRESHOLD Memr[P2R($1+32)] # threshold in divisor for model
+
+define PM_NORMFACTOR Memr[P2R($1+34)] # the normalization factor
+
+#define PM_PRATIO Memr[P2R($1+24)] # power ration threshold
+#define PM_XSHIFTS Memi[$1+26] # pointer to x shifts
+#define PM_YSHIFTS Memi[$1+27] # pointer to y shifts
+
+define PM_REFFFT Memi[$1+35] # pointer to reference fft
+define PM_IMFFT Memi[$1+36] # pointer to image fft
+define PM_FFT Memi[$1+37] # pointer to unfiltered fft
+define PM_CONV Memi[$1+38] # pointer to kernel
+define PM_ASFFT Memi[$1+39] # pointer to power spectrum
+define PM_NXFFT Memi[$1+40] # x dimension of FFT
+define PM_NYFFT Memi[$1+41] # y dimension of FFT
+
+define PM_BSTRING Memc[P2C($1+42)] # background string
+define PM_CSTRING Memc[P2C($1+42+SZ_FNAME+1)] # convolution string
+define PM_FSTRING Memc[P2C($1+42+2*SZ_FNAME+2)] # convolution string
+
+define PM_IMAGE Memc[P2C($1+42+4*SZ_FNAME+4)] # input image
+define PM_REFIMAGE Memc[P2C($1+42+5*SZ_FNAME+5)] # reference image
+define PM_PSFDATA Memc[P2C($1+42+6*SZ_FNAME+6)] # psf data
+define PM_PSFIMAGE Memc[P2C($1+42+7*SZ_FNAME+7)] # psf image if any
+define PM_OBJLIST Memc[P2C($1+42+8*SZ_FNAME+8)] # object list if any
+define PM_KERNEL Memc[P2C($1+42+9*SZ_FNAME+9)] # kernel image
+define PM_OUTIMAGE Memc[P2C($1+42+10*SZ_FNAME+10)] # output convolved image
+
+# Define the paramerter ids
+
+define RC1 1
+define RC2 2
+define RL1 3
+define RL2 4
+define RZERO 5
+define RXSLOPE 6
+define RYSLOPE 7
+define NREGIONS 8
+define CNREGION 9
+
+define CENTER 10
+define BACKGRD 11
+define BVALUER 12
+define BVALUE 13
+define LOREJECT 15
+define HIREJECT 16
+define APODIZE 17
+
+define CONVOLUTION 18
+define DNX 19
+define DNY 20
+define PNX 21
+define PNY 22
+define KNX 23
+define KNY 24
+define POWER 25
+
+#define XSHIFTS 20
+#define YSHIFTS 21
+
+define REFFFT 26
+define IMFFT 27
+define FFT 28
+define CONV 29
+define ASFFT 30
+define NXFFT 31
+define NYFFT 32
+
+define UFLUXRATIO 33
+define FLUXRATIO 34
+define FILTER 35
+define SXINNER 36
+define SXOUTER 37
+define SYINNER 38
+define SYOUTER 39
+define RADSYM 40
+define THRESHOLD 41
+
+define NORMFACTOR 43
+
+#define PRATIO 34
+
+define BSTRING 44
+define CSTRING 45
+define FSTRING 46
+
+define REFIMAGE 48
+define IMAGE 49
+define PSFDATA 50
+define PSFIMAGE 51
+define OBJLIST 52
+define KERNEL 53
+define OUTIMAGE 54
+
+# Define the default parameter values
+
+define DEF_CENTER YES
+define DEF_BACKGRD PM_BMEDIAN
+define DEF_LOREJECT INDEFR
+define DEF_HIREJECT INDEFR
+
+define DEF_CONVOLUTION PM_CONIMAGE
+define DEF_DNX 63
+define DEF_DNY 63
+define DEF_PNX 31
+define DEF_PNY 31
+define DEF_POWER NO
+
+define DEF_FILTER PM_FREPLACE
+define DEF_SXINNER INDEFR
+define DEF_SXOUTER INDEFR
+define DEF_SYINNER INDEFR
+define DEF_SYOUTER INDEFR
+define DEF_RADSYM NO
+define DEF_THRESHOLD 0.0
+
+#define DEF_PRATIO 0.0
+
+define DEF_NORMFACTOR 1.0
+define DEF_UFLUXRATIO INDEFR
+
+# Define the background fitting techniques
+
+define PM_BNONE 1
+define PM_BMEAN 2
+define PM_BMEDIAN 3
+define PM_BSLOPE 4
+define PM_BNUMBER 5
+
+define PM_BTYPES "|none|mean|median|plane|"
+
+# Define the convolution computation options
+
+define PM_CONIMAGE 1
+define PM_CONPSF 2
+define PM_CONKERNEL 3
+
+define PM_CTYPES "|image|psf|kernel|"
+
+# Define the filtering options
+
+define PM_FNONE 1
+define PM_FCOSBELL 2
+define PM_FREPLACE 3
+define PM_FMODEL 4
+
+define PM_FTYPES "|none|cosbell|replace|model|"
+
+# Define the normalization options
+
+define PM_UNIT 1
+define PM_RATIO 2
+define PM_NUMBER 3
+
+define PM_NTYPES "|unit|ratio|"
+
+# Miscellaneous
+
+define MAX_NREGIONS 100
+
+# Commands
+
+define PMCMDS "|input|reference|psfdata|psfimage|kernel|output|dnx|dny|\
+pnx|pny|center|background|loreject|hireject|apodize|convolution|fluxratio|\
+filter|sx1|sx2|sy1|sy2|radsym|threshold|normfactor|show|mark|"
+
+define PMCMD_IMAGE 1
+define PMCMD_REFIMAGE 2
+define PMCMD_PSFDATA 3
+define PMCMD_PSFIMAGE 4
+define PMCMD_KERNEL 5
+define PMCMD_OUTIMAGE 6
+
+define PMCMD_DNX 7
+define PMCMD_DNY 8
+define PMCMD_PNX 9
+define PMCMD_PNY 10
+
+define PMCMD_CENTER 11
+define PMCMD_BACKGRD 12
+define PMCMD_LOREJECT 13
+define PMCMD_HIREJECT 14
+define PMCMD_APODIZE 15
+
+define PMCMD_CONVOLUTION 16
+define PMCMD_UFLUXRATIO 17
+define PMCMD_FILTER 18
+define PMCMD_SXINNER 19
+define PMCMD_SXOUTER 20
+define PMCMD_SYINNER 21
+define PMCMD_SYOUTER 22
+define PMCMD_RADSYM 23
+define PMCMD_THRESHOLD 24
+
+define PMCMD_NORMFACTOR 25
+
+define PMCMD_SHOW 26
+define PMCMD_MARK 27
+
+# Keywords
+
+define KY_IMAGE "input"
+define KY_REFIMAGE "reference"
+define KY_PSFDATA "psfdata"
+define KY_PSFIMAGE "psfimage"
+define KY_KERNEL "kernel"
+define KY_OUTIMAGE "output"
+
+define KY_DNX "dnx"
+define KY_DNY "dny"
+define KY_PNX "pnx"
+define KY_PNY "pny"
+
+define KY_CENTER "center"
+define KY_BACKGRD "background"
+define KY_LOREJECT "loreject"
+define KY_HIREJECT "hireject"
+define KY_APODIZE "apodize"
+
+define KY_CONVOLUTION "convolution"
+
+define KY_UFLUXRATIO "fluxratio"
+define KY_FILTER "filter"
+define KY_SXINNER "sx1"
+define KY_SXOUTER "sx2"
+define KY_SYINNER "sy1"
+define KY_SYOUTER "sy2"
+define KY_RADSYM "radsym"
+define KY_THRESHOLD "threshold"
+
+define KY_NORMFACTOR "normfactor"
+
diff --git a/pkg/images/immatch/src/psfmatch/psfmatch.key b/pkg/images/immatch/src/psfmatch/psfmatch.key
new file mode 100644
index 00000000..57ef3b2e
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/psfmatch.key
@@ -0,0 +1,50 @@
+ Interactive Keystroke Commands
+
+
+? Print help
+: Colon commands
+k Draw a contour plot of the psf matching kernel
+p Draw a contour plot of the psf matching kernel power spectrum
+x Draw a column plot of the psf matching kernel / power spectrum
+y Draw a line plot of the psf matching kernel / power spectrum
+r Redraw the current plot
+f Recompute the psf matching kernel
+w Update the task parameters
+q Exit
+
+
+ Colon Commands
+
+
+:mark [file] Mark objects on the display
+:show Show current values of the parameters
+
+
+ Show/Set Parameters
+
+:input [string] Show/set the current input image name
+:reference [string] Show/set the current reference image/psf name
+:psf [file/string] Show/set the objects/input psf list
+:psfimage [string] Show/set the current input psf name
+:kernel [string] Show/set the current psf matching kernel name
+:output [string] Show/set the current output image name
+
+:dnx [value] Show/set x width of data region(s) to extract
+:dny [value] Show/set y width of data region(s) to extract
+:pnx [value] Show/set x width of psf matching kernel
+:pny [value] Show/set y width of psf matching kernel
+:center [yes/no] Show/set the centering switch
+:background [string] Show/set the background fitting function
+:loreject [value] Show/set low side k-sigma rejection parameter
+:hireject [value] Show/set high side k-sigma rejection parameter
+:apodize [value] Show/set percent of endpoints to apodize
+
+:filter [string] Show/set the filtering algorithm
+:fluxratio [value] Show/set the reference/input psf flux ratio
+:sx1 [value] Show/set inner x frequency for cosbell filter
+:sx2 [value] Show/set outer x frequency for cosbell filter
+:sy1 [value] Show/set inner y frequency for cosbell filter
+:sy2 [value] Show/set outer y frequency for cosbell filter
+:radsym [yes/no] Show/set radial symmetry for cosbell filter
+:threshold [value] Show/set %threshold for replace/modeling filter
+:normfactor [value] Show/set the kernel normalization factor
diff --git a/pkg/images/immatch/src/psfmatch/rgpbckgrd.x b/pkg/images/immatch/src/psfmatch/rgpbckgrd.x
new file mode 100644
index 00000000..1670b943
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpbckgrd.x
@@ -0,0 +1,70 @@
+include <math.h>
+include <math/gsurfit.h>
+include "psfmatch.h"
+
+# RG_PSCALE -- Compute the background offset and x and y slope.
+
+procedure rg_pscale (pm, data, npts, nx, ny, pnx, pny, offset, coeff)
+
+pointer pm #I pointer to the psfmatch structure
+real data[ARB] #I the input data
+int npts #I the number of points
+int nx, ny #I the dimensions of the original subraster
+int pnx, pny #I the dimensions of the data region
+real offset #I the input offset
+real coeff[ARB] #O the output coefficients
+
+int wxborder, wyborder
+pointer gs
+real loreject, hireject, zero
+int rg_pstati(), rg_znsum(), rg_znmedian(), rg_slope()
+real rg_pstatr()
+
+begin
+ loreject = rg_pstatr (pm, LOREJECT)
+ hireject = rg_pstatr (pm, HIREJECT)
+
+ switch (rg_pstati (pm, BACKGRD)) {
+ case PM_BNONE:
+ coeff[1] = 0.0
+ coeff[2] = 0.0
+ coeff[3] = 0.0
+ case PM_BNUMBER:
+ coeff[1] = offset
+ coeff[2] = 0.0
+ coeff[3] = 0.0
+ case PM_BMEAN:
+ if (rg_znsum (data, npts, zero, loreject, hireject) <= 0)
+ zero = 0.0
+ coeff[1] = zero
+ coeff[2] = 0.0
+ coeff[3] = 0.0
+ case PM_BMEDIAN:
+ if (rg_znmedian (data, npts, zero, loreject, hireject) <= 0)
+ zero = 0.0
+ coeff[1] = zero
+ coeff[2] = 0.0
+ coeff[3] = 0.0
+ case PM_BSLOPE:
+ call gsinit (gs, GS_POLYNOMIAL, 2, 2, GS_XNONE, 1.0, real (nx), 1.0,
+ real (ny))
+ wxborder = (nx - pnx) / 2
+ wyborder = (ny - pny) / 2
+ if (rg_slope (gs, data, npts, nx, ny, wxborder, wyborder, loreject,
+ hireject) == ERR) {
+ coeff[1] = 0.0
+ coeff[2] = 0.0
+ coeff[3] = 0.0
+ } else {
+ call gssave (gs, coeff)
+ coeff[1] = coeff[GS_SAVECOEFF+1]
+ coeff[2] = coeff[GS_SAVECOEFF+2]
+ coeff[3] = coeff[GS_SAVECOEFF+3]
+ }
+ call gsfree (gs)
+ default:
+ coeff[1] = 0.0
+ coeff[2] = 0.0
+ coeff[3] = 0.0
+ }
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgpcolon.x b/pkg/images/immatch/src/psfmatch/rgpcolon.x
new file mode 100644
index 00000000..8eefb22d
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpcolon.x
@@ -0,0 +1,501 @@
+include <imhdr.h>
+include <imset.h>
+include <error.h>
+include "psfmatch.h"
+
+# RG_PCOLON -- Show/set the psfmatch task algorithm parameters.
+
+procedure rg_pcolon (gd, pm, imr, reglist, impsf, im1, imk, imfourier, im2,
+ cmdstr, newref, newdata, newfourier, newfilter)
+
+pointer gd #I pointer to the graphics stream
+pointer pm #I pointer to psfmatch structure
+pointer imr #I pointer to the reference image
+int reglist #I the regions / psf list descriptor
+pointer impsf #I pointer to the regions list
+pointer im1 #I pointer to the input image
+pointer imk #I pointer to kernel image
+pointer imfourier #I pointer to fourier spectrum image
+pointer im2 #I pointer to the output image
+char cmdstr[ARB] #I command string
+int newref #I/O new reference image
+int newdata #I/O new input image
+int newfourier #I/O new FFT
+int newfilter #I/O new filter
+
+bool bval
+int ncmd, ival, stat, fd, ip
+pointer sp, cmd, str
+real rval
+bool itob()
+bool streq()
+int strdic(), nscan(), rg_pstati(), btoi(), rg_pregions(), fntopnb()
+int access(), rg_pmkregions(), open(), ctor()
+pointer immap()
+real rg_pstatr()
+errchk immap(), fntopnb()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Get the command.
+ call sscan (cmdstr)
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (Memc[cmd] == EOS) {
+ call sfree (sp)
+ return
+ }
+
+ # Process the command.
+ ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PMCMDS)
+ switch (ncmd) {
+ case PMCMD_REFIMAGE:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ if (imr == NULL || Memc[cmd] == EOS || streq (Memc[cmd],
+ Memc[str])) {
+ call printf ("%s: %s\n")
+ call pargstr (KY_REFIMAGE)
+ call pargstr (Memc[str])
+ } else {
+ if (imr != NULL) {
+ call imunmap (imr)
+ imr = NULL
+ }
+ iferr {
+ imr = immap (Memc[cmd], READ_ONLY, 0)
+ } then {
+ call erract (EA_WARN)
+ imr = immap (Memc[str], READ_ONLY, 0)
+ } else if (IM_NDIM(imr) > 2 || IM_NDIM(imr) != IM_NDIM(im1)) {
+ call printf (
+ "Reference image has the wrong number of dimensions\n")
+ call imunmap (imr)
+ imr = immap (Memc[str], READ_ONLY, 0)
+ } else {
+ call rg_psets (pm, REFIMAGE, Memc[cmd])
+ newref = YES; newdata = YES
+ newfourier = YES; newfilter = YES
+ }
+ }
+
+ case PMCMD_IMAGE:
+
+ call gargwrd (Memc[cmd], SZ_LINE)
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) {
+ call printf ("%s: %s\n")
+ call pargstr (KY_IMAGE)
+ call pargstr (Memc[str])
+ } else {
+ if (im1 != NULL) {
+ call imunmap (im1)
+ im1 = NULL
+ }
+ iferr {
+ im1 = immap (Memc[cmd], READ_ONLY, 0)
+ call imseti (im1, IM_TYBNDRY, BT_NEAREST)
+ if (IM_NDIM(im1) == 1)
+ call imseti (im1, IM_NBNDRYPIX, IM_LEN(im1,1))
+ else
+ call imseti (im1, IM_NBNDRYPIX,
+ max (IM_LEN(im1,1), IM_LEN(im1,2)))
+ } then {
+ call erract (EA_WARN)
+ im1 = immap (Memc[str], READ_ONLY, 0)
+ call imseti (im1, IM_TYBNDRY, BT_NEAREST)
+ if (IM_NDIM(im1) == 1)
+ call imseti (im1, IM_NBNDRYPIX, IM_LEN(im1,1))
+ else
+ call imseti (im1, IM_NBNDRYPIX,
+ max (IM_LEN(im1,1), IM_LEN(im1,2)))
+ } else if (IM_NDIM(im1) > 2 || IM_NDIM(im1) != IM_NDIM(imr)) {
+ call printf (
+ "Reference image has the wrong number of dimensions\n")
+ call imunmap (im1)
+ im1 = immap (Memc[str], READ_ONLY, 0)
+ call imseti (im1, IM_TYBNDRY, BT_NEAREST)
+ if (IM_NDIM(im1) == 1)
+ call imseti (im1, IM_NBNDRYPIX, IM_LEN(im1,1))
+ else
+ call imseti (im1, IM_NBNDRYPIX,
+ max (IM_LEN(im1,1), IM_LEN(im1,2)))
+ } else {
+ call rg_psets (pm, IMAGE, Memc[cmd])
+ newdata = YES; newref = YES
+ newfourier = YES; newfilter = YES
+ }
+ }
+
+ case PMCMD_PSFDATA:
+
+ call gargwrd (Memc[cmd], SZ_LINE)
+ call rg_pstats (pm, PSFDATA, Memc[str], SZ_FNAME)
+ if (reglist == NULL || nscan() == 1 || (streq (Memc[cmd],
+ Memc[str]) && Memc[cmd] != EOS)) {
+ call printf ("%s [string/file]: %s\n")
+ call pargstr (KY_PSFDATA)
+ call pargstr (Memc[str])
+ } else if (rg_pstati(pm, CONVOLUTION) == PM_CONIMAGE) {
+ call fntclsb (reglist)
+ iferr {
+ reglist = fntopnb (Memc[cmd], NO)
+ } then {
+ reglist = fntopnb (Memc[str], NO)
+ } else {
+ if (rg_pregions (reglist, imr, pm, 1, NO) > 0)
+ ;
+ call rg_psets (pm, PSFDATA, Memc[cmd])
+ newdata = YES; newref = YES
+ newfourier = YES; newfilter = YES
+ }
+ }
+
+ case PMCMD_PSFIMAGE:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ call rg_pstats (pm, PSFIMAGE, Memc[str], SZ_FNAME)
+ if (impsf == NULL || Memc[cmd] == EOS || streq (Memc[cmd],
+ Memc[str])) {
+ call printf ("%s: %s\n")
+ call pargstr (KY_PSFIMAGE)
+ call pargstr (Memc[str])
+ } else {
+ if (impsf != NULL) {
+ call imunmap (impsf)
+ impsf = NULL
+ }
+ iferr {
+ impsf = immap (Memc[cmd], READ_ONLY, 0)
+ } then {
+ call erract (EA_WARN)
+ impsf = immap (Memc[str], READ_ONLY, 0)
+ } else if (IM_NDIM(impsf) > 2 || IM_NDIM(impsf) !=
+ IM_NDIM(imr)) {
+ call printf (
+ "PSF image has the wrong number of dimensions\n")
+ call imunmap (impsf)
+ impsf = immap (Memc[str], READ_ONLY, 0)
+ } else {
+ call rg_psets (pm, PSFIMAGE, Memc[cmd])
+ newref = YES; newdata = YES
+ newfourier = YES; newfilter = YES
+ }
+ }
+
+ case PMCMD_KERNEL:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ call rg_pstats (pm, KERNEL, Memc[str], SZ_FNAME)
+ if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) {
+ call printf ("%s: %s\n")
+ call pargstr (KY_KERNEL)
+ call pargstr (Memc[str])
+ } else {
+ if (imk != NULL) {
+ call imunmap (imk)
+ call imdelete (Memc[str])
+ imk = NULL
+ }
+ iferr {
+ imk = immap (Memc[cmd], NEW_IMAGE, 0)
+ } then {
+ call erract (EA_WARN)
+ imk = NULL
+ call rg_psets (pm, KERNEL, "")
+ } else
+ call rg_psets (pm, KERNEL, Memc[cmd])
+ }
+
+
+ case PMCMD_OUTIMAGE:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ call rg_pstats (pm, OUTIMAGE, Memc[str], SZ_FNAME)
+ if (im2 == NULL || Memc[cmd] == EOS || streq (Memc[cmd],
+ Memc[str])) {
+ call printf ("%s: %s\n")
+ call pargstr (KY_OUTIMAGE)
+ call pargstr (Memc[str])
+ } else {
+ if (im2 != NULL) {
+ call imunmap (im2)
+ im2 = NULL
+ }
+ iferr {
+ im2 = immap (Memc[cmd], NEW_COPY, im1)
+ } then {
+ call erract (EA_WARN)
+ im2 = immap (Memc[str], NEW_COPY, im1)
+ } else {
+ call rg_psets (pm, OUTIMAGE, Memc[cmd])
+ }
+ }
+
+ case PMCMD_DNX:
+ call gargi (ival)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_DNX)
+ call pargi (rg_pstati (pm, DNX))
+ } else {
+ if (mod (ival, 2) == 0)
+ ival = ival + 1
+ call rg_pseti (pm, DNX, ival)
+ newref = YES; newdata = YES; newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_DNY:
+ call gargi (ival)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_DNY)
+ call pargi (rg_pstati (pm, DNY))
+ } else {
+ if (mod (ival, 2) == 0)
+ ival = ival + 1
+ call rg_pseti (pm, DNY, ival)
+ newref = YES; newdata = YES; newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_PNX:
+ call gargi (ival)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_PNX)
+ call pargi (rg_pstati (pm, PNX))
+ } else {
+ if (mod (ival, 2) == 0)
+ ival = ival + 1
+ call rg_pseti (pm, PNX, min (ival, rg_pstati (pm, DNX)))
+ newref = YES; newdata = YES; newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_PNY:
+ call gargi (ival)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_PNY)
+ call pargi (rg_pstati (pm, PNY))
+ } else {
+ if (mod (ival, 2) == 0)
+ ival = ival + 1
+ call rg_pseti (pm, PNY, min (ival, rg_pstati(pm, DNY)))
+ newref = YES; newdata = YES; newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_CENTER:
+ call gargb (bval)
+ if (nscan() == 1) {
+ call printf ("%s = %b\n")
+ call pargstr (KY_CENTER)
+ call pargb (itob (rg_pstati (pm, CENTER)))
+ } else {
+ call rg_pseti (pm, CENTER, btoi (bval))
+ newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_BACKGRD:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (Memc[cmd] == EOS) {
+ call rg_pstats (pm, BSTRING, Memc[str], SZ_FNAME)
+ call printf ("%s: %s\n")
+ call pargstr (KY_BACKGRD)
+ call pargstr (Memc[str])
+ } else {
+ stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PM_BTYPES)
+ ip = 1
+ if (stat > 0) {
+ call rg_pseti (pm, BACKGRD, stat)
+ call rg_psets (pm, BSTRING, Memc[cmd])
+ newfourier = YES; newfilter = YES
+ } else if (ctor (str, ip, rval) > 0) {
+ call rg_psetr (pm, BVALUE, rval)
+ if (ctor (str, ip, rval) > 0) {
+ call rg_psetr (pm, BVALUER, rval)
+ call strcpy (str, PM_BSTRING(pm), SZ_FNAME)
+ call rg_pseti (pm, BACKGRD, PM_NUMBER)
+ } else {
+ call rg_psetr (pm, BVALUE, 0.0)
+ call rg_psetr (pm, BVALUER, 0.0)
+ }
+ }
+ }
+
+ case PMCMD_LOREJECT:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_LOREJECT)
+ call pargr (rg_pstatr (pm, LOREJECT))
+ } else {
+ call rg_psetr (pm, LOREJECT, rval)
+ newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_HIREJECT:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_HIREJECT)
+ call pargr (rg_pstatr (pm, HIREJECT))
+ } else {
+ call rg_psetr (pm, HIREJECT, rval)
+ newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_APODIZE:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_APODIZE)
+ call pargr (rg_pstatr (pm, APODIZE))
+ } else {
+ call rg_psetr (pm, APODIZE, rval)
+ newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_CONVOLUTION:
+ if (Memc[cmd] == EOS) {
+ call rg_pstats (pm, CSTRING, Memc[str], SZ_LINE)
+ call printf ("%s: %s\n")
+ call pargstr (KY_CONVOLUTION)
+ call pargstr (Memc[str])
+ }
+
+ case PMCMD_UFLUXRATIO:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_UFLUXRATIO)
+ call pargr (rg_pstatr (pm, UFLUXRATIO))
+ } else {
+ call rg_psetr (pm, UFLUXRATIO, rval)
+ newfourier = YES; newfilter = YES
+ }
+
+ case PMCMD_FILTER:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (Memc[cmd] == EOS) {
+ call rg_pstats (pm, FSTRING, Memc[str], SZ_LINE)
+ call printf ("%s: %s\n")
+ call pargstr (KY_FILTER)
+ call pargstr (Memc[str])
+ } else {
+ stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PM_FTYPES)
+ if (stat > 0) {
+ call rg_pseti (pm, FILTER, stat)
+ call rg_psets (pm, FSTRING, Memc[cmd])
+ }
+ newfilter = YES
+ }
+
+ case PMCMD_SXINNER:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_SXINNER)
+ call pargr (rg_pstatr (pm, SXINNER))
+ } else {
+ call rg_psetr (pm, SXINNER, rval)
+ newfilter = YES
+ }
+
+ case PMCMD_SXOUTER:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_SXOUTER)
+ call pargr (rg_pstatr (pm, SXOUTER))
+ } else {
+ call rg_psetr (pm, SXOUTER, rval)
+ newfilter = YES
+ }
+
+ case PMCMD_SYINNER:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_SYINNER)
+ call pargr (rg_pstatr (pm, SYINNER))
+ } else {
+ call rg_psetr (pm, SYINNER, rval)
+ newfilter = YES
+ }
+
+ case PMCMD_SYOUTER:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_SYOUTER)
+ call pargr (rg_pstatr (pm, SYOUTER))
+ } else {
+ call rg_psetr (pm, SYOUTER, rval)
+ newfilter = YES
+ }
+
+ case PMCMD_RADSYM:
+ call gargb (bval)
+ if (nscan() == 1) {
+ call printf ("%s = %b\n")
+ call pargstr (KY_RADSYM)
+ call pargb (itob (rg_pstati (pm, RADSYM)))
+ } else {
+ call rg_pseti (pm, RADSYM, btoi (bval))
+ newfilter = YES
+ }
+
+ case PMCMD_THRESHOLD:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_THRESHOLD)
+ call pargr (rg_pstatr (pm, THRESHOLD))
+ } else {
+ call rg_psetr (pm, THRESHOLD, rval)
+ newfilter = YES
+ }
+
+ case PMCMD_NORMFACTOR:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("%s = %g\n")
+ call pargstr (KY_NORMFACTOR)
+ call pargr (rg_pstatr (pm, NORMFACTOR))
+ } else {
+ call rg_psetr (pm, NORMFACTOR, rval)
+ newfilter = YES
+ }
+
+ case PMCMD_SHOW:
+ call gdeactivate (gd, 0)
+ call rg_pshow (pm)
+ call greactivate (gd, 0)
+
+ case PMCMD_MARK:
+ call gdeactivate (gd, 0)
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (Memc[cmd] == EOS) {
+ fd = NULL
+ } else if (access (Memc[cmd], 0, 0) == YES) {
+ call printf ("Warning: file %s already exists\n")
+ call pargstr (Memc[cmd])
+ fd = NULL
+ } else {
+ fd = open (Memc[cmd], NEW_FILE, TEXT_FILE)
+ }
+ call printf ("\n")
+ if (rg_pmkregions (fd, imr, pm, 1, MAX_NREGIONS) <= 0)
+ call printf ("The regions list is empty\n")
+ newdata = YES; newref = YES
+ newfourier = YES; newfilter = YES
+ call printf ("\n")
+ if (fd != NULL)
+ call close (fd)
+ call greactivate (gd, 0)
+
+ default:
+ call printf ("Unknown or ambiguous colon command\7\n")
+ }
+
+ call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgpconvolve.x b/pkg/images/immatch/src/psfmatch/rgpconvolve.x
new file mode 100644
index 00000000..6b516a95
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpconvolve.x
@@ -0,0 +1,106 @@
+include <error.h>
+include <imhdr.h>
+include <imset.h>
+
+# RG_PCONVOLVE -- Convolve an image with an nxk by nyk kernel. The kernel
+# dimensions are assumed to be odd.
+
+procedure rg_pconvolve (im1, im2, kernel, nxk, nyk, boundary, constant)
+
+pointer im1 # pointer to the input image
+pointer im2 # pointer to the output image
+real kernel[nxk,nyk] # the convolution kernel
+int nxk, nyk # dimensions of the kernel
+int boundary # type of boundary extension
+real constant # constant for constant boundary extension
+
+int i, ncols, nlines, col1, col2, nincols, inline, outline
+pointer sp, lineptrs, linebuf, outbuf, nkern
+pointer imgs2r(), impl2r()
+errchk imgs2r, impl2r
+
+begin
+ # Set up an array of line pointers.
+ call smark (sp)
+ call salloc (lineptrs, nyk, TY_POINTER)
+ call salloc (nkern, nxk * nyk, TY_REAL)
+
+ # Set the number of image buffers.
+ call imseti (im1, IM_NBUFS, nyk)
+
+ # Set the input image boundary conditions.
+ call imseti (im1, IM_TYBNDRY, boundary)
+ call imseti (im1, IM_NBNDRYPIX, max (nxk / 2 + 1, nyk / 2 + 1))
+ if (boundary == BT_CONSTANT)
+ call imsetr (im1, IM_BNDRYPIXVAL, constant)
+
+ # Define the number of output image lines and columns.
+ ncols = IM_LEN(im2,1)
+ if (IM_NDIM(im2) == 1)
+ nlines = 1
+ else
+ nlines = IM_LEN(im2,2)
+
+ # Set the input image column limits.
+ col1 = 1 - nxk / 2
+ col2 = IM_LEN(im1,1) + nxk / 2
+ nincols = col2 - col1 + 1
+
+ # Flip the kernel
+ call rg_pflip (kernel, Memr[nkern], nxk, nyk)
+
+ # Initialise the line buffers.
+ inline = 1 - nyk / 2
+ do i = 1 , nyk - 1 {
+ Memi[lineptrs+i] = imgs2r (im1, col1, col2, inline, inline)
+ inline = inline + 1
+ }
+
+ # Generate the output image line by line
+ call salloc (linebuf, nincols, TY_REAL)
+ do outline = 1, nlines {
+
+ # Scroll the input buffers
+ do i = 1, nyk - 1
+ Memi[lineptrs+i-1] = Memi[lineptrs+i]
+
+ # Read in new image line
+ Memi[lineptrs+nyk-1] = imgs2r (im1, col1, col2, inline,
+ inline)
+
+ # Get output image line
+ outbuf = impl2r (im2, outline)
+ if (outbuf == EOF)
+ call error (0, "Error writing output image.")
+
+ # Generate output image line
+ call aclrr (Memr[outbuf], ncols)
+ do i = 1, nyk
+ call acnvr (Memr[Memi[lineptrs+i-1]], Memr[outbuf], ncols,
+ Memr[nkern+(i-1)*nxk], nxk)
+
+ inline = inline + 1
+ }
+
+ # Free the image buffer pointers
+ call sfree (sp)
+end
+
+
+# RG_PFLIP -- Flip the kernel in preparation for convolution.
+
+procedure rg_pflip (inkern, outkern, nxk, nyk)
+
+real inkern[nxk,nyk] # the input kernel
+real outkern[nxk,nyk] # the output kernel
+int nxk, nyk # the kernel dimensions
+
+int i, j
+
+begin
+ do j = 1, nyk {
+ do i = 1, nxk {
+ outkern[i,j] = inkern[nxk+1-i,nyk+1-j]
+ }
+ }
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgpfft.x b/pkg/images/immatch/src/psfmatch/rgpfft.x
new file mode 100644
index 00000000..b5f36375
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpfft.x
@@ -0,0 +1,443 @@
+
+# RG_PG10F -- Fetch the 0 component of the fft.
+
+real procedure rg_pg10f (fft, nxfft, nyfft)
+
+real fft[nxfft,nyfft] #I array containing 2 real ffts
+int nxfft #I x dimension of complex array
+int nyfft #I y dimension of complex array
+
+int xcen, ycen
+
+begin
+ xcen = nxfft / 2 + 1
+ ycen = nyfft / 2 + 1
+
+ return (fft[xcen,ycen])
+end
+
+
+# RG_PG1NORM -- Estimate the normalization factor by computing the amplitude
+# of the best fitting Gaussian. This routine may eventually be replaced by
+# on which does a complete Gaussian fit. The Gaussian is assumed to be
+# of the form g = a * exp (b * r * r). The input array is a 2D real array
+# storing 1 fft of dimension nxfft by nyfft in complex order with the
+# zero frequency in the center.
+
+real procedure rg_pg1norm (fft, nxfft, nyfft)
+
+real fft[nxfft,nyfft] #I array containing 2 real ffts
+int nxfft #I x dimension of complex array
+int nyfft #I y dimension of complex array
+
+int xcen, ycen
+real ln1, ln2, cx, cy
+
+begin
+ xcen = nxfft / 2 + 1
+ ycen = nyfft / 2 + 1
+
+ if (nxfft >= 8) {
+ ln1 = log (sqrt (fft[xcen-2,ycen] ** 2 + fft[xcen-1,ycen] ** 2))
+ ln2 = log (sqrt (fft[xcen-4,ycen] ** 2 + fft[xcen-3,ycen] ** 2))
+ cx = exp ((4.0 * ln1 - ln2) / 3.0)
+ } else
+ cx = 0.0
+
+ if (nyfft >= 4) {
+ ln1 = log (sqrt (fft[xcen,ycen-1] ** 2 + fft[xcen+1,ycen-1] ** 2))
+ ln2 = log (sqrt (fft[xcen,ycen-2] ** 2 + fft[xcen+1,ycen-2] ** 2))
+ cy = exp ((4.0 * ln1 - ln2) / 3.0)
+ } else
+ cy = 0.0
+
+ if (cx <= 0.0)
+ return (cy)
+ else if (cy <= 0.0)
+ return (cx)
+ else
+ return (0.5 * (cx + cy))
+end
+
+
+# RG_PG20F -- Fetch the 0 component of the fft.
+
+real procedure rg_pg20f (fft, nxfft, nyfft)
+
+real fft[nxfft,nyfft] #I array containing 2 real ffts
+int nxfft #I x dimension of complex array
+int nyfft #I y dimension of complex array
+
+int xcen, ycen
+
+begin
+ xcen = nxfft / 2 + 1
+ ycen = nyfft / 2 + 1
+
+ return (fft[xcen,ycen] / fft[xcen+1,ycen])
+end
+
+
+# RG_PG2NORM -- Estimate the normalization factor by computing the amplitude
+# of the best fitting Gaussian. This routine may eventually be replaced by
+# on which does a complete Gaussian fit. The Gaussian is assumed to be
+# of the form g = a * exp (b * r * r). The input array is a 2D real array
+# storing 2 2D ffts of dimension nxfft by nyfft in complex order with the
+# zero frequency in the center.
+
+real procedure rg_pg2norm (fft, nxfft, nyfft)
+
+real fft[nxfft,nyfft] #I array containing 2 real ffts
+int nxfft #I x dimension of complex array
+int nyfft #I y dimension of complex array
+
+int xcen, ycen
+real fftr, ffti, ln1r, ln2r, ln1i, ln2i, cxr, cyr, cxi, cyi, ampr, ampi
+
+begin
+
+ xcen = nxfft / 2 + 1
+ ycen = nyfft / 2 + 1
+
+ # Compute the x amplitude for the first fft.
+ if (nxfft >= 8) {
+
+ fftr = 0.5 * (fft[xcen+2,ycen] + fft[xcen-2,ycen])
+ ffti = 0.5 * (fft[xcen+3,ycen] - fft[xcen-1,ycen])
+ ln1r = log (sqrt (fftr ** 2 + ffti ** 2))
+ fftr = 0.5 * (fft[xcen+4,ycen] + fft[xcen-4,ycen])
+ ffti = 0.5 * (fft[xcen+5,ycen] - fft[xcen-3,ycen])
+ ln2r = log (sqrt (fftr ** 2 + ffti ** 2))
+
+ fftr = 0.5 * (fft[xcen+3,ycen] + fft[xcen-1,ycen])
+ ffti = -0.5 * (fft[xcen+2,ycen] - fft[xcen-2,ycen])
+ ln1i = log (sqrt (fftr ** 2 + ffti ** 2))
+ fftr = 0.5 * (fft[xcen+5,ycen] + fft[xcen-3,ycen])
+ ffti = -0.5 * (fft[xcen+4,ycen] - fft[xcen-4,ycen])
+ ln2i = log (sqrt (fftr ** 2 + ffti ** 2))
+
+ cxr = exp ((4.0 * ln1r - ln2r) / 3.0)
+ cxi = exp ((4.0 * ln1i - ln2i) / 3.0)
+
+ } else {
+
+ cxr = 0.0
+ cxi = 0.0
+
+ }
+
+ # Compute the y ratio.
+ if (nyfft >= 4) {
+
+ fftr = 0.5 * (fft[xcen,ycen+1] + fft[xcen,ycen-1])
+ ffti = 0.5 * (fft[xcen+1,ycen+1] - fft[xcen+1,ycen-1])
+ ln1r = log (sqrt (fftr ** 2 + ffti ** 2))
+ fftr = 0.5 * (fft[xcen,ycen+2] + fft[xcen,ycen-2])
+ ffti = 0.5 * (fft[xcen+1,ycen+2] - fft[xcen+1,ycen-2])
+ ln2r = log (sqrt (fftr ** 2 + ffti ** 2))
+
+ fftr = 0.5 * (fft[xcen+1,ycen+1] + fft[xcen+1,ycen-1])
+ ffti = -0.5 * (fft[xcen,ycen+1] - fft[xcen,ycen-1])
+ ln1i = log (sqrt (fftr ** 2 + ffti ** 2))
+ fftr = 0.5 * (fft[xcen+1,ycen+2] + fft[xcen+1,ycen-2])
+ ffti = -0.5 * (fft[xcen,ycen+2] - fft[xcen,ycen-2])
+ ln2i = log (sqrt (fftr ** 2 + ffti ** 2))
+
+ cyr = exp ((4.0 * ln1r - ln2r) / 3.0)
+ cyi = exp ((4.0 * ln1i - ln2i) / 3.0)
+
+ } else {
+
+ cyr = 0.0
+ cyi = 0.0
+
+ }
+
+ if (cxr <= 0.0)
+ ampr = cyr
+ else if (cyr <= 0.0)
+ ampr = cxr
+ else
+ ampr = 0.5 * (cxr + cyr)
+
+ if (cxi <= 0.0)
+ ampi = cyi
+ else if (cyi <= 0.0)
+ ampi = cxi
+ else
+ ampi = 0.5 * (cxi + cyi)
+
+ if (ampi <= 0.0)
+ return (INDEFR)
+ else
+ return (ampr /ampi)
+end
+
+
+# RG_PDIVFFT -- Unpack the two fft's, save the first fft, and compute the
+# quotient of the two ffts.
+
+procedure rg_pdivfft (fft1, fftnum, fftdenom, fft2, nxfft, nyfft)
+
+real fft1[nxfft,nyfft] # array containing 2 ffts of 2 real functions
+real fftnum[nxfft,nyfft] # the numerator fft
+real fftdenom[nxfft,nyfft] # the denominator fft
+real fft2[nxfft,nyfft] # fft of psf matching function
+int nxfft, nyfft # dimensions of fft
+
+int i, j, xcen, ycen, nxp2, nxp3, nyp2
+real c1, c2, h1r, h1i, h2r, h2i, denom
+
+begin
+ c1 = 0.5
+ c2 = -0.5
+ xcen = nxfft / 2 + 1
+ ycen = nyfft / 2 + 1
+ nxp2 = nxfft + 2
+ nxp3 = nxfft + 3
+ nyp2 = nyfft + 2
+
+ # Compute the 0 frequency point.
+ h1r = fft1[xcen,ycen]
+ h1i = 0.0
+ h2r = fft1[xcen+1,ycen]
+ h2i = 0.0
+ fftnum[xcen,ycen] = h1r
+ fftnum[xcen+1,ycen] = 0.0
+ fftdenom[xcen,ycen] = h2r
+ fftdenom[xcen+1,ycen] = 0.0
+ fft2[xcen,ycen] = h1r / h2r
+ fft2[xcen+1,ycen] = 0.0
+
+ #call eprintf ("fft11=%g fft21=%g\n")
+ #call pargr (fft1[1,1])
+ #call pargr (fft1[2,1])
+
+ # Compute the first point.
+ h1r = c1 * (fft1[1,1] + fft1[1,1])
+ h1i = 0.0
+ h2r = -c2 * (fft1[2,1] + fft1[2,1])
+ h2i = 0.0
+
+ fftnum[1,1] = h1r
+ fftnum[2,1] = h1i
+ fftdenom[1,1] = h2r
+ fftdenom[2,1] = h2i
+ denom = h2r * h2r + h2i * h2i
+ if (denom == 0.0) {
+ fft2[1,1] = 1.0
+ fft2[2,1] = 0.0
+ } else {
+ fft2[1,1] = (h1r * h2r + h1i * h2i) / denom
+ fft2[2,1] = (h1i * h2r - h2i * h1r) / denom
+ }
+
+ # Compute the x symmetry axis points.
+ do i = 3, xcen - 1, 2 {
+
+ h1r = c1 * (fft1[i,ycen] + fft1[nxp2-i,ycen])
+ h1i = c1 * (fft1[i+1,ycen] - fft1[nxp3-i,ycen])
+ h2r = -c2 * (fft1[i+1,ycen] + fft1[nxp3-i,ycen])
+ h2i = c2 * (fft1[i,ycen] - fft1[nxp2-i,ycen])
+
+ fftnum[i,ycen] = h1r
+ fftnum[i+1,ycen] = h1i
+ fftnum[nxp2-i,ycen] = h1r
+ fftnum[nxp3-i,ycen] = -h1i
+
+ fftdenom[i,ycen] = h2r
+ fftdenom[i+1,ycen] = h2i
+ fftdenom[nxp2-i,ycen] = h2r
+ fftdenom[nxp3-i,ycen] = -h2i
+
+ denom = h2r * h2r + h2i * h2i
+ if (denom == 0.0) {
+ fft2[i,ycen] = 1.0
+ fft2[i+1,ycen] = 0.0
+ } else {
+ fft2[i,ycen] = (h1r * h2r + h1i * h2i) / denom
+ fft2[i+1,ycen] = (h1i * h2r - h2i * h1r) / denom
+ }
+ fft2[nxp2-i,ycen] = fft2[i,ycen]
+ fft2[nxp3-i,ycen] = -fft2[i+1,ycen]
+
+ }
+
+ # Quit if the transform is 1D.
+ if (nyfft < 2)
+ return
+
+ # Compute the x axis points.
+ do i = 3, xcen + 1, 2 {
+
+ h1r = c1 * (fft1[i,1] + fft1[nxp2-i,1])
+ h1i = c1 * (fft1[i+1,1] - fft1[nxp3-i,1])
+ h2r = -c2 * (fft1[i+1,1] + fft1[nxp3-i,1])
+ h2i = c2 * (fft1[i,1] - fft1[nxp2-i,1])
+
+ fftnum[i,1] = h1r
+ fftnum[i+1,1] = h1i
+ fftnum[nxp2-i,1] = h1r
+ fftnum[nxp3-i,1] = -h1i
+
+ fftdenom[i,1] = h2r
+ fftdenom[i+1,1] = h2i
+ fftdenom[nxp2-i,1] = h2r
+ fftdenom[nxp3-i,1] = -h2i
+
+ denom = h2r * h2r + h2i * h2i
+ if (denom == 0) {
+ fft2[i,1] = 1.0
+ fft2[i+1,1] = 0.0
+ } else {
+ fft2[i,1] = (h1r * h2r + h1i * h2i) / denom
+ fft2[i+1,1] = (h1i * h2r - h2i * h1r) / denom
+ }
+ fft2[nxp2-i,1] = fft2[i,1]
+ fft2[nxp3-i,1] = -fft2[i+1,1]
+ }
+
+ # Compute the y symmetry axis points.
+ do i = 2, ycen - 1 {
+
+ h1r = c1 * (fft1[xcen,i] + fft1[xcen, nyp2-i])
+ h1i = c1 * (fft1[xcen+1,i] - fft1[xcen+1,nyp2-i])
+ h2r = -c2 * (fft1[xcen+1,i] + fft1[xcen+1,nyp2-i])
+ h2i = c2 * (fft1[xcen,i] - fft1[xcen,nyp2-i])
+
+ fftnum[xcen,i] = h1r
+ fftnum[xcen+1,i] = h1i
+ fftnum[xcen,nyp2-i] = h1r
+ fftnum[xcen+1,nyp2-i] = -h1i
+
+ fftdenom[xcen,i] = h2r
+ fftdenom[xcen+1,i] = h2i
+ fftdenom[xcen,nyp2-i] = h2r
+ fftdenom[xcen+1,nyp2-i] = -h2i
+
+ denom = h2r * h2r + h2i * h2i
+ if (denom == 0.0) {
+ fft2[xcen,i] = 1.0
+ fft2[xcen+1,i] = 0.0
+ } else {
+ fft2[xcen,i] = (h1r * h2r + h1i * h2i) / denom
+ fft2[xcen+1,i] = (h1i * h2r - h2i * h1r) / denom
+ }
+ fft2[xcen,nyp2-i] = fft2[xcen,i]
+ fft2[xcen+1,nyp2-i] = -fft2[xcen+1,i]
+
+ }
+
+ # Compute the y axis points.
+ do i = 2, ycen {
+
+ h1r = c1 * (fft1[1,i] + fft1[1,nyp2-i])
+ h1i = c1 * (fft1[2,i] - fft1[2,nyp2-i])
+ h2r = -c2 * (fft1[2,i] + fft1[2,nyp2-i])
+ h2i = c2 * (fft1[1,i] - fft1[1,nyp2-i])
+
+ fftnum[1,i] = h1r
+ fftnum[2,i] = h1i
+ fftnum[1,nyp2-i] = h1r
+ fftnum[2,nyp2-i] = -h1i
+
+ fftdenom[1,i] = h2r
+ fftdenom[2,i] = h2i
+ fftdenom[1,nyp2-i] = h2r
+ fftdenom[2,nyp2-i] = -h2i
+
+ denom = h2r * h2r + h2i * h2i
+ if (denom == 0.0) {
+ fft2[1,i] = 1.0
+ fft2[2,i] = 0.0
+ } else {
+ fft2[1,i] = (h1r * h2r + h1i * h2i) / denom
+ fft2[2,i] = (h1i * h2r - h2i * h1r) / denom
+ }
+ fft2[1,nyp2-i] = fft2[1,i]
+ fft2[2,nyp2-i] = -fft2[2,i]
+ }
+
+ # Compute the remainder of the transform.
+ do j = 2, ycen - 1 {
+
+ do i = 3, xcen - 1, 2 {
+
+ h1r = c1 * (fft1[i,j] + fft1[nxp2-i, nyp2-j])
+ h1i = c1 * (fft1[i+1,j] - fft1[nxp3-i,nyp2-j])
+ h2r = -c2 * (fft1[i+1,j] + fft1[nxp3-i,nyp2-j])
+ h2i = c2 * (fft1[i,j] - fft1[nxp2-i,nyp2-j])
+
+ fftnum[i,j] = h1r
+ fftnum[i+1,j] = h1i
+ fftnum[nxp2-i,nyp2-j] = h1r
+ fftnum[nxp3-i,nyp2-j] = -h1i
+
+ fftdenom[i,j] = h2r
+ fftdenom[i+1,j] = h2i
+ fftdenom[nxp2-i,nyp2-j] = h2r
+ fftdenom[nxp3-i,nyp2-j] = -h2i
+
+ denom = h2r * h2r + h2i * h2i
+ if (denom == 0.0) {
+ fft2[i,j] = 1.0
+ fft2[i+1,j] = 0.0
+ } else {
+ fft2[i,j] = (h1r * h2r + h1i * h2i) / denom
+ fft2[i+1,j] = (h1i * h2r - h2i * h1r) / denom
+ }
+ fft2[nxp2-i,nyp2-j] = fft2[i,j]
+ fft2[nxp3-i,nyp2-j] = - fft2[i+1,j]
+ }
+
+ do i = xcen + 2, nxfft, 2 {
+
+ h1r = c1 * (fft1[i,j] + fft1[nxp2-i, nyp2-j])
+ h1i = c1 * (fft1[i+1,j] - fft1[nxp3-i,nyp2-j])
+ h2r = -c2 * (fft1[i+1,j] + fft1[nxp3-i,nyp2-j])
+ h2i = c2 * (fft1[i,j] - fft1[nxp2-i,nyp2-j])
+
+ fftnum[i,j] = h1r
+ fftnum[i+1,j] = h1i
+ fftnum[nxp2-i,nyp2-j] = h1r
+ fftnum[nxp3-i,nyp2-j] = -h1i
+
+ fftdenom[i,j] = h2r
+ fftdenom[i+1,j] = h2i
+ fftdenom[nxp2-i,nyp2-j] = h2r
+ fftdenom[nxp3-i,nyp2-j] = -h2i
+
+ denom = h2r * h2r + h2i * h2i
+ if (denom == 0.0) {
+ fft2[i,j] = 1.0
+ fft2[i+1,j] = 0.0
+ } else {
+ fft2[i,j] = (h1r * h2r + h1i * h2i) / denom
+ fft2[i+1,j] = (h1i * h2r - h2i * h1r) / denom
+ }
+ fft2[nxp2-i,nyp2-j] = fft2[i,j]
+ fft2[nxp3-i,nyp2-j] = - fft2[i+1,j]
+
+ }
+ }
+end
+
+
+# RG_PNORM -- Insert the normalization value into the 0 frequency of the
+# fft. The fft is a 2D fft stored in a real array in complex order.
+# The fft is assumed to be centered.
+
+procedure rg_pnorm (fft, nxfft, nyfft, norm)
+
+real fft[ARB] #I the input fft
+int nxfft #I the x dimension of fft (complex storage)
+int nyfft #I the y dimension of the fft
+real norm #I the flux ratio
+
+int index
+
+begin
+ index = nxfft + 1 + 2 * (nyfft / 2) * nxfft
+ fft[index] = norm
+ fft[index+1] = 0.0
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgpfilter.x b/pkg/images/immatch/src/psfmatch/rgpfilter.x
new file mode 100644
index 00000000..63040b63
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpfilter.x
@@ -0,0 +1,502 @@
+include <math.h>
+
+# RG_PCOSBELL -- Apply a cosine bell function to the data.
+
+procedure rg_pcosbell (fft, nxfft, nyfft, sx1, sx2, sy1, sy2, radsym)
+
+real fft[ARB] #I/O the ifft to be filtered
+int nxfft #I the x dimension of the fft
+int nyfft #I the y dimension of the fft
+real sx1 #I inner x radius of the cosine bell filter
+real sx2 #I outer x radius of the cosine bell filter
+real sy1 #I inner y radius of the cosine bell filter
+real sy2 #I outer y radius of the cosine bell filter
+int radsym #I radial symmetry ?
+
+int i, j, index, xcen, ycen
+real factorx, factory, r1, r2, r, rj, cos2
+
+begin
+ # Compute the center of the fft.
+ xcen = (nxfft / 2) + 1
+ ycen = (nyfft / 2) + 1
+
+ if (radsym == NO) {
+
+ # Filter in the y direction independently.
+ if (IS_INDEFR(sy1))
+ r1 = 0.0
+ else
+ r1 = sy1
+ if (IS_INDEFR(sy2))
+ r2 = nyfft - ycen + 1
+ else
+ r2 = sy2
+ factory = HALFPI / (r2 - r1)
+ index = 1
+ do j = 1, nyfft {
+ r = abs (ycen - j)
+ if (r >= r2)
+ cos2 = 0.0
+ else if (r <= r1)
+ cos2 = 1.0
+ else
+ cos2 = cos ((r - r1) * factory) ** 2
+ call amulkr (fft[index], cos2, fft[index], 2 * nxfft)
+ index = index + 2 * nxfft
+ }
+
+ # Filter in the x direction independently.
+ if (IS_INDEFR(sx1))
+ r1 = 0.0
+ else
+ r1 = sx1
+ if (IS_INDEFR(sx2))
+ r2 = nxfft - xcen + 1
+ else
+ r2 = sx2
+ factorx = HALFPI / (r2 - r1)
+
+ do i = 1, nxfft {
+ r = abs (xcen - i)
+ if (r >= r2)
+ cos2 = 0.0
+ else if (r <= r1)
+ cos2 = 1.0
+ else
+ cos2 = cos ((r - r1) * factorx) ** 2
+ do j = 2 * i - 1, 2 * nxfft * nyfft, 2 * nxfft {
+ fft[j] = fft[j] * cos2
+ fft[j+1] = fft[j+1] * cos2
+ }
+ }
+
+ } else {
+
+ if (IS_INDEFR(sx1) && IS_INDEFR(sy1))
+ r1 = 0.0
+ else if (IS_INDEFR(sx1))
+ r1 = sy1
+ else if (IS_INDEFR(sy1))
+ r1 = sx1
+ else
+ r1 = (sx1 + sy1) / 2.0
+ if (IS_INDEFR(sx2) && IS_INDEFR(sy2))
+ r2 = (nxfft - xcen + 1 + nyfft - ycen + 1) / 2.0
+ else if (IS_INDEFR(sx2))
+ r2 = sy2
+ else if (IS_INDEFR(sy2))
+ r2 = sx2
+ else
+ r2 = (sx2 + sy2) / 2.0
+ factorx = HALFPI / (r2 - r1)
+
+ index = 0
+ do j = 1, nyfft {
+ rj = (ycen - j) ** 2
+ do i = 1, nxfft {
+ r = sqrt ((i - xcen) ** 2 + rj)
+ if (r >= r2) {
+ fft[index+2*i-1] = 0.0
+ fft[index+2*i] = 0.0
+ } else if (r > r1) {
+ fft[index+2*i-1] = fft[index+2*i-1] * cos ((r - r1) *
+ factorx) ** 2
+ fft[index+2*i] = fft[index+2*i] * cos ((r - r1) *
+ factorx) ** 2
+ }
+ }
+ index = index + 2 * nxfft
+ }
+ }
+end
+
+
+# RG_PREPLACE -- Replace low valued regions in the kernel fft with a Gaussian
+# extension.
+
+procedure rg_preplace (fft, fftdiv, nxfft, nyfft, pthreshold, norm)
+
+real fft[ARB] #I/O the fft of the kernel
+real fftdiv[ARB] #I the divisor fft
+int nxfft #I x dimension of the fft (complex storage)
+int nyfft #I y dimension of the fft
+real pthreshold #I the minimum percent amplitude in the divisor
+real norm #I the normalization value
+
+pointer sp, params
+int xcen, ycen, i, j, ri, rj, index
+real divpeak, a1, a2, a3, u, v, divisor, absv, phi
+
+begin
+ call smark (sp)
+ call salloc (params, 5, TY_REAL)
+
+ # Compute the central amplitude peak.
+ xcen = nxfft / 2 + 1
+ ycen = nyfft / 2 + 1
+ divpeak = pthreshold * fftdiv[1+nxfft+2*(ycen-1)*nxfft]
+
+ # Fit the parameters.
+ call rg_pgaussfit (fft, fftdiv, nxfft, nyfft, divpeak, norm,
+ Memr[params])
+
+ # Store the parameters in temporary variables.
+ a1 = Memr[params]
+ a2 = Memr[params+1]
+ a3 = Memr[params+2]
+ u = Memr[params+3]
+ v = Memr[params+4]
+
+ # Perform the extension.
+ index = 0
+ do j = 1, nyfft {
+ rj = j - ycen
+ do i = 1, nxfft {
+ ri = i - xcen
+ divisor = sqrt (fftdiv[index+2*i-1] ** 2 +
+ fftdiv[index+2*i] ** 2)
+ if (divisor < divpeak) {
+ absv = norm * exp (a1 * ri * ri + a2 * ri * rj + a3 *
+ rj * rj)
+ phi = u * ri + v * rj
+ fft[index+2*i-1] = absv * cos (phi)
+ fft[index+2*i] = absv * sin (phi)
+ }
+ }
+ index = index + 2 * nxfft
+ }
+
+ # Correct the first row.
+ do i = 1, 2 * nxfft, 2 {
+ fft[i] = sqrt (fft[i] ** 2 + fft[i+1] ** 2)
+ fft[i+1] = 0.0
+ }
+
+ # Correct the first column.
+ index = 1
+ do j = 2, nyfft {
+ fft[index] = sqrt (fft[index] ** 2 + fft[index+1] ** 2)
+ fft[index+1] = 0.0
+ index = index + 2 * nxfft
+ }
+
+ call sfree (sp)
+end
+
+
+# RG_PGMODEL -- Replace low values with a Gaussian mode.
+
+procedure rg_pgmodel (fft, fftdiv, nxfft, nyfft, pthreshold, norm)
+
+real fft[ARB] #I/O the fft of the kernel
+real fftdiv[ARB] #I the divisor fft
+int nxfft #I the x dimension of the fft
+int nyfft #I the y dimension of the fft
+real pthreshold #I the minimum percent amplitude in the divisor
+real norm #I the normalization factor
+
+pointer sp, params
+int xcen, ycen, i, j, index
+real divpeak, a1, a2, a3, u, v, absv, phi, ri, rj
+
+begin
+ call smark (sp)
+ call salloc (params, 5, TY_REAL)
+
+ # Compute the central amplitude peak.
+ xcen = nxfft / 2 + 1
+ ycen = nyfft / 2 + 1
+ divpeak = pthreshold * fftdiv[1+nxfft+2*(ycen-1)*nxfft]
+
+ # Fit the parameters.
+ call rg_pgaussfit (fft, fftdiv, nxfft, nyfft, divpeak, norm,
+ Memr[params])
+
+ # Store the parameters in temporary variables
+ a1 = Memr[params]
+ a2 = Memr[params+1]
+ a3 = Memr[params+2]
+ u = Memr[params+3]
+ v = Memr[params+4]
+
+ # Perform the extension.
+ index = 0
+ do j = 1, nyfft {
+ rj = j - ycen
+ do i = 1, nxfft {
+ ri = i - xcen
+ absv = norm * exp (a1 * ri * ri + a2 * ri * rj + a3 * rj * rj)
+ phi = u * ri + v * rj
+ fft[index+2*i-1] = absv * cos (phi)
+ fft[index+2*i] = absv * sin (phi)
+ }
+ index = index + 2 * nxfft
+ }
+
+ # Correct the first row.
+ do i = 1, 2 * nxfft, 2 {
+ fft[i] = sqrt (fft[i] ** 2 + fft[i+1] ** 2)
+ fft[i+1] = 0.0
+ }
+
+ # Correct the first column.
+ index = 1
+ do j = 2, nyfft {
+ fft[index] = sqrt (fft[index] ** 2 + fft[index+1] ** 2)
+ fft[index+1] = 0.0
+ index = index + 2 * nxfft
+ }
+
+ call sfree (sp)
+end
+
+
+# RG_PGAUSSFIT -- Procedure to compute the Gaussian parameters
+
+procedure rg_pgaussfit (fft, fftdiv, nxfft, nyfft, divpeak, norm, param)
+
+real fft[ARB] #I the fft of the kernel
+real fftdiv[ARB] #I the divisor fft
+int nxfft #I the x dimension of the fft
+int nyfft #I the y dimension of the fft
+real divpeak #I the minimum value in the divisor
+real norm #I the normalization value norm value
+real param[ARB] #O the output fitted parameters
+
+int i, j, yj, xcen, ycen
+double x, y, x2, xy, y2, z, wt, x2w, y2w, xyw, zw, xzw, yzw
+double sxxxx, sxxxy, sxxyy, sxyyy, syyyy, sxxz, sxyz, syyz, sxx, sxy
+double syy, sxz, syz
+pointer sp, mat
+real divisor
+
+begin
+ # Allocate temporary space.
+ call smark (sp)
+ call salloc (mat, 12, TY_DOUBLE)
+
+ # Define the center of the fft.
+ xcen = nxfft / 2 + 1
+ ycen = nyfft / 2 + 1
+
+ # Initialize.
+ sxxxx = 0.0d0
+ sxxxy = 0.0d0
+ sxxyy = 0.0d0
+ sxyyy = 0.0d0
+ syyyy = 0.0d0
+ sxxz = 0.0d0
+ sxyz = 0.0d0
+ syyz = 0.0d0
+ sxx = 0.0d0
+ sxy = 0.0d0
+ syy = 0.0d0
+ sxz = 0.0d0
+ syz = 0.0d0
+
+ do i = 1, nxfft {
+ x = i - xcen
+ yj = - ycen
+ do j = 2 * i - 1, 2 * nxfft * nyfft, 2 * nxfft {
+ yj = yj + 1
+ y = yj
+
+ # Skip low points in the fit.
+ divisor = sqrt (fftdiv[j] ** 2 + fftdiv[j+1] ** 2)
+ if (divisor < divpeak)
+ next
+ if (i == xcen || yj == ycen)
+ next
+
+ # Accumulate the intermediate products.
+ divisor = sqrt (fft[j] ** 2 + fft[j+1] ** 2)
+ if (divisor <= 0.0)
+ next
+ z = log (divisor / norm)
+ x2 = x * x
+ y2 = y * y
+ wt = 1.0 / sqrt (x2 + y2)
+ xy = x * y
+ x2w = x2 * wt
+ y2w = y2 * wt
+ xyw = xy * wt
+ zw = z * wt
+ xzw = x * zw
+ yzw = y * zw
+
+ # Accumulate the sums for the Gaussian.
+ sxxxx = sxxxx + x2 * x2w
+ sxxxy = sxxxy + x2 * xyw
+ sxxyy = sxxyy + x2 * y2w
+ sxyyy = sxyyy + xy * y2w
+ syyyy = syyyy + y2 * y2w
+ sxxz = sxxz + x * xzw
+ sxyz = sxyz + x * yzw
+ syyz = syyz + y * yzw
+
+ # New weight and z point.
+ wt = sqrt (fft[j] ** 2 + fft[j+1] ** 2) / norm
+ z = atan2 (fft[j+1], fft[j])
+
+ # Accumulate the sums for the shift determinantion.
+ sxx = sxx + x2 * wt
+ sxy = sxy + xy * wt
+ syy = syy + y2 * wt
+ sxz = sxz + x * z * wt
+ syz = syz + y * z * wt
+ }
+ }
+
+ # Solve for the gaussian.
+ Memd[mat] = sxxxx
+ Memd[mat+1] = sxxxy
+ Memd[mat+2] = sxxyy
+ Memd[mat+3] = sxxz
+ Memd[mat+4] = sxxxy
+ Memd[mat+5] = sxxyy
+ Memd[mat+6] = sxyyy
+ Memd[mat+7] = sxyz
+ Memd[mat+8] = sxxyy
+ Memd[mat+9] = sxyyy
+ Memd[mat+10] = syyyy
+ Memd[mat+11] = syyz
+ call rg_pgelim (Memd[mat], 3)
+ param[1] = Memd[mat+3]
+ param[2] = Memd[mat+7]
+ param[3] = Memd[mat+11]
+
+ # Solve for the shift.
+ Memd[mat] = sxx
+ Memd[mat+1] = sxy
+ Memd[mat+2] = sxz
+ Memd[mat+3] = sxy
+ Memd[mat+4] = syy
+ Memd[mat+5] = syz
+ call rg_pgelim (Memd[mat], 2)
+ param[4] = Memd[mat+2]
+ param[5] = Memd[mat+5]
+
+ call sfree (sp)
+end
+
+
+# RG_PGELIM -- Solve a matrix using Gaussian elimination.
+
+procedure rg_pgelim (a, n)
+
+double a[n+1,n] #I/O matrix to be solved
+int n #I number of variables
+
+int i, j, k
+double den, hold
+
+begin
+ do k = 1, n {
+
+ den = a[k,k]
+ if (den == 0.0d0) { # look for non-zero switch
+ do j = k + 1, n {
+ if (a[k,k] != 0.0d0) {
+ do i = k, n + 1 {
+ hold = a[i,j]
+ a[i,j] = a[i,k]
+ a[i,k] = hold
+ }
+ den = a[k,k]
+ }
+ }
+ if (den == 0.0d0) # if still zero, skip
+ next
+ }
+
+ do i = k, n + 1
+ a[i,k] = a[i,k] / den
+ do j = 1, n {
+ if (j != k) {
+ den = a[k,j]
+ do i = k, n + 1
+ a[i,j] = a[i,j] - a[i,k] * den
+ }
+ }
+ }
+end
+
+
+# RG_PNORMFILT -- Filter out any values greater than the normalization
+# from the kernel fft.
+
+procedure rg_pnormfilt (fft, nxfft, nyfft, norm)
+
+real fft[ARB] #I/O the input fft
+int nxfft #I the x length of the fft
+int nyfft #I the y length of the fft
+real norm #I the normalization factor
+
+int j, i_index
+
+begin
+ do j = 1, nyfft {
+ i_index = 1 + 2 * (j - 1) * nxfft
+ call rg_pnreplace (fft[i_index], nxfft, norm)
+ }
+end
+
+
+# RG_PFOURIER -- Compute the fourier spectrum of the convolution kernel.
+
+procedure rg_pfourier (fft, psfft, nxfft, nyfft)
+
+real fft[ARB] # the input fft
+real psfft[ARB] # fourier spectrum of the fft
+int nxfft # the x dimension of the fft
+int nyfft # the y dimension of the fft
+
+int j, i_index, o_index
+
+begin
+ do j = 1, nyfft {
+ i_index = 1 + 2 * (j - 1) * nxfft
+ o_index = 1 + (j - 1) * nxfft
+ call rg_pvfourier (fft[i_index], psfft[o_index], nxfft)
+ }
+end
+
+
+# RG_PVFOURIER -- Procedure to compute the fourier spectrum of a vector.
+
+procedure rg_pvfourier (a, b, nxfft)
+
+real a[ARB] # input vector in complex storage order
+real b[ARB] # output vector in real storage order
+int nxfft # length of vector
+
+int i
+
+begin
+ do i = 1, nxfft
+ b[i] = sqrt (a[2*i-1] ** 2 + a[2*i] ** 2)
+end
+
+
+# RG_PNREPLACE -- Replace values whose absolute value is greater than the
+# flux ratio.
+
+procedure rg_pnreplace (a, nxfft, norm)
+
+real a[ARB] #I/O ithe nput vector in complex storage order
+int nxfft #I the length of the vector
+real norm #I the flux ratio
+
+int i
+real val
+
+begin
+ do i = 1, 2 * nxfft, 2 {
+ val = sqrt (a[i] ** 2 + a[i+1] ** 2)
+ if (val > norm) {
+ a[i] = a[i] / val * norm
+ a[i+1] = a[i+1] / val * norm
+ }
+ }
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgpisfm.x b/pkg/images/immatch/src/psfmatch/rgpisfm.x
new file mode 100644
index 00000000..24df8fd7
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpisfm.x
@@ -0,0 +1,556 @@
+include <imhdr.h>
+include <ctype.h>
+include <gset.h>
+include "psfmatch.h"
+
+define HELPFILE "immatch$src/psfmatch/psfmatch.key"
+
+# Define the plot functions
+
+define PM_PPOWER 1
+define PM_PKERNEL 2
+
+# Define the plot types
+
+define PM_PCONTOUR 1
+define PM_PLINE 2
+define PM_PCOL 3
+
+# RG_PISFM -- Procedure to compute the shifts interactively.
+
+int procedure rg_pisfm (pm, imr, reglist, impsf, im1, imk, imp, im2, gd, id)
+
+pointer pm #I pointer to the psfmatch structure
+pointer imr #I/O pointer to the reference image/psf
+pointer reglist #I/O pointer to the regions list
+pointer impsf #I/O pointer to the input psf
+pointer im1 #I/O pointer to the input image
+pointer imp #I/O pointer to the fourier spectrum image
+pointer imk #I/O pointer to the kernel image
+pointer im2 #I/O pointer to the output image
+pointer gd #I graphics stream pointer
+pointer id #I display stream pointer
+
+int newref, newimage, newfourier, newfilter, plotfunc, plottype, wcs, key
+int newplot, ncolr, nliner, ip
+pointer sp, cmd
+real wx, wy
+int rg_pstati(), rg_psfm(), clgcur(), rg_pgqverify(), rg_pgtverify()
+int ctoi(), rg_pregions()
+pointer rg_pstatp()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ newref = YES
+ newimage = YES
+ newfourier = YES
+ newfilter = YES
+ ncolr = INDEFI
+ nliner = INDEFI
+ plotfunc = PM_PKERNEL
+ plottype = PM_PCONTOUR
+
+ # Compute the convolution kernel for the current image.
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONIMAGE && rg_pstati (pm,
+ NREGIONS) <= 0) {
+ call gclear (gd)
+ call gflush (gd)
+ call printf ("The objects list is empty\n")
+ } else {
+ if (rg_psfm (pm, imr, im1, impsf, imk, newref) == OK) {
+ call rg_pplot (gd, pm, ncolr, nliner, plotfunc, plottype)
+ newref = NO
+ newimage = NO
+ newfourier = NO
+ newfilter = NO
+ } else {
+ call gclear (gd)
+ call gflush (gd)
+ call rg_pstats (pm, IMAGE, Memc[cmd], SZ_FNAME)
+ call printf ("Error computing kernel for image %s\n")
+ call pargstr (Memc[cmd])
+ }
+ }
+ newplot = NO
+
+ # Loop over the cursor commands.
+ while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+ EOF) {
+
+ switch (key) {
+
+ # Print the help page.
+ case '?':
+ call gpagefile (gd, HELPFILE, "")
+
+ # Quit the task gracefully.
+ case 'q':
+ if (rg_pgqverify ("psfmatch", pm, imk, key) == YES) {
+ call sfree (sp)
+ return (rg_pgtverify (key))
+ }
+
+ # Process colon commands.
+ case ':':
+ for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1)
+ ;
+ switch (Memc[cmd+ip-1]) {
+
+ case 'x':
+ if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') {
+ call rg_pcolon (gd, pm, imr, reglist, impsf, im1, imk,
+ NULL, im2, Memc[cmd], newref, newimage,
+ newfourier, newfilter)
+ } else {
+ ip = ip + 1
+ if (ctoi (Memc[cmd], ip, ncolr) <= 0) {
+ switch (plotfunc) {
+ case PM_PPOWER:
+ ncolr = rg_pstati (pm, NXFFT) / 2 + 1
+ case PM_PKERNEL:
+ ncolr = rg_pstati (pm, KNX) / 2 + 1
+ default:
+ ncolr = rg_pstati (pm, KNX) / 2 + 1
+ }
+ }
+ plottype = PM_PCOL
+ newplot = YES
+ }
+
+ case 'y':
+ if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') {
+ call rg_pcolon (gd, pm, imr, reglist, impsf, im1, imk,
+ NULL, im2, Memc[cmd], newref, newimage,
+ newfourier, newfilter)
+ } else {
+ ip = ip + 1
+ if (ctoi (Memc[cmd], ip, nliner) <= 0) {
+ switch (plotfunc) {
+ case PM_PPOWER:
+ nliner = rg_pstati (pm, NYFFT) / 2 + 1
+ case PM_PKERNEL:
+ nliner = rg_pstati (pm, KNY) / 2 + 1
+ default:
+ nliner = rg_pstati (pm, KNY) / 2 + 1
+ }
+ }
+ plottype = PM_PLINE
+ newplot = YES
+ }
+
+
+ default:
+ call rg_pcolon (gd, pm, imr, reglist, impsf, im1, imk, NULL,
+ im2, Memc[cmd], newref, newimage, newfourier,
+ newfilter)
+ }
+
+ # Write the parameters to the parameter file.
+ case 'w':
+ call rg_pppars (pm)
+
+ # Recompute the convolution kernel function.
+ case 'f':
+
+ if (rg_pstati(pm,CONVOLUTION) == PM_CONIMAGE) {
+ if (newref == YES)
+ if (rg_pregions (reglist, imr, pm, 1, YES) > 0)
+ ;
+ else if (newimage == YES)
+ call rg_pindefr (pm)
+ }
+
+ if (rg_pstati (pm, NREGIONS) > 0 || rg_pstati (pm,
+ CONVOLUTION) != PM_CONIMAGE) {
+
+ if (newfourier == YES) {
+ call printf (
+ "\nRecomputing convolution kernel ...\n")
+ if (rg_psfm (pm, imr, im1, impsf, imk,
+ newref) == OK) {
+ ncolr = INDEFI
+ nliner = INDEFI
+ call rg_pplot (gd, pm, ncolr, nliner, plotfunc,
+ plottype)
+ newref = NO
+ newimage = NO
+ newfourier = NO
+ newfilter = NO
+ newplot = NO
+ } else
+ call printf (
+ "\nError computing new kernel ...\n")
+ }
+
+ if (newfilter == YES) {
+ if (Memr[rg_pstatp(pm,FFT)] != NULL) {
+ call rg_pfilter (pm)
+ ncolr = INDEFI
+ nliner = INDEFI
+ call rg_pplot (gd, pm, ncolr, nliner, plotfunc,
+ plottype)
+ newfilter = NO
+ newplot = NO
+ } else
+ call printf (
+ "The kernel fourier spectrum is undefined\n")
+ }
+
+ } else
+ call printf ("The objects list is empty\n")
+
+ # Draw a contour plot of the kernel.
+ case 'k':
+ if (plotfunc != PM_PKERNEL)
+ newplot = YES
+ if (plottype != PM_PCONTOUR)
+ newplot = YES
+ plotfunc = PM_PKERNEL
+ plottype = PM_PCONTOUR
+ ncolr = (1 + rg_pstati (pm, KNX)) / 2
+ nliner = (1 + rg_pstati (pm, KNY)) / 2
+
+ # Draw a contour plot of the fourier spectrum.
+ case 'p':
+ if (plotfunc != PM_PPOWER)
+ newplot = YES
+ if (plottype != PM_PCONTOUR)
+ newplot = YES
+ plotfunc = PM_PPOWER
+ plottype = PM_PCONTOUR
+ ncolr = (1 + rg_pstati (pm, NXFFT)) / 2
+ nliner = (1 + rg_pstati (pm, NYFFT)) / 2
+
+ # Plot a line of the current plot.
+ case 'x':
+ if (plottype != PM_PCOL)
+ newplot = YES
+ if (plottype == PM_PCONTOUR) {
+ ncolr = nint (wx)
+ nliner = nint (wy)
+ } else if (plottype == PM_PLINE) {
+ ncolr = nint (wx)
+ }
+ plottype = PM_PCOL
+
+ # Plot a line of the current plot.
+ case 'y':
+ if (plottype != PM_PLINE)
+ newplot = YES
+ if (plottype == PM_PCONTOUR) {
+ ncolr = nint (wx)
+ nliner = nint (wy)
+ } else if (plottype == PM_PCOL) {
+ ncolr = nint (wx)
+ }
+ plottype = PM_PLINE
+
+ # Redraw the current plot.
+ case 'r':
+ newplot = YES
+
+ # Do nothing gracefully.
+ default:
+ ;
+
+ }
+
+ if (newplot == YES) {
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONIMAGE &&
+ rg_pstati (pm, NREGIONS) <= 0) {
+ call printf ("Warning: The objects list is empty\n")
+ } else if (newref == YES || newimage == YES ||
+ newfourier == YES || newfilter == YES) {
+ call printf (
+ "Warning: Convolution kernel should be refit\n")
+ } else if (rg_pstatp (pm, CONV) != NULL) {
+ call rg_pplot (gd, pm, ncolr, nliner, plotfunc, plottype)
+ newplot = NO
+ } else {
+ call printf (
+ "Warning: The convolution kernel is undefined\n")
+ }
+ }
+
+ }
+
+ call sfree (sp)
+end
+
+
+define QUERY "[Hit return to continue, n next image, q quit, w quit and update parameters]"
+
+# RG_PGQVERIFY -- Print a message in the status line asking the user if they
+# really want to quit, returning YES if they really want to quit, NO otherwise.
+
+int procedure rg_pgqverify (task, pm, imk, ch)
+
+char task[ARB] # task name
+pointer pm # pointer to psfmatch structure
+pointer imk # pointer to kernel image
+int ch # character keystroke command
+
+int wcs, stat
+pointer sp, cmd
+real wx, wy
+bool streq()
+int clgcur(), rg_pstati()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ # Print the status line query in reverse video and get the keystroke.
+ call printf (QUERY)
+ if (clgcur ("gcommands", wx, wy, wcs, ch, Memc[cmd], SZ_LINE) == EOF)
+ ;
+
+ # Process the command.
+ if (ch == 'q') {
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL)
+ call rg_pwrite (pm, imk, NULL)
+ stat = YES
+ } else if (ch == 'w') {
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL)
+ call rg_pwrite (pm, imk, NULL)
+ if (streq ("psfmatch", task))
+ call rg_pppars (pm)
+ stat = YES
+ } else if (ch == 'n') {
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL)
+ call rg_pwrite (pm, imk, NULL)
+ stat = YES
+ } else {
+ stat = NO
+ }
+
+ call sfree (sp)
+
+ return (stat)
+end
+
+
+# RG_PGTVERIFY -- Verify whether or not the user truly wishes to quit the
+# task.
+
+int procedure rg_pgtverify (ch)
+
+int ch #I the input keystroke command
+
+begin
+ if (ch == 'q') {
+ return (YES)
+ } else if (ch == 'w') {
+ return (YES)
+ } else if (ch == 'n') {
+ return (NO)
+ } else {
+ return (NO)
+ }
+end
+
+
+# RG_PPLOT -- Draw the default plot of the kernel fourier spectrum or the
+# kernel itself.
+
+procedure rg_pplot (gd, pm, col, line, plotfunc, plottype)
+
+pointer gd #I pointer to the graphics stream
+pointer pm #I pointer to the psfmatch structure
+int col #I column of cross-correlation function to plot
+int line #I line of cross-correlation function to plot
+int plotfunc #I the default plot function type
+int plottype #I the default plot type
+
+int nx, ny
+pointer sp, title, str, data
+int rg_pstati(), strlen()
+pointer rg_pstatp()
+
+begin
+ if (gd == NULL)
+ return
+
+ # Allocate working space.
+ call smark (sp)
+ call salloc (title, SZ_LINE, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Initialize the plot title and data.
+ switch (plotfunc) {
+ case PM_PPOWER:
+ call sprintf (Memc[title], SZ_LINE,
+ "Fourier Spectrum for Reference: %s Image: %s")
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ data = rg_pstatp (pm, ASFFT)
+ nx = rg_pstati (pm, NXFFT)
+ ny = rg_pstati (pm, NYFFT)
+ case PM_PKERNEL:
+ call sprintf (Memc[title], SZ_LINE,
+ "Convolution Kernel for Reference: %s Image: %s")
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ data = rg_pstatp (pm, CONV)
+ nx = rg_pstati (pm, KNX)
+ ny = rg_pstati (pm, KNY)
+ default:
+ call sprintf (Memc[title], SZ_LINE,
+ "Convolution Kernel for Reference: %s Image: %s")
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ data = rg_pstatp (pm, CONV)
+ nx = rg_pstati (pm, KNX)
+ nx = rg_pstati (pm, KNY)
+ }
+ if (IS_INDEFI(col))
+ col = 1 + nx / 2
+ if (IS_INDEFI(line))
+ line = 1 + ny / 2
+
+ # Draw the plot.
+ if (ny == 1) {
+ switch (plotfunc) {
+ case PM_PPOWER:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nLine %d")
+ call pargi (1)
+ call rg_pcpline (gd, Memc[title], Memr[rg_pstatp(pm,ASFFT)],
+ nx, ny, 1)
+ case PM_PKERNEL:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nLine %d")
+ call pargi (1)
+ call rg_pcpline (gd, Memc[title], Memr[rg_pstatp(pm,CONV)],
+ nx, ny, 1)
+ default:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nLine %d")
+ call pargi (1)
+ call rg_pcpline (gd, Memc[title], Memr[rg_pstatp(pm,CONV)],
+ nx, ny, 1)
+ }
+ } else {
+ switch (plottype) {
+ case PM_PCONTOUR:
+ call rg_contour (gd, Memc[title], "", Memr[data], nx, ny)
+ case PM_PLINE:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nLine %d")
+ call pargi (line)
+ call rg_pcpline (gd, Memc[title], Memr[data], nx, ny, line)
+ case PM_PCOL:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nColumn %d")
+ call pargi (col)
+ call rg_pcpcol (gd, Memc[title], Memr[data], nx, ny, col)
+ default:
+ call rg_contour (gd, Memc[title], "", Memr[data], nx, ny)
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# RG_PCPLINE -- Plot a line of a 2D function.
+
+procedure rg_pcpline (gd, title, data, nx, ny, nline)
+
+pointer gd #I pointer to the graphics stream
+char title[ARB] #I title for the plot
+real data[nx,ARB] #I the input data array
+int nx, ny #I dimensions of the input data array
+int nline #I the line number
+
+int i
+pointer sp, str, x
+real ymin, ymax
+
+begin
+ # Return if no graphics stream.
+ if (gd == NULL)
+ return
+
+ # Check for valid line number.
+ if (nline < 1 || nline > ny)
+ return
+
+ # Allocate some working space.
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (x, nx, TY_REAL)
+
+ # Initialize the data.
+ do i = 1, nx
+ Memr[x+i-1] = i
+ call alimr (data[1,nline], nx, ymin, ymax)
+
+ # Set up the labels and the axes.
+ call gclear (gd)
+ call gswind (gd, 1.0, real (nx), ymin, ymax)
+ call glabax (gd, title, "X Lag", "X-Correlation Function")
+
+ # Plot the line profile.
+ call gseti (gd, G_PLTYPE, GL_SOLID)
+ call gpline (gd, Memr[x], data[1,nline], nx)
+ call gflush (gd)
+
+ call sfree (sp)
+end
+
+
+# RG_PCPCOL -- Plot a column of the cross-correlation function.
+
+procedure rg_pcpcol (gd, title, data, nx, ny, ncol)
+
+pointer gd #I pointer to the graphics stream
+char title[ARB] #I title of the column plot
+real data[nx,ARB] #I the input data array
+int nx, ny #I the dimensions of the input data array
+int ncol #I line number
+
+int i
+pointer sp, x, y
+real ymin, ymax
+
+begin
+ # Return if no graphics stream.
+ if (gd == NULL)
+ return
+
+ # Check for valid column number.
+ if (ncol < 1 || ncol > nx)
+ return
+
+ # Initialize.
+ call smark (sp)
+ call salloc (x, ny, TY_REAL)
+ call salloc (y, ny, TY_REAL)
+
+ # Get the data to be plotted.
+ do i = 1, ny {
+ Memr[x+i-1] = i
+ Memr[y+i-1] = data[ncol,i]
+ }
+ call alimr (Memr[y], ny, ymin, ymax)
+
+ # Set up the labels and the axes.
+ call gclear (gd)
+ call gswind (gd, 1.0, real (ny), ymin, ymax)
+ call glabax (gd, title, "Y Lag", "X-Correlation Function")
+
+ # Plot the profile.
+ call gseti (gd, G_PLTYPE, GL_SOLID)
+ call gpline (gd, Memr[x], Memr[y], ny)
+
+ call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgppars.x b/pkg/images/immatch/src/psfmatch/rgppars.x
new file mode 100644
index 00000000..c8d49baa
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgppars.x
@@ -0,0 +1,124 @@
+include "psfmatch.h"
+
+# RG_PGPARS -- Read in the psf matching algorithm parameters.
+
+procedure rg_pgpars (pm)
+
+pointer pm #I pointer to psfmatch structure
+
+int ival
+pointer sp, str
+bool clgetb()
+int clgwrd(), clgeti(), btoi()
+real clgetr()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Initialize the psf matching structure.
+ call rg_pinit (pm, clgwrd ("convolution", Memc[str], SZ_LINE,
+ PM_CTYPES))
+
+ # Define the data and kernel sizes.
+ ival = clgeti ("dnx")
+ if (mod (ival, 2) == 0)
+ ival = ival + 1
+ call rg_pseti (pm, DNX, ival)
+ ival = clgeti ("dny")
+ if (mod (ival, 2) == 0)
+ ival = ival + 1
+ call rg_pseti (pm, DNY, ival)
+ ival = clgeti ("pnx")
+ if (mod (ival, 2) == 0)
+ ival = ival + 1
+ call rg_pseti (pm, PNX, ival)
+ ival = clgeti ("pny")
+ if (mod (ival, 2) == 0)
+ ival = ival + 1
+ call rg_pseti (pm, PNY, ival)
+
+ # Centering parameters.
+ call rg_pseti (pm, CENTER, btoi (clgetb ("center")))
+
+ # Background value computation.
+ call clgstr ("background", Memc[str], SZ_LINE)
+ call rg_psets (pm, BSTRING, Memc[str])
+ call rg_psetr (pm, LOREJECT, clgetr ("loreject"))
+ call rg_psetr (pm, HIREJECT, clgetr ("hireject"))
+ call rg_psetr (pm, APODIZE, clgetr ("apodize"))
+
+ # Filtering parameters.
+ call rg_psetr (pm, UFLUXRATIO, clgetr ("fluxratio"))
+ call clgstr ("filter", Memc[str], SZ_LINE)
+ call rg_psets (pm, FSTRING, Memc[str])
+ call rg_psetr (pm, SXINNER, clgetr ("sx1"))
+ call rg_psetr (pm, SXOUTER, clgetr ("sx2"))
+ call rg_psetr (pm, SYINNER, clgetr ("sy1"))
+ call rg_psetr (pm, SYOUTER, clgetr ("sy2"))
+ call rg_pseti (pm, RADSYM, btoi (clgetb ("radsym")))
+ call rg_psetr (pm, THRESHOLD, (clgetr ("threshold")))
+
+ # Normalization parameter.
+ call rg_psetr (pm, NORMFACTOR, clgetr ("normfactor"))
+
+ #call rg_psetr (pm, PRATIO, clgetr ("pratio"))
+
+ call sfree (sp)
+end
+
+
+# RG_PPPARS -- Put the parameters required for the psf matching from
+# the cl to the parameter file.
+
+procedure rg_pppars (pm)
+
+pointer pm #I pointer to the psf matching structure
+
+pointer sp, str
+bool itob()
+int rg_pstati()
+real rg_pstatr()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Store the psf data string.
+ call rg_pstats (pm, PSFDATA, Memc[str], SZ_LINE)
+ call clpstr ("psf", Memc[str])
+
+ # Store the size parameters.
+ call clputi ("dnx", rg_pstati (pm, DNX))
+ call clputi ("dny", rg_pstati (pm, DNY))
+ call clputi ("pnx", rg_pstati (pm, PNX))
+ call clputi ("pny", rg_pstati (pm, PNY))
+
+ # Store the centering parameters.
+ call clputb ("center", itob (rg_pstati (pm, CENTER)))
+
+ # Store the background fitting parameters.
+ call rg_pstats (pm, BSTRING, Memc[str], SZ_LINE)
+ call clpstr ("background", Memc[str])
+ call clputr ("loreject", rg_pstatr (pm, LOREJECT))
+ call clputr ("hireject", rg_pstatr (pm, HIREJECT))
+ call clputr ("apodize", rg_pstatr (pm, APODIZE))
+
+ # Store the filtering parameters.
+ call clputr ("fluxratio", rg_pstatr(pm, UFLUXRATIO))
+ call rg_pstats (pm, FSTRING, Memc[str], SZ_LINE)
+ call clpstr ("filter", Memc[str])
+ call clputr ("sx1", rg_pstatr (pm, SXINNER))
+ call clputr ("sx2", rg_pstatr (pm, SXOUTER))
+ call clputr ("sy1", rg_pstatr (pm, SYINNER))
+ call clputr ("sy2", rg_pstatr (pm, SYOUTER))
+ call clputb ("radsym", itob (rg_pstati (pm, RADSYM)))
+ call clputr ("threshold", rg_pstatr (pm, THRESHOLD))
+
+ # Store the normalization parameters.
+ call clputr ("normfactor", rg_pstatr (pm, NORMFACTOR))
+
+ call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgpregions.x b/pkg/images/immatch/src/psfmatch/rgpregions.x
new file mode 100644
index 00000000..c04dcf97
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpregions.x
@@ -0,0 +1,464 @@
+include <fset.h>
+include <imhdr.h>
+include "psfmatch.h"
+
+# RG_PREGIONS -- Decoode the regions specification. If the sections
+# string is NULL then a default region dnx by dny pixels wide centered
+# on the reference image is used. Otherwise the section centers are
+# read from the regions string or from the objects list.
+
+int procedure rg_pregions (list, im, pm, rp, reread)
+
+int list #I pointer to regions file list
+pointer im #I pointer to the image
+pointer pm #I pointer to the psfmatch structure
+int rp #I region pointer
+int reread #I reread the current file
+
+char fname[SZ_FNAME]
+int nregions, fd
+int open(), rg_prregions(), rg_pgregions(), fntgfnb()
+int rg_pstati()
+data fname[1] /EOS/
+errchk open(), fntgfnb(), close()
+
+begin
+ if (rp < 1 || rp > MAX_NREGIONS) {
+ nregions = 0
+ } else if (rg_pgregions (im, pm, rp, MAX_NREGIONS) > 0) {
+ nregions = rg_pstati (pm, NREGIONS)
+ } else if (list != NULL) {
+ if (reread == NO) {
+ iferr {
+ if (fntgfnb (list, fname, SZ_FNAME) != EOF) {
+ fd = open (fname, READ_ONLY, TEXT_FILE)
+ nregions= rg_prregions (fd, im, pm, rp, MAX_NREGIONS)
+ call close (fd)
+ }
+ } then
+ nregions = 0
+ } else if (fname[1] != EOS) {
+ iferr {
+ fd = open (fname, READ_ONLY, TEXT_FILE)
+ nregions= rg_prregions (fd, im, pm, rp, MAX_NREGIONS)
+ call close (fd)
+ } then
+ nregions = 0
+ }
+ } else
+ nregions = 0
+
+ return (nregions)
+end
+
+
+# RG_PMKREGIONS -- Create a list of psf objects by selecting objects with
+# the image display cursor.
+
+int procedure rg_pmkregions (fd, im, pm, rp, max_nregions)
+
+int fd #I the output coordinates file descriptor
+pointer im #I pointer to the image
+pointer pm #I pointer to the psf matching structure
+int rp #I pointer to current region
+int max_nregions #I maximum number of regions
+
+int nregions, wcs, key, x1, x2, y1, y2
+pointer sp, region, cmd
+real x, y, xc, yc
+int clgcur(), rg_pstati()
+pointer rg_pstatp()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (region, SZ_FNAME, TY_CHAR)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ # Allocate the arrays to hold the regions information,
+ call rg_prealloc (pm, max_nregions)
+
+ nregions = min (rp-1, rg_pstati (pm, NREGIONS))
+ while (nregions < max_nregions) {
+
+ # Identify the object.
+ call printf ("Mark object %d [any key=mark,q=quit]:\n")
+ call pargi (nregions + 1)
+ if (clgcur ("icommands", x, y, wcs, key, Memc[cmd], SZ_LINE) == EOF)
+ break
+ if (key == 'q')
+ break
+
+ # Center the object.
+ if (rg_pstati (pm, CENTER) == YES) {
+ call rg_pcntr (im, x, y, max (rg_pstati(pm, PNX),
+ rg_pstati(pm, PNY)), xc, yc)
+ } else {
+ xc = x
+ yc = y
+ }
+
+ # Compute the data section.
+ x1 = xc - rg_pstati (pm, DNX) / 2
+ x2 = x1 + rg_pstati (pm, DNX) - 1
+ y1 = yc - rg_pstati (pm, DNY) / 2
+ y2 = y1 + rg_pstati (pm, DNY) - 1
+
+ # Make sure that the region is on the image.
+ if (x1 < 1 || x2 > IM_LEN(im,1) || y1 < 1 || y2 >
+ IM_LEN(im,2))
+ next
+
+ if (fd != NULL) {
+ call fprintf (fd, "%0.3f %0.3f\n")
+ call pargr (xc)
+ call pargr (yc)
+ }
+
+ Memi[rg_pstatp(pm,RC1)+nregions] = x1
+ Memi[rg_pstatp(pm,RC2)+nregions] = x2
+ Memi[rg_pstatp(pm,RL1)+nregions] = y1
+ Memi[rg_pstatp(pm,RL2)+nregions] = y2
+ Memr[rg_pstatp(pm,RZERO)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RXSLOPE)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RYSLOPE)+nregions] = INDEFR
+ nregions = nregions + 1
+
+ }
+
+ # Reallocate the correct amount of space.
+ call rg_pseti (pm, NREGIONS, nregions)
+ if (nregions > 0) {
+ call rg_prealloc (pm, nregions)
+ if (fd != NULL) {
+ call fstats (fd, F_FILENAME, Memc[region], SZ_FNAME)
+ call rg_psets (pm, PSFDATA, Memc[region])
+ } else
+ call rg_psets (pm, PSFDATA, "")
+ } else {
+ call rg_prfree (pm)
+ call rg_psets (pm, PSFDATA, "")
+ }
+
+ call sfree (sp)
+ return (nregions)
+end
+
+
+# RG_PRREGIONS -- Procedure to read the regions from a file.
+
+int procedure rg_prregions (fd, im, pm, rp, max_nregions)
+
+int fd #I regions file descriptor
+pointer im #I pointer to the image
+pointer pm #I pointer to psf matching structure
+int rp #I pointer to current region
+int max_nregions #I maximum number of regions
+
+int nregions, x1, y1, x2, y2
+pointer sp, line
+real x, y, xc, yc
+int rg_pstati(), getline()
+pointer rg_pstatp()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (line, SZ_LINE, TY_CHAR)
+
+ # Allocate the arrays to hold the regions information,
+ call rg_prealloc (pm, max_nregions)
+
+ # Decode the regions string.
+ nregions = min (rp - 1, rg_pstati (pm, NREGIONS))
+ while (getline (fd, Memc[line]) != EOF) {
+
+ if (nregions >= max_nregions)
+ break
+
+ call sscan (Memc[line])
+ call gargr (x)
+ call gargr (y)
+ if (rg_pstati (pm, CENTER) == YES) {
+ call rg_pcntr (im, x, y, max (rg_pstati(pm, PNX),
+ rg_pstati(pm, PNY)), xc, yc)
+ } else {
+ xc = x
+ yc = y
+ }
+
+ # Compute the data section.
+ x1 = xc - rg_pstati (pm, DNX) / 2
+ x2 = x1 + rg_pstati (pm, DNX) - 1
+ if (IM_NDIM(im) == 1) {
+ y1 = 1
+ y2 = 1
+ } else {
+ y1 = yc - rg_pstati (pm, DNY) / 2
+ y2 = y1 + rg_pstati (pm, DNY) - 1
+ }
+
+ # Make sure that the region is on the image.
+ if (x1 < 1 || x2 > IM_LEN(im,1) || y1 < 1 || y2 >
+ IM_LEN(im,2))
+ next
+
+ # Add the new region to the list.
+ Memi[rg_pstatp(pm,RC1)+nregions] = x1
+ Memi[rg_pstatp(pm,RC2)+nregions] = x2
+ Memi[rg_pstatp(pm,RL1)+nregions] = y1
+ Memi[rg_pstatp(pm,RL2)+nregions] = y2
+ Memr[rg_pstatp(pm,RZERO)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RXSLOPE)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RYSLOPE)+nregions] = INDEFR
+ nregions = nregions + 1
+ }
+
+ call rg_pseti (pm, NREGIONS, nregions)
+ if (nregions > 0)
+ call rg_prealloc (pm, nregions)
+ else
+ call rg_prfree (pm)
+
+ call sfree (sp)
+ return (nregions)
+end
+
+
+# RG_PGREGIONS -- Procedure to compute the column and line limits given
+# an x and y position and a default size.
+
+int procedure rg_pgregions (im, pm, rp, max_nregions)
+
+pointer im #I pointer to the image
+pointer pm #I pointer to psf matching structure
+int rp #I pointer to the current region
+int max_nregions #I maximum number of regions
+
+int ncols, nlines, nregions
+int x1, x2, y1, y2
+pointer sp, region
+real x, y, xc, yc
+int rg_pstati(), nscan()
+pointer rg_pstatp()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (region, SZ_LINE, TY_CHAR)
+
+ # Allocate the arrays to hold the regions information.
+ call rg_prealloc (pm, max_nregions)
+
+ # Get the constants.
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+
+ # Decode the center.
+ call rg_pstats (pm, PSFDATA, Memc[region], SZ_LINE)
+ nregions = min (rp - 1, rg_pstati (pm, NREGIONS))
+ call sscan (Memc[region])
+ call gargr (x)
+ call gargr (y)
+
+ # Compute the data region.
+ if (nscan() >= 2) {
+
+ # Compute a more accurate center.
+ if (rg_pstati (pm, CENTER) == YES) {
+ call rg_pcntr (im, x, y, max (rg_pstati(pm, PNX),
+ rg_pstati(pm, PNY)), xc, yc)
+ } else {
+ xc = x
+ yc = y
+ }
+
+ # Compute the data section.
+ x1 = xc - rg_pstati (pm, DNX) / 2
+ x2 = x1 + rg_pstati (pm, DNX) - 1
+ if (IM_NDIM(im) == 1) {
+ y1 = 1
+ y2 = 1
+ } else {
+ y1 = yc - rg_pstati (pm, DNY) / 2
+ y2 = y1 + rg_pstati (pm, DNY) - 1
+ }
+
+ # Make sure that the region is on the image.
+ if (x1 >= 1 && x2 <= IM_LEN(im,1) && y1 >= 1 &&
+ y2 <= IM_LEN(im,2)) {
+ Memi[rg_pstatp(pm,RC1)+nregions] = x1
+ Memi[rg_pstatp(pm,RC2)+nregions] = x2
+ Memi[rg_pstatp(pm,RL1)+nregions] = y1
+ Memi[rg_pstatp(pm,RL2)+nregions] = y2
+ Memr[rg_pstatp(pm,RZERO)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RXSLOPE)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RYSLOPE)+nregions] = INDEFR
+ nregions = nregions + 1
+ }
+ }
+
+
+ # Reallocate the correct amount of space.
+ call rg_pseti (pm, NREGIONS, nregions)
+ if (nregions > 0)
+ call rg_prealloc (pm, nregions)
+ else
+ call rg_prfree (pm)
+
+ call sfree (sp)
+
+ return (nregions)
+end
+
+
+# RG_PCNTR -- Compute star center using MPC algorithm.
+
+procedure rg_pcntr (im, xstart, ystart, boxsize, xcntr, ycntr)
+
+pointer im #I pointer to the input image
+real xstart, ystart #I initial position
+int boxsize #I width of the centering box
+real xcntr, ycntr #O computed center
+
+int x1, x2, y1, y2, half_box
+int ncols, nrows, nx, ny, try
+real xinit, yinit
+pointer bufptr, sp, x_vect, y_vect
+int imgs2r()
+
+begin
+ # Inialize.
+ half_box = (boxsize - 1) / 2
+ xinit = xstart
+ ncols = IM_LEN (im, 1)
+ if (IM_NDIM(im) == 1) {
+ yinit = 1
+ nrows = 1
+ } else {
+ yinit = ystart
+ nrows = IM_LEN (im, 2)
+ }
+ try = 0
+
+ # Iterate until pixel shifts are less than one.
+ repeat {
+
+ # Define region to extract.
+ x1 = max (xinit - half_box, 1.0) +0.5
+ x2 = min (xinit + half_box, real(ncols)) +0.5
+ y1 = max (yinit - half_box, 1.0) +0.5
+ y2 = min (yinit + half_box, real(nrows)) +0.5
+ nx = x2 - x1 + 1
+ ny = y2 - y1 + 1
+
+ # Extract region around center
+ bufptr = imgs2r (im, x1, x2, y1, y2)
+
+ # Compute the new center.
+ call smark (sp)
+ if (IM_NDIM(im) == 1) {
+ call salloc (x_vect, nx, TY_REAL)
+ call aclrr (Memr[x_vect], nx)
+ call rg_prowsum (Memr[bufptr], Memr[x_vect], nx, ny)
+ call rg_pcenter (Memr[x_vect], nx, xcntr)
+ ycntr = 1
+ } else {
+ call salloc (x_vect, nx, TY_REAL)
+ call salloc (y_vect, ny, TY_REAL)
+ call aclrr (Memr[x_vect], nx)
+ call aclrr (Memr[y_vect], ny)
+ call rg_prowsum (Memr[bufptr], Memr[x_vect], nx, ny)
+ call rg_pcolsum (Memr[bufptr], Memr[y_vect], nx, ny)
+ call rg_pcenter (Memr[x_vect], nx, xcntr)
+ call rg_pcenter (Memr[y_vect], ny, ycntr)
+ }
+ call sfree (sp)
+
+ # Check for INDEF centers.
+ if (IS_INDEFR(xcntr) || IS_INDEFR(ycntr)) {
+ xcntr = xinit
+ ycntr = yinit
+ break
+ }
+
+ # Add in offsets
+ xcntr = xcntr + x1
+ ycntr = ycntr + y1
+
+ try = try + 1
+ if (try == 1) {
+ if ((abs(xcntr-xinit) > 1.0) || (abs(ycntr-yinit) > 1.0)) {
+ xinit = xcntr
+ yinit = ycntr
+ }
+ } else
+ break
+ }
+end
+
+
+# RG_PROWSUM -- Sum all rows in a raster.
+
+procedure rg_prowsum (v, row, nx, ny)
+
+real v[nx,ny] #I the input subraster
+real row[ARB] #O the output row sum
+int nx, ny #I the dimensions of the subraster
+
+int i, j
+
+begin
+ do i = 1, ny
+ do j = 1, nx
+ row[j] = row[j] + v[j,i]
+end
+
+
+# RG_PCOLSUM -- Sum all columns in a raster.
+
+procedure rg_pcolsum (v, col, nx, ny)
+
+real v[nx,ny] #I the input subraster
+real col[ARB] #O the output column sum
+int nx, ny #I the dimensions of the subraster
+
+int i, j
+
+begin
+ do i = 1, ny
+ do j = 1, nx
+ col[j] = col[j] + v[i,j]
+end
+
+
+# RG_PCENTER -- Compute center of gravity of array.
+
+procedure rg_pcenter (v, nv, vc)
+
+real v[ARB] #I the input vector
+int nv #I the length of the vector
+real vc #O the output center
+
+int i
+real sum1, sum2, sigma, cont
+
+begin
+ # Compute first moment
+ sum1 = 0.0
+ sum2 = 0.0
+
+ call aavgr (v, nv, cont, sigma)
+
+ do i = 1, nv
+ if (v[i] > cont) {
+ sum1 = sum1 + (i-1) * (v[i] - cont)
+ sum2 = sum2 + (v[i] - cont)
+ }
+
+ # Determine center
+ if (sum2 == 0.0)
+ vc = INDEFR
+ else
+ vc = sum1 / sum2
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgpsfm.x b/pkg/images/immatch/src/psfmatch/rgpsfm.x
new file mode 100644
index 00000000..493d48c9
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpsfm.x
@@ -0,0 +1,815 @@
+include <imhdr.h>
+include <math/gsurfit.h>
+include "psfmatch.h"
+
+# RG_PSFM -- Procedure to match the psf functions of two images.
+
+int procedure rg_psfm (pm, imr, im1, impsf, imk, newref)
+
+pointer pm #I pointer to psf matching structure
+pointer imr #I pointer to reference image
+pointer im1 #I pointer to input image
+pointer impsf #I pointer to the psf image
+pointer imk #I pointer to kernel image
+int newref #I new reference image ?
+
+int stat
+int rg_pstati(), rg_pfget(), rg_psfget(), rg_kget()
+pointer rg_pstatp()
+
+begin
+ # Compute the convolution kernel.
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL) {
+
+ # Compute the kernel using raw image data or the psf image.
+ if (rg_pstati (pm,CONVOLUTION) == PM_CONIMAGE) {
+
+ # Set the kernel size to the user specified kernel size.
+ call rg_pseti (pm, KNX, rg_pstati (pm, PNX))
+ if (IM_NDIM(imr) == 1)
+ call rg_pseti (pm, KNY, 1)
+ else
+ call rg_pseti (pm, KNY, rg_pstati (pm, PNY))
+
+ # Compute the FFTS of the input and reference image.
+ stat = rg_pfget (pm, imr, im1, newref)
+
+ } else {
+
+ # Set the kernel size to the psf image size
+ call rg_pseti (pm, KNX, IM_LEN (impsf,1))
+ if (IM_NDIM(imr) == 1)
+ call rg_pseti (pm, KNY, 1)
+ else
+ call rg_pseti (pm, KNY, IM_LEN(impsf,2))
+
+ # Compute the FFTS of the input and reference psf images.
+ stat = rg_psfget (pm, imr, impsf, newref)
+ }
+
+ # Delete working arrays if an error occurs.
+ if (stat == ERR) {
+ if (rg_pstatp (pm, REFFFT) != NULL)
+ call mfree (rg_pstatp (pm, REFFFT), TY_REAL)
+ call rg_psetp (pm, REFFFT, NULL)
+ if (rg_pstatp (pm, IMFFT) != NULL)
+ call mfree (rg_pstatp (pm, IMFFT), TY_REAL)
+ call rg_psetp (pm, IMFFT, NULL)
+ if (rg_pstatp (pm, FFT) != NULL)
+ call mfree (rg_pstatp (pm, FFT), TY_REAL)
+ call rg_psetp (pm, FFT, NULL)
+ if (rg_pstatp (pm, CONV) != NULL)
+ call mfree (rg_pstatp (pm, CONV), TY_REAL)
+ call rg_psetp (pm, CONV, NULL)
+ if (rg_pstatp (pm, ASFFT) != NULL)
+ call mfree (rg_pstatp (pm, ASFFT), TY_REAL)
+ call rg_psetp (pm, ASFFT, NULL)
+ }
+
+ # Do the filtering in frequency space.
+ if (rg_pstatp (pm, FFT) != NULL)
+ call rg_pfilter (pm)
+
+ } else {
+
+ # Set the kernel size.
+ call rg_pseti (pm, KNX, IM_LEN(imk,1))
+ if (IM_NDIM(im1) == 1)
+ call rg_pseti (pm, KNY, 1)
+ else
+ call rg_pseti (pm, KNY, IM_LEN(imk,2))
+
+ # Read in the convolution kernel.
+ stat = rg_kget (pm, imk)
+
+ # Delete working arrays if an error occurs.
+ if (stat == ERR) {
+ if (rg_pstatp (pm, REFFFT) != NULL)
+ call mfree (rg_pstatp (pm, REFFFT), TY_REAL)
+ call rg_psetp (pm, REFFFT, NULL)
+ if (rg_pstatp (pm, IMFFT) != NULL)
+ call mfree (rg_pstatp (pm, IMFFT), TY_REAL)
+ call rg_psetp (pm, IMFFT, NULL)
+ if (rg_pstatp (pm, FFT) != NULL)
+ call mfree (rg_pstatp (pm, FFT), TY_REAL)
+ call rg_psetp (pm, FFT, NULL)
+ if (rg_pstatp (pm, CONV) != NULL)
+ call mfree (rg_pstatp (pm, CONV), TY_REAL)
+ call rg_psetp (pm, CONV, NULL)
+ if (rg_pstatp (pm, ASFFT) != NULL)
+ call mfree (rg_pstatp (pm, ASFFT), TY_REAL)
+ call rg_psetp (pm, ASFFT, NULL)
+ }
+ }
+
+ return (stat)
+end
+
+
+# RG_PFGET -- Compute the psfmatching function using Fourier techniques.
+
+int procedure rg_pfget (pm, imr, im1, newref)
+
+pointer pm #I pointer to psfmatch structure
+pointer imr #I pointer to reference image
+pointer im1 #I pointer to input image
+int newref #I new reference image ?
+
+int i, nregions, nrimcols, nrimlines, nrcols, nrlines, nrpcols, nrplines
+int nborder, stat, rc1, rc2, rl1, rl2, nxfft, nyfft
+pointer sp, str, coeff, dim, rbuf, ibuf, rsum, isum, border
+pointer prc1, prc2, prl1, prl2, przero, prxslope, pryslope, reffft, imfft, fft
+real rwtsum, iwtsum, rscale, iscale, rnscale, inscale
+bool fp_equalr()
+int rg_pstati(), rg_border(), rg_szfft()
+pointer rg_pstatp(), rg_pgdata()
+real rg_pstatr(), rg_pnsum(), rg_pg1norm(), rg_pg2norm()
+real rg_pg10f(), rg_pg20f()
+
+define nextimage_ 11
+
+begin
+ # Assemble the PSF data by looping over the regions list.
+ nregions = rg_pstati (pm, NREGIONS)
+ if (nregions <= 0)
+ return (ERR)
+
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (coeff, max (GS_SAVECOEFF+6, 9), TY_REAL)
+ call salloc (dim, 2, TY_INT)
+
+ # Get the reference region pointers.
+ prc1 = rg_pstatp (pm, RC1)
+ prc2 = rg_pstatp (pm, RC2)
+ prl1 = rg_pstatp (pm, RL1)
+ prl2 = rg_pstatp (pm, RL2)
+ przero = rg_pstatp (pm, RZERO)
+ prxslope = rg_pstatp (pm, RXSLOPE)
+ pryslope = rg_pstatp (pm, RYSLOPE)
+
+ # Check to see if the reference / input images are 1D.
+ nrimcols = IM_LEN(imr,1)
+ nrpcols = rg_pstati (pm, PNX)
+ if (IM_NDIM(imr) == 1) {
+ nrimlines = 1
+ nrplines = 1
+ } else {
+ nrimlines = IM_LEN(imr,2)
+ nrplines = rg_pstati (pm, PNY)
+ }
+
+ # Initialize
+ rwtsum = 0.0
+ iwtsum = 0.0
+ rnscale = INDEFR
+ inscale = INDEFR
+ rbuf = NULL
+ ibuf = NULL
+ stat = OK
+ if (newref == YES)
+ call calloc (rsum, rg_pstati (pm, DNX) * rg_pstati (pm, DNY),
+ TY_REAL)
+ call calloc (isum, rg_pstati (pm, DNX) * rg_pstati (pm, DNY),
+ TY_REAL)
+
+ do i = 1, nregions {
+
+ # Get the reference subraster regions.
+ rc1 = max (1, min (nrimcols, Memi[prc1+i-1]))
+ rc2 = min (nrimcols, max (1, Memi[prc2+i-1]))
+ rl1 = max (1, min (nrimlines, Memi[prl1+i-1]))
+ rl2 = min (nrimlines, max (1, Memi[prl2+i-1]))
+ nrcols = rc2 - rc1 + 1
+ nrlines = rl2 - rl1 + 1
+
+ # Go to next object if reference region is off the image.
+ if (nrcols < rg_pstati (pm, DNX) || (IM_NDIM(imr) == 2 &&
+ nrlines < rg_pstati(pm, DNY))) {
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_LINE)
+ call eprintf (
+ "Reference object %d: %s[%d:%d,%d:%d] is off image.\n")
+ call pargi (i)
+ call pargstr (Memc[str])
+ call pargi (rc1)
+ call pargi (rc2)
+ call pargi (rl1)
+ call pargi (rl2)
+ next
+ }
+
+ if (newref == YES) {
+
+ # Get the reference data.
+ rbuf = rg_pgdata (imr, rc1, rc2, rl1, rl2)
+
+ # Do the reference image background subtraction.
+ border = NULL
+ nborder = rg_border (Memr[rbuf], nrcols, nrlines, nrpcols,
+ nrplines, border)
+ call rg_pscale (pm, Memr[border], nborder, nrcols,
+ nrlines, nrpcols, nrplines, rg_pstatr (pm, BVALUER),
+ Memr[coeff])
+ if (border != NULL)
+ call mfree (border, TY_REAL)
+
+ # Save the coefficients.
+ Memr[przero+i-1] = Memr[coeff]
+ Memr[prxslope+i-1] = Memr[coeff+1]
+ Memr[pryslope+i-1] = Memr[coeff+2]
+
+ # Subtract the reference background.
+ call rg_subtract (Memr[rbuf], nrcols, nrlines,
+ Memr[przero+i-1], Memr[prxslope+i-1], Memr[pryslope+i-1])
+
+ # Apodize the reference image data.
+ if (rg_pstatr (pm, APODIZE) > 0.0)
+ call rg_apodize (Memr[rbuf], nrcols, nrlines,
+ rg_pstatr (pm, APODIZE), YES)
+
+ # Compute the scale factors and accumulate the weighted sums.
+ rscale = rg_pnsum (Memr[rbuf], nrcols, nrlines, nrpcols,
+ nrplines)
+ if (! IS_INDEFR(rscale)) {
+ if (IS_INDEFR(rnscale))
+ rnscale = 1.0 / rscale
+ }
+ if (IS_INDEFR(rscale))
+ rscale = 1.0
+ else
+ rscale = rscale / rnscale
+
+ call amulkr (Memr[rbuf], rscale, Memr[rbuf], nrcols *
+ nrlines)
+ rwtsum = rwtsum + rscale
+ call aaddr (Memr[rsum], Memr[rbuf], Memr[rsum], nrcols *
+ nrlines)
+
+ call mfree (rbuf, TY_REAL)
+ }
+
+ # Get the input image data
+ ibuf = rg_pgdata (im1, rc1, rc2, rl1, rl2)
+
+ # Compute the zero point, and the x and y slopes of input image.
+ border = NULL
+ nborder = rg_border (Memr[ibuf], nrcols, nrlines, nrpcols,
+ nrplines, border)
+ call rg_pscale (pm, Memr[border], nborder, nrcols, nrlines,
+ nrpcols, nrplines, rg_pstatr (pm, BVALUE), Memr[coeff])
+ if (border != NULL)
+ call mfree (border, TY_REAL)
+
+ # Subtract the background from the input image.
+ call rg_subtract (Memr[ibuf], nrcols, nrlines, Memr[coeff],
+ Memr[coeff+1], Memr[coeff+2])
+
+ # Apodize the data.
+ if (rg_pstatr (pm, APODIZE) > 0.0)
+ call rg_apodize (Memr[ibuf], nrcols, nrlines, rg_pstatr (pm,
+ APODIZE), YES)
+
+ # Compute the scale factors and accumulate the weighted sums for
+ # input image.
+ iscale = rg_pnsum (Memr[ibuf], nrcols, nrlines, nrpcols, nrplines)
+ if (! IS_INDEFR(iscale)) {
+ if (IS_INDEFR(inscale))
+ inscale = 1.0 / iscale
+ }
+ if (IS_INDEFR(iscale))
+ iscale = 1.0
+ else
+ iscale = iscale / inscale
+
+ call amulkr (Memr[ibuf], iscale, Memr[ibuf], nrcols * nrlines)
+ iwtsum = iwtsum + iscale
+ call aaddr (Memr[isum], Memr[ibuf], Memr[isum], nrcols * nrlines)
+
+ # Free the individual image buffers.
+ call mfree (ibuf, TY_REAL)
+ }
+
+ # Check to see if any data was read.
+ if (iwtsum <= 0.0) {
+ stat = ERR
+ goto nextimage_
+ }
+
+ # Normalize the summed buffers by the weights.
+ if (newref == YES) {
+ if (! fp_equalr (rwtsum, 0.0))
+ call adivkr (Memr[rsum], rwtsum, Memr[rsum], nrcols * nrlines)
+ }
+ if (! fp_equalr (iwtsum, 0.0))
+ call adivkr (Memr[isum], iwtsum, Memr[isum], nrcols * nrlines)
+
+ # Figure out how big the Fourier transform has to be, given
+ # the size of the reference subraster, the window size and
+ # the fact that the FFT must be a power of 2.
+
+ nxfft = rg_szfft (nrcols, 0)
+ if (nrlines == 1)
+ nyfft = 1
+ else
+ nyfft = rg_szfft (nrlines, 0)
+ call rg_pseti (pm, NXFFT, nxfft)
+ call rg_pseti (pm, NYFFT, nyfft)
+
+ imfft = rg_pstatp (pm, IMFFT)
+ if (imfft != NULL)
+ call mfree (imfft, TY_REAL)
+ call calloc (imfft, 2 * nxfft * nyfft, TY_REAL)
+ call rg_psetp (pm, IMFFT, imfft)
+
+ # Allocate space for the fft.
+ fft = rg_pstatp (pm, FFT)
+ if (fft != NULL)
+ call mfree (fft, TY_REAL)
+ call calloc (fft, 2 * nxfft * nyfft, TY_REAL)
+ call rg_psetp (pm, FFT, fft)
+
+ # Allocate space for the reference and input image ffts
+ if (newref == YES) {
+
+ reffft = rg_pstatp (pm, REFFFT)
+ if (reffft != NULL)
+ call mfree (reffft, TY_REAL)
+ call calloc (reffft, 2 * nxfft * nyfft, TY_REAL)
+ call rg_psetp (pm, REFFFT, reffft)
+
+ # Load the reference image FFT.
+ call rg_rload (Memr[rsum], nrcols, nrlines, Memr[fft], nxfft,
+ nyfft)
+ call mfree (rsum, TY_REAL)
+ rsum = NULL
+
+ # Load the input image FFT.
+ call rg_iload (Memr[isum], nrcols, nrlines, Memr[fft], nxfft,
+ nyfft)
+ call mfree (isum, TY_REAL)
+ isum = NULL
+
+ # Shift the data for easy of filtering.
+ call rg_fshift (Memr[fft], Memr[fft], 2 * nxfft, nyfft)
+
+ # Compute the Fourier Transform of the reference and input image
+ # data.
+ Memi[dim] = nxfft
+ Memi[dim+1] = nyfft
+ if (Memi[dim+1] == 1)
+ call rg_fourn (Memr[fft], Memi[dim], 1, 1)
+ else
+ call rg_fourn (Memr[fft], Memi[dim], 2, 1)
+
+ # Compute the flux ratio between the two data sets.
+ if (IS_INDEFR(rg_pstatr(pm, UFLUXRATIO))) {
+ if (rg_pstati (pm, BACKGRD) == PM_BNONE)
+ call rg_psetr (pm, FLUXRATIO, rg_pg2norm (Memr[fft],
+ 2 * nxfft, nyfft))
+ else
+ call rg_psetr (pm, FLUXRATIO, rg_pg20f (Memr[fft],
+ 2 * nxfft, nyfft))
+ } else
+ call rg_psetr (pm, FLUXRATIO, rg_pstatr (pm, UFLUXRATIO))
+
+ # Separate the two transforms and compute the division.
+ call rg_pdivfft (Memr[fft], Memr[reffft], Memr[imfft], Memr[fft],
+ 2 * nxfft, nyfft)
+
+ } else {
+
+
+ # Get the reference image FFT.
+ reffft = rg_pstatp (pm, REFFFT)
+
+ # Load the input image FFT.
+ call rg_rload (Memr[isum], nrcols, nrlines, Memr[imfft], nxfft,
+ nyfft)
+ call mfree (isum, TY_REAL)
+ isum = NULL
+
+ # Shift the data for easy of filtering.
+ call rg_fshift (Memr[imfft], Memr[imfft], 2 * nxfft, nyfft)
+
+ # Compute the Fourier Transform of the input image data.
+ Memi[dim] = nxfft
+ Memi[dim+1] = nyfft
+ if (Memi[dim+1] == 1)
+ call rg_fourn (Memr[imfft], Memi[dim], 1, 1)
+ else
+ call rg_fourn (Memr[imfft], Memi[dim], 2, 1)
+
+ # Compute the flux ratio between the two data sets.
+ if (IS_INDEFR(rg_pstatr(pm, UFLUXRATIO))) {
+ if (rg_pstati (pm, BACKGRD) == PM_BNONE)
+ call rg_psetr (pm, FLUXRATIO, rg_pg1norm (Memr[reffft],
+ 2 * nxfft, nyfft) / rg_pg1norm (Memr[imfft], 2 * nxfft,
+ nyfft))
+ else
+ call rg_psetr (pm, FLUXRATIO, rg_pg10f (Memr[reffft],
+ 2 * nxfft, nyfft) / rg_pg10f (Memr[imfft], 2 * nxfft,
+ nyfft))
+ } else
+ call rg_psetr (pm, FLUXRATIO, rg_pstatr (pm, UFLUXRATIO))
+
+ # Divide the two functions.
+ call adivx (Memr[reffft], Memr[imfft], Memr[fft], nxfft * nyfft)
+ }
+
+ # Normalize the FFT.
+ call rg_pnorm (Memr[fft], nxfft, nyfft, rg_pstatr (pm, FLUXRATIO))
+
+
+nextimage_
+
+ if (rsum != NULL)
+ call mfree (rsum, TY_REAL)
+ if (isum != NULL)
+ call mfree (isum, TY_REAL)
+ call sfree (sp)
+ if (stat == ERR)
+ return (ERR)
+ else
+ return (OK)
+end
+
+
+# RG_PSFGET -- Compute the psfmatching function using Fourier techniques.
+
+int procedure rg_psfget (pm, imr, impsf, newref)
+
+pointer pm #I pointer to the psfmatch structure
+pointer imr #I pointer to the reference psf
+pointer impsf #I pointer to the input image psf
+int newref #I new reference image
+
+int nrcols, nrlines, nxfft, nyfft
+pointer sp, dim, rbuf, ibuf, imfft, fft, reffft
+int rg_szfft()
+pointer rg_pgdata(), rg_pstatp()
+real rg_pstatr(), rg_pg2norm(), rg_pg1norm()
+
+begin
+ call smark (sp)
+ call salloc (dim, 2, TY_INT)
+
+ nrcols = IM_LEN(imr,1)
+ if (IM_NDIM(imr) == 1)
+ nrlines = 1
+ else
+ nrlines = IM_LEN(imr,2)
+
+ # Get the psf data.
+ rbuf = NULL
+ ibuf = NULL
+ if (newref == YES) {
+ call calloc (rbuf, nrcols * nrlines, TY_REAL)
+ rbuf = rg_pgdata (imr, 1, nrcols, 1, nrlines)
+ }
+ call calloc (ibuf, nrcols * nrlines, TY_REAL)
+ ibuf = rg_pgdata (impsf, 1, nrcols, 1, nrlines)
+
+ # Compute the size for the FFT buffers.
+ nxfft = rg_szfft (nrcols, 0)
+ if (nrlines == 1)
+ nyfft = 1
+ else
+ nyfft = rg_szfft (nrlines, 0)
+ call rg_pseti (pm, NXFFT, nxfft)
+ call rg_pseti (pm, NYFFT, nyfft)
+
+ imfft = rg_pstatp (pm, IMFFT)
+ if (imfft != NULL)
+ call mfree (imfft, TY_REAL)
+ call calloc (imfft, 2 * nxfft * nyfft, TY_REAL)
+ call rg_psetp (pm, IMFFT, imfft)
+
+ # Allocate space for the fft.
+ fft = rg_pstatp (pm, FFT)
+ if (fft != NULL)
+ call mfree (fft, TY_REAL)
+ call calloc (fft, 2 * nxfft * nyfft, TY_REAL)
+ call rg_psetp (pm, FFT, fft)
+
+ if (newref == YES) {
+
+ reffft = rg_pstatp (pm, REFFFT)
+ if (reffft != NULL)
+ call mfree (reffft, TY_REAL)
+ call calloc (reffft, 2 * nxfft * nyfft, TY_REAL)
+ call rg_psetp (pm, REFFFT, reffft)
+
+ # Load the reference image FFT.
+ call rg_rload (Memr[rbuf], nrcols, nrlines, Memr[fft], nxfft,
+ nyfft)
+
+ # Load the input image FFT.
+ call rg_iload (Memr[ibuf], nrcols, nrlines, Memr[fft], nxfft,
+ nyfft)
+
+ # Shift the data for easy of filtering.
+ call rg_fshift (Memr[fft], Memr[fft], 2 * nxfft, nyfft)
+
+ # Compute the Fourier Transform of the reference and input image
+ # data.
+ Memi[dim] = nxfft
+ Memi[dim+1] = nyfft
+ if (Memi[dim+1] == 1)
+ call rg_fourn (Memr[fft], Memi[dim], 1, 1)
+ else
+ call rg_fourn (Memr[fft], Memi[dim], 2, 1)
+
+ # Compute the flux ratio between the two data sets.
+ if (IS_INDEFR(rg_pstatr(pm, UFLUXRATIO)))
+ call rg_psetr (pm, FLUXRATIO, rg_pg2norm (Memr[fft],
+ 2 * nxfft, nyfft))
+ else
+ call rg_psetr (pm, FLUXRATIO, rg_pstatr(pm, UFLUXRATIO))
+
+ # Separate the two transforms and compute the division.
+ call rg_pdivfft (Memr[fft], Memr[reffft], Memr[imfft], Memr[fft],
+ 2 * nxfft, nyfft)
+
+ } else {
+
+ # Get the reference image FFT.
+ reffft = rg_pstatp (pm, REFFFT)
+
+ # Load the input image FFT.
+ call rg_rload (Memr[ibuf], nrcols, nrlines, Memr[imfft], nxfft,
+ nyfft)
+
+ # Shift the data for easy of filtering.
+ call rg_fshift (Memr[imfft], Memr[imfft], 2 * nxfft, nyfft)
+
+ # Compute the Fourier Transform of the input image data.
+ Memi[dim] = nxfft
+ Memi[dim+1] = nyfft
+ if (Memi[dim+1] == 1)
+ call rg_fourn (Memr[imfft], Memi[dim], 1, 1)
+ else
+ call rg_fourn (Memr[imfft], Memi[dim], 2, 1)
+
+ # Compute the flux ratio between the two data sets.
+ if (IS_INDEFR(rg_pstatr(pm, UFLUXRATIO)))
+ call rg_psetr (pm, FLUXRATIO, rg_pg1norm (Memr[reffft],
+ 2 * nxfft, nyfft) / rg_pg1norm (Memr[imfft], 2 * nxfft,
+ nyfft))
+ else
+ call rg_psetr (pm, FLUXRATIO, rg_pstatr(pm, UFLUXRATIO))
+
+ # Divide the two functions.
+ call adivx (Memr[reffft], Memr[imfft], Memr[fft], nxfft * nyfft)
+
+ }
+
+ # Normalize the FFT.
+ call rg_pnorm (Memr[fft], nxfft, nyfft, rg_pstatr (pm, FLUXRATIO))
+
+ # Free the data buffers.
+ if (rbuf != NULL)
+ call mfree (rbuf, TY_REAL)
+ if (ibuf != NULL)
+ call mfree (ibuf, TY_REAL)
+
+ call sfree (sp)
+
+ return (OK)
+end
+
+
+# RG_KGET -- Read in the convolution kernel.
+
+int procedure rg_kget (pm, imk)
+
+pointer pm #I pointer to the psfmatch structure
+pointer imk #I pointer to the kernel image
+
+int nrlines
+pointer conv
+pointer rg_pstatp(), rg_pgdata()
+
+begin
+ if (IM_NDIM(imk) == 1)
+ nrlines = 1
+ else
+ nrlines = IM_LEN(imk,2)
+ conv = rg_pstatp (pm, CONV)
+ if (conv != NULL)
+ call mfree (conv, TY_REAL)
+ conv = rg_pgdata (imk, 1, int(IM_LEN(imk,1)), 1, nrlines)
+ call rg_psetp (pm, CONV, conv)
+
+ return (OK)
+end
+
+
+# RG_PFILTER -- Procedure to filter the FFT in frequency space.
+
+procedure rg_pfilter (pm)
+
+pointer pm #I pointer to the psf matching structure
+
+pointer sp, dim, psfft, conv
+real nfactor
+int rg_pstati()
+pointer rg_pstatp()
+real rg_pstatr(), asumr()
+
+begin
+ call smark (sp)
+ call salloc (dim, 2, TY_INT)
+
+ # Allocate space for the fourier spectrum.
+ if (rg_pstatp (pm, ASFFT) != NULL)
+ call mfree (rg_pstatp (pm, ASFFT), TY_REAL)
+ call calloc (psfft, rg_pstati (pm, NXFFT) * rg_pstati (pm, NYFFT),
+ TY_REAL)
+ call rg_psetp (pm, ASFFT, psfft)
+
+ # Allocate space for the convolution kernel.
+ if (rg_pstatp (pm, CONV) != NULL)
+ call mfree (rg_pstatp (pm, CONV), TY_REAL)
+ call malloc (conv, 2 * rg_pstati (pm, NXFFT) * rg_pstati (pm, NYFFT),
+ TY_REAL)
+ call rg_psetp (pm, CONV, conv)
+ call amovr (Memr[rg_pstatp(pm,FFT)], Memr[rg_pstatp(pm,CONV)],
+ 2 * rg_pstati (pm, NXFFT) * rg_pstati (pm, NYFFT))
+
+# # Compute the zextend parameter.
+# call rg_psetr (pm, THRESHOLD, rg_pstatr (pm, PRATIO) *
+# rg_gnorm (Memr[rg_pstatp(pm,IMFFT)], rg_pstati(pm,NXFFT),
+# rg_pstati(pm,NYFFT)))
+
+ # Filter the frequency spectrum.
+ switch (rg_pstati(pm,FILTER)) {
+ case PM_FCOSBELL:
+ call rg_pcosbell (Memr[rg_pstatp(pm,CONV)], rg_pstati (pm, NXFFT),
+ rg_pstati (pm, NYFFT), rg_pstatr (pm, SXINNER), rg_pstatr (pm,
+ SXOUTER), rg_pstatr (pm, SYINNER), rg_pstatr (pm, SYOUTER),
+ rg_pstati (pm, RADSYM))
+ case PM_FREPLACE:
+ call rg_preplace (Memr[rg_pstatp(pm,CONV)], Memr[rg_pstatp(pm,
+ IMFFT)], rg_pstati (pm, NXFFT), rg_pstati (pm, NYFFT),
+ rg_pstatr (pm,THRESHOLD), rg_pstatr (pm,FLUXRATIO))
+ case PM_FMODEL:
+ call rg_pgmodel (Memr[rg_pstatp(pm,CONV)], Memr[rg_pstatp(pm,
+ IMFFT)], rg_pstati (pm, NXFFT), rg_pstati (pm, NYFFT),
+ rg_pstatr (pm, THRESHOLD), rg_pstatr (pm, FLUXRATIO))
+ default:
+ ;
+ }
+
+ # Filter out any values greater than the normalization.
+ call rg_pnormfilt (Memr[rg_pstatp(pm,CONV)], rg_pstati(pm,NXFFT),
+ rg_pstati(pm,NYFFT), rg_pstatr (pm, FLUXRATIO))
+
+ # Compute the fourier spectrum.
+ call rg_pfourier (Memr[rg_pstatp(pm,CONV)], Memr[rg_pstatp(pm,ASFFT)],
+ rg_pstati(pm,NXFFT), rg_pstati(pm,NYFFT))
+
+ Memi[dim] = rg_pstati (pm, NXFFT)
+ Memi[dim+1] = rg_pstati (pm, NYFFT)
+ call rg_fshift (Memr[rg_pstatp(pm,CONV)], Memr[rg_pstatp(pm,CONV)],
+ 2 * rg_pstati(pm, NXFFT), rg_pstati(pm, NYFFT))
+ call rg_fourn (Memr[rg_pstatp(pm,CONV)], Memi[dim], 2, -1)
+ call rg_fshift (Memr[rg_pstatp(pm,CONV)], Memr[rg_pstatp(pm,CONV)],
+ 2 * rg_pstati(pm, NXFFT), rg_pstati(pm, NYFFT))
+ call adivkr (Memr[rg_pstatp(pm,CONV)], real (rg_pstati(pm,NXFFT) *
+ rg_pstati(pm,NYFFT)), Memr[rg_pstatp(pm,CONV)], 2 * rg_pstati(pm,
+ NXFFT) * rg_pstati(pm,NYFFT))
+
+ # Unpack the convolution kernel.
+ call rg_movexr (Memr[rg_pstatp(pm,CONV)], rg_pstati(pm,NXFFT),
+ rg_pstati(pm,NYFFT), Memr[rg_pstatp(pm,CONV)], rg_pstati(pm,KNX),
+ rg_pstati(pm,KNY))
+
+ # Normalize the kernel.
+ if (! IS_INDEFR(rg_pstatr (pm, NORMFACTOR))) {
+ nfactor = rg_pstatr (pm, NORMFACTOR) / asumr (Memr[rg_pstatp(pm,
+ CONV)], rg_pstati (pm, KNX) * rg_pstati(pm,KNY))
+ call amulkr (Memr[rg_pstatp (pm,CONV)], nfactor,
+ Memr[rg_pstatp(pm, CONV)], rg_pstati (pm, KNX) *
+ rg_pstati (pm, KNY))
+ }
+
+ # Reallocate the convolution kernel array
+ #conv = rg_pstatp (pm, CONV)
+ #if (conv != NULL) {
+ #call realloc (conv, rg_pstati(pm, KNX) * rg_pstati(pm, KNY),
+ #TY_REAL)
+ #call rg_psetp (pm, CONV, conv)
+ #}
+
+ call sfree (sp)
+end
+
+
+# RG_PGDATA -- Fill a buffer from a specified region of the image.
+
+pointer procedure rg_pgdata (im, c1, c2, l1, l2)
+
+pointer im #I pointer to the iraf image
+int c1, c2 #I column limits in the input image
+int l1, l2 #I line limits in the input image
+
+int i, ncols, nlines, npts
+pointer ptr, index, buf
+pointer imgs1r(), imgs2r()
+
+begin
+ ncols = c2 - c1 + 1
+ nlines = l2 - l1 + 1
+ npts = ncols * nlines
+ call malloc (ptr, npts, TY_REAL)
+
+ index = ptr
+ do i = l1, l2 {
+ if (IM_NDIM(im) == 1)
+ buf = imgs1r (im, c1, c2)
+ else
+ buf = imgs2r (im, c1, c2, i, i)
+ call amovr (Memr[buf], Memr[index], ncols)
+ index = index + ncols
+ }
+
+ return (ptr)
+end
+
+
+# RG_PNSUM -- Compute the total intensity in the subtracted subraster.
+
+real procedure rg_pnsum (data, ncols, nlines, nxdata, nydata)
+
+real data[ncols,nlines] #I the input data subraster
+int ncols, nlines #I the size of the input subraster
+int nxdata, nydata #I the size of the data region
+
+int j, wxborder, wyborder, npts
+real sum
+bool fp_equalr()
+real asumr()
+
+begin
+ wxborder = (ncols - nxdata) / 2
+ wyborder = (nlines - nydata) / 2
+
+ sum = 0.0
+ npts = 0
+ do j = 1 + wyborder, nlines - wyborder {
+ sum = sum + asumr (data[1+wxborder,j], nxdata)
+ npts = npts + nxdata
+ }
+ if (npts <= 0 || fp_equalr (sum, 0.0))
+ return (INDEFR)
+ else
+ return (sum)
+end
+
+
+# RG_PWRITE -- Save the convolution kernel and the fourier spectrum of the
+# convolution kernel in an image.
+
+procedure rg_pwrite (pm, imk, imf)
+
+pointer pm #I pointer to psf matching structure
+pointer imk #I pointer to kernel image
+pointer imf #I pointer to fourier spectrum image
+
+int nx, ny
+pointer buf
+int rg_pstati()
+pointer rg_pstatp(), imps2r()
+
+begin
+ # Write out the kernel image.
+ if (imk != NULL && rg_pstatp(pm, CONV) != NULL) {
+ nx = rg_pstati (pm, KNX)
+ ny = rg_pstati (pm, KNY)
+ IM_NDIM(imk) = 2
+ IM_LEN(imk,1) = nx
+ IM_LEN(imk,2) = ny
+ IM_PIXTYPE(imk) = TY_REAL
+ buf = imps2r (imk, 1, nx, 1, ny)
+ if (rg_pstatp (pm, CONV) != NULL)
+ call amovr (Memr[rg_pstatp(pm,CONV)], Memr[buf], nx * ny)
+ else
+ call amovkr (0.0, Memr[buf], nx * ny)
+ }
+
+ # Write out the fourier spectrum.
+ if (imf != NULL && rg_pstatp(pm,ASFFT) != NULL) {
+ nx = rg_pstati (pm, NXFFT)
+ ny = rg_pstati (pm, NYFFT)
+ IM_NDIM(imf) = 2
+ IM_LEN(imf,1) = nx
+ IM_LEN(imf,2) = ny
+ IM_PIXTYPE(imf) = TY_REAL
+ buf = imps2r (imf, 1, nx, 1, ny)
+ if (rg_pstatp (pm, CONV) != NULL)
+ call amovr (Memr[rg_pstatp(pm,ASFFT)], Memr[buf], nx * ny)
+ else
+ call amovkr (0.0, Memr[buf], nx * ny)
+ }
+end
+
diff --git a/pkg/images/immatch/src/psfmatch/rgpshow.x b/pkg/images/immatch/src/psfmatch/rgpshow.x
new file mode 100644
index 00000000..c94349a6
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpshow.x
@@ -0,0 +1,116 @@
+include "psfmatch.h"
+
+# RG_PSHOW -- Print the PSFMATCH task parameters.
+
+procedure rg_pshow (pm)
+
+pointer pm #I pointer to psfmatch structure
+
+pointer sp, str
+bool itob()
+int rg_pstati()
+real rg_pstatr()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ call rg_pstats (pm, CSTRING, Memc[str], SZ_FNAME)
+ call printf ("\nConvolution: %s\n")
+ call pargstr (Memc[str])
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONIMAGE) {
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call printf (" %s: %s\n")
+ call pargstr (KY_IMAGE)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ call printf (" %s: %s\n")
+ call pargstr (KY_REFIMAGE)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, PSFDATA, Memc[str], SZ_FNAME)
+ call printf (" %s: %s\n")
+ call pargstr (KY_PSFDATA)
+ call pargstr (Memc[str])
+ } else if (rg_pstati (pm, CONVOLUTION) == PM_CONPSF) {
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call printf (" %s: %s\n")
+ call pargstr (KY_IMAGE)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, PSFIMAGE, Memc[str], SZ_FNAME)
+ call printf (" input psf: %s\n")
+ call pargstr (Memc[str])
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ call printf (" reference psf: %s\n")
+ call pargstr (Memc[str])
+ } else {
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call printf (" %s: %s\n")
+ call pargstr (KY_IMAGE)
+ call pargstr (Memc[str])
+ }
+
+ call rg_pstats (pm, KERNEL, Memc[str], SZ_FNAME)
+ call printf (" %s: %s\n")
+ call pargstr (KY_KERNEL)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, OUTIMAGE, Memc[str], SZ_FNAME)
+ if (Memc[str] != EOS) {
+ call printf (" %s: %s\n")
+ call pargstr (KY_OUTIMAGE)
+ call pargstr (Memc[str])
+ }
+
+ call printf ("Centering and background fitting\n")
+ call printf (" %s: %b\n")
+ call pargstr (KY_CENTER)
+ call pargb (itob(rg_pstati(pm,CENTER)))
+ call rg_pstats (pm, BSTRING, Memc[str], SZ_LINE)
+ call printf (" %s: %s\n")
+ call pargstr (KY_BACKGRD)
+ call pargstr (Memc[str])
+ call printf (" %s = %g %s = %g\n")
+ call pargstr (KY_LOREJECT)
+ call pargr (rg_pstatr (pm, LOREJECT))
+ call pargstr (KY_HIREJECT)
+ call pargr (rg_pstatr (pm, HIREJECT))
+ call printf (" %s = %g\n")
+ call pargstr (KY_APODIZE)
+ call pargr (rg_pstatr (pm, APODIZE))
+
+ call printf ("Filtering:\n")
+ call rg_pstats (pm, FSTRING, Memc[str], SZ_LINE)
+ call printf (" %s: %s\n")
+ call pargstr (KY_FILTER)
+ call pargstr (Memc[str])
+ if (rg_pstati(pm,FILTER) == PM_FCOSBELL) {
+ call printf (" %s: %g %s: %g\n")
+ call pargstr (KY_SXINNER)
+ call pargr (rg_pstatr (pm, SXINNER))
+ call pargstr (KY_SXOUTER)
+ call pargr (rg_pstatr (pm, SXOUTER))
+ call printf (" %s: %g %s: %g\n")
+ call pargstr (KY_SYINNER)
+ call pargr (rg_pstatr (pm, SYINNER))
+ call pargstr (KY_SYOUTER)
+ call pargr (rg_pstatr (pm, SYOUTER))
+ call printf (" %s: %b\n")
+ call pargstr (KY_RADSYM)
+ call pargb (itob(rg_pstati(pm,RADSYM)))
+ } else {
+ call printf (" %s: %g\n")
+ call pargstr (KY_UFLUXRATIO)
+ call pargr (rg_pstatr (pm, UFLUXRATIO))
+ call printf (" %s: %g\n")
+ call pargstr (KY_THRESHOLD)
+ call pargr (rg_pstatr(pm,THRESHOLD))
+ }
+
+ call printf ("Normalization\n")
+ call printf (" %s: %g\n")
+ call pargstr (KY_NORMFACTOR)
+ call pargr (rg_pstatr (pm, NORMFACTOR))
+
+ call printf ("\n")
+
+ call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/psfmatch/rgptools.x b/pkg/images/immatch/src/psfmatch/rgptools.x
new file mode 100644
index 00000000..df36c166
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgptools.x
@@ -0,0 +1,641 @@
+include "psfmatch.h"
+
+# RG_PINIT -- Initialize the main psfmatch data structure.
+
+procedure rg_pinit (pm, cfunc)
+
+pointer pm #O pointer to psfmatch structure
+int cfunc #I mode of computing the convolution function
+
+begin
+ call malloc (pm, LEN_PSFSTRUCT, TY_STRUCT)
+
+ # Initialize the pointers.
+ PM_RC1(pm) = NULL
+ PM_RC2(pm) = NULL
+ PM_RL1(pm) = NULL
+ PM_RL2(pm) = NULL
+ PM_RZERO(pm) = NULL
+ PM_RXSLOPE(pm) = NULL
+ PM_RYSLOPE(pm) = NULL
+ PM_NREGIONS(pm) = 0
+ PM_CNREGION(pm) = 1
+
+ # Define the background fitting parameters.
+ PM_CENTER(pm) = DEF_CENTER
+ PM_BACKGRD(pm) = DEF_BACKGRD
+ PM_BVALUER(pm) = 0.0
+ PM_BVALUE(pm) = 0.0
+ call strcpy ("median", PM_BSTRING(pm), SZ_FNAME)
+ PM_LOREJECT(pm) = DEF_LOREJECT
+ PM_HIREJECT(pm) = DEF_HIREJECT
+ PM_APODIZE(pm) = 0.0
+
+ PM_UFLUXRATIO(pm) = DEF_UFLUXRATIO
+ PM_FILTER(pm) = DEF_FILTER
+ call strcpy ("replace", PM_FSTRING(pm), SZ_FNAME)
+ PM_SXINNER(pm) = DEF_SXINNER
+ PM_SXOUTER(pm) = DEF_SXOUTER
+ PM_SYINNER(pm) = DEF_SYINNER
+ PM_SYOUTER(pm) = DEF_SYOUTER
+ PM_RADSYM(pm) = DEF_RADSYM
+ PM_THRESHOLD(pm) = DEF_THRESHOLD
+
+ PM_NORMFACTOR(pm) = DEF_NORMFACTOR
+
+ PM_CONVOLUTION(pm) = cfunc
+ switch (cfunc) {
+ case PM_CONIMAGE:
+ PM_CONVOLUTION(pm) = PM_CONIMAGE
+ call strcpy ("image", PM_CSTRING(pm), SZ_FNAME)
+ case PM_CONPSF:
+ PM_CONVOLUTION(pm) = PM_CONPSF
+ call strcpy ("psf", PM_CSTRING(pm), SZ_FNAME)
+ case PM_CONKERNEL:
+ PM_CONVOLUTION(pm) = PM_CONKERNEL
+ call strcpy ("kernel", PM_CSTRING(pm), SZ_FNAME)
+ default:
+ PM_CONVOLUTION(pm) = PM_CONIMAGE
+ call strcpy ("image", PM_CSTRING(pm), SZ_FNAME)
+ }
+ PM_DNX(pm) = DEF_DNX
+ PM_DNY(pm) = DEF_DNY
+ PM_PNX(pm) = DEF_PNX
+ PM_PNY(pm) = DEF_PNY
+ PM_KNX(pm) = 0
+ PM_KNY(pm) = 0
+ PM_POWER(pm) = DEF_POWER
+
+ PM_REFFFT(pm) = NULL
+ PM_IMFFT(pm) = NULL
+ PM_FFT(pm) = NULL
+ PM_CONV(pm) = NULL
+ PM_ASFFT(pm) = NULL
+ PM_NXFFT(pm) = 0
+ PM_NYFFT(pm) = 0
+
+ # Initialize the strings.
+ PM_IMAGE(pm) = EOS
+ PM_REFIMAGE(pm) = EOS
+ PM_PSFDATA(pm) = EOS
+ PM_PSFIMAGE(pm) = EOS
+ PM_OBJLIST(pm) = EOS
+ PM_KERNEL(pm) = EOS
+ PM_OUTIMAGE(pm) = EOS
+
+ # Initialize the buffers.
+ call rg_prinit (pm)
+end
+
+
+# RG_PRINIT -- Initialize the regions definition portion of the psf matching
+# code fitting structure.
+
+procedure rg_prinit (pm)
+
+pointer pm #I pointer to psfmatch structure
+
+begin
+ call rg_prfree (pm)
+
+ PM_NREGIONS(pm) = 0
+ PM_CNREGION(pm) = 1
+
+ call malloc (PM_RC1(pm), MAX_NREGIONS, TY_INT)
+ call malloc (PM_RC2(pm), MAX_NREGIONS, TY_INT)
+ call malloc (PM_RL1(pm), MAX_NREGIONS, TY_INT)
+ call malloc (PM_RL2(pm), MAX_NREGIONS, TY_INT)
+ call malloc (PM_RZERO(pm), MAX_NREGIONS, TY_REAL)
+ call malloc (PM_RXSLOPE(pm), MAX_NREGIONS, TY_REAL)
+ call malloc (PM_RYSLOPE(pm), MAX_NREGIONS, TY_REAL)
+
+ call amovki (INDEFI, Memi[PM_RC1(pm)], MAX_NREGIONS)
+ call amovki (INDEFI, Memi[PM_RC2(pm)], MAX_NREGIONS)
+ call amovki (INDEFI, Memi[PM_RL1(pm)], MAX_NREGIONS)
+ call amovki (INDEFI, Memi[PM_RL2(pm)], MAX_NREGIONS)
+ call amovkr (INDEFR, Memr[PM_RZERO(pm)], MAX_NREGIONS)
+ call amovkr (INDEFR, Memr[PM_RXSLOPE(pm)], MAX_NREGIONS)
+ call amovkr (INDEFR, Memr[PM_RYSLOPE(pm)], MAX_NREGIONS)
+end
+
+
+# RG_PINDEFR -- Re-initialize the background and answers regions portion of
+# the psf-matching structure.
+
+procedure rg_pindefr (pm)
+
+pointer pm #I pointer to the psfmatch structure
+
+int nregions
+int rg_pstati ()
+
+begin
+ nregions = rg_pstati (pm, NREGIONS)
+
+ if (nregions > 0) {
+ call amovkr (INDEFR, Memr[PM_RZERO(pm)], nregions)
+ call amovkr (INDEFR, Memr[PM_RXSLOPE(pm)], nregions)
+ call amovkr (INDEFR, Memr[PM_RYSLOPE(pm)], nregions)
+ }
+end
+
+
+# RG_PREALLOC -- Reallocate the regions buffers and initialize if necessary.
+
+procedure rg_prealloc (pm, nregions)
+
+pointer pm #I pointer to psfmatch structure
+int nregions #I number of regions
+
+int nr
+int rg_pstati()
+
+begin
+ nr = rg_pstati (pm, NREGIONS)
+
+ call realloc (PM_RC1(pm), nregions, TY_INT)
+ call realloc (PM_RC2(pm), nregions, TY_INT)
+ call realloc (PM_RL1(pm), nregions, TY_INT)
+ call realloc (PM_RL2(pm), nregions, TY_INT)
+ call realloc (PM_RZERO(pm), nregions, TY_REAL)
+ call realloc (PM_RXSLOPE(pm), nregions, TY_REAL)
+ call realloc (PM_RYSLOPE(pm), nregions, TY_REAL)
+
+ call amovki (INDEFI, Memi[PM_RC1(pm)+nr], nregions - nr)
+ call amovki (INDEFI, Memi[PM_RC2(pm)+nr], nregions - nr)
+ call amovki (INDEFI, Memi[PM_RL1(pm)+nr], nregions - nr)
+ call amovki (INDEFI, Memi[PM_RL2(pm)+nr], nregions - nr)
+ call amovkr (INDEFR, Memr[PM_RZERO(pm)+nr], nregions - nr)
+ call amovkr (INDEFR, Memr[PM_RXSLOPE(pm)+nr], nregions - nr)
+ call amovkr (INDEFR, Memr[PM_RYSLOPE(pm)+nr], nregions - nr)
+ #call amovkr (INDEFR, Memr[PM_XSHIFTS(pm)+nr], nregions - nr)
+ #call amovkr (INDEFR, Memr[PM_YSHIFTS(pm)+nr], nregions - nr)
+end
+
+
+# RG_PRFREE -- Free the regions portion of the psfmatch structure.
+
+procedure rg_prfree (pm)
+
+pointer pm #I/O pointer to psfmatch structure
+
+begin
+ call rg_pseti (pm, NREGIONS, 0)
+ if (PM_RC1(pm) != NULL)
+ call mfree (PM_RC1(pm), TY_INT)
+ PM_RC1(pm) = NULL
+ if (PM_RC2(pm) != NULL)
+ call mfree (PM_RC2(pm), TY_INT)
+ PM_RC2(pm) = NULL
+ if (PM_RL1(pm) != NULL)
+ call mfree (PM_RL1(pm), TY_INT)
+ PM_RL1(pm) = NULL
+ if (PM_RL2(pm) != NULL)
+ call mfree (PM_RL2(pm), TY_INT)
+ PM_RL2(pm) = NULL
+ if (PM_RZERO(pm) != NULL)
+ call mfree (PM_RZERO(pm), TY_REAL)
+ PM_RZERO(pm) = NULL
+ if (PM_RXSLOPE(pm) != NULL)
+ call mfree (PM_RXSLOPE(pm), TY_REAL)
+ PM_RXSLOPE(pm) = NULL
+ if (PM_RYSLOPE(pm) != NULL)
+ call mfree (PM_RYSLOPE(pm), TY_REAL)
+ PM_RYSLOPE(pm) = NULL
+end
+
+
+# RG_PFREE -- Free the psfmatch structure.
+
+procedure rg_pfree (pm)
+
+pointer pm #I pointer to psfmatch structure
+
+begin
+ # Free the region descriptors
+ call rg_prfree (pm)
+
+ if (PM_REFFFT(pm) != NULL)
+ call mfree (PM_REFFFT(pm), TY_REAL)
+ if (PM_IMFFT(pm) != NULL)
+ call mfree (PM_IMFFT(pm), TY_REAL)
+ if (PM_FFT(pm) != NULL)
+ call mfree (PM_FFT(pm), TY_REAL)
+ if (PM_CONV(pm) != NULL)
+ call mfree (PM_CONV(pm), TY_REAL)
+ if (PM_ASFFT(pm) != NULL)
+ call mfree (PM_ASFFT(pm), TY_REAL)
+
+ call mfree (pm, TY_STRUCT)
+end
+
+
+# RG_PSTATI -- Fetch the value of a psfmatch task integer parameter.
+
+int procedure rg_pstati (pm, param)
+
+pointer pm # pointer to psfmatch structure
+int param # parameter to be fetched
+
+begin
+ switch (param) {
+ case NREGIONS:
+ return (PM_NREGIONS(pm))
+ case CNREGION:
+ return (PM_CNREGION(pm))
+ case CENTER:
+ return (PM_CENTER(pm))
+ case BACKGRD:
+ return (PM_BACKGRD(pm))
+ case CONVOLUTION:
+ return (PM_CONVOLUTION(pm))
+ case DNX:
+ return (PM_DNX(pm))
+ case DNY:
+ return (PM_DNY(pm))
+ case PNX:
+ return (PM_PNX(pm))
+ case PNY:
+ return (PM_PNY(pm))
+ case KNX:
+ return (PM_KNX(pm))
+ case KNY:
+ return (PM_KNY(pm))
+ case POWER:
+ return (PM_POWER(pm))
+
+ case FILTER:
+ return (PM_FILTER(pm))
+ case RADSYM:
+ return (PM_RADSYM(pm))
+
+ case NXFFT:
+ return (PM_NXFFT(pm))
+ case NYFFT:
+ return (PM_NYFFT(pm))
+
+ default:
+ call error (0, "RG_PSTATI: Unknown integer parameter.")
+ }
+end
+
+
+# RG_PSTATP -- Fetch the value of a psfmatch task pointer parameter.
+
+pointer procedure rg_pstatp (pm, param)
+
+pointer pm # pointer to psfmatch structure
+int param # parameter to be fetched
+
+begin
+ switch (param) {
+ case RC1:
+ return (PM_RC1(pm))
+ case RC2:
+ return (PM_RC2(pm))
+ case RL1:
+ return (PM_RL1(pm))
+ case RL2:
+ return (PM_RL2(pm))
+ case RZERO:
+ return (PM_RZERO(pm))
+ case RXSLOPE:
+ return (PM_RXSLOPE(pm))
+ case RYSLOPE:
+ return (PM_RYSLOPE(pm))
+ case REFFFT:
+ return (PM_REFFFT(pm))
+ case IMFFT:
+ return (PM_IMFFT(pm))
+ case FFT:
+ return (PM_FFT(pm))
+ case CONV:
+ return (PM_CONV(pm))
+ case ASFFT:
+ return (PM_ASFFT(pm))
+ default:
+ call error (0, "RG_PSTATP: Unknown pointer parameter.")
+ }
+end
+
+
+# RG_PSTATR -- Fetch the value of a psfmath task real parameter.
+
+real procedure rg_pstatr (pm, param)
+
+pointer pm # pointer to psfmatch structure
+int param # parameter to be fetched
+
+begin
+ switch (param) {
+ case BVALUER:
+ return (PM_BVALUER(pm))
+ case BVALUE:
+ return (PM_BVALUE(pm))
+ case APODIZE:
+ return (PM_APODIZE(pm))
+ case LOREJECT:
+ return (PM_LOREJECT(pm))
+ case HIREJECT:
+ return (PM_HIREJECT(pm))
+ case UFLUXRATIO:
+ return (PM_UFLUXRATIO(pm))
+ case FLUXRATIO:
+ return (PM_FLUXRATIO(pm))
+ case SXINNER:
+ return (PM_SXINNER(pm))
+ case SXOUTER:
+ return (PM_SXOUTER(pm))
+ case SYINNER:
+ return (PM_SYINNER(pm))
+ case SYOUTER:
+ return (PM_SYOUTER(pm))
+ case THRESHOLD:
+ return (PM_THRESHOLD(pm))
+ case NORMFACTOR:
+ return (PM_NORMFACTOR(pm))
+ default:
+ call error (0, "RG_PSTATR: Unknown real parameter.")
+ }
+end
+
+
+# RG_PSTATS -- Fetch the value of a psfmatch string string parameter.
+
+procedure rg_pstats (pm, param, str, maxch)
+
+pointer pm # pointer to psfmatch structure
+int param # parameter to be fetched
+char str[ARB] # output string
+int maxch # maximum number of characters
+
+begin
+ switch (param) {
+ case BSTRING:
+ call strcpy (PM_BSTRING(pm), str, maxch)
+ case CSTRING:
+ call strcpy (PM_CSTRING(pm), str, maxch)
+ case FSTRING:
+ call strcpy (PM_FSTRING(pm), str, maxch)
+ case IMAGE:
+ call strcpy (PM_IMAGE(pm), str, maxch)
+ case REFIMAGE:
+ call strcpy (PM_REFIMAGE(pm), str, maxch)
+ case PSFDATA:
+ call strcpy (PM_PSFDATA(pm), str, maxch)
+ case PSFIMAGE:
+ call strcpy (PM_PSFIMAGE(pm), str, maxch)
+ case OBJLIST:
+ call strcpy (PM_OBJLIST(pm), str, maxch)
+ case KERNEL:
+ call strcpy (PM_KERNEL(pm), str, maxch)
+ case OUTIMAGE:
+ call strcpy (PM_OUTIMAGE(pm), str, maxch)
+ default:
+ call error (0, "RG_PSTATS: Unknown string parameter.")
+ }
+end
+
+
+# RG_PSETI -- Set the value of a psfmatch task integer parameter.
+
+procedure rg_pseti (pm, param, value)
+
+pointer pm # pointer to psfmatch structure
+int param # parameter to be fetched
+int value # value of the integer parameter
+
+begin
+ switch (param) {
+ case NREGIONS:
+ PM_NREGIONS(pm) = value
+ case CNREGION:
+ PM_CNREGION(pm) = value
+ case CENTER:
+ PM_CENTER(pm) = value
+ case BACKGRD:
+ PM_BACKGRD(pm) = value
+ switch (value) {
+ case PM_BNONE:
+ call strcpy ("none", PM_BSTRING(pm), SZ_FNAME)
+ case PM_BMEAN:
+ call strcpy ("mean", PM_BSTRING(pm), SZ_FNAME)
+ case PM_BMEDIAN:
+ call strcpy ("median", PM_BSTRING(pm), SZ_FNAME)
+ case PM_BSLOPE:
+ call strcpy ("plane", PM_BSTRING(pm), SZ_FNAME)
+ case PM_BNUMBER:
+ ;
+ default:
+ call strcpy ("none", PM_BSTRING(pm), SZ_FNAME)
+ }
+ case CONVOLUTION:
+ PM_CONVOLUTION(pm) = value
+ switch (value) {
+ case PM_CONIMAGE:
+ call strcpy ("image", PM_CSTRING(pm), SZ_FNAME)
+ case PM_CONPSF:
+ call strcpy ("psf", PM_CSTRING(pm), SZ_FNAME)
+ case PM_CONKERNEL:
+ call strcpy ("kernel", PM_CSTRING(pm), SZ_FNAME)
+ default:
+ call strcpy ("image", PM_CSTRING(pm), SZ_FNAME)
+ }
+ case DNX:
+ PM_DNX(pm) = value
+ case DNY:
+ PM_DNY(pm) = value
+ case PNX:
+ PM_PNX(pm) = value
+ case PNY:
+ PM_PNY(pm) = value
+ case KNX:
+ PM_KNX(pm) = value
+ case KNY:
+ PM_KNY(pm) = value
+ case POWER:
+ PM_POWER(pm) = value
+ case RADSYM:
+ PM_RADSYM(pm) = value
+ case NXFFT:
+ PM_NXFFT(pm) = value
+ case NYFFT:
+ PM_NYFFT(pm) = value
+ case FILTER:
+ PM_FILTER(pm) = value
+ switch (value) {
+ case PM_FNONE:
+ call strcpy ("none", PM_FSTRING(pm), SZ_FNAME)
+ case PM_FCOSBELL:
+ call strcpy ("cosbell", PM_FSTRING(pm), SZ_FNAME)
+ case PM_FREPLACE:
+ call strcpy ("replace", PM_FSTRING(pm), SZ_FNAME)
+ case PM_FMODEL:
+ call strcpy ("model", PM_FSTRING(pm), SZ_FNAME)
+ default:
+ call strcpy ("none", PM_FSTRING(pm), SZ_FNAME)
+ }
+ default:
+ call error (0, "RG_PSETI: Unknown integer parameter.")
+ }
+end
+
+
+# RG_PSETP -- Set the value of a psfmatch task pointer parameter.
+
+procedure rg_psetp (pm, param, value)
+
+pointer pm # pointer to psfmatch structure
+int param # parameter to be fetched
+pointer value # value of the pointer parameter
+
+begin
+ switch (param) {
+ case RC1:
+ PM_RC1(pm) = value
+ case RC2:
+ PM_RC2(pm) = value
+ case RL1:
+ PM_RL1(pm) = value
+ case RL2:
+ PM_RL2(pm) = value
+ case RZERO:
+ PM_RZERO(pm) = value
+ case RXSLOPE:
+ PM_RXSLOPE(pm) = value
+ case RYSLOPE:
+ PM_RYSLOPE(pm) = value
+ case REFFFT:
+ PM_REFFFT(pm) = value
+ case IMFFT:
+ PM_IMFFT(pm) = value
+ case FFT:
+ PM_FFT(pm) = value
+ case CONV:
+ PM_CONV(pm) = value
+ case ASFFT:
+ PM_ASFFT(pm) = value
+
+ default:
+ call error (0, "RG_PSETP: Unknown pointer parameter.")
+ }
+end
+
+
+# RG_PSETR -- Set the value of a psfmatch task real parameter.
+
+procedure rg_psetr (pm, param, value)
+
+pointer pm # pointer to psfmatch structure
+int param # parameter to be fetched
+real value # real parameter
+
+begin
+ switch (param) {
+ case BVALUER:
+ PM_BVALUER(pm) = value
+ case BVALUE:
+ PM_BVALUE(pm) = value
+ case LOREJECT:
+ PM_LOREJECT(pm) = value
+ case HIREJECT:
+ PM_HIREJECT(pm) = value
+ case APODIZE:
+ PM_APODIZE(pm) = value
+ case UFLUXRATIO:
+ PM_UFLUXRATIO(pm) = value
+ case FLUXRATIO:
+ PM_FLUXRATIO(pm) = value
+ case SXINNER:
+ PM_SXINNER(pm) = value
+ case SXOUTER:
+ PM_SXOUTER(pm) = value
+ case SYINNER:
+ PM_SYINNER(pm) = value
+ case SYOUTER:
+ PM_SYOUTER(pm) = value
+ case THRESHOLD:
+ PM_THRESHOLD(pm) = value
+ case NORMFACTOR:
+ PM_NORMFACTOR(pm) = value
+ default:
+ call error (0, "RG_PSETR: Unknown real parameter.")
+ }
+end
+
+
+# RG_PSETS -- Procedure to set the value of a string parameter.
+
+procedure rg_psets (pm, param, str)
+
+pointer pm # pointer to psfmatch structure
+int param # parameter to be fetched
+char str[ARB] # output string
+
+int index, ip
+pointer sp, temp
+real rval
+int strdic(), fnldir(), ctor()
+
+begin
+ call smark (sp)
+ call salloc (temp, SZ_LINE, TY_CHAR)
+
+ switch (param) {
+ case BSTRING:
+ ip = 1
+ index = strdic (str, str, SZ_LINE, PM_BTYPES)
+ if (index > 0) {
+ call strcpy (str, PM_BSTRING(pm), SZ_FNAME)
+ call rg_pseti (pm, BACKGRD, index)
+ } else if (ctor (str, ip, rval) > 0) {
+ call rg_psetr (pm, BVALUE, rval)
+ if (ctor (str, ip, rval) > 0) {
+ call rg_psetr (pm, BVALUER, rval)
+ call strcpy (str, PM_BSTRING(pm), SZ_FNAME)
+ call rg_pseti (pm, BACKGRD, PM_NUMBER)
+ } else {
+ call rg_psetr (pm, BVALUE, 0.0)
+ call rg_psetr (pm, BVALUER, 0.0)
+ }
+ }
+ case CSTRING:
+ index = strdic (str, str, SZ_LINE, PM_CTYPES)
+ if (index > 0) {
+ call strcpy (str, PM_CSTRING(pm), SZ_FNAME)
+ call rg_pseti (pm, CONVOLUTION, index)
+ }
+ case FSTRING:
+ index = strdic (str, str, SZ_LINE, PM_FTYPES)
+ if (index > 0) {
+ call strcpy (str, PM_FSTRING(pm), SZ_FNAME)
+ call rg_pseti (pm, FILTER, index)
+ }
+ case IMAGE:
+ call imgcluster (str, Memc[temp], SZ_FNAME)
+ index = fnldir (Memc[temp], PM_IMAGE(pm), SZ_FNAME)
+ call strcpy (Memc[temp+index], PM_IMAGE(pm), SZ_FNAME)
+ case REFIMAGE:
+ call imgcluster (str, Memc[temp], SZ_FNAME)
+ index = fnldir (Memc[temp], PM_REFIMAGE(pm), SZ_FNAME)
+ call strcpy (Memc[temp+index], PM_REFIMAGE(pm), SZ_FNAME)
+ case PSFDATA:
+ call strcpy (str, PM_PSFDATA(pm), SZ_FNAME)
+ case PSFIMAGE:
+ call imgcluster (str, Memc[temp], SZ_FNAME)
+ index = fnldir (Memc[temp], PM_PSFIMAGE(pm), SZ_FNAME)
+ call strcpy (Memc[temp+index], PM_PSFIMAGE(pm), SZ_FNAME)
+ case OBJLIST:
+ call strcpy (str, PM_OBJLIST(pm), SZ_FNAME)
+ case KERNEL:
+ call imgcluster (str, Memc[temp], SZ_FNAME)
+ index = fnldir (Memc[temp], PM_KERNEL(pm), SZ_FNAME)
+ call strcpy (Memc[temp+index], PM_KERNEL(pm), SZ_FNAME)
+ case OUTIMAGE:
+ call strcpy (str, PM_OUTIMAGE(pm), SZ_FNAME)
+ default:
+ call error (0, "RG_PSETS: Unknown string parameter.")
+ }
+
+ call sfree (sp)
+end
diff --git a/pkg/images/immatch/src/psfmatch/t_psfmatch.x b/pkg/images/immatch/src/psfmatch/t_psfmatch.x
new file mode 100644
index 00000000..182ac286
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/t_psfmatch.x
@@ -0,0 +1,365 @@
+include <fset.h>
+include <imhdr.h>
+include "psfmatch.h"
+
+# T_PSFMATCH -- Match the resolution of an image to that of a reference
+# image.
+
+procedure t_psfmatch ()
+
+pointer image1 # pointer to the input image name
+pointer imager # pointer to the reference image name
+pointer fpsflist # pointer to the regions list
+pointer image2 # pointer to the output image name
+pointer kernel # pointer to the kernel image name
+pointer pspectra # pointer to the fourier spectra image name
+int interactive # interactive mode ?
+int verbose # verbose mode ?
+int boundary # boundary extension type
+real constant # constant for boundary extension
+
+int list1, listr, psflist, listk, list2
+int nregions, newref, stat
+pointer sp, imtemp, str, pm, gd, id, imr, im1, impsf, imk, im2
+bool clgetb()
+int imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb(), clgwrd(), btoi()
+int rg_pstati(), rg_ptmpimage(), rg_pregions(), rg_psfm(), rg_pisfm()
+pointer gopen(), immap(), rg_pstatp()
+real clgetr()
+errchk fntopnb(), fntclsb()
+
+begin
+ call fseti (STDOUT, F_FLUSHNL, YES)
+
+ # Allocate temporary space.
+ call smark (sp)
+ call salloc (image1, SZ_FNAME, TY_CHAR)
+ call salloc (imager, SZ_FNAME, TY_CHAR)
+ call salloc (fpsflist, SZ_LINE, TY_CHAR)
+ call salloc (kernel, SZ_FNAME, TY_CHAR)
+ call salloc (image2, SZ_FNAME, TY_CHAR)
+ call salloc (pspectra, SZ_FNAME, TY_CHAR)
+ call salloc (imtemp, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Get task parameters.
+ call clgstr ("input", Memc[str], SZ_LINE)
+ list1 = imtopen (Memc[str])
+ call clgstr ("reference", Memc[str], SZ_LINE)
+ listr = imtopen (Memc[str])
+ call clgstr ("psfdata", Memc[fpsflist], SZ_LINE)
+ call clgstr ("kernel", Memc[str], SZ_LINE)
+ listk = imtopen (Memc[str])
+ call clgstr ("output", Memc[str], SZ_LINE)
+ list2 = imtopen (Memc[str])
+
+ # Open the psf matching fitting structure.
+ call rg_pgpars (pm)
+
+ # Will the task run in interactive mode?
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONKERNEL)
+ interactive = NO
+ else
+ interactive = btoi (clgetb ("interactive"))
+
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL) {
+ if (imtlen (listr) <= 0)
+ call error (0, "The reference image list is empty.")
+ if (imtlen (listr) > 1 && imtlen (listr) != imtlen (list1))
+ call error (0,
+ "The number of reference and input images is not the same.")
+ if (interactive == NO && Memc[fpsflist] == EOS) {
+ call error (0, "The objects list is empty.")
+ } else if (rg_pstati (pm, CONVOLUTION) == PM_CONIMAGE) {
+ psflist = fntopnb (Memc[fpsflist], NO)
+ if (fntlenb(psflist) > 0 && imtlen (listr) != fntlenb (psflist))
+ call error (0,
+ "The number of reference images and objects lists is not the same")
+ } else {
+ psflist = imtopen (Memc[fpsflist])
+ if (imtlen (list1) != imtlen (psflist))
+ call error (0,
+ "The number of input and psf images is not the same")
+ }
+ call rg_psets (pm, PSFDATA, Memc[fpsflist])
+ } else {
+ call imtclose (listr)
+ listr = NULL
+ psflist = NULL
+ call rg_psets (pm, PSFDATA, "")
+ }
+
+ # Compare the lengths of the input and output lists.
+ if (imtlen(listk) <= 0) {
+ call imtclose (listk)
+ listk = NULL
+ } else if (imtlen (list1) != imtlen (listk))
+ call error (0,
+ "The number of input and kernel images is not the same.")
+
+ if (imtlen (list2) <= 0) {
+ call imtclose (list2)
+ list2 = NULL
+ } else if (imtlen (list1) != imtlen (list2))
+ call error (0,
+ "The number of input and output images are not the same.")
+
+ # Get the boundary extension parameters for the image convolution.
+ boundary = clgwrd ("boundary", Memc[str], SZ_LINE,
+ "|constant|nearest|reflect|wrap|")
+ constant = clgetr ("constant")
+
+ if (interactive == YES) {
+ call clgstr ("graphics", Memc[str], SZ_FNAME)
+ iferr (gd = gopen (Memc[str], NEW_FILE, STDGRAPH))
+ gd = NULL
+ call clgstr ("display", Memc[str], SZ_FNAME)
+ iferr (id = gopen (Memc[str], APPEND, STDIMAGE))
+ id = NULL
+ verbose = YES
+ } else {
+ gd = NULL
+ id = NULL
+ verbose = btoi (clgetb ("verbose"))
+ }
+
+ imr = NULL
+ impsf = NULL
+
+ # Do each set of input and output images.
+ while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF)) {
+
+ # Open reference image and the associated objects file
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL) {
+ if (imtgetim (listr, Memc[imager], SZ_FNAME) != EOF) {
+ if (imr != NULL)
+ call imunmap (imr)
+ imr = immap (Memc[imager], READ_ONLY, 0)
+ if (IM_NDIM(imr) > 2)
+ call error (0, "Reference psf/image must be 1D or 2D")
+ call rg_psets (pm, REFIMAGE, Memc[imager])
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONIMAGE) {
+ nregions = rg_pregions (psflist, imr, pm, 1, NO)
+ if (nregions <= 0 && interactive == NO)
+ call error (0, "The objects list is empty.")
+ call rg_psets (pm, PSFIMAGE, "")
+ }
+ newref = YES
+ }
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONPSF) {
+ if (imtgetim (psflist, Memc[str], SZ_FNAME) != EOF) {
+ impsf = immap (Memc[str], READ_ONLY, 0)
+ if (IM_NDIM(impsf) != IM_NDIM(imr))
+ call error (0,
+ "Image and reference psf must have same dimensionality")
+ if (IM_LEN(impsf,1) != IM_LEN(imr,1))
+ call error (0,
+ "Image and reference psf are not the same size")
+ if (IM_NDIM(impsf) == 2 && (IM_LEN(impsf,2) !=
+ IM_LEN(imr,2)))
+ call error (0,
+ "Image and reference psf are not the same size")
+ call rg_psets (pm, PSFIMAGE, Memc[str])
+ newref = YES
+ }
+ }
+ } else {
+ imr = NULL
+ impsf = NULL
+ call rg_psets (pm, REFIMAGE, "")
+ call rg_psets (pm, PSFIMAGE, "")
+ call rg_psets (pm, OBJLIST, "")
+ newref = NO
+ }
+
+ # Open the input image.
+ im1 = immap (Memc[image1], READ_ONLY, 0)
+ if (IM_NDIM(im1) > 2) {
+ call error (0, "Input image must be 1D or 2D")
+ } else if (imr != NULL) {
+ if (IM_NDIM(im1) != IM_NDIM(imr))
+ call error (0,
+ "Input and reference images must have same dimensionality")
+ }
+ call rg_psets (pm, IMAGE, Memc[image1])
+
+ # Open the kernel image name.
+ if (listk != NULL) {
+ if (imtgetim (listk, Memc[kernel], SZ_FNAME) != EOF)
+ ;
+ } else {
+ if (rg_ptmpimage (Memc[image1], "ker", "ker", Memc[kernel],
+ SZ_FNAME) == NO)
+ ;
+ }
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL)
+ imk = immap (Memc[kernel], NEW_IMAGE, 0)
+ else
+ imk = immap (Memc[kernel], READ_ONLY, 0)
+ call rg_psets (pm, KERNEL, Memc[kernel])
+
+
+ # Construct the output image name.
+ if (list2 == NULL) {
+ im2 = NULL
+ Memc[image2] = NULL
+ } else if (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF) {
+ call xt_mkimtemp (Memc[image1], Memc[image2], Memc[imtemp],
+ SZ_FNAME)
+ im2 = immap (Memc[image2], NEW_COPY, im1)
+ } else {
+ im2 = NULL
+ Memc[image2] = NULL
+ }
+ call rg_psets (pm, OUTIMAGE, Memc[image2])
+
+ # Compute the the psf matching kernel.
+ if (interactive == YES) {
+ stat = rg_pisfm (pm, imr, psflist, impsf, im1, imk, NULL, im2,
+ gd, id)
+ } else {
+ if (rg_psfm (pm, imr, im1, impsf, imk, newref) == OK) {
+ if (verbose == YES) {
+ call printf (
+ "Completed computing/reading kernel %s for image %s\n")
+ call pargstr (Memc[kernel])
+ call pargstr (Memc[image1])
+ if (rg_pstati(pm, CONVOLUTION) != PM_CONKERNEL)
+ call rg_pwrite (pm, imk, NULL)
+ }
+ } else {
+ if (verbose == YES) {
+ call printf (
+ "Error computing/reading kernel %s for image %s\n")
+ call pargstr (Memc[kernel])
+ call pargstr (Memc[image1])
+ }
+ }
+ stat = NO
+ }
+
+ # Convolve the image.
+ if (im2 != NULL && stat == NO) {
+ if (verbose == YES) {
+ if (rg_pstatp(pm, CONV) != NULL)
+ call printf (
+ "\tComputing matched image %s ...\n")
+ else
+ call printf (
+ "\tComputing matched image %s ...\n")
+ call pargstr (Memc[imtemp])
+ call pargstr (Memc[kernel])
+ }
+ if (rg_pstatp(pm, CONV) != NULL)
+ call rg_pconvolve (im1, im2, Memr[rg_pstatp(pm,CONV)],
+ rg_pstati(pm,KNX), rg_pstati(pm,KNY), boundary,
+ constant)
+ }
+
+ # Close up the images.
+ if (im2 != NULL) {
+ call imunmap (im2)
+ if (rg_pstatp(pm, CONV) == NULL)
+ call imdelete (Memc[image2])
+ else
+ call xt_delimtemp (Memc[image2], Memc[imtemp])
+ }
+ if (impsf != NULL)
+ call imunmap (impsf)
+ if (imk != NULL) {
+ call imunmap (imk)
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL &&
+ rg_pstatp(pm, CONV) == NULL)
+ call imdelete (Memc[kernel])
+ }
+ call imunmap (im1)
+
+ if (stat == YES)
+ break
+ newref = NO
+ }
+
+ # Close up the lists.
+ if (imr != NULL)
+ call imunmap (imr)
+
+ if (list2 != NULL)
+ call imtclose (list2)
+ if (listk != NULL)
+ call imtclose (listk)
+ if (psflist != NULL) {
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONIMAGE)
+ call fntclsb (psflist)
+ else
+ call imtclose (psflist)
+ }
+ if (listr != NULL)
+ call imtclose (listr)
+ call imtclose (list1)
+
+ call rg_pfree (pm)
+
+ # Close up te graphics and the display.
+ if (gd != NULL)
+ call gclose (gd)
+ if (id != NULL)
+ call gclose (id)
+
+ call sfree (sp)
+end
+
+
+# RG_PTMPIMAGE -- Generate either a permanent image name using a user specified
+# prefix or temporary image name using a default prefix. Return NO if the
+# image is temporary or YES if it is permanent.
+
+int procedure rg_ptmpimage (image, prefix, tmp, name, maxch)
+
+char image[ARB] #I image name
+char prefix[ARB] #I user supplied prefix
+char tmp[ARB] #I user supplied temporary root
+char name[ARB] #O output name
+int maxch #I max number of chars
+
+int npref, ndir
+int fnldir(), rg_pimroot(), strlen()
+
+begin
+ npref = strlen (prefix)
+ ndir = fnldir (prefix, name, maxch)
+ if (npref == ndir) {
+ call mktemp (tmp, name[ndir+1], maxch)
+ return (NO)
+ } else {
+ call strcpy (prefix, name, npref)
+ if (rg_pimroot (image, name[npref+1], maxch) <= 0)
+ ;
+ return (YES)
+ }
+end
+
+
+# RG_PIMROOT -- Fetch the root image name minus the directory specification
+# and the section notation. The length of the root name is returned.
+
+int procedure rg_pimroot (image, root, maxch)
+
+char image[ARB] #I image specification
+char root[ARB] #O rootname
+int maxch #I maximum number of characters
+
+int nchars
+pointer sp, str
+int fnldir(), strlen()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ call imgimage (image, root, maxch)
+ nchars = fnldir (root, Memc[str], maxch)
+ call strcpy (root[nchars+1], root, maxch)
+
+ call sfree (sp)
+ return (strlen (root))
+end
generated by cgit (git 2.34.1) at 2025-09-20 08:55:13 -0400