From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- pkg/images/immatch/src/xregister/mkpkg | 25 + pkg/images/immatch/src/xregister/oxregister.key | 33 + pkg/images/immatch/src/xregister/rgxbckgrd.x | 63 ++ pkg/images/immatch/src/xregister/rgxcolon.x | 508 +++++++++++ pkg/images/immatch/src/xregister/rgxcorr.x | 1034 +++++++++++++++++++++++ pkg/images/immatch/src/xregister/rgxdbio.x | 290 +++++++ pkg/images/immatch/src/xregister/rgxfft.x | 179 ++++ pkg/images/immatch/src/xregister/rgxfit.x | 814 ++++++++++++++++++ pkg/images/immatch/src/xregister/rgxgpars.x | 68 ++ pkg/images/immatch/src/xregister/rgxicorr.x | 583 +++++++++++++ pkg/images/immatch/src/xregister/rgximshift.x | 391 +++++++++ pkg/images/immatch/src/xregister/rgxplot.x | 317 +++++++ pkg/images/immatch/src/xregister/rgxppars.x | 49 ++ pkg/images/immatch/src/xregister/rgxregions.x | 459 ++++++++++ pkg/images/immatch/src/xregister/rgxshow.x | 172 ++++ pkg/images/immatch/src/xregister/rgxtools.x | 685 +++++++++++++++ pkg/images/immatch/src/xregister/rgxtransform.x | 446 ++++++++++ pkg/images/immatch/src/xregister/t_xregister.x | 440 ++++++++++ pkg/images/immatch/src/xregister/xregister.h | 250 ++++++ pkg/images/immatch/src/xregister/xregister.key | 47 ++ 20 files changed, 6853 insertions(+) create mode 100644 pkg/images/immatch/src/xregister/mkpkg create mode 100644 pkg/images/immatch/src/xregister/oxregister.key create mode 100644 pkg/images/immatch/src/xregister/rgxbckgrd.x create mode 100644 pkg/images/immatch/src/xregister/rgxcolon.x create mode 100644 pkg/images/immatch/src/xregister/rgxcorr.x create mode 100644 pkg/images/immatch/src/xregister/rgxdbio.x create mode 100644 pkg/images/immatch/src/xregister/rgxfft.x create mode 100644 pkg/images/immatch/src/xregister/rgxfit.x create mode 100644 pkg/images/immatch/src/xregister/rgxgpars.x create mode 100644 pkg/images/immatch/src/xregister/rgxicorr.x create mode 100644 pkg/images/immatch/src/xregister/rgximshift.x create mode 100644 pkg/images/immatch/src/xregister/rgxplot.x create mode 100644 pkg/images/immatch/src/xregister/rgxppars.x create mode 100644 pkg/images/immatch/src/xregister/rgxregions.x create mode 100644 pkg/images/immatch/src/xregister/rgxshow.x create mode 100644 pkg/images/immatch/src/xregister/rgxtools.x create mode 100644 pkg/images/immatch/src/xregister/rgxtransform.x create mode 100644 pkg/images/immatch/src/xregister/t_xregister.x create mode 100644 pkg/images/immatch/src/xregister/xregister.h create mode 100644 pkg/images/immatch/src/xregister/xregister.key (limited to 'pkg/images/immatch/src/xregister') diff --git a/pkg/images/immatch/src/xregister/mkpkg b/pkg/images/immatch/src/xregister/mkpkg new file mode 100644 index 00000000..262b721d --- /dev/null +++ b/pkg/images/immatch/src/xregister/mkpkg @@ -0,0 +1,25 @@ +# Make the XREGISTER task + +$checkout libpkg.a ../../../ +$update libpkg.a +$checkin libpkg.a ../../../ +$exit + +libpkg.a: + rgxbckgrd.x "xregister.h" + rgxcolon.x "xregister.h" + rgxcorr.x "xregister.h" + rgxdbio.x "xregister.h" + rgxfft.x + rgxfit.x "xregister.h" + rgxgpars.x "xregister.h" + rgxicorr.x "xregister.h" + rgximshift.x + rgxplot.x + rgxppars.x "xregister.h" + rgxregions.x "xregister.h" + rgxshow.x "xregister.h" + rgxtools.x "xregister.h" + rgxtransform.x "xregister.h" + t_xregister.x "xregister.h" + ; diff --git a/pkg/images/immatch/src/xregister/oxregister.key b/pkg/images/immatch/src/xregister/oxregister.key new file mode 100644 index 00000000..91064ff8 --- /dev/null +++ b/pkg/images/immatch/src/xregister/oxregister.key @@ -0,0 +1,33 @@ + Xregister Image Overlay Sub-menu + + +? Print help +c Overlay the marked column of the reference image + with the same column of the input image +l Overlay the marked line of the reference image + with the sname line of the input image +x Overlay the marked column of the reference image + with the x and y lagged column of the input image +y Overlay the marked line of the reference image + with the x and y lagged line of the input image +v Overlay the marked column of the reference image + with the x and y shifted column of the input image +h Overlay the marked line of the reference image + with the x and y shifted line of the input image +q Quit + + + Image Overlay Sub-menu Colon Commands + +:c [m] [n] Overlay the middle [mth] column of the reference image + with the mth [nth] column of the input image +:l [m] [n] Overlay the middle [mth] line of the reference image + with the mth [nth] line of the input image +:x [m] Overlay the middle [mth] column of the reference image + with the x and y lagged column of the input image +:y [m] Overlay the middle [mth] line of the reference image + with the x and y lagged line of the input image +:v [m] Overlay the middle [mth] column of the reference image + with the x and y shifted column of the input image +:h [m] Overlay the middle [mth] line of the reference image + with the x and y shifted line of the input image diff --git a/pkg/images/immatch/src/xregister/rgxbckgrd.x b/pkg/images/immatch/src/xregister/rgxbckgrd.x new file mode 100644 index 00000000..c9747ee6 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxbckgrd.x @@ -0,0 +1,63 @@ +include +include "xregister.h" + +# RG_XSCALE -- Compute the background offset and x and y slope. + +procedure rg_xscale (xc, data, npts, nx, ny, offset, coeff) + +pointer xc #I pointer to the cross-correlation function +real data[ARB] #I the input data +int npts #I the number of points +int nx, ny #I the dimensions of the original subraster +real offset #I the input offset +real coeff[ARB] #O the output coefficients + +int wborder +pointer gs +real loreject, hireject, zero +int rg_xstati(), rg_znsum(), rg_znmedian(), rg_slope() +real rg_xstatr() + +begin + loreject = rg_xstatr (xc, LOREJECT) + hireject = rg_xstatr (xc, HIREJECT) + wborder = rg_xstati (xc, BORDER) + + switch (rg_xstati (xc, BACKGRD)) { + case XC_BNONE: + coeff[1] = offset + coeff[2] = 0.0 + coeff[3] = 0.0 + case XC_MEAN: + if (rg_znsum (data, npts, zero, loreject, hireject) <= 0) + zero = 0.0 + coeff[1] = zero + coeff[2] = 0.0 + coeff[3] = 0.0 + case XC_MEDIAN: + if (rg_znmedian (data, npts, zero, loreject, hireject) <= 0) + zero = 0.0 + coeff[1] = zero + coeff[2] = 0.0 + coeff[3] = 0.0 + case XC_SLOPE: + call gsinit (gs, GS_POLYNOMIAL, 2, 2, GS_XNONE, 1.0, real (nx), 1.0, + real (ny)) + if (rg_slope (gs, data, npts, nx, ny, wborder, wborder, 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] = offset + coeff[2] = 0.0 + coeff[3] = 0.0 + } +end diff --git a/pkg/images/immatch/src/xregister/rgxcolon.x b/pkg/images/immatch/src/xregister/rgxcolon.x new file mode 100644 index 00000000..cb007473 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxcolon.x @@ -0,0 +1,508 @@ +include +include +include +include "xregister.h" + +# RG_XCOLON-- Procedure to process colon commands for setting the cross- +# correlation parameters. + +procedure rg_xcolon (gd, xc, imr, im1, im2, db, dformat, tfd, reglist, cmdstr, + newdata, newcross, newcenter) + +pointer gd #I pointer to the graphics stream +pointer xc #I pointer to cross-correlation structure +pointer imr #I/O pointer to the reference image +pointer im1 #I/O pointer to the input image +pointer im2 #I/O pointer to the output image +pointer db #I/O pointer to the shifts database file +int dformat #I is the shifts file in database format +int tfd #I/O the transformations file descriptor +pointer reglist #I/O pointer to the regions list +char cmdstr[ARB] #I input command string +int newdata #I/O new input data +int newcross #I/O new cross-correlation function flag +int newcenter #I/O new cross-correlation peak flag + +bool streq() +int ncmd, creg, nreg, ival, stat +pointer sp, cmd, str +real rval +int strdic(), open(), nscan(), rg_xstati(), fntopnb() +int rg_xregions(), rg_xmkregions(), strlen() +pointer immap(), dtmap(), rg_xstatp() +real rg_xstatr() +errchk immap(), dtmap(), open(), 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, XCMDS) + switch (ncmd) { + case XCMD_REFIMAGE: + call gargwrd (Memc[cmd], SZ_LINE) + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_FNAME) + if (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 ( + "Image has the wrong number of dimensions\n") + call imunmap (imr) + imr = immap (Memc[str], READ_ONLY, 0) + } else { + call rg_xsets (xc, REFIMAGE, Memc[cmd]) + newdata = YES; newcross = YES; newcenter = YES + } + } + + case XCMD_IMAGE: + call gargwrd (Memc[cmd], SZ_LINE) + call rg_xstats (xc, 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 ( + "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_xsets (xc, IMAGE, Memc[cmd]) + newdata = YES; newcross = YES; newcenter = YES + } + } + + case XCMD_OUTIMAGE: + call gargwrd (Memc[cmd], SZ_LINE) + call rg_xstats (xc, 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_xsets (xc, OUTIMAGE, Memc[cmd]) + } + } + + case XCMD_DATABASE: + call gargwrd (Memc[cmd], SZ_LINE) + call rg_xstats (xc, DATABASE, Memc[str], SZ_FNAME) + if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) { + call printf ("%s: %s\n") + call pargstr (KY_DATABASE) + call pargstr (Memc[str]) + } else { + if (db != NULL) { + if (dformat == YES) + call dtunmap (db) + else + call close (db) + db = NULL + } + iferr { + if (dformat == YES) + db = dtmap (Memc[cmd], APPEND) + else + db = open (Memc[cmd], NEW_FILE, TEXT_FILE) + } then { + call erract (EA_WARN) + if (dformat == YES) + db = dtmap (Memc[str], APPEND) + else + db = open (Memc[str], APPEND, TEXT_FILE) + } else { + call rg_xsets (xc, DATABASE, Memc[cmd]) + } + } + + CASE XCMD_RECORD: + call gargstr (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call rg_xstats (xc, RECORD, Memc[str], SZ_FNAME) + call printf ("%s: %s\n") + call pargstr (KY_RECORD) + call pargstr (Memc[str]) + } else + call rg_xsets (xc, RECORD, Memc[cmd]) + + case XCMD_CREGION: + + call gargi (nreg) + creg = rg_xstati (xc, CREGION) + + if (nscan() == 1 || (nreg == creg)) { + call printf ("%s: %d/%d") + call pargstr (KY_CREGION) + call pargi (creg) + call pargi (rg_xstati (xc, NREGIONS)) + call printf (" [%d:%d,%d:%d]\n") + call pargi (Memi[rg_xstatp (xc,RC1)+creg-1]) + call pargi (Memi[rg_xstatp (xc,RC2)+creg-1]) + call pargi (Memi[rg_xstatp (xc,RL1)+creg-1]) + call pargi (Memi[rg_xstatp (xc,RL2)+creg-1]) + + } else { + if (nreg < 1 || nreg > rg_xstati (xc,NREGIONS)) { + call printf ("Region %d is out of range\n") + call pargi (nreg) + } else { + call printf ( + "Setting current region to %d: [%d:%d,%d:%d]\n") + call pargi (nreg) + call pargi (Memi[rg_xstatp (xc,RC1)+nreg-1]) + call pargi (Memi[rg_xstatp (xc,RC2)+nreg-1]) + call pargi (Memi[rg_xstatp (xc,RL1)+nreg-1]) + call pargi (Memi[rg_xstatp (xc,RL2)+nreg-1]) + call rg_xseti (xc, CREGION, nreg) + newdata = YES; newcross = YES; newcenter = YES + } + + } + + case XCMD_REGIONS: + + call gargwrd (Memc[cmd], SZ_LINE) + call rg_xstats (xc, REGIONS, Memc[str], SZ_FNAME) + if (nscan() == 1 || streq (Memc[cmd], Memc[str]) || Memc[cmd] == + EOS) { + call printf ("%s [string/file]: %s\n") + call pargstr (KY_REGIONS) + call pargstr (Memc[str]) + } else { + if (reglist != NULL) { + call fntclsb (reglist) + reglist = NULL + } + iferr (reglist = fntopnb (Memc[cmd], NO)) + reglist = NULL + if (rg_xregions (reglist, imr, xc, 1) > 0) { + call rg_xseti (xc, CREGION, 1) + call rg_xsets (xc, REGIONS, Memc[cmd]) + newdata = YES; newcross = YES; newcenter = YES + } else { + if (reglist != NULL) { + call fntclsb (reglist) + reglist = NULL + } + iferr (reglist = fntopnb (Memc[str], NO)) + reglist = NULL + if (rg_xregions (reglist, imr, xc, 1) > 0) + ; + call rg_xsets (xc, REGIONS, Memc[str]) + call rg_xseti (xc, CREGION, 1) + } + } + + case XCMD_REFFILE: + + call gargwrd (Memc[cmd], SZ_LINE) + call rg_xstats (xc, REFFILE, Memc[str], SZ_FNAME) + if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) { + call printf ("%s: %s\n") + call pargstr (KY_REFFILE) + call pargstr (Memc[str]) + } else { + if (tfd != NULL) { + call close (tfd) + tfd = NULL + } + iferr { + tfd = open (Memc[cmd], READ_ONLY, TEXT_FILE) + } then { + tfd = NULL + call erract (EA_WARN) + call rg_xsets (xc, REFFILE, "") + call printf ("Coords file is undefined.\n") + } else + call rg_xsets (xc, REFFILE, Memc[cmd]) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_XLAG: + call gargi (ival) + if (nscan () == 1) { + call printf ("%s = %d\n") + call pargstr (KY_XLAG) + call pargi (rg_xstati (xc, XLAG)) + } else { + call rg_xseti (xc, XLAG, ival) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_YLAG: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_YLAG) + call pargi (rg_xstati (xc, YLAG)) + } else { + call rg_xseti (xc, YLAG, ival) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_DXLAG: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_DXLAG) + call pargi (rg_xstati (xc, DXLAG)) + } else { + call rg_xseti (xc, DXLAG, ival) + } + + case XCMD_DYLAG: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_DYLAG) + call pargi (rg_xstati (xc, DYLAG)) + } else { + call rg_xseti (xc, DYLAG, ival) + } + + case XCMD_BACKGROUND: + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] != EOS) + call strcat (" ", Memc[cmd], SZ_LINE) + call gargwrd (Memc[cmd+strlen(Memc[cmd])], SZ_LINE) + if (Memc[cmd] == EOS) { + call rg_xstats (xc, BSTRING, Memc[str], SZ_FNAME) + call printf ("%s: %s\n") + call pargstr (KY_BACKGROUND) + call pargstr (Memc[str]) + } else { + call rg_xsets (xc, BSTRING, Memc[cmd]) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_BORDER: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_BORDER) + call pargi (rg_xstati (xc, BORDER)) + } else { + call rg_xseti (xc, BORDER, ival) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_LOREJECT: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_LOREJECT) + call pargr (rg_xstatr (xc, LOREJECT)) + } else { + call rg_xsetr (xc, LOREJECT, rval) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_HIREJECT: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_HIREJECT) + call pargr (rg_xstatr (xc, HIREJECT)) + } else { + call rg_xsetr (xc, HIREJECT, rval) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_APODIZE: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_APODIZE) + call pargr (rg_xstatr (xc, APODIZE)) + } else { + call rg_xsetr (xc, APODIZE, max (0.0, min (rval, 0.50))) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_CORRELATION: + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call rg_xstats (xc, CSTRING, Memc[str], SZ_FNAME) + call printf ("%s = %s\n") + call pargstr (KY_CORRELATION) + call pargstr (Memc[str]) + } else { + stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, XC_CTYPES) + if (stat > 0) { + call rg_xseti (xc, CFUNC, stat) + call rg_xsets (xc, CSTRING, Memc[cmd]) + newcross = YES; newcenter = YES + } + } + + case XCMD_XWINDOW: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_XWINDOW) + call pargi (rg_xstati (xc, XWINDOW)) + } else { + call rg_xseti (xc, XWINDOW, ival) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_YWINDOW: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_YWINDOW) + call pargi (rg_xstati (xc, YWINDOW)) + } else { + call rg_xseti (xc, YWINDOW, ival) + newdata = YES; newcross = YES; newcenter = YES + } + + case XCMD_PEAKCENTER: + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call rg_xstats (xc, PSTRING, Memc[str], SZ_FNAME) + call printf ("%s: %s\n") + call pargstr (KY_PEAKCENTER) + call pargstr (Memc[str]) + } else { + stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, XC_PTYPES) + if (stat > 0) { + call rg_xseti (xc, PFUNC, stat) + call rg_xsets (xc, PSTRING, Memc[cmd]) + newcenter = YES + } + } + + case XCMD_XCBOX: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_XCBOX) + call pargi (rg_xstati (xc, XCBOX)) + } else { + if (mod (ival, 2) == 0) + ival = ival + 1 + call rg_xseti (xc, XCBOX, ival) + newcenter = YES + } + + case XCMD_YCBOX: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_YCBOX) + call pargi (rg_xstati (xc, YCBOX)) + } else { + if (mod (ival, 2) == 0) + ival = ival + 1 + call rg_xseti (xc, YCBOX, ival) + newcenter = YES + } + + case XCMD_SHOW: + call gdeactivate (gd, 0) + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, XSHOW) + switch (ncmd) { + case XSHOW_DATA: + call rg_xnshow (xc) + case XSHOW_BACKGROUND: + call rg_xbshow (xc) + case XSHOW_CORRELATION: + call rg_xxshow (xc) + case XSHOW_PEAKCENTER: + call rg_xpshow (xc) + default: + call rg_xshow (xc) + } + call greactivate (gd, 0) + + case XCMD_MARK: + call gdeactivate (gd, 0) + if (reglist != NULL) { + call fntclsb (reglist) + reglist = NULL + } + if (rg_xmkregions (imr, xc, 1, MAX_NREGIONS, Memc[str], + SZ_LINE) <= 0) { + call rg_xstats (xc, REGIONS, Memc[str], SZ_LINE) + iferr (reglist = fntopnb (Memc[str], NO)) + reglist = NULL + if (rg_xregions (reglist, imr, xc, 1) > 0) + ; + call rg_xsets (xc, REGIONS, Memc[str]) + call rg_xseti (xc, CREGION, 1) + } else { + call rg_xseti (xc, CREGION, 1) + call rg_xsets (xc, REGIONS, Memc[str]) + newdata = YES; newcross = YES; newcenter = YES + } + 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/xregister/rgxcorr.x b/pkg/images/immatch/src/xregister/rgxcorr.x new file mode 100644 index 00000000..a708bf7a --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxcorr.x @@ -0,0 +1,1034 @@ +include +include +include +include "xregister.h" + +# RG_XCORR -- Compute the shift shift for an image relative to a reference +# image using cross-correlation techniques. + +int procedure rg_xcorr (imr, im1, db, dformat, xc) + +pointer imr #I pointer to the reference image +pointer im1 #I pointer to the input image +pointer db #I pointer to the shifts database +int dformat #I write shifts file in database format ? +pointer xc #I pointer to the cross-correlation structure + +pointer sp, image, imname +real xshift, yshift +bool streq() +int rg_xstati(), fscan(), nscan() +errchk rg_cross(), rg_xfile() + +begin + # Allocate working space. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (imname, SZ_FNAME, TY_CHAR) + call rg_xstats (xc, IMAGE, Memc[image], SZ_FNAME) + + # Initialize. + xshift = 0.0 + yshift = 0.0 + + # Compute the average shift for the image. + switch (rg_xstati (xc, CFUNC)) { + case XC_DISCRETE, XC_DIFFERENCE, XC_FOURIER: + + # Write out the parameters. + if (dformat == YES) + call rg_xdbparams (db, xc) + + # Compute the cross-correlation function. + call rg_cross (imr, im1, xc, NULL, xshift, yshift) + call rg_xsetr (xc, TXSHIFT, xshift) + call rg_xsetr (xc, TYSHIFT, yshift) + + # Write out the results for the individual regions. + if (dformat == YES) + call rg_xwreg (db, xc) + + # Write out the total shifts. + if (dformat == YES) + call rg_xdbshift (db, xc) + else { + call fprintf (db, "%s %g %g\n") + call pargstr (Memc[image]) + call pargr (xshift) + call pargr (yshift) + } + + # Set the x and y lags for the next picture. + if (rg_xstati (xc, NREFPTS) > 0) { + call rg_xseti (xc, XLAG, 0) + call rg_xseti (xc, YLAG, 0) + } else if (IS_INDEFI (rg_xstati (xc, DXLAG)) || + IS_INDEFI (rg_xstati (xc, DYLAG))) { + call rg_xseti (xc, XLAG, nint (-xshift)) + call rg_xseti (xc, YLAG, nint (-yshift)) + } else { + call rg_xseti (xc, XLAG, rg_xstati (xc, XLAG) + rg_xstati (xc, + DXLAG)) + call rg_xseti (xc, YLAG, rg_xstati (xc, YLAG) + rg_xstati (xc, + DYLAG)) + } + + case XC_FILE: + if (dformat == YES) + call rg_xfile (db, xc, xshift, yshift) + else { + if (fscan (db) != EOF) { + call gargwrd (Memc[imname], SZ_FNAME) + call gargr (xshift) + call gargr (yshift) + if (! streq (Memc[imname], Memc[image]) || nscan() != 3) { + xshift = 0.0 + yshift = 0.0 + } + } else { + xshift = 0.0 + yshift = 0.0 + } + } + call rg_xsetr (xc, TXSHIFT, xshift) + call rg_xsetr (xc, TYSHIFT, yshift) + + default: + call error (0, "The correlation function is undefined.") + } + + call sfree (sp) + + return (NO) +end + + +# RG_CROSS -- Compute the cross-correlation function for all the regions +# using discrete, fourier, or difference techniques and compute the position +# of its peak using one of several centering algorithms. + +procedure rg_cross (imr, im1, xc, gd, xavshift, yavshift) + +pointer imr #I pointer to the reference image +pointer im1 #I pointer to the input image +pointer xc #I pointer to the cross correlation structure +pointer gd #I pointer to graphics stream +real xavshift #O x coord shift +real yavshift #O y coord shift + +int i, nregions, ngood +pointer pxshift, pyshift +real xshift, yshift +int rg_xstati(), rg_xcget(), rg_xfget() +pointer rg_xstatp() + +begin + # Get the pointers. + pxshift = rg_xstatp (xc, XSHIFTS) + pyshift = rg_xstatp (xc, YSHIFTS) + nregions = rg_xstati (xc, NREGIONS) + + # Loop over the regions. + xavshift = 0.0 + yavshift = 0.0 + ngood = 0 + do i = 1, nregions { + + # Compute the cross_correlation function. + switch (rg_xstati (xc, CFUNC)) { + case XC_DISCRETE, XC_DIFFERENCE: + if (rg_xcget (xc, imr, im1, i) == ERR) { + Memr[pxshift+i-1] = INDEFR + Memr[pyshift+i-1] = INDEFR + if (rg_xstatp (xc, XCOR) != NULL) + call mfree (rg_xstatp (xc, XCOR), TY_REAL) + call rg_xsetp (xc, XCOR, NULL) + next + } + case XC_FOURIER: + if (rg_xfget (xc, imr, im1, i) == ERR) { + Memr[pxshift+i-1] = INDEFR + Memr[pyshift+i-1] = INDEFR + if (rg_xstatp (xc, XCOR) != NULL) + call mfree (rg_xstatp (xc, XCOR), TY_REAL) + call rg_xsetp (xc, XCOR, NULL) + next + } + default: + call error (0, "The correlation function is undefined") + } + + # Find the peak of the cross-correlation function. + call rg_fit (xc, i, gd, xshift, yshift) + + # Accumulate the shifts. + xavshift = xavshift + xshift + yavshift = yavshift + yshift + ngood = ngood + 1 + } + + # Compute the average shift. + if (ngood > 0) { + xavshift = xavshift / ngood + yavshift = yavshift / ngood + } +end + + +# RG_XFILE -- Read the average x and y shifts from the shifts database. + +procedure rg_xfile (db, xc, xshift, yshift) + +pointer db #I pointer to the database +pointer xc #I pointer to the cross correlation structure +real xshift #O shift in x +real yshift #O shift in y + +int rec +pointer sp, str +int dtlocate() +real dtgetr() +errchk dtlocate(), dtgetr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + call rg_xstats (xc, RECORD, Memc[str], SZ_LINE) + iferr { + rec = dtlocate (db, Memc[str]) + xshift = dtgetr (db, rec, "xshift") + yshift = dtgetr (db, rec, "yshift") + } then { + xshift = 0.0 + yshift = 0.0 + } + + call sfree (sp) +end + + +# RG_ICROSS -- Compute the cross-correlation function for a given region. + +int procedure rg_icross (xc, imr, im1, nreg) + +pointer xc #I pointer to the cross-correlation structure +pointer imr #I pointer to the reference image +pointer im1 #I pointer to the input image +int nreg #I the index of the current region + +int stat +pointer pxshift, pyshift +int rg_xstati(), rg_xcget(), rg_xfget() +pointer rg_xstatp() + +begin + pxshift = rg_xstatp (xc, XSHIFTS) + pyshift = rg_xstatp (xc, YSHIFTS) + + switch (rg_xstati (xc, CFUNC)) { + case XC_DISCRETE, XC_DIFFERENCE: + stat = rg_xcget (xc, imr, im1, nreg) + if (stat == ERR) { + Memr[pxshift+nreg-1] = INDEFR + Memr[pyshift+nreg-1] = INDEFR + if (rg_xstatp (xc, XCOR) != NULL) + call mfree (rg_xstatp (xc, XCOR), TY_REAL) + call rg_xsetp (xc, XCOR, NULL) + } + case XC_FOURIER: + stat = rg_xfget (xc, imr, im1, nreg) + if (stat == ERR) { + Memr[pxshift+nreg-1] = INDEFR + Memr[pyshift+nreg-1] = INDEFR + if (rg_xstatp (xc, XCOR) != NULL) + call mfree (rg_xstatp (xc, XCOR), TY_REAL) + call rg_xsetp (xc, XCOR, NULL) + } + case XC_FILE: + stat = OK + } + + return (stat) +end + + +# RG_XCGET -- Compute the convolution using the discrete or difference +# correlation functions. + +int procedure rg_xcget (xc, imr, im1, i) + +pointer xc #I pointer to the cross-correlation structure +pointer imr #I pointer to the reference image +pointer im1 #I pointer to input image image +int i #I index of region + +int stat, xwindow, ywindow, nrimcols, nrimlines, nimcols, nimlines +int nrcols, nrlines, ncols, nlines +int xlag, ylag, nborder, rc1, rc2, rl1, rl2, c1, c2, l1, l2 +pointer sp, str, coeff, rbuf, ibuf, xcor +pointer prc1, prc2, prl1, prl2, przero, prxslope, pryslope, border +real rxlag, rylag +int rg_xstati(), rg_border() +pointer rg_xstatp(), rg_ximget() +real rg_xstatr() + +define nextregion_ 10 + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (coeff, max (GS_SAVECOEFF + 6, 9), TY_REAL) + rbuf = NULL + ibuf = NULL + + # Check for regions. + if (i > rg_xstati (xc, NREGIONS)) { + stat = ERR + goto nextregion_ + } + + # Get the image sizes. + nrimcols = IM_LEN(imr,1) + if (IM_NDIM(imr) == 1) + nrimlines = 1 + else + nrimlines = IM_LEN(imr,2) + nimcols = IM_LEN(im1,1) + if (IM_NDIM(im1) == 1) + nimlines = 1 + else + nimlines = IM_LEN(im1,2) + + # Get the reference region pointers. + prc1 = rg_xstatp (xc, RC1) + prc2 = rg_xstatp (xc, RC2) + prl1 = rg_xstatp (xc, RL1) + prl2 = rg_xstatp (xc, RL2) + przero = rg_xstatp (xc, RZERO) + prxslope = rg_xstatp (xc, RXSLOPE) + pryslope = rg_xstatp (xc, RYSLOPE) + + # Compute the reference region limits. + rc1 = max (1, min (int (nrimcols), Memi[prc1+i-1])) + rc2 = min (int (nrimcols), max (1, Memi[prc2+i-1])) + rl1 = max (1, min (int (nrimlines), Memi[prl1+i-1])) + rl2 = min (int (nrimlines), max (1, Memi[prl2+i-1])) + nrcols = rc2 - rc1 + 1 + nrlines = rl2 - rl1 + 1 + + # Move to the next reference region if current region is off the image. + if (rc1 > nrimcols || rc2 < 1 || rl1 > nrimlines || rl2 < 1) { + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE) + call eprintf ( + "Reference section: %s[%d:%d,%d:%d] is off image.\n") + call pargstr (Memc[str]) + call pargi (rc1) + call pargi (rc2) + call pargi (rl1) + call pargi (rl2) + stat = ERR + goto nextregion_ + } + + # Check the window sizes. + xwindow = rg_xstati (xc, XWINDOW) + if (nrlines == 1) + ywindow = 1 + else + ywindow = rg_xstati (xc, YWINDOW) + + # Move to next ref regions if current region is too small. + if (nrcols < xwindow || (IM_NDIM(imr) == 2 && nrlines < ywindow)) { + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE) + call eprintf ( + "Reference section: %s[%d:%d,%d:%d] has too few points.\n") + call pargstr (Memc[str]) + call pargi (rc1) + call pargi (rc2) + call pargi (rl1) + call pargi (rl2) + stat = ERR + goto nextregion_ + } + + # Apply the transformation if defined or lag to the ref regions. + if (rg_xstati (xc, NREFPTS) > 0) { + call rg_etransform (xc, (rc1 + rc2) / 2.0, (rl1 + rl2) / 2.0, + rxlag, rylag) + xlag = rxlag - (rc1 + rc2) / 2.0 + if (ywindow == 1) + ylag = 0 + else + ylag = rylag - (rl1 + rl2) / 2.0 + } else { + xlag = rg_xstati (xc, XLAG) + if (ywindow == 1) + ylag = 0 + else + ylag = rg_xstati (xc, YLAG) + } + + # Get the input image limits. + c1 = rc1 + xlag - xwindow / 2 + c2 = rc2 + xlag + xwindow / 2 + l1 = rl1 + ylag - ywindow / 2 + l2 = rl2 + ylag + ywindow / 2 + ncols = c2 - c1 + 1 + nlines = l2 - l1 + 1 + + # Move to the next ref region if input region is off image. + if (c1 > nimcols || c2 < 1 || l1 > nimlines || l2 < 1) { + call rg_xstats (xc, IMAGE, Memc[str], SZ_LINE) + call eprintf ( + "Image section: %s[%d:%d,%d:%d] is off image.\n") + call pargstr (Memc[str]) + call pargi (c1) + call pargi (c2) + call pargi (l1) + call pargi (l2) + stat = ERR + goto nextregion_ + } + + # Move to the next ref region if input region is less than 3 by 3. + if ((ncols < xwindow) || (IM_NDIM(im1) == 2 && nlines < ywindow)) { + call rg_xstats (xc, IMAGE, Memc[str], SZ_LINE) + call eprintf ( + "Image section: %s[%d:%d,%d:%d] has too few points.\n") + call pargstr (Memc[str]) + call pargi (c1) + call pargi (c2) + call pargi (l1) + call pargi (l2) + stat = ERR + goto nextregion_ + } + + # Get the input reference and input image data. + rbuf = rg_ximget (imr, rc1, rc2, rl1, rl2) + if (rbuf == NULL) { + stat = ERR + goto nextregion_ + } + ibuf = rg_ximget (im1, c1, c2, l1, l2) + if (ibuf == NULL) { + stat = ERR + goto nextregion_ + } + + # Do the background subtraction. + + # Compute the zero point, x slope and y slope of ref image. + if (IS_INDEFR(Memr[przero+i-1]) || IS_INDEFR(Memr[prxslope+i- 1]) || + IS_INDEFR(Memr[pryslope+i-1])) { + if (IS_INDEFI(rg_xstati (xc, BORDER))) { + call rg_xscale (xc, Memr[rbuf], nrcols * nrlines, nrcols, + nrlines, rg_xstatr (xc, BVALUER), Memr[coeff]) + } else { + border = NULL + nborder = rg_border (Memr[rbuf], nrcols, nrlines, + max (0, nrcols - 2 * rg_xstati (xc, BORDER)), + max (0, nrlines - 2 * rg_xstati (xc, BORDER)), + border) + call rg_xscale (xc, Memr[border], nborder, nrcols, + nrlines, rg_xstatr (xc, 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] + } + + call rg_subtract (Memr[rbuf], nrcols, nrlines, Memr[przero+i-1], + Memr[prxslope+i-1], Memr[pryslope+i-1]) + + # Compute the zero point, and the x and y slopes of input image. + if (IS_INDEFI(rg_xstati (xc, BORDER))) { + call rg_xscale (xc, Memr[ibuf], ncols * nlines, ncols, + nlines, rg_xstatr (xc, BVALUE), Memr[coeff]) + } else { + border = NULL + nborder = rg_border (Memr[ibuf], ncols, nlines, + max (0, ncols - 2 * rg_xstati (xc, BORDER)), + max (0, nlines - 2 * rg_xstati (xc, BORDER)), + border) + call rg_xscale (xc, Memr[border], nborder, ncols, nlines, + rg_xstatr (xc, BVALUE), Memr[coeff]) + if (border != NULL) + call mfree (border, TY_REAL) + } + + # Subtract the baseline. + call rg_subtract (Memr[ibuf], ncols, nlines, Memr[coeff], + Memr[coeff+1], Memr[coeff+2]) + + # Apodize the data. + if (rg_xstatr (xc, APODIZE) > 0.0) { + call rg_apodize (Memr[rbuf], nrcols, nrlines, rg_xstatr (xc, + APODIZE), YES) + call rg_apodize (Memr[ibuf], ncols, nlines, rg_xstatr (xc, + APODIZE), YES) + } + + # Spatially filter the data with a Laplacian. + switch (rg_xstati (xc, FILTER)) { + case XC_LAPLACE: + call rg_xlaplace (Memr[rbuf], nrcols, nrlines, 1.0) + call rg_xlaplace (Memr[ibuf], ncols, nlines, 1.0) + default: + ; + } + + # Allocate space for the cross-correlation function. + if (rg_xstatp (xc, XCOR) == NULL) { + call malloc (xcor, xwindow * ywindow, TY_REAL) + call rg_xsetp (xc, XCOR, xcor) + } else { + xcor = rg_xstatp (xc, XCOR) + call realloc (xcor, xwindow * ywindow, TY_REAL) + call rg_xsetp (xc, XCOR, xcor) + } + + # Clear the correlation function. + call aclrr (Memr[xcor], xwindow * ywindow) + + # Compute the cross-correlation function. + if (rg_xstati (xc, CFUNC) == XC_DISCRETE) { + call rg_xconv (Memr[rbuf], nrcols, nrlines, Memr[ibuf], ncols, + nlines, Memr[xcor], xwindow, ywindow) + } else { + call rg_xdiff (Memr[rbuf], nrcols, nrlines, Memr[ibuf], ncols, + nlines, Memr[xcor], xwindow, ywindow) + } + + stat = OK + +nextregion_ + + # Free memory. + call sfree (sp) + if (rbuf != NULL) + call mfree (rbuf, TY_REAL) + if (ibuf != NULL) + call mfree (ibuf, TY_REAL) + if (stat == ERR) + return (ERR) + else + return (OK) +end + + +# RG_XFGET -- Compute the cross-correlation function using Fourier techniques. + +int procedure rg_xfget (xc, imr, im1, i) + +pointer xc #I pointer to the cross-correlation structure +pointer imr #I pointer to the reference image +pointer im1 #I pointer to the input image +int i #I index of the current region + +int rc1, rc2, rl1, rl2, nrcols, nrlines, c1, c2, l1, l2, ncols, nlines +int nrimcols, nrimlines, nimcols, nimlines +int xwindow, ywindow, xlag, nxfft, nyfft, ylag, stat, nborder +pointer sp, str, coeff, xcor, rbuf, ibuf, fft, border +pointer prc1, prc2, prl1, prl2, przero, prxslope, pryslope +real rxlag, rylag +int rg_xstati(), rg_border(), rg_szfft() +pointer rg_xstatp(), rg_ximget() +real rg_xstatr() + +define nextregion_ 11 + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (coeff, max (GS_SAVECOEFF+6, 9), TY_REAL) + + # Check for number of regions. + if (i > rg_xstati (xc, NREGIONS)) { + stat = ERR + goto nextregion_ + } + + # Allocate space for the cross-correlation function. + nrimcols = IM_LEN(imr,1) + if (IM_NDIM(imr) == 1) + nrimlines = 1 + else + nrimlines = IM_LEN(imr,2) + nimcols = IM_LEN(im1,1) + if (IM_NDIM(im1) == 1) + nimlines = 1 + else + nimlines = IM_LEN(im1,2) + + # Get the regions pointers. + prc1 = rg_xstatp (xc, RC1) + prc2 = rg_xstatp (xc, RC2) + prl1 = rg_xstatp (xc, RL1) + prl2 = rg_xstatp (xc, RL2) + przero = rg_xstatp (xc, RZERO) + prxslope = rg_xstatp (xc, RXSLOPE) + pryslope = rg_xstatp (xc, RYSLOPE) + + # Get the reference subraster region. + rc1 = max (1, min (int (nrimcols), Memi[prc1+i-1])) + rc2 = min (int (nrimcols), max (1, Memi[prc2+i-1])) + rl1 = max (1, min (int (nrimlines), Memi[prl1+i-1])) + rl2 = min (int (nrimlines), max (1, Memi[prl2+i-1])) + nrcols = rc2 - rc1 + 1 + nrlines = rl2 - rl1 + 1 + + # Go to next region if the reference region is off the image. + if (rc1 > nrimcols || rc2 < 1 || rl1 > nrimlines || rl2 < 1) { + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE) + call eprintf ( + "Reference section: %s[%d:%d,%d:%d] is off image.\n") + call pargstr (Memc[str]) + call pargi (rc1) + call pargi (rc2) + call pargi (rl1) + call pargi (rl2) + stat = ERR + goto nextregion_ + } + + # Check the window sizes. + xwindow = rg_xstati (xc, XWINDOW) + if (nrlines == 1) + ywindow = 1 + else + ywindow = rg_xstati (xc, YWINDOW) + + # Go to the next region if the reference region has too few points. + if ((nrcols < xwindow) || (IM_NDIM(im1) == 2 && nrlines < ywindow)) { + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE) + call eprintf ( + "Reference section: %s[%d:%d,%d:%d] has too few points.\n") + call pargstr (Memc[str]) + call pargi (rc1) + call pargi (rc2) + call pargi (rl1) + call pargi (rl2) + stat = ERR + goto nextregion_ + } + + # Apply the transformation if defined or the lag. + if (rg_xstati (xc, NREFPTS) > 0) { + call rg_etransform (xc, (rc1 + rc2) / 2.0, (rl1 + rl2) / 2.0, + rxlag, rylag) + xlag = rxlag - (rc1 + rc2) / 2.0 + if (ywindow == 1) + ylag = 0 + else + ylag = rylag - (rl1 + rl2) / 2.0 + } else { + xlag = rg_xstati (xc, XLAG) + if (ywindow == 1) + ylag = 0 + else + ylag = rg_xstati (xc, YLAG) + } + + # Get the input image subraster regions. + c1 = rc1 + xlag + c2 = rc2 + xlag + l1 = rl1 + ylag + l2 = rl2 + ylag + ncols = c2 - c1 + 1 + nlines = l2 - l1 + 1 + + # Go to next region if region is off the image. + if (c1 > nimcols || c2 < 1 || l1 > nimlines || l2 < 1) { + call rg_xstats (xc, IMAGE, Memc[str], SZ_LINE) + call eprintf ( + "Image section: %s[%d:%d,%d:%d] is off image.\n") + call pargstr (Memc[str]) + call pargi (c1) + call pargi (c2) + call pargi (l1) + call pargi (l2) + stat = ERR + goto nextregion_ + } + + # Go to next region if region has too few points. + if ((ncols < xwindow) || (IM_NDIM(im1) == 2 && nlines < ywindow)) { + call rg_xstats (xc, IMAGE, Memc[str], SZ_LINE) + call eprintf ( + "Image section: %s[%d:%d,%d:%d] has too few points.\n") + call pargstr (Memc[str]) + call pargi (c1) + call pargi (c2) + call pargi (l1) + call pargi (l2) + stat = ERR + goto nextregion_ + } + + # 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, xwindow) + if (ywindow == 1) + nyfft = 1 + else + nyfft = rg_szfft (nrlines, ywindow) + call calloc (fft, 2 * nxfft * nyfft, TY_REAL) + + # Get the input reference and input image data. + rbuf = NULL + rbuf = rg_ximget (imr, rc1, rc2, rl1, rl2) + if (rbuf == NULL) { + stat = ERR + goto nextregion_ + } + + # Do the background subtraction. + + # Compute the zero point, x slope and y slope of ref image. + if (IS_INDEFR(Memr[przero+i-1]) || IS_INDEFR(Memr[prxslope+i- 1]) || + IS_INDEFR(Memr[pryslope+i-1])) { + if (IS_INDEFI(rg_xstati (xc, BORDER))) { + call rg_xscale (xc, Memr[rbuf], nrcols * nrlines, nrcols, + nrlines, rg_xstatr (xc, BVALUER), Memr[coeff]) + } else { + border = NULL + nborder = rg_border (Memr[rbuf], nrcols, nrlines, + max (0, nrcols - 2 * rg_xstati (xc, BORDER)), + max (0, nrlines - 2 * rg_xstati (xc, BORDER)), + border) + call rg_xscale (xc, Memr[border], nborder, nrcols, + nrlines, rg_xstatr (xc, 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] + } + + call rg_subtract (Memr[rbuf], nrcols, nrlines, Memr[przero+i-1], + Memr[prxslope+i-1], Memr[pryslope+i-1]) + + # Apodize the data. + if (rg_xstatr (xc, APODIZE) > 0.0) + call rg_apodize (Memr[rbuf], nrcols, nrlines, rg_xstatr (xc, + APODIZE), YES) + + # Spatially filter the data with a Laplacian. + switch (rg_xstati (xc, FILTER)) { + case XC_LAPLACE: + call rg_xlaplace (Memr[rbuf], nrcols, nrlines, 1.0) + default: + ; + } + + # Load the reference data into the FFT. + call rg_rload (Memr[rbuf], nrcols, nrlines, Memr[fft], nxfft, nyfft) + call mfree (rbuf, TY_REAL) + + ibuf = NULL + ibuf = rg_ximget (im1, c1, c2, l1, l2) + if (ibuf == NULL) { + stat = ERR + goto nextregion_ + } + + # Compute the zero point, and the x and y slopes of input image. + if (IS_INDEFI(rg_xstati (xc, BORDER))) { + call rg_xscale (xc, Memr[ibuf], ncols * nlines, ncols, + nlines, rg_xstatr (xc, BVALUE), Memr[coeff]) + } else { + border = NULL + nborder = rg_border (Memr[ibuf], ncols, nlines, + max (0, ncols - 2 * rg_xstati (xc, BORDER)), + max (0, nlines - 2 * rg_xstati (xc, BORDER)), + border) + call rg_xscale (xc, Memr[border], nborder, ncols, nlines, + rg_xstatr (xc, BVALUE), Memr[coeff]) + if (border != NULL) + call mfree (border, TY_REAL) + } + + # Subtract the baseline. + call rg_subtract (Memr[ibuf], ncols, nlines, Memr[coeff], + Memr[coeff+1], Memr[coeff+2]) + + # Apodize the data. + if (rg_xstatr (xc, APODIZE) > 0.0) + call rg_apodize (Memr[ibuf], ncols, nlines, rg_xstatr (xc, + APODIZE), YES) + + # Spatially filter the data with a Laplacian. + switch (rg_xstati (xc, FILTER)) { + case XC_LAPLACE: + call rg_xlaplace (Memr[ibuf], ncols, nlines, 1.0) + default: + ; + } + + # Load the image data into the FFT. + call rg_iload (Memr[ibuf], ncols, nlines, Memr[fft], nxfft, nyfft) + call mfree (ibuf, TY_REAL) + + # Normalize the data. + call rg_fnorm (Memr[fft], nrcols, nrlines, nxfft, nyfft) + + # Compute the cross-correlation function. + call rg_fftcor (Memr[fft], nxfft, nyfft) + + # Allocate space for the correlation function. + if (rg_xstatp (xc, XCOR) == NULL) { + call malloc (xcor, xwindow * ywindow, TY_REAL) + call rg_xsetp (xc, XCOR, xcor) + } else { + xcor = rg_xstatp (xc, XCOR) + call realloc (xcor, xwindow * ywindow, TY_REAL) + call rg_xsetp (xc, XCOR, xcor) + } + + # Move the valid lags into the crosscorrelation array + call rg_movexr (Memr[fft], nxfft, nyfft, Memr[xcor], xwindow, ywindow) + + # Free space. + call mfree (fft, TY_REAL) + + stat = OK + +nextregion_ + + call sfree (sp) + if (stat == ERR) + return (ERR) + else + return (OK) +end + + +# RG_XIMGET -- Fill a buffer from a specified region of the image. + +pointer procedure rg_ximget (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_XLAPLACE -- Compute the Laplacian of an image subraster in place. + +procedure rg_xlaplace (data, nx, ny, rho) + +real data[nx,ARB] #I the input array +int nx, ny #I the size of the input/output data array +real rho #I the pixel to pixel correlation factor + +int i, inline, outline, nxk, nyk, nxc +pointer sp, lineptrs, ptr +real rhosq, kernel[3,3] +data nxk /3/, nyk /3/ + +begin + # Define the kernel. + rhosq = rho * rho + kernel[1,1] = rhosq + kernel[2,1] = -rho * (1.0 + rhosq) + kernel[3,1] = rhosq + kernel[1,2] = -rho * (1.0 + rhosq) + kernel[2,2] = (1.0 + rhosq) * (1 + rhosq) + kernel[3,2] = -rho * (1.0 + rhosq) + kernel[1,3] = rhosq + kernel[2,3] = -rho * (1.0 + rhosq) + kernel[3,3] = rhosq + + # Set up an array of line pointers. + call smark (sp) + call salloc (lineptrs, nyk, TY_POINTER) + + # Allocate working space. + nxc = nx + 2 * (nxk / 2) + do i = 1, nyk + call salloc (Memi[lineptrs+i-1], nxc, TY_REAL) + + inline = 1 - nyk / 2 + do i = 1, nyk - 1 { + if (inline < 1) { + call amovr (data[1,1], Memr[Memi[lineptrs+i]+nxk/2], nx) + Memr[Memi[lineptrs+i]] = data[1,1] + Memr[Memi[lineptrs+i]+nxc-1] = data[nx,1] + } else { + call amovr (data[1,i-1], Memr[Memi[lineptrs+i]+nxk/2], nx) + Memr[Memi[lineptrs+i]] = data[1,i-1] + Memr[Memi[lineptrs+i]+nxc-1] = data[nx,i-1] + } + inline = inline + 1 + } + + # Generate the output image line by line + do outline = 1, ny { + + # Scroll the input buffers + ptr = Memi[lineptrs] + do i = 1, nyk - 1 + Memi[lineptrs+i-1] = Memi[lineptrs+i] + Memi[lineptrs+nyk-1] = ptr + + # Read in new image line + if (inline > ny) { + call amovr (data[1,ny], Memr[Memi[lineptrs+nyk-1]+nxk/2], + nx) + Memr[Memi[lineptrs+nyk-1]] = data[1,ny] + Memr[Memi[lineptrs+nyk-1]+nxc-1] = data[nx,ny] + } else { + call amovr (data[1,inline], Memr[Memi[lineptrs+nyk-1]+nxk/2], + nx) + Memr[Memi[lineptrs+nyk-1]] = data[1,inline] + Memr[Memi[lineptrs+nyk-1]+nxc-1] = data[nx,inline] + } + + # Generate output image line + call aclrr (data[1,outline], nx) + do i = 1, nyk + call acnvr (Memr[Memi[lineptrs+i-1]], data[1,outline], nx, + kernel[1,i], nxk) + + inline = inline + 1 + } + + # Free the image buffer pointers + call sfree (sp) +end + + +# RG_XCONV -- Compute the cross-correlation function directly in the spatial +# domain. + +procedure rg_xconv (ref, nrcols, nrlines, image, ncols, nlines, xcor, xwindow, + ywindow) + +real ref[nrcols,nrlines] #I the input reference subraster +int nrcols, nrlines #I size of the reference subraster +real image[ncols,nlines] #I the input image subraster +int ncols, nlines #I size of the image subraster +real xcor[xwindow,ywindow] #O the output cross-correlation function +int xwindow, ywindow #I size of the cross-correlation function + +int lagx, lagy, i, j +real meanr, facr, meani, faci, sum +real asumr() +#real cxmin, cxmax + +begin + meanr = asumr (ref, nrcols * nrlines) / (nrcols * nrlines) + facr = 0.0 + do j = 1, nrlines { + do i = 1, nrcols + facr = facr + (ref[i,j] - meanr) ** 2 + } + if (facr <= 0.0) + facr = 1.0 + else + facr = sqrt (facr) + + do lagy = 1, ywindow { + do lagx = 1, xwindow { + meani = 0.0 + do j = 1, nrlines { + do i = 1, nrcols + meani = meani + image[i+lagx-1,j+lagy-1] + } + meani = meani / (nrcols * nrlines) + faci = 0.0 + sum = 0.0 + do j = 1, nrlines { + do i = 1, nrcols { + faci = faci + (image[i+lagx-1,j+lagy-1] - meani) ** 2 + sum = sum + (ref[i,j] - meanr) * + (image[i+lagx-1,j+lagy-1] - meani) + } + } + if (faci <= 0.0) + faci = 1.0 + else + faci = sqrt (faci) + xcor[lagx,lagy] = sum / facr / faci + } + } +end + + +# RG_XDIFF -- Compute the error function at each of several templates. + +procedure rg_xdiff (ref, nrcols, nrlines, image, ncols, nlines, xcor, xwindow, + ywindow) + +real ref[nrcols,nrlines] #I reference subraste +int nrcols, nrlines #I size of the reference subraster +real image[ncols,nlines] #I image subraster +int ncols, nlines #I size of image subraster +real xcor[xwindow,ywindow] #O crosscorrelation function +int xwindow, ywindow #I size of correlation function + +int lagx, lagy, i, j +real meanr, meani, sum, cormin, cormax +real asumr() + + +begin + meanr = asumr (ref, nrcols * nrlines) / (nrcols * nrlines) + do lagy = 1, ywindow { + do lagx = 1, xwindow { + meani = 0.0 + do j = 1, nrlines { + do i = 1, nrcols + meani = meani + image[i+lagx-1,j+lagy-1] + } + meani = meani / (nrcols * nrlines) + sum = 0.0 + do j = 1, nrlines { + do i = 1, nrcols { + sum = sum + abs ((ref[i,j] - meanr) - + (image[i+lagx-1,j+lagy-1] - meani)) + } + } + xcor[lagx,lagy] = sum + } + } + + call alimr (xcor, xwindow * ywindow, cormin, cormax) + call adivkr (xcor, cormax, xcor, xwindow * ywindow) + call asubkr (xcor, 1.0, xcor, xwindow * ywindow) + call anegr (xcor, xcor, xwindow * ywindow) +end + diff --git a/pkg/images/immatch/src/xregister/rgxdbio.x b/pkg/images/immatch/src/xregister/rgxdbio.x new file mode 100644 index 00000000..3e197636 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxdbio.x @@ -0,0 +1,290 @@ +include "xregister.h" + +# RG_XWREC -- Procedure to write out the whole record. + +procedure rg_xwrec (db, dformat, xc) + +pointer db #I pointer to the database file +int dformat #I is the shifts file in database format +pointer xc #I pointer to the cross correlation structure + +int i, nregions, ngood, c1, c2, l1, l2, xlag, ylag +pointer sp, image, prc1, prc2, prl1, prl2, pxshift, pyshift +real xin, yin, xout, yout, xavshift, yavshift +int rg_xstati() +pointer rg_xstatp() + +begin + # Allocate working space. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + + # Write the header record. + if (dformat == YES) + call rg_xdbparams (db, xc) + + # Fetch the pointers to the columns. + prc1 = rg_xstatp (xc, RC1) + prc2 = rg_xstatp (xc, RC2) + prl1 = rg_xstatp (xc, RL1) + prl2 = rg_xstatp (xc, RL2) + pxshift = rg_xstatp (xc, XSHIFTS) + pyshift = rg_xstatp (xc, YSHIFTS) + nregions = rg_xstati (xc, NREGIONS) + + xavshift = 0.0 + yavshift = 0.0 + ngood = 0 + do i = 1, nregions { + + xin = (Memi[prc1+i-1] + Memi[prc2+i-1]) / 2.0 + yin = (Memi[prl1+i-1] + Memi[prl2+i-1]) / 2.0 + if (rg_xstati (xc, NREFPTS) > 0) { + call rg_etransform (xc, xin, yin, xout, yout) + xlag = xout - xin + ylag = yout - yin + } else { + xlag = rg_xstati (xc, XLAG) + ylag = rg_xstati (xc, YLAG) + } + c1 = Memi[prc1+i-1] + xlag + c2 = Memi[prc2+i-1] + xlag + l1 = Memi[prl1+i-1] + ylag + l2 = Memi[prl2+i-1] + ylag + + if (IS_INDEFR(Memr[pxshift+i-1]) || IS_INDEFR(Memr[pyshift+i-1])) { + if (dformat == YES) + call rg_xdbshiftr (db, Memi[prc1+i-1], Memi[prc2+i-1], + Memi[prl1+i-1], Memi[prl2+i-1], c1, c2, l1, l2, + INDEFR, INDEFR) + } else { + if (dformat == YES) + call rg_xdbshiftr (db, Memi[prc1+i-1], Memi[prc2+i-1], + Memi[prl1+i-1], Memi[prl2+i-1], c1, c2, l1, l2, + Memr[pxshift+i-1], Memr[pyshift+i-1]) + ngood = ngood + 1 + xavshift = xavshift + Memr[pxshift+i-1] + yavshift = yavshift + Memr[pyshift+i-1] + } + } + + # Compute the average shift. + if (ngood <= 0) { + xavshift = 0.0 + yavshift = 0.0 + } else { + xavshift = xavshift / ngood + yavshift = yavshift / ngood + } + call rg_xsetr (xc, TXSHIFT, xavshift) + call rg_xsetr (xc, TYSHIFT, yavshift) + + if (dformat == YES) + call rg_xdbshift (db, xc) + else { + call rg_xstats (xc, IMAGE, Memc[image], SZ_FNAME) + call fprintf (db, "%s %g %g\n") + call pargstr (Memc[image]) + call pargr (xavshift) + call pargr (yavshift) + } + + call sfree (sp) +end + + +# RG_XDBPARAMS -- Write the cross-correlation parameters to the database file. + +procedure rg_xdbparams (db, xc) + +pointer db #I pointer to the database file +pointer xc #I pointer to the cross-correlation structure + +pointer sp, str +int rg_xstati() +#real rg_xstatr() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Write out the time record was written. + call dtput (db, "\n") + call dtptime (db) + + # Write out the record name. + call rg_xstats (xc, RECORD, Memc[str], SZ_FNAME) + call dtput (db, "begin\t%s\n") + call pargstr (Memc[str]) + + # Write the image names. + call rg_xstats (xc, IMAGE, Memc[str], SZ_FNAME) + call dtput (db, "\t%s\t\t%s\n") + call pargstr (KY_IMAGE) + call pargstr (Memc[str]) + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_FNAME) + call dtput (db, "\t%s\t%s\n") + call pargstr (KY_REFIMAGE) + call pargstr (Memc[str]) + + call dtput (db, "\t%s\t%d\n") + call pargstr (KY_NREGIONS) + call pargi (rg_xstati (xc, NREGIONS)) + + call sfree (sp) +end + + +# RG_XWREG -- Write out the results for each region individually into +# the shifts file. + +procedure rg_xwreg (db, xc) + +pointer db #I pointer to the database file +pointer xc #I pointer to the cross-correlation structure + +int i, nregions, c1, c2, l1, l2, xlag, ylag +pointer prc1, prc2, prl1, prl2, pxshift, pyshift +real xin, yin, xout, yout +int rg_xstati() +pointer rg_xstatp() + +begin + # Fetch the regions pointers. + prc1 = rg_xstatp (xc, RC1) + prc2 = rg_xstatp (xc, RC2) + prl1 = rg_xstatp (xc, RL1) + prl2 = rg_xstatp (xc, RL2) + pxshift = rg_xstatp (xc, XSHIFTS) + pyshift = rg_xstatp (xc, YSHIFTS) + nregions = rg_xstati (xc, NREGIONS) + + # Write out the reference image region(s) and the equivalent + # input image regions. + do i = 1, nregions { + + xin = (Memi[prc1+i-1] + Memi[prc2+i-1]) / 2.0 + yin = (Memi[prl1+i-1] + Memi[prl2+i-1]) / 2.0 + if (rg_xstati (xc, NREFPTS) > 0) { + call rg_etransform (xc, xin, yin, xout, yout) + xlag = xout - xin + ylag = yout - yin + } else { + xlag = rg_xstati (xc, XLAG) + ylag = rg_xstati (xc, YLAG) + } + c1 = Memi[prc1+i-1] + xlag + c2 = Memi[prc2+i-1] + xlag + l1 = Memi[prl1+i-1] + ylag + l2 = Memi[prl2+i-1] + ylag + + if (IS_INDEFR(Memr[pxshift+i-1]) || IS_INDEFR(Memr[pyshift+i-1])) + call rg_xdbshiftr (db, Memi[prc1+i-1], Memi[prc2+i-1], + Memi[prl1+i-1], Memi[prl2+i-1], c1, c2, l1, l2, + INDEFR, INDEFR) + else + call rg_xdbshiftr (db, Memi[prc1+i-1], Memi[prc2+i-1], + Memi[prl1+i-1], Memi[prl2+i-1], c1, c2, l1, l2, + Memr[pxshift+i-1], Memr[pyshift+i-1]) + } +end + + +# RG_XDBSHIFTR -- Write out the reference image section, input image +# section and x and y shifts for each region. + +procedure rg_xdbshiftr (db, rc1, rc2, rl1, rl2, c1, c2, l1, l2, xshift, yshift) + +pointer db #I pointer to the database file +int rc1, rc2 #I reference region column limits +int rl1, rl2 #I reference region line limits +int c1, c2 #I image region column limits +int l1, l2 #I image region line limits +real xshift #I x shift +real yshift #I y shift + +begin + call dtput (db,"\t[%d:%d,%d:%d]\t[%d:%d,%d:%d]\t%g\t%g\n") + call pargi (rc1) + call pargi (rc2) + call pargi (rl1) + call pargi (rl2) + call pargi (c1) + call pargi (c2) + call pargi (l1) + call pargi (l2) + call pargr (xshift) + call pargr (yshift) +end + + +# RG_XDBSHIFT -- Write the average shifts to the shifts database. + +procedure rg_xdbshift (db, xc) + +pointer db #I pointer to text database file +pointer xc #I pointer to the cross-correlation structure + +real rg_xstatr() + +begin + call dtput (db, "\t%s\t\t%g\n") + call pargstr (KY_TXSHIFT) + call pargr (rg_xstatr (xc, TXSHIFT)) + call dtput (db, "\t%s\t\t%g\n") + call pargstr (KY_TYSHIFT) + call pargr (rg_xstatr (xc, TYSHIFT)) +end + + +# RG_XPWREC -- Print the computed shift for a region. + +procedure rg_xpwrec (xc, i) + +pointer xc #I pointer to the cross-correlation structure +int i #I the current region + +int xlag, ylag, c1, c2, l1, l2 +pointer prc1, prc2, prl1, prl2 +real xin, yin, rxlag, rylag +int rg_xstati() +pointer rg_xstatp() + +begin + # Fetch the pointers to the reference regions. + prc1 = rg_xstatp (xc, RC1) + prc2 = rg_xstatp (xc, RC2) + prl1 = rg_xstatp (xc, RL1) + prl2 = rg_xstatp (xc, RL2) + + # Transform the reference region to the input region. + xin = (Memi[prc1+i-1] + Memi[prc2+i-1]) / 2.0 + yin = (Memi[prl1+i-1] + Memi[prl2+i-1]) / 2.0 + if (rg_xstati (xc, NREFPTS) > 0) { + call rg_etransform (xc, xin, yin, rxlag, rylag) + xlag = rxlag - xin + ylag = rylag - yin + } else { + xlag = rg_xstati (xc, XLAG) + ylag = rg_xstati (xc, YLAG) + } + + c1 = Memi[prc1+i-1] + xlag + c2 = Memi[prc2+i-1] + xlag + l1 = Memi[prl1+i-1] + ylag + l2 = Memi[prl2+i-1] + ylag + + # Print the results. + call printf ("Region %d: [%d:%d,%d:%d] [%d:%d,%d:%d] %g %g\n") + call pargi (i) + call pargi (Memi[prc1+i-1]) + call pargi (Memi[prc2+i-1]) + call pargi (Memi[prl1+i-1]) + call pargi (Memi[prl2+i-1]) + call pargi (c1) + call pargi (c2) + call pargi (l1) + call pargi (l2) + call pargr (Memr[rg_xstatp(xc,XSHIFTS)+i-1]) + call pargr (Memr[rg_xstatp(xc,YSHIFTS)+i-1]) +end diff --git a/pkg/images/immatch/src/xregister/rgxfft.x b/pkg/images/immatch/src/xregister/rgxfft.x new file mode 100644 index 00000000..8847cf56 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxfft.x @@ -0,0 +1,179 @@ +# RG_FFTCOR -- Compute the FFT of the reference and image data, take their +# product, and compute the inverse transform to get the cross-correlation +# function. The reference and input image are loaded into alternate memory +# locations. + +procedure rg_fftcor (fft, nxfft nyfft) + +real fft[ARB] #I/O array to be fft'd +int nxfft, nyfft #I dimensions of the fft + +pointer sp, dim + +begin + call smark (sp) + call salloc (dim, 2, TY_INT) + + # Fourier transform the two arrays. + Memi[dim] = nxfft + Memi[dim+1] = nyfft + if (Memi[dim+1] == 1) + call rg_fourn (fft, Memi[dim], 1, 1) + else + call rg_fourn (fft, Memi[dim], 2, 1) + + # Compute the product of the two transforms. + call rg_mulfft (fft, fft, 2 * nxfft, nyfft) + + # Shift the array to center the transform. + call rg_fshift (fft, fft, 2 * nxfft, nyfft) + + # Normalize the transform. + call adivkr (fft, real (nxfft * nyfft), fft, 2 * nxfft * nyfft) + + # Compute the inverse transform. + if (Memi[dim+1] == 1) + call rg_fourn (fft, Memi[dim], 1, -1) + else + call rg_fourn (fft, Memi[dim], 2, -1) + + call sfree (sp) +end + + +# RG_MULFFT -- Unpack the two individual ffts and compute their product. + +procedure rg_mulfft (fft1, fft2, nxfft, nyfft) + +real fft1[nxfft,nyfft] #I array containing 2 ffts of 2 real functions +real fft2[nxfft,nyfft] #O fft of correlation function +int nxfft, nyfft #I dimensions of fft + +int i,j, nxd2p2, nxp2, nxp3, nyd2p1, nyp2 +real c1, c2, h1r, h1i, h2r, h2i + +begin + c1 = 0.5 + c2 = -0.5 + + nxd2p2 = nxfft / 2 + 2 + nxp2 = nxfft + 2 + nxp3 = nxfft + 3 + nyd2p1 = nyfft / 2 + 1 + nyp2 = nyfft + 2 + + # Compute the 0 frequency point. + h1r = fft1[1,1] + h1i = 0.0 + h2r = fft1[2,1] + h2i = 0.0 + fft2[1,1] = h1r * h2r + fft2[2,1] = 0.0 + + # Compute the x axis points. + do i = 3, nxd2p2, 2 { + h2r = c1 * (fft1[i,1] + fft1[nxp2-i,1]) + h2i = c1 * (fft1[i+1,1] - fft1[nxp3-i,1]) + h1r = -c2 * (fft1[i+1,1] + fft1[nxp3-i,1]) + h1i = c2 * (fft1[i,1] - fft1[nxp2-i,1]) + fft2[i,1] = (h1r * h2r + h1i * h2i) + fft2[i+1,1] = (h1i * h2r - h2i * h1r) + fft2[nxp2-i,1] = fft2[i,1] + fft2[nxp3-i,1] = - fft2[i+1,1] + } + + # Quit if the transform is 1D. + if (nyfft < 2) + return + + # Compute the y axis points. + do i = 2, nyd2p1 { + h2r = c1 * (fft1[1,i] + fft1[1, nyp2-i]) + h2i = c1 * (fft1[2,i] - fft1[2,nyp2-i]) + h1r = -c2 * (fft1[2,i] + fft1[2,nyp2-i]) + h1i = c2 * (fft1[1,i] - fft1[1,nyp2-i]) + fft2[1,i] = (h1r * h2r + h1i * h2i) + fft2[2,i] = (h1i * h2r - h2i * h1r) + fft2[1,nyp2-i] = fft2[1,i] + fft2[2,nyp2-i] = - fft2[2,i] + } + + # Compute along the axis of symmetry. + do i = 3, nxd2p2, 2 { + h2r = c1 * (fft1[i,nyd2p1] + fft1[nxp2-i, nyd2p1]) + h2i = c1 * (fft1[i+1,nyd2p1] - fft1[nxp3-i,nyd2p1]) + h1r = -c2 * (fft1[i+1,nyd2p1] + fft1[nxp3-i,nyd2p1]) + h1i = c2 * (fft1[i,nyd2p1] - fft1[nxp2-i,nyd2p1]) + fft2[i,nyd2p1] = (h1r * h2r + h1i * h2i) + fft2[i+1,nyd2p1] = (h1i * h2r - h2i * h1r) + fft2[nxp2-i,nyd2p1] = fft2[i,nyd2p1] + fft2[nxp3-i,nyd2p1] = - fft2[i+1,nyd2p1] + } + + # Compute the remainder of the transform. + do j = 2, nyd2p1 - 1 { + do i = 3, nxfft, 2 { + h2r = c1 * (fft1[i,j] + fft1[nxp2-i, nyp2-j]) + h2i = c1 * (fft1[i+1,j] - fft1[nxp3-i,nyp2-j]) + h1r = -c2 * (fft1[i+1,j] + fft1[nxp3-i,nyp2-j]) + h1i = c2 * (fft1[i,j] - fft1[nxp2-i,nyp2-j]) + fft2[i,j] = (h1r * h2r + h1i * h2i) + fft2[i+1,j] = (h1i * h2r - h2i * h1r) + fft2[nxp2-i,nyp2-j] = fft2[i,j] + fft2[nxp3-i,nyp2-j] = - fft2[i+1,j] + } + } +end + + +# RG_FNORM -- Normalize the reference and image data before computing +# the fft's. + +procedure rg_fnorm (array, ncols, nlines, nxfft, nyfft) + +real array[ARB] #I/O the input/output data array +int ncols, nlines #I dimensions of the input data array +int nxfft, nyfft #I dimensions of the fft + +int i, j, index +real sumr, sumi, meanr, meani + +begin + # Compute the mean. + sumr = 0.0 + sumi = 0.0 + index = 0 + do j = 1, nlines { + do i = 1, ncols { + sumr = sumr + array[index+2*i-1] + sumi = sumi + array[index+2*i] + } + index = index + 2 * nxfft + } + meanr = sumr / (ncols * nlines) + meani = sumi / (ncols * nlines) + + # Compute the sigma. + sumr = 0.0 + sumi = 0.0 + index = 0 + do j = 1, nlines { + do i = 1, ncols { + sumr = sumr + (array[index+2*i-1] - meanr) ** 2 + sumi = sumi + (array[index+2*i] - meani) ** 2 + } + index = index + 2 * nxfft + } + sumr = sqrt (sumr) + sumi = sqrt (sumi) + + # Normalize the data. + index = 0 + do j = 1, nlines { + do i = 1, ncols { + array[index+2*i-1] = (array[index+2*i-1] - meanr) / sumr + array[index+2*i] = (array[index+2*i] - meani) / sumi + } + index = index + 2 * nxfft + } +end diff --git a/pkg/images/immatch/src/xregister/rgxfit.x b/pkg/images/immatch/src/xregister/rgxfit.x new file mode 100644 index 00000000..34e6398c --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxfit.x @@ -0,0 +1,814 @@ +include +include +include +include "xregister.h" + +define NL_MAXITER 10 +define NL_TOL 0.001 + +# RG_FIT -- Fit the peak of the cross-correlation function using one of the +# fitting functions. + +procedure rg_fit (xc, nreg, gd, xshift, yshift) + +pointer xc #I the pointer to the cross-corrrelation structure +int nreg #I the current region +pointer gd #I the pointer to the graphics stream +real xshift, yshift #O the computed shifts + +int nrlines, xwindow, ywindow, xcbox, ycbox, xlag, ylag +real xin, yin, xout, yout +int rg_xstati() +pointer rg_xstatp() + +begin + # Check the window and centering box sizes. + nrlines = Memi[rg_xstatp(xc,RL2)+nreg-1] - + Memi[rg_xstatp(xc,RL1)+nreg-1] + 1 + xwindow = rg_xstati (xc, XWINDOW) + if (nrlines == 1) + ywindow = 1 + else + ywindow = rg_xstati (xc, YWINDOW) + xcbox = rg_xstati (xc, XCBOX) + if (nrlines == 1) + ycbox = 1 + else + ycbox = rg_xstati (xc, YCBOX) + + # Do the centering. + switch (rg_xstati (xc, PFUNC)) { + case XC_PNONE: + call rg_maxmin (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow, + xshift, yshift) + case XC_CENTROID: + call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow, + ywindow, xcbox, ycbox, xshift, yshift) + case XC_SAWTOOTH: + call rg_sawtooth (Memr[rg_xstatp(xc,XCOR)], xwindow, + ywindow, xcbox, ycbox, xshift, yshift) + case XC_PARABOLA: + call rg_iparabolic (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow, + xcbox, ycbox, xshift, yshift) + case XC_MARK: + if (gd == NULL) + call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow, + ywindow, xcbox, ycbox, xshift, yshift) + else + call rg_xmkpeak (gd, xwindow, ywindow, xshift, yshift) + default: + call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow, + xcbox, ycbox, xshift, yshift) + } + + # Store the shifts. + if (rg_xstati (xc, NREFPTS) > 0) { + xin = (Memi[rg_xstatp(xc,RC1)+nreg-1] + + Memi[rg_xstatp(xc,RC2)+nreg-1]) / 2.0 + yin = (Memi[rg_xstatp(xc,RL1)+nreg-1] + + Memi[rg_xstatp(xc,RL2)+nreg-1]) / 2.0 + call rg_etransform (xc, xin, yin, xout, yout) + xlag = xout - xin + ylag = yout - yin + } else { + xlag = rg_xstati (xc, XLAG) + ylag = rg_xstati (xc, YLAG) + } + xshift = - (xshift + xlag) + yshift = - (yshift + ylag) + Memr[rg_xstatp(xc,XSHIFTS)+nreg-1] = xshift + Memr[rg_xstatp(xc,YSHIFTS)+nreg-1] = yshift +end + + +# RG_MAXMIN -- Procedure to compute the peak of the cross-correlation function +# by determining the maximum point. + +procedure rg_maxmin (xcor, xwindow, ywindow, xshift, yshift) + +real xcor[xwindow,ywindow] #I the cross-correlation function +int xwindow, ywindow #I dimensions of cross-correlation function +real xshift, yshift #O x and shift of the peak + +int xindex, yindex + +begin + # Locate the maximum point. + call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex) + xshift = xindex - (1.0 + xwindow) / 2.0 + yshift = yindex - (1.0 + ywindow) / 2.0 +end + + +# RG_IMEAN -- Compute the peak of the cross-correlation function using the +# intensity weighted mean of the marginal distributions in x and y. + +procedure rg_imean (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift) + +real xcor[xwindow,ARB] #I the cross-correlation function +int xwindow, ywindow #I dimensions of the cross-correlation function +int xcbox, ycbox #I dimensions of the centering box +real xshift, yshift #O x and y shift of cross-correlation function + +int xindex, yindex, xlo, xhi, ylo, yhi, nx, ny +pointer sp, xmarg, ymarg + +begin + call smark (sp) + call salloc (xmarg, xcbox, TY_REAL) + call salloc (ymarg, ycbox, TY_REAL) + + # Locate the maximum point and normalize. + call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex) + + # Compute the limits of the centering box. + xlo = max (1, xindex - xcbox / 2) + xhi = min (xwindow, xindex + xcbox / 2) + nx = xhi - xlo + 1 + ylo = max (1, yindex - ycbox / 2) + yhi = min (ywindow, yindex + ycbox / 2) + ny = yhi - ylo + 1 + + # Accumulate the marginals. + call rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi, + Memr[xmarg], Memr[ymarg]) + + # Compute the shifts. + call rg_centroid (Memr[xmarg], nx, xshift) + xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0 + call rg_centroid (Memr[ymarg], ny, yshift) + yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0 + + call sfree (sp) +end + + +# RG_IPARABOLIC -- Computer the peak of the cross-correlation function by +# doing parabolic interpolation around the peak. + +procedure rg_iparabolic (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift) + +real xcor[xwindow,ARB] #I the cross-correlation function +int xwindow, ywindow #I dimensions of the cross-correlation fucntion +int xcbox, ycbox #I the dimensions of the centering box +real xshift, yshift #O the x and y shift of the peak + +int i, j, xindex, yindex, xlo, xhi, nx, ylo, yhi, ny +pointer sp, x, y, c, xfit, yfit + +begin + # Allocate working space. + call smark (sp) + call salloc (x, 3, TY_REAL) + call salloc (y, 3, TY_REAL) + call salloc (c, 3, TY_REAL) + call salloc (xfit, 3, TY_REAL) + call salloc (yfit, 3, TY_REAL) + + # Locate the maximum point. + call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex) + + xlo = max (1, xindex - 1) + xhi = min (xwindow, xindex + 1) + nx = xhi - xlo + 1 + ylo = max (1, yindex - 1) + yhi = min (ywindow, yindex + 1) + ny = yhi - ylo + 1 + + # Initialize. + do i = 1, 3 + Memr[x+i-1] = i + + # Fit the x shift. + if (nx >= 3) { + do j = ylo, yhi { + do i = xlo, xhi + Memr[y+i-xlo] = xcor[i,j] + call rg_iparab (Memr[x], Memr[y], Memr[c]) + Memr[xfit+j-ylo] = - Memr[c+1] / (2.0 * Memr[c+2]) + Memr[yfit+j-ylo] = Memr[c] + Memr[c+1] * Memr[xfit+j-ylo] + + Memr[c+2] * Memr[xfit+j-ylo] ** 2 + } + if (ny >= 3) + call rg_iparab (Memr[xfit], Memr[yfit], Memr[c]) + xshift = - Memr[c+1] / (2.0 * Memr[c+2]) + } else + xshift = xindex - xlo + 1 + + # Fit the y shift. + if (ny >= 3) { + do i = xlo, xhi { + do j = ylo, yhi + Memr[y+j-ylo] = xcor[i,j] + call rg_iparab (Memr[x], Memr[y], Memr[c]) + Memr[xfit+i-xlo] = - Memr[c+1] / (2.0 * Memr[c+2]) + Memr[yfit+i-xlo] = Memr[c] + Memr[c+1] * Memr[xfit+i-xlo] + + Memr[c+2] * Memr[xfit+i-xlo] ** 2 + } + call rg_iparab (Memr[xfit], Memr[yfit], Memr[c]) + yshift = - Memr[c+1] / (2.0 * Memr[c+2]) + } else + yshift = yindex - ylo + 1 + + xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0 + yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0 + + call sfree (sp) +end + + +define NPARS_PARABOLA 3 + +# RG_PARABOLIC -- Compute the peak of the cross-correlation function by fitting +# a parabola to the peak. + +procedure rg_parabolic (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift) + +real xcor[xwindow,ARB] #I the cross-correlation function +int xwindow, ywindow #I dimensions of the cross-correlation fucntion +int xcbox, ycbox #I the dimensions of the centering box +real xshift, yshift #O the x and y shift of the peak + +extern rg_polyfit, rg_dpolyfit +int i, xindex, yindex, xlo, xhi, ylo, yhi, nx, ny, npar, ier +pointer sp, x, w, xmarg, ymarg, params, eparams, list, nl +int locpr() + +begin + call smark (sp) + call salloc (x, max (xwindow, ywindow), TY_REAL) + call salloc (w, max (xwindow, ywindow), TY_REAL) + call salloc (xmarg, max (xwindow, ywindow), TY_REAL) + call salloc (ymarg, max (xwindow, ywindow), TY_REAL) + call salloc (params, NPARS_PARABOLA, TY_REAL) + call salloc (eparams, NPARS_PARABOLA, TY_REAL) + call salloc (list, NPARS_PARABOLA, TY_INT) + + # Locate the maximum point. + call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex) + + xlo = max (1, xindex - xcbox / 2) + xhi = min (xwindow, xindex + xcbox / 2) + nx = xhi - xlo + 1 + ylo = max (1, yindex - ycbox / 2) + yhi = min (ywindow, yindex + ycbox / 2) + ny = yhi - ylo + 1 + + # Accumulate the marginals. + call rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi, + Memr[xmarg], Memr[ymarg]) + + # Compute the x shift. + if (nx >= 3) { + do i = 1, nx + Memr[x+i-1] = i + do i = 1, nx + Memr[w+i-1] = Memr[xmarg+i-1] + call rg_iparab (Memr[x+xindex-xlo-1], Memr[xmarg+xindex-xlo-1], + Memr[params]) + xshift = - Memr[params+1] / (2.0 * Memr[params+2]) + call eprintf ("\txshift=%g\n") + call pargr (xshift) + call aclrr (Memr[eparams], NPARS_PARABOLA) + do i = 1, NPARS_PARABOLA + Memi[list+i-1] = i + call nlinitr (nl, locpr (rg_polyfit), locpr (rg_dpolyfit), + Memr[params], Memr[eparams], NPARS_PARABOLA, Memi[list], + NPARS_PARABOLA, .0001, NL_MAXITER) + call nlfitr (nl, Memr[x], Memr[xmarg], Memr[w], nx, 1, WTS_USER, + ier) + call nlvectorr (nl, Memr[x], Memr[w], nx, 1) + do i = 1, nx { + call eprintf ("x=%g y=%g yfit=%g\n") + call pargr (Memr[x+i-1]) + call pargr (Memr[xmarg+i-1]) + call pargr (Memr[w+i-1]) + } + if (ier != NO_DEG_FREEDOM) { + call nlpgetr (nl, Memr[params], npar) + if (Memr[params+2] != 0) + xshift = - Memr[params+1] / (2.0 * Memr[params+2]) + else + xshift = xindex - xlo + 1 + } else + xshift = xindex - xlo + 1 + call nlfreer (nl) + } else + xshift = xindex - xlo + 1 + + # Compute the y shift. + if (ny >= 3) { + do i = 1, ny + Memr[x+i-1] = i + do i = 1, ny + Memr[w+i-1] = Memr[ymarg+i-1] + call rg_iparab (Memr[x+yindex-ylo-1], Memr[ymarg+yindex-ylo-1], + Memr[params]) + yshift = - Memr[params+1] / (2.0 * Memr[params+2]) + call eprintf ("\tyshift=%g\n") + call pargr (yshift) + call aclrr (Memr[eparams], NPARS_PARABOLA) + do i = 1, NPARS_PARABOLA + Memi[list+i-1] = i + call nlinitr (nl, locpr (rg_polyfit), locpr (rg_dpolyfit), + Memr[params], Memr[eparams], NPARS_PARABOLA, Memi[list], + NPARS_PARABOLA, 0.0001, NL_MAXITER) + call nlfitr (nl, Memr[x], Memr[ymarg], Memr[w], ny, 1, WTS_USER, + ier) + call nlvectorr (nl, Memr[x], Memr[w], ny, 1) + do i = 1, ny { + call eprintf ("x=%g y=%g yfit=%g\n") + call pargr (Memr[x+i-1]) + call pargr (Memr[ymarg+i-1]) + call pargr (Memr[w+i-1]) + } + if (ier != NO_DEG_FREEDOM) { + call nlpgetr (nl, Memr[params], npar) + if (Memr[params+2] != 0) + yshift = -Memr[params+1] / (2.0 * Memr[params+2]) + else + yshift = yindex - ylo + 1 + } else + yshift = yindex - ylo + 1 + call nlfreer (nl) + } else + yshift = yindex - ylo + 1 + + xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0 + yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0 + + call sfree (sp) +end + +define EMISSION 1 # emission features +define ABSORPTION 2 # emission features + +# RG_SAWTOOTH -- Compute the the x and y centers using a sawtooth +# convolution function. + +procedure rg_sawtooth (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift) + +real xcor[xwindow,ARB] #I the cross-correlation function +int xwindow, ywindow #I the dimensions of the cross-correlation +int xcbox, ycbox #I the dimensions of the centering box +real xshift, yshift #O the x and y shifts + +int i, j, xindex, yindex, xlo, xhi, ylo, yhi, nx, ny +pointer sp, data, xfit, yfit, yclean +real ic + +begin + call smark (sp) + call salloc (data, max (xwindow, ywindow), TY_REAL) + call salloc (xfit, max (xwindow, ywindow), TY_REAL) + call salloc (yfit, max (xwindow, ywindow), TY_REAL) + call salloc (yclean, max (xwindow, ywindow), TY_REAL) + + # Locate the maximum point and normalize. + call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex) + + xlo = max (1, xindex - xcbox) + xhi = min (xwindow, xindex + xcbox) + nx = xhi - xlo + 1 + ylo = max (1, yindex - ycbox) + yhi = min (ywindow, yindex + ycbox) + ny = yhi - ylo + 1 + + # Compute the y shift. + if (ny >= 3) { + do j = ylo, yhi { + do i = xlo, xhi + Memr[data+i-xlo] = xcor[i,j] + call rg_x1dcenter (real (xindex - xlo + 1), Memr[data], nx, + Memr[xfit+j-ylo], Memr[yfit+j-ylo], real (nx / 2.0), + EMISSION, real (nx / 2.0), 0.0) + } + call arbpix (Memr[yfit], Memr[yclean], ny, II_SPLINE3, + II_BOUNDARYEXT) + call rg_x1dcenter (real (yindex - ylo + 1), Memr[yclean], ny, + yshift, ic, real (ny / 2.0), EMISSION, real (ny / 2.0), 0.0) + if (IS_INDEFR(yshift)) + yshift = yindex - ylo + 1 + } else + yshift = yindex - ylo + 1 + yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0 + + # Compute the x shift. + if (nx >= 3) { + if (ny >= 3) { + do i = xlo, xhi { + do j = ylo, yhi + Memr[data+j-ylo] = xcor[i,j] + call rg_x1dcenter (real (yindex - ylo + 1), Memr[data], ny, + Memr[xfit+i-xlo], Memr[yfit+i-xlo], real (ny / 2.0), + EMISSION, real (ny / 2.0), 0.0) + } + call arbpix (Memr[yfit], Memr[yclean], nx, II_SPLINE3, + II_BOUNDARYEXT) + call rg_x1dcenter (real (xindex - xlo + 1), Memr[yclean], nx, + xshift, ic, real (nx / 2.0), EMISSION, real (nx / 2.0), 0.0) + } else { + call rg_x1dcenter (real (xindex - xlo + 1), xcor[xlo,1], nx, + xshift, ic, real (nx / 2.0), EMISSION, real (nx / 2.0), 0.0) + } + if (IS_INDEFR(xshift)) + xshift = xindex - xlo + 1 + } else + xshift = xindex - xlo + 1 + xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0 + + call sfree (sp) +end + + +# RG_ALIM2R -- Determine the pixel position of the data maximum. + +procedure rg_alim2r (data, nx, ny, i, j) + +real data[nx,ARB] #I the input data +int nx, ny #I the dimensions of the input array +int i, j #O the indices of the maximum pixel + +int ii, jj +real datamax + +begin + datamax = -MAX_REAL + do jj = 1, ny { + do ii = 1, nx { + if (data[ii,jj] > datamax) { + datamax = data[ii,jj] + i = ii + j = jj + } + } + } +end + + +# RG_XMKMARG -- Acumulate the marginal arrays in x and y. + +procedure rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi, xmarg, + ymarg) + +real xcor[xwindow,ARB] #I the cross-correlation function +int xwindow, ywindow #I dimensions of cross-correlation function +int xlo, xhi #I the x limits for centering +int ylo, yhi #I the y limits for centering +real xmarg[ARB] #O the output x marginal array +real ymarg[ARB] #O the output y marginal array + +int i, j, index, nx, ny + +begin + nx = xhi - xlo + 1 + ny = yhi - ylo + 1 + + # Compute the x marginal. + index = 1 - xlo + do i = xlo, xhi { + xmarg[index+i] = 0.0 + do j = ylo, yhi + xmarg[index+i] = xmarg[index+i] + xcor[i,j] + } + + # Normalize the x marginal. + call adivkr (xmarg, real (ny), xmarg, nx) + + # Compute the y marginal. + index = 1 - ylo + do j = ylo, yhi { + ymarg[index+j] = 0.0 + do i = xlo, xhi + ymarg[index+j] = ymarg[index+j] + xcor[i,j] + } + + # Normalize the ymarginal. + call adivkr (ymarg, real (nx), ymarg, ny) +end + + +# RG_CENTROID -- Compute the intensity weighted maximum of an array. + +procedure rg_centroid (a, npts, shift) + +real a[ARB] #I the input array +int npts #I the number of points +real shift #O the position of the maximum + +int i +real mean, dif, sumi, sumix +bool fp_equalr() +real asumr() + +begin + sumi = 0.0 + sumix = 0.0 + mean = asumr (a, npts) / npts + + do i = 1, npts { + dif = a[i] + dif = a[i] - mean + if (dif < 0.0) + next + sumi = sumi + dif + sumix = sumix + i * dif + } + + if (fp_equalr (sumi, 0.0)) + shift = (1.0 + npts) / 2.0 + else + shift = sumix / sumi +end + + +define MIN_WIDTH 3. # minimum centering width +define EPSILON 0.001 # accuracy of centering +define EPSILON1 0.005 # tolerance for convergence check +define ITERATIONS 100 # maximum number of iterations +define MAX_DXCHECK 3 # look back for failed convergence +define INTERPTYPE II_SPLINE3 # image interpolation type + + +# RG_X1DCENTER -- Locate the center of a one dimensional feature. +# A value of INDEF is returned in the centering fails for any reason. +# This procedure just sets up the data and adjusts for emission or +# absorption features. The actual centering is done by C1D_CENTER. + +procedure rg_x1dcenter (x, data, npts, xc, ic, width, type, radius, threshold) + +real x #I initial guess +real data[npts] #I data points +int npts #I number of data points +real xc #O computed center +real ic #O intensity at computed center +real width #I feature width +int type #I feature type +real radius #I centering radius +real threshold #I minimum range in feature + +int x1, x2, nx +real a, b, rad, wid +pointer sp, data1 + +begin + # Check starting value. + if (IS_INDEF(x) || (x < 1) || (x > npts)) { + xc = INDEF + ic = INDEF + return + } + + # Set minimum width and error radius. The minimum in the error radius + # is for defining the data window. The user error radius is used to + # check for an error in the derived center at the end of the centering. + + wid = max (width, MIN_WIDTH) + rad = max (2., radius) + + # Determine the pixel value range around the initial center, including + # the width and error radius buffer. Check for a minimum range. + + x1 = max (1., x - wid / 2 - rad - wid) + x2 = min (real (npts), x + wid / 2 + rad + wid + 1) + nx = x2 - x1 + 1 + call alimr (data[x1], nx, a, b) + if (b - a < threshold) { + xc = INDEF + ic = INDEF + return + } + + # Allocate memory for the continuum subtracted data vector. The X + # range is just large enough to include the error radius and the + # half width. + + x1 = max (1., x - wid / 2 - rad) + x2 = min (real (npts), x + wid / 2 + rad + 1) + nx = x2 - x1 + 1 + + call smark (sp) + call salloc (data1, nx, TY_REAL) + call amovr (data[x1], Memr[data1], nx) + + # Make the centering data positive, subtract the continuum, and + # apply a threshold to eliminate noise spikes. + + switch (type) { + case EMISSION: + a = 0. + call asubkr (data[x1], a + threshold, Memr[data1], nx) + call amaxkr (Memr[data1], 0., Memr[data1], nx) + case ABSORPTION: + call anegr (data[x1], Memr[data1], nx) + call asubkr (Memr[data1], threshold - b, Memr[data1], nx) + call amaxkr (Memr[data1], 0., Memr[data1], nx) + default: + call error (0, "Unknown feature type") + } + + # Determine the center. + call rg_xcenter (x - x1 + 1, Memr[data1], nx, xc, ic, wid) + + # Check user centering error radius. + if (!IS_INDEF(xc)) { + xc = xc + x1 - 1 + if (abs (x - xc) > radius) { + xc = INDEF + ic = INDEF + } + } + + # Free memory and return the center position. + call sfree (sp) +end + + +# RG_XCENTER -- One dimensional centering algorithm. + +procedure rg_xcenter (x, data, npts, xc, ic, width) + +real x #I starting guess +int npts #I number of points in data vector +real data[npts] #I data vector +real xc #O computed xc +real ic #O computed intensity at xc +real width #I centering width + +int i, j, iteration, dxcheck +real hwidth, dx, dxabs, dxlast +real a, b, sum1, sum2, intgrl1, intgrl2 +pointer asi1, asi2, sp, data1 + +real asigrl(), asieval() + +define done_ 99 + +begin + # Find the nearest local maxima as the starting point. + # This is required because the threshold limit may have set + # large regions of the data to zero and without a gradient + # the centering will fail. + + i = x + for (i=x+.5; (i1) && (data[j]<=data[j-1]); j=j-1) + ; + + if (i-x < x-j) + xc = i + else + xc = j + + # Check data range. + hwidth = width / 2 + if ((xc - hwidth < 1) || (xc + hwidth > npts)) { + xc = INDEF + ic = INDEF + return + } + + # Set interpolation functions. + call asiinit (asi1, INTERPTYPE) + call asiinit (asi2, INTERPTYPE) + call asifit (asi1, data, npts) + + # Allocate, compute, and interpolate the x*y values. + call smark (sp) + call salloc (data1, npts, TY_REAL) + do i = 1, npts + Memr[data1+i-1] = data[i] * i + call asifit (asi2, Memr[data1], npts) + call sfree (sp) + + # Iterate to find center. This loop exits when 1) the maximum + # number of iterations is reached, 2) the delta is less than + # the required accuracy (criterion for finding a center), 3) + # there is a problem in the computation, 4) successive steps + # continue to exceed the minimum delta. + + dxlast = 1. + do iteration = 1, ITERATIONS { + + # Triangle centering function. + a = xc - hwidth + b = xc - hwidth / 2 + intgrl1 = asigrl (asi1, a, b) + intgrl2 = asigrl (asi2, a, b) + sum1 = (xc - hwidth) * intgrl1 - intgrl2 + sum2 = -intgrl1 + a = b + b = xc + hwidth / 2 + intgrl1 = asigrl (asi1, a, b) + intgrl2 = asigrl (asi2, a, b) + sum1 = sum1 - xc * intgrl1 + intgrl2 + sum2 = sum2 + intgrl1 + a = b + b = xc + hwidth + intgrl1 = asigrl (asi1, a, b) + intgrl2 = asigrl (asi2, a, b) + sum1 = sum1 + (xc + hwidth) * intgrl1 - intgrl2 + sum2 = sum2 - intgrl1 + + # Return no center if sum2 is zero. + if (sum2 == 0.) + break + + # Limit dx change in one iteration to 1 pixel. + dx = max (-1., min (1., sum1 / abs (sum2))) + dxabs = abs (dx) + xc = xc + dx + ic = asieval (asi1, xc) + + # Check data range. Return no center if at edge of data. + if ((xc - hwidth < 1) || (xc + hwidth > npts)) + break + + # Convergence tests. + if (dxabs < EPSILON) + goto done_ + if (dxabs > dxlast + EPSILON1) { + dxcheck = dxcheck + 1 + if (dxcheck > MAX_DXCHECK) + break + } else { + dxcheck = 0 + dxlast = dxabs + } + } + + # If we get here then no center was found. + xc = INDEF + ic = INDEF + +done_ call asifree (asi1) + call asifree (asi2) +end + + +# RG_IPARAB -- Compute the coefficients of the parabola through three +# evenly spaced points. + +procedure rg_iparab (x, y, c) + +real x[NPARS_PARABOLA] #I input x values +real y[NPARS_PARABOLA] #I input y values +real c[NPARS_PARABOLA] #O computed coefficients + +begin + c[3] = (y[1]-y[2]) * (x[2]-x[3]) / (x[1]-x[2]) - (y[2]-y[3]) + c[3] = c[3] / ((x[1]**2-x[2]**2) * (x[2]-x[3]) / (x[1]-x[2]) - + (x[2]**2-x[3]**2)) + + c[2] = (y[1] - y[2]) - c[3] * (x[1]**2 - x[2]**2) + c[2] = c[2] / (x[1] - x[2]) + + c[1] = y[1] - c[2] * x[1] - c[3] * x[1]**2 +end + + +# RG_POLYFIT -- Evaluate an nth order polynomial. + +procedure rg_polyfit (x, nvars, p, np, z) + +real x #I position coordinate +int nvars #I number of variables +real p[ARB] #I coefficients of polynomial +int np #I number of parameters +real z #O function return + +int i +real r + +begin + r = 0.0 + do i = 2, np + r = r + x**(i-1) * p[i] + z = p[1] + r +end + + +# RG_DPOLYFIT -- Evaluate an nth order polynomial and its derivatives. + +procedure rg_dpolyfit (x, nvars, p, dp, np, z, der) + +real x #I position coordinate +int nvars #I number of variables +real p[ARB] #I coefficients of polynomial +real dp[ARB] #I parameter derivative increments +int np #I number of parameters +real z #O function value +real der[ARB] #O derivatives + +int i + +begin + der[1] = 1.0 + z = 0.0 + do i = 2, np { + der[i] = x ** (i-1) + z = z + x**(i-1) * p[i] + } + z = p[1] + z +end diff --git a/pkg/images/immatch/src/xregister/rgxgpars.x b/pkg/images/immatch/src/xregister/rgxgpars.x new file mode 100644 index 00000000..82943730 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxgpars.x @@ -0,0 +1,68 @@ +include "xregister.h" + +# RG_XGPARS -- Read in the XREGISTER task algorithm parameters. + +procedure rg_xgpars (xc) + +pointer xc #I pointer to the main structure + +int xlag, ylag, xwindow, ywindow, xcbox, ycbox +pointer sp, str +int clgwrd(), clgeti() +real clgetr() + +begin + # Allocate working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Initialize the correlation structure. + call rg_xinit (xc, clgwrd ("correlation", Memc[str], SZ_LINE, + XC_CTYPES)) + + # Fetch the initial shift information. + xlag = clgeti ("xlag") + ylag = clgeti ("ylag") + call rg_xseti (xc, IXLAG, xlag) + call rg_xseti (xc, IYLAG, ylag) + call rg_xseti (xc, XLAG, xlag) + call rg_xseti (xc, YLAG, ylag) + call rg_xseti (xc, DXLAG, clgeti ("dxlag")) + call rg_xseti (xc, DYLAG, clgeti ("dylag")) + + # Get the background value computation parameters. + call rg_xseti (xc, BACKGRD, clgwrd ("background", Memc[str], SZ_LINE, + XC_BTYPES)) + call rg_xsets (xc, BSTRING, Memc[str]) + call rg_xseti (xc, BORDER, clgeti ("border")) + call rg_xsetr (xc, LOREJECT, clgetr ("loreject")) + call rg_xsetr (xc, HIREJECT, clgetr ("hireject")) + call rg_xsetr (xc, APODIZE, clgetr ("apodize")) + call rg_xseti (xc, FILTER, clgwrd ("filter", Memc[str], SZ_LINE, + XC_FTYPES)) + call rg_xsets (xc, FSTRING, Memc[str]) + + # Get the window parameters and force the window size to be odd. + xwindow = clgeti ("xwindow") + if (mod (xwindow,2) == 0) + xwindow = xwindow + 1 + call rg_xseti (xc, XWINDOW, xwindow) + ywindow = clgeti ("ywindow") + if (mod (ywindow,2) == 0) + ywindow = ywindow + 1 + call rg_xseti (xc, YWINDOW, ywindow) + + # Get the peak fitting parameters. + call rg_xseti (xc, PFUNC, clgwrd ("function", Memc[str], SZ_LINE, + XC_PTYPES)) + xcbox = clgeti ("xcbox") + if (mod (xcbox,2) == 0) + xcbox = xcbox + 1 + call rg_xseti (xc, XCBOX, xcbox) + ycbox = clgeti ("ycbox") + if (mod (ycbox,2) == 0) + ycbox = ycbox + 1 + call rg_xseti (xc, YCBOX, ycbox) + + call sfree (sp) +end diff --git a/pkg/images/immatch/src/xregister/rgxicorr.x b/pkg/images/immatch/src/xregister/rgxicorr.x new file mode 100644 index 00000000..e96c6dec --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxicorr.x @@ -0,0 +1,583 @@ +include +include +include +include "xregister.h" + +define HELPFILE "immatch$src/xregister/xregister.key" +define OHELPFILE "immatch$src/xregister/oxregister.key" + +define XC_PCONTOUR 1 +define XC_PLINE 2 +define XC_PCOL 3 + + +# RG_XICORR -- Compute the shifts for each image interactively using +# cross-correlation techniques. + +int procedure rg_xicorr (imr, im1, im2, db, dformat, reglist, tfd, xc, gd, id) + +pointer imr #I/O pointer to the reference image +pointer im1 #I/O pointer to the input image +pointer im2 #I/O pointer to the output image +pointer db #I/O pointer to the shifts database file +int dformat #I is the shifts file in database format +int reglist #I/O the regions list descriptor +int tfd #I/O the transform file descriptor +pointer xc #I pointer to the cross-corrrelation structure +pointer gd #I the graphics stream pointer +pointer id #I the display stream pointer + +int newdata, newcross, newcenter, wcs, key, cplottype, newplot +int ip, ncolr, nliner +pointer sp, cmd +real xshift, yshift, wx, wy +int rg_xstati(), rg_icross(), clgcur(), rg_xgtverify(), rg_xgqverify() +int ctoi() +pointer rg_xstatp() + +begin + # Allocate working space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Initialize. + newdata = YES + newcross = YES + newcenter = YES + ncolr = (1 + rg_xstati (xc, XWINDOW)) / 2 + nliner = (1 + rg_xstati (xc, YWINDOW)) / 2 + cplottype = XC_PCONTOUR + newplot = YES + xshift = 0.0 + yshift = 0.0 + + # Compute the cross-correlation function for the first region + # and print the results. + if (rg_xstati (xc, NREGIONS) <= 0) { + call gclear (gd) + call printf ("The regions list is empty\n") + } else if (rg_icross (xc, imr, im1, rg_xstati (xc, CREGION)) != ERR) { + call rg_xcplot (xc, gd, ncolr, nliner, cplottype) + call rg_fit (xc, rg_xstati (xc, CREGION), gd, xshift, yshift) + call rg_xpwrec (xc, rg_xstati (xc, CREGION)) + newdata = NO + newcross = NO + newcenter = NO + newplot = NO + } else { + call gclear (gd) + call printf ( + "Error computing X-correlation function for region %d\n") + call pargi (rg_xstati (xc, CREGION)) + } + + + # 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, "") + + # Redraw the current plot. + case 'r': + newplot = YES + + # Draw a contour plot of the cross-correlation function. + case 'c': + if (cplottype != XC_PCONTOUR) + newplot = YES + ncolr = (rg_xstati (xc, XWINDOW) + 1) / 2 + nliner = (rg_xstati (xc, YWINDOW) + 1) / 2 + cplottype = XC_PCONTOUR + + # Plot a column of the cross-correlation function. + case 'x': + if (cplottype != XC_PCOL) + newplot = YES + if (cplottype == XC_PCONTOUR) { + ncolr = nint (wx) + nliner = nint (wy) + } else if (cplottype == XC_PLINE) { + ncolr = nint (wx) + } + cplottype = XC_PCOL + + # Plot a line of the cross-correlation function. + case 'y': + if (cplottype != XC_PLINE) + newplot = YES + if (cplottype == XC_PCONTOUR) { + ncolr = nint (wx) + nliner = nint (wy) + } else if (cplottype == XC_PCOL) { + ncolr = nint (wx) + } + cplottype = XC_PLINE + + # Quit the task gracefully. + case 'q': + if (rg_xgqverify ("xregister", db, dformat, xc, key) == YES) { + call sfree (sp) + return (rg_xgtverify (key)) + } + + # The Data overlay menu. + case 'o': + #call gdeactivate (gd, 0) + call rg_xoverlay (gd, xc, rg_xstati (xc, CREGION), imr, im1) + #call greactivate (gd, 0) + newplot = YES + + # 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_xcolon (gd, xc, imr, im1, im2, db, dformat, + tfd, reglist, Memc[cmd], newdata, newcross, + newcenter) + } else { + ip = ip + 1 + if (ctoi (Memc[cmd], ip, ncolr) <= 0) + ncolr = (1 + rg_xstati (xc, XWINDOW)) / 2 + cplottype = XC_PCOL + newplot = YES + } + case 'y': + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call rg_xcolon (gd, xc, imr, im1, im2, db, dformat, + tfd, reglist, Memc[cmd], newdata, newcross, + newcenter) + } else { + ip = ip + 1 + if (ctoi (Memc[cmd], ip, nliner) <= 0) + nliner = (1 + rg_xstati (xc, YWINDOW)) / 2 + cplottype = XC_PLINE + newplot = YES + } + default: + call rg_xcolon (gd, xc, imr, im1, im2, db, dformat, tfd, + reglist, Memc[cmd], newdata, newcross, newcenter) + } + + # Compute an image lag interactively. + case 't': + call gdeactivate (gd, 0) + call rg_itransform (xc, imr, im1, id) + newdata = YES; newcross = YES; newcenter = YES + call greactivate (gd, 0) + + # Write the parameters to the parameter file. + case 'w': + call rg_pxpars (xc) + + case 'f': + + if (rg_xstati (xc, NREGIONS) > 0) { + + if (newdata == YES) { + call rg_xcindefr (xc, rg_xstati(xc,CREGION)) + newdata = NO + } + + if (newcross == YES) { + call printf ( + "Recomputing X-correlation function ...\n") + if (rg_icross (xc, imr, im1, rg_xstati (xc, + CREGION)) != ERR) { + ncolr = (1 + rg_xstati (xc, XWINDOW)) / 2 + if (IM_NDIM(imr) == 1) + nliner = 1 + else + nliner = (1 + rg_xstati (xc, YWINDOW)) / 2 + call rg_xcplot (xc, gd, ncolr, nliner, cplottype) + call rg_fit (xc, rg_xstati (xc, CREGION), gd, + xshift, yshift) + call rg_xpwrec (xc, rg_xstati (xc, CREGION)) + newcross = NO + newcenter = NO + newplot = NO + } else { + call printf ( + "Error computing X-correlation function for region %d\n") + call pargi (rg_xstati (xc, CREGION)) + } + } + + if (newcenter == YES) { + call rg_fit (xc, rg_xstati (xc, CREGION), gd, + xshift, yshift) + call rg_xpwrec (xc, rg_xstati (xc, CREGION)) + newcenter = NO + } + + } else + call printf ("The regions list is empty\n") + + + + # Do nothing gracefully. + default: + call printf ("Unknown or ambiguous keystroke command\n") + } + + # Replot the correlation function. + if (newplot == YES) { + if (newdata == YES) { + call printf ( + "Warning: X-correlation function should be refit\n") + } else if (newcross == YES) { + call printf ( + "Warning: X-correlation function should be refit\n") + } else if (newcenter == YES) { + call printf ( + "Warning: X-correlation function should be refit\n") + } else if (rg_xstatp (xc, XCOR) != NULL) { + call rg_xcplot (xc, gd, ncolr, nliner, cplottype) + call rg_xpwrec (xc, rg_xstati (xc, CREGION)) + newplot = NO + } else { + call printf ( + "Warning: X-correlation function is undefined\n") + } + } + } + + call sfree (sp) +end + + +# RG_XOVERLAY -- The image overlay plot menu. + +procedure rg_xoverlay (gd, xc, nreg, imr, im1) + +pointer gd #I graphics stream pointer +pointer xc #I pointer to the crosscor structure +int nreg #I the current region number +pointer imr #I pointer to the reference image +pointer im1 #I pointer to the input image + +int ip, wcs, key, ixlag, iylag, ixshift, iyshift +int nrimcols, nrimlines, nimcols, nimlines, ncolr, ncoli, nliner, nlinei +pointer sp, cmd +real wx, wy, rxlag, rylag, xshift, yshift +int clgcur(), ctoi(), rg_xstati() +pointer rg_xstatp() + +begin + if (gd == NULL) + return + + nrimcols = IM_LEN(imr,1) + if (IM_NDIM(imr) == 1) + nrimlines = 1 + else + nrimlines = IM_LEN(imr,2) + nimcols = IM_LEN(im1,1) + if (IM_NDIM(im1) == 1) + nimlines = 1 + else + nimlines = IM_LEN(im1,2) + if (rg_xstati (xc, NREFPTS) > 0) { + wx = (1. + nrimcols) / 2.0 + wy = (1. + nrimlines) / 2.0 + call rg_etransform (xc, wx, wy, rxlag, rylag) + ixlag = rxlag - wx + iylag = rylag - wy + } else { + ixlag = rg_xstati (xc, XLAG) + iylag = rg_xstati (xc, YLAG) + } + xshift = -Memr[rg_xstatp(xc,XSHIFTS)+nreg-1] + yshift = -Memr[rg_xstatp(xc,YSHIFTS)+nreg-1] + + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + + switch (key) { + + # Print the help menu. + case '?': + call gdeactivate (gd, 0) + call pagefile (OHELPFILE, "") + call greactivate (gd, 0) + + # Quit. + case 'q': + break + + # Plot the same line of the reference and input image. + case 'l': + call rg_xpline (gd, imr, im1, nint (wy), 0, 0) + + # Plot the same column of the reference and input image + case 'c': + call rg_xpcol (gd, imr, im1, nint (wx), 0, 0) + + case 'y': + call rg_xpline (gd, imr, im1, nint (wy), ixlag, iylag) + + case 'x': + call rg_xpcol (gd, imr, im1, nint (wx), ixlag, iylag) + + case 'h': + call rg_xpline (gd, imr, im1, nint (wy), nint (xshift), + nint (yshift)) + + case 'v': + call rg_xpcol (gd, imr, im1, nint (wx), nint (xshift), + nint (yshift)) + + case ':': + ip = 1 + call rg_cokeys (Memc[cmd], ip, SZ_LINE, key) + switch (key) { + case 'l': + ixshift = 0 + if (ctoi (Memc[cmd], ip, nliner) <= 0) + nliner = (1 + nrimlines) / 2 + nliner = max (1, min (nliner, nrimlines)) + if (ctoi (Memc[cmd], ip, nlinei) <= 0) + nlinei = nliner + iyshift = nlinei - nliner + call rg_xpline (gd, imr, im1, nliner, ixshift, iyshift) + + case 'c': + if (ctoi (Memc[cmd], ip, ncolr) <= 0) + ncolr = (1 + nrimcols) / 2 + ncolr = max (1, min (ncolr, nrimcols)) + if (ctoi (Memc[cmd], ip, ncoli) <= 0) + ncoli = ncolr + ncoli = max (1, min (ncoli, nimcols)) + ixshift = ncoli - ncolr + iyshift = 0 + call rg_xpcol (gd, imr, im1, ncolr, ixshift, iyshift) + + case 'y': + if (ctoi (Memc[cmd], ip, nliner) <= 0) + nliner = (1 + nrimlines) / 2 + nliner = max (1, min (nliner, nrimlines)) + call rg_xpline (gd, imr, im1, nliner, ixlag, iylag) + + case 'x': + if (ctoi (Memc[cmd], ip, ncolr) <= 0) + ncolr = (1 + nrimcols) / 2 + ncolr = max (1, min (ncolr, nrimcols)) + call rg_xpcol (gd, imr, im1, ncolr, ixlag, iylag) + + case 'h': + if (ctoi (Memc[cmd], ip, nliner) <= 0) + nliner = (1 + nrimlines) / 2 + nliner = max (1, min (nliner, nrimlines)) + call rg_xpline (gd, imr, im1, nliner, nint (xshift), + nint (yshift)) + + case 'v': + if (ctoi (Memc[cmd], ip, ncolr) <= 0) + ncolr = (1 + nrimcols) / 2 + ncolr = max (1, min (ncolr, nrimcols)) + call rg_xpcol (gd, imr, im1, ncolr, nint (xshift), + nint (yshift)) + default: + call printf ("Ambiguous or unknown overlay menu command\n") + } + case 'g': + while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + if (key == 'q') + break + } + default: + call printf ("Ambiguous or unknown overlay menu command\n") + } + + } + + call sfree (sp) +end + + +# RG_XCPLOT -- Draw the default plot of the cross-correlation function. + +procedure rg_xcplot (xc, gd, col, line, plottype) + +pointer xc #I pointer to cross-correlation structure +pointer gd #I pointer to the graphics stream +int col #I column of cross-correlation function to plot +int line #I line of cross-correlation function to plot +int plottype #I the default plot type + +int nreg, xwindow, ywindow +pointer sp, title, str, prc1, prc2, prl1, prl2 +int rg_xstati(), strlen() +pointer rg_xstatp() + +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) + + # Get the regions. + nreg = rg_xstati (xc, CREGION) + prc1 = rg_xstatp (xc, RC1) + prc2 = rg_xstatp (xc, RC2) + prl1 = rg_xstatp (xc, RL1) + prl2 = rg_xstatp (xc, RL2) + + # Initialize the window size. + xwindow = rg_xstati (xc, XWINDOW) + if ((Memi[prl2+nreg-1] - Memi[prl1+nreg-1] + 1) == 1) + ywindow = 1 + else + ywindow = rg_xstati (xc, YWINDOW) + + # Construct a title. + call sprintf (Memc[title], SZ_LINE, + "Reference: %s Image: %s Region: [%d:%d,%d:%d]") + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_FNAME) + call pargstr (Memc[str]) + call rg_xstats (xc, IMAGE, Memc[str], SZ_FNAME) + call pargstr (Memc[str]) + call pargi (Memi[prc1+nreg-1]) + call pargi (Memi[prc2+nreg-1]) + call pargi (Memi[prl1+nreg-1]) + call pargi (Memi[prl2+nreg-1]) + + # Draw the plot. + if (ywindow == 1) { + call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE, + "\nX-Correlation Function: line %d") + call pargi (1) + call rg_xcpline (gd, Memc[title], Memr[rg_xstatp(xc,XCOR)], + xwindow, ywindow, 1) + } else { + switch (plottype) { + case XC_PCONTOUR: + call rg_contour (gd, "X-Correlation Function", Memc[title], + Memr[rg_xstatp (xc, XCOR)], xwindow, ywindow) + case XC_PLINE: + call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE, + "\nX-Correlation Function: line %d") + call pargi (line) + call rg_xcpline (gd, Memc[title], Memr[rg_xstatp(xc,XCOR)], + xwindow, ywindow, line) + case XC_PCOL: + call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE, + "\nX-Correlation Function: column %d") + call pargi (col) + call rg_xcpcol (gd, Memc[title], Memr[rg_xstatp(xc,XCOR)], + xwindow, ywindow, col) + default: + call rg_contour (gd, "X-Correlation Function", Memc[title], + Memr[rg_xstatp (xc, XCOR)], xwindow, ywindow) + } + } + + call sfree (sp) +end + + +# RG_COKEYS -- Fetch the first keystroke of a colon command. + +procedure rg_cokeys (cmd, ip, maxch, key) + +char cmd[ARB] #I the command string +int ip #I/O pointer into the command string +int maxch #I maximum number of characters +int key #O the keystroke + +begin + ip = 1 + while (IS_WHITE(cmd[ip]) && cmd[ip] != EOS && ip <= maxch) + ip = ip + 1 + + if (cmd[ip] == EOS && ip > maxch) + key = EOS + else { + key = cmd[ip] + ip = ip + 1 + } +end + + +define QUERY "Hit [return=continue, n=next image, q=quit, w=quit and update parameters]: " + +# RG_XGQVERIFY -- Print a message on 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_xgqverify (task, db, dformat, rg, ch) + +char task[ARB] #I the calling task name +pointer db #I pointer to the shifts database file +int dformat #I is the shifts file in database format +pointer rg #I pointer to the task structure +int ch #I the input keystroke command + +int wcs, stat +pointer sp, cmd +real wx, wy +bool streq() +int clgcur() + +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) + #call flush (STDOUT) + if (clgcur ("gcommands", wx, wy, wcs, ch, Memc[cmd], SZ_LINE) == EOF) + ; + + # Process the command. + if (ch == 'q') { + call rg_xwrec (db, dformat, rg) + stat = YES + } else if (ch == 'w') { + call rg_xwrec (db, dformat, rg) + if (streq ("xregister", task)) + call rg_pxpars (rg) + stat = YES + } else if (ch == 'n') { + call rg_xwrec (db, dformat, rg) + stat = YES + } else { + stat = NO + } + + call sfree (sp) + return (stat) +end + + +# RG_XGTVERIFY -- Verify whether or not the user truly wishes to quit the +# task. + +int procedure rg_xgtverify (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 diff --git a/pkg/images/immatch/src/xregister/rgximshift.x b/pkg/images/immatch/src/xregister/rgximshift.x new file mode 100644 index 00000000..08cb3f62 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgximshift.x @@ -0,0 +1,391 @@ +include +include +include + +define NYOUT 16 # number of lines output at once +define NMARGIN 3 # number of boundary pixels required +define NMARGIN_SPLINE3 16 # number of spline boundary pixels required + + +# RG_XSHIFTIM - Shift a 1 or 2D image by a fractional pixel amount +# x and y + +procedure rg_xshiftim (im1, im2, xshift, yshift, interpstr, boundary_type, + constant) + +pointer im1 #I pointer to input image +pointer im2 #I pointer to output image +real xshift #I shift in x direction +real yshift #I shift in y direction +char interpstr[ARB] #I type of interpolant +int boundary_type #I type of boundary extension +real constant #I value of constant for boundary extension + +int interp_type +pointer sp, str +bool fp_equalr() +int strdic() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + interp_type = strdic (interpstr, Memc[str], SZ_FNAME, II_BFUNCTIONS) + + if (interp_type == II_NEAREST) + call rg_xishiftim (im1, im2, nint (xshift), nint (yshift), + interp_type, boundary_type, constant) + else if (fp_equalr (xshift, real (int (xshift))) && fp_equalr (yshift, + real (int (xshift)))) + call rg_xishiftim (im1, im2, int (xshift), int (yshift), + interp_type, boundary_type, constant) + else + call rg_xfshiftim (im1, im2, xshift, yshift, interpstr, + boundary_type, constant) + call sfree (sp) +end + + +# RG_XISHIFTIM -- Shift a 2-D image by integral pixels in x and y. + +procedure rg_xishiftim (im1, im2, nxshift, nyshift, interp_type, boundary_type, + constant) + +pointer im1 #I pointer to the input image +pointer im2 #I pointer to the output image +int nxshift, nyshift #I shift in x and y +int interp_type #I type of interpolant +int boundary_type #I type of boundary extension +real constant #I constant for boundary extension + +int ixshift, iyshift +pointer buf1, buf2 +long v[IM_MAXDIM] +int ncols, nlines, nbpix +int i, x1col, x2col, yline + +int impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx() +pointer imgs2s(), imgs2i(), imgs2l(), imgs2r(), imgs2d(), imgs2x() +errchk impnls, impnli, impnll, impnlr, impnld, impnlx +errchk imgs2s, imgs2i, imgs2l, imgs2r, imgs2d, imgs2x +string wrerr "ISHIFTXY: Error writing in image." + +begin + ixshift = nxshift + iyshift = nyshift + + ncols = IM_LEN(im1,1) + nlines = IM_LEN(im1,2) + + # Cannot shift off image. + if (ixshift < -ncols || ixshift > ncols) + call error (3, "ISHIFTXY: X shift out of bounds.") + if (iyshift < -nlines || iyshift > nlines) + call error (4, "ISHIFTXY: Y shift out of bounds.") + + # Calculate the shift. + switch (boundary_type) { + case BT_CONSTANT,BT_REFLECT,BT_NEAREST: + ixshift = min (ncols, max (-ncols, ixshift)) + iyshift = min (nlines, max (-nlines, iyshift)) + case BT_WRAP: + ixshift = mod (ixshift, ncols) + iyshift = mod (iyshift, nlines) + } + + # Set the boundary extension values. + nbpix = max (abs (ixshift), abs (iyshift)) + call imseti (im1, IM_NBNDRYPIX, nbpix) + call imseti (im1, IM_TYBNDRY, boundary_type) + if (boundary_type == BT_CONSTANT) + call imsetr (im1, IM_BNDRYPIXVAL, constant) + + # Get column boundaries in the input image. + x1col = max (-ncols + 1, - ixshift + 1) + x2col = min (2 * ncols, ncols - ixshift) + + call amovkl (long (1), v, IM_MAXDIM) + + # Shift the image using the appropriate data type operators. + switch (IM_PIXTYPE(im1)) { + case TY_SHORT: + do i = 1, nlines { + if (impnls (im2, buf2, v) == EOF) + call error (5, wrerr) + yline = i - iyshift + buf1 = imgs2s (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (5, wrerr) + call amovs (Mems[buf1], Mems[buf2], ncols) + } + case TY_INT: + do i = 1, nlines { + if (impnli (im2, buf2, v) == EOF) + call error (5, wrerr) + yline = i - iyshift + buf1 = imgs2i (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (5, wrerr) + call amovi (Memi[buf1], Memi[buf2], ncols) + } + case TY_USHORT, TY_LONG: + do i = 1, nlines { + if (impnll (im2, buf2, v) == EOF) + call error (5, wrerr) + yline = i - iyshift + buf1 = imgs2l (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (5, wrerr) + call amovl (Meml[buf1], Meml[buf2], ncols) + } + case TY_REAL: + do i = 1, nlines { + if (impnlr (im2, buf2, v) == EOF) + call error (5, wrerr) + yline = i - iyshift + buf1 = imgs2r (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (5, wrerr) + call amovr (Memr[buf1], Memr[buf2], ncols) + } + case TY_DOUBLE: + do i = 1, nlines { + if (impnld (im2, buf2, v) == EOF) + call error (0, wrerr) + yline = i - iyshift + buf1 = imgs2d (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (0, wrerr) + call amovd (Memd[buf1], Memd[buf2], ncols) + } + case TY_COMPLEX: + do i = 1, nlines { + if (impnlx (im2, buf2, v) == EOF) + call error (0, wrerr) + yline = i - iyshift + buf1 = imgs2x (im1, x1col, x2col, yline, yline) + if (buf1 == EOF) + call error (0, wrerr) + call amovx (Memx[buf1], Memx[buf2], ncols) + } + default: + call error (6, "ISHIFTXY: Unknown IRAF type.") + } +end + + + +# RG_XFSHIFTIM -- Shift a 1 or 2D image by a fractional pixel amount +# in x and y. + +procedure rg_xfshiftim (im1, im2, xshift, yshift, interpstr, boundary_type, + constant) + +pointer im1 #I pointer to input image +pointer im2 #I pointer to output image +real xshift #I shift in x direction +real yshift #I shift in y direction +char interpstr[ARB] #I type of interpolant +int boundary_type #I type of boundary extension +real constant #I value of constant for boundary extension + +int i, interp_type, nsinc, nincr +int ncols, nlines, nbpix, fstline, lstline, nxymargin +int cin1, cin2, nxin, lin1, lin2, nyin +int lout1, lout2, nyout +real xshft, yshft, deltax, deltay, dx, dy, cx, ly +pointer sp, x, y, msi, sinbuf, soutbuf +bool fp_equalr() +int msigeti() +pointer imps2r() + +errchk imgs2r, imps2r +errchk msiinit, msifree, msifit, msigrid +errchk smark, salloc, sfree + +begin + ncols = IM_LEN(im1,1) + nlines = IM_LEN(im1,2) + + # Check for out of bounds shift. + if (xshift < -ncols || xshift > ncols) + call error (0, "XC_SHIFTIM: X shift out of bounds.") + if (yshift < -nlines || yshift > nlines) + call error (0, "XC_SHIFTIM: Y shift out of bounds.") + + # Get the real shift. + if (boundary_type == BT_WRAP) { + xshft = mod (xshift, real (ncols)) + yshft = mod (yshift, real (nlines)) + } else { + xshft = xshift + yshft = yshift + } + + # Allocate temporary space. + call smark (sp) + call salloc (x, 2 * ncols, TY_REAL) + call salloc (y, 2 * nlines, TY_REAL) + sinbuf = NULL + + # Define the x and y interpolation coordinates. + dx = abs (xshft - int (xshft)) + if (fp_equalr (dx, 0.0)) + deltax = 0.0 + else if (xshft > 0.) + deltax = 1. - dx + else + deltax = dx + dy = abs (yshft - int (yshft)) + if (fp_equalr (dy, 0.0)) + deltay = 0.0 + else if (yshft > 0.) + deltay = 1. - dy + else + deltay = dy + + # Initialize the 2-D interpolation routines. + call msitype (interpstr, interp_type, nsinc, nincr, cx) + if (interp_type == II_BILSINC || interp_type == II_BISINC) + call msisinit (msi, interp_type, nsinc, 1, 1, + deltax - nint (deltax), deltay - nint (deltay), 0.0) + else + call msisinit (msi, interp_type, nsinc, 1, 1, cx, cx, 0.0) + + # Set boundary extension parameters. + if (interp_type == II_BISPLINE3) + nxymargin = NMARGIN_SPLINE3 + else if (interp_type == II_BISINC || interp_type == II_BILSINC) + nxymargin = msigeti (msi, II_MSINSINC) + else + nxymargin = NMARGIN + nbpix = max (int (abs(xshft)+1.0), int (abs(yshft)+1.0)) + nxymargin + call imseti (im1, IM_NBNDRYPIX, nbpix) + call imseti (im1, IM_TYBNDRY, boundary_type) + if (boundary_type == BT_CONSTANT) + call imsetr (im1, IM_BNDRYPIXVAL, constant) + + # Define the x interpolation coordinates. + deltax = deltax + nxymargin + if (interp_type == II_BIDRIZZLE) { + do i = 1, ncols { + Memr[x+2*i-2] = i + deltax - 0.5 + Memr[x+2*i-1] = i + deltax + 0.5 + } + } else { + do i = 1, ncols + Memr[x+i-1] = i + deltax + } + + # Define the y interpolation coordinates. + deltay = deltay + nxymargin + if (interp_type == II_BIDRIZZLE) { + do i = 1, NYOUT { + Memr[y+2*i-2] = i + deltay - 0.5 + Memr[y+2*i-1] = i + deltay + 0.5 + } + } else { + do i = 1, NYOUT + Memr[y+i-1] = i + deltay + } + + # Define column range in the input image. + cx = 1. - nxymargin - xshft + if ((cx <= 0.) && (! fp_equalr (dx, 0.0))) + cin1 = int (cx) - 1 + else + cin1 = int (cx) + cin2 = ncols - xshft + nxymargin + 1 + nxin = cin2 - cin1 + 1 + + # Loop over output sections. + for (lout1 = 1; lout1 <= nlines; lout1 = lout1 + NYOUT) { + + # Define range of output lines. + lout2 = min (lout1 + NYOUT - 1, nlines) + nyout = lout2 - lout1 + 1 + + # Define correspoding range of input lines. + ly = lout1 - nxymargin - yshft + if ((ly <= 0) && (! fp_equalr (dy, 0.0))) + lin1 = int (ly) - 1 + else + lin1 = int (ly) + lin2 = lout2 - yshft + nxymargin + 1 + nyin = lin2 - lin1 + 1 + + # Get appropriate input image section and compute the coefficients. + if ((sinbuf == NULL) || (lin1 < fstline) || (lin2 > lstline)) { + fstline = lin1 + lstline = lin2 + call rg_buf (im1, cin1, cin2, lin1, lin2, sinbuf) + call msifit (msi, Memr[sinbuf], nxin, nyin, nxin) + } + + # Output the image section. + soutbuf = imps2r (im2, 1, ncols, lout1, lout2) + if (soutbuf == EOF) + call error (0, "GSHIFTXY: Error writing output image.") + + # Evaluate the interpolant. + call msigrid (msi, Memr[x], Memr[y], Memr[soutbuf], ncols, nyout, + ncols) + } + + call msifree (msi) + call sfree (sp) +end + + +# RG_BUF -- Procedure to provide a buffer of image lines with minimum reads + +procedure rg_buf (im, col1, col2, line1, line2, buf) + +pointer im #I pointer to input image +int col1, col2 #I column range of input buffer +int line1, line2 #I line range of input buffer +pointer buf #I buffer + +int i, ncols, nlines, nclast, llast1, llast2, nllast +pointer buf1, buf2 + +pointer imgs2r() + +begin + ncols = col2 - col1 + 1 + nlines = line2 - line1 + 1 + + if (buf == NULL) { + call malloc (buf, ncols * nlines, TY_REAL) + llast1 = line1 - nlines + llast2 = line2 - nlines + } else if ((nlines != nllast) || (ncols != nclast)) { + call realloc (buf, ncols * nlines, TY_REAL) + llast1 = line1 - nlines + llast2 = line2 - nlines + } + + if (line1 < llast1) { + do i = line2, line1, -1 { + if (i > llast1) + buf1 = buf + (i - llast1) * ncols + else + buf1 = imgs2r (im, col1, col2, i, i) + buf2 = buf + (i - line1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } else if (line2 > llast2) { + do i = line1, line2 { + if (i < llast2) + buf1 = buf + (i - llast1) * ncols + else + buf1 = imgs2r (im, col1, col2, i, i) + buf2 = buf + (i - line1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } + + llast1 = line1 + llast2 = line2 + nclast = ncols + nllast = nlines +end diff --git a/pkg/images/immatch/src/xregister/rgxplot.x b/pkg/images/immatch/src/xregister/rgxplot.x new file mode 100644 index 00000000..8b347ab5 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxplot.x @@ -0,0 +1,317 @@ +include +include + +# RG_XPLINE -- Plot a line of reference and input image. + +procedure rg_xpline (gd, imr, im, nliner, xshift, yshift) + +pointer gd #I pointer to the graphics stream +pointer imr #I pointer to the reference image +pointer im #I pointer to the image +int nliner #I the reference line +int xshift #I x shift +int yshift #I y shift + +int i, rncols, rnlines, incols, inlines +pointer sp, title, xr, xi, ptrr, ptri +real ymin, ymax, tymin, tymax +int strlen() +pointer imgl1r(), imgl2r() + +begin + # Return if no graphics stream. + if (gd == NULL) + return + + # Check for valid line number. + rncols = IM_LEN(imr,1) + if (IM_NDIM(imr) == 1) + rnlines = 1 + else + rnlines = IM_LEN(imr,2) + incols = IM_LEN(im,1) + if (IM_NDIM(im) == 1) + inlines = 1 + else + inlines = IM_LEN(im,2) + if ((nliner < 1) || (nliner > rnlines)) + return + if (((nliner + yshift) < 1) || ((nliner + yshift) > inlines)) + return + + # Allocate working space. + call smark (sp) + call salloc (title, SZ_LINE, TY_CHAR) + call salloc (xr, rncols, TY_REAL) + call salloc (xi, rncols, TY_REAL) + + # Initialize the x data data. + do i = 1, rncols { + Memr[xr+i-1] = i + Memr[xi+i-1] = i - xshift + } + + # Initalize the y data. + if (IM_NDIM(imr) == 1) + ptrr = imgl1r (imr) + else + ptrr = imgl2r (imr, nliner) + if (IM_NDIM(im) == 1) + ptri = imgl1r (im) + else + ptri = imgl2r (im, nliner + yshift) + call alimr (Memr[ptrr], rncols, ymin, ymax) + call alimr (Memr[ptri], incols, tymin, tymax) + ymin = min (ymin, tymin) + ymax = max (ymax, tymax) + + # Construct the title. + call sprintf (Memc[title], SZ_LINE, + "Refimage: %s Image: %s\n") + call pargstr (IM_HDRFILE(imr)) + call pargstr (IM_HDRFILE(im)) + call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE, + "Refline (solid): %d Inline (dashed): %d Xlag: %d Ylag: %d") + call pargi (nliner) + call pargi (nliner + yshift) + call pargi (xshift) + call pargi (yshift) + + # Set up the axes labels and window. + call gclear (gd) + call gswind (gd, 1.0, real(rncols), ymin, ymax) + call glabax (gd, Memc[title], "Column Number", "Counts") + + # Plot the two lines. + call gseti (gd, G_PLTYPE, GL_SOLID) + call gpline (gd, Memr[xr], Memr[ptrr], rncols) + call gseti (gd, G_PLTYPE, GL_DASHED) + call gpline (gd, Memr[xi], Memr[ptri], incols) + call gflush (gd) + + call sfree (sp) +end + + +# RG_XPCOL -- Plot a column in the reference and input image. + +procedure rg_xpcol (gd, imr, im, ncolr, xshift, yshift) + +pointer gd #I pointer to the graphics stream +pointer imr #I pointer to the reference image +pointer im #I pointer to the image +int ncolr #I the line number +int xshift #I xshift to be applied +int yshift #I yshift to be applied + +int i, rncols, rnlines, incols, inlines +pointer sp, title, xr, xi, ptrr, ptri +real ymin, ymax, tymin, tymax +int strlen() +pointer imgs1r(), imgs2r() + +begin + # Return if no graphics stream. + if (gd == NULL) + return + + # Check for valid column number. + rncols = IM_LEN(imr,1) + if (IM_NDIM(imr) == 1) + rnlines = 1 + else + rnlines = IM_LEN(imr,2) + incols = IM_LEN(im,1) + if (IM_NDIM(im) == 1) + inlines = 1 + else + inlines = IM_LEN(im,2) + if ((ncolr < 1) || (ncolr > rncols)) + return + if (((ncolr - xshift) < 1) || ((ncolr - xshift) > incols)) + return + + # Allocate valid working space. + call smark (sp) + call salloc (title, SZ_LINE, TY_CHAR) + call salloc (xr, rnlines, TY_REAL) + call salloc (xi, inlines, TY_REAL) + + # Initialize the data. + do i = 1, rnlines { + Memr[xr+i-1] = i + Memr[xi+i-1] = i - yshift + } + if (IM_NDIM(imr) == 1) + ptrr = imgs1r (imr, ncolr, ncolr) + else + ptrr = imgs2r (imr, ncolr, ncolr, 1, rnlines) + if (IM_NDIM(im) == 1) + ptri = imgs1r (im, ncolr + xshift, ncolr + xshift) + else + ptri = imgs2r (im, ncolr + xshift, ncolr + xshift, 1, inlines) + call alimr (Memr[ptrr], rnlines, ymin, ymax) + call alimr (Memr[ptri], inlines, tymin, tymax) + ymin = min (ymin, tymin) + ymax = max (ymax, tymax) + + # Construct the title. + call sprintf (Memc[title], SZ_LINE, "Refimage: %s Image: %s\n") + call pargstr (IM_HDRFILE(imr)) + call pargstr (IM_HDRFILE(im)) + call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE, + "Refcol (solid): %d Imcol (dashed): %d Xlag: %d Ylag: %d") + call pargi (ncolr) + call pargi (ncolr + xshift) + call pargi (xshift) + call pargi (yshift) + + # Set up the labels and the axes. + call gclear (gd) + call gswind (gd, 1.0, real (rnlines), ymin, ymax) + call glabax (gd, Memc[title], "Line Number", "Counts") + + # Plot the profile. + call gseti (gd, G_PLTYPE, GL_SOLID) + call gpline (gd, Memr[xr], Memr[ptrr], rnlines) + call gseti (gd, G_PLTYPE, GL_DASHED) + call gpline (gd, Memr[xi], Memr[ptri], rnlines) + call gflush (gd) + + call sfree (sp) +end + + +# RG_XCPLINE -- Plot a line of the 2D correlation function. + +procedure rg_xcpline (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_XCPCOL -- Plot a column of the cross-correlation function. + +procedure rg_xcpcol (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 gflush (gd) + + call sfree (sp) +end + + +# RG_XMKPEAK -- Procedure to mark the peak from a correlation function +# contour plot. + +procedure rg_xmkpeak (gd, xwindow, ywindow, xshift, yshift) + +pointer gd #I pointer to the graphics stream +int xwindow #I x dimension of correlation function +int ywindow #I y dimension of correlation function +real xshift #O x shift +real yshift #O y shift + +int wcs, key +pointer sp, cmd +real wx, wy +int clgcur() + +begin + if (gd == NULL) + return + + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call printf ("Mark peak of the cross correlation function\n") + if (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF) + ; + if (wx < 1.0 || wx > real (xwindow) || wy < 1.0 || wy > + real (ywindow)) { + xshift = 0.0 + yshift = 0.0 + } else { + xshift = wx - (1 + xwindow) / 2 + yshift = wy - (1 + ywindow) / 2 + } + + call sfree (sp) +end diff --git a/pkg/images/immatch/src/xregister/rgxppars.x b/pkg/images/immatch/src/xregister/rgxppars.x new file mode 100644 index 00000000..2dc6aafd --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxppars.x @@ -0,0 +1,49 @@ +include "xregister.h" + +# RG_PXPARS -- Update the cross-correlation algorithm parameters. + +procedure rg_pxpars (xc) + +pointer xc #I pointer to the cross-correlation structure + +pointer sp, str +int rg_xstati() +real rg_xstatr() + +begin + # Allocate working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Define the regions. + call rg_xstats (xc, REGIONS, Memc[str], SZ_LINE) + call clpstr ("regions", Memc[str]) + call clputi ("xlag", rg_xstati (xc, XLAG)) + call clputi ("ylag", rg_xstati (xc, YLAG)) + call clputi ("dxlag", rg_xstati (xc, DXLAG)) + call clputi ("dylag", rg_xstati (xc, DYLAG)) + + # Store the background fitting parameters. + call rg_xstats (xc, BSTRING, Memc[str], SZ_LINE) + call clpstr ("background", Memc[str]) + call clputi ("border", rg_xstati (xc, BORDER)) + call clputr ("loreject", rg_xstatr (xc, LOREJECT)) + call clputr ("hireject", rg_xstatr (xc, HIREJECT)) + call clputr ("apodize", rg_xstatr (xc, APODIZE)) + call rg_xstats (xc, FSTRING, Memc[str], SZ_LINE) + call clpstr ("filter", Memc[str]) + + # Store the cross-correlation parameters. + call rg_xstats (xc, CSTRING, Memc[str], SZ_LINE) + call clpstr ("correlation", Memc[str]) + call clputi ("xwindow", rg_xstati (xc, XWINDOW)) + call clputi ("ywindow", rg_xstati (xc, YWINDOW)) + + # Store the peak centering parameters. + call rg_xstats (xc, PSTRING, Memc[str], SZ_LINE) + call clpstr ("function", Memc[str]) + call clputi ("xcbox", rg_xstati (xc, XCBOX)) + call clputi ("ycbox", rg_xstati (xc, YCBOX)) + + call sfree (sp) +end diff --git a/pkg/images/immatch/src/xregister/rgxregions.x b/pkg/images/immatch/src/xregister/rgxregions.x new file mode 100644 index 00000000..ed682f61 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxregions.x @@ -0,0 +1,459 @@ +include +include +include +include "xregister.h" + +# RG_XREGIONS -- Decode the image sections into regions. If the sections string +# is NULL then the regions list is initially empty and depending on the mode +# of the task, XREGISTER will or will not complain.Otherwise the image +# sections specified in the sections string or file are decoded into a +# regions list. + +int procedure rg_xregions (list, im, xc, rp) + +int list #I pointer to the regions list +pointer im #I pointer to the reference image +pointer xc #I pointer to the cross-correlation structure +int rp #I index of the current region + +int fd, nregions +pointer sp, fname, regions +int rg_xgrid(), rg_xgregions(), rg_xrregions(), rg_xstati(), fntgfnb() +int open() +errchk fntgfnb(), open(), close() + +begin + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (regions, SZ_LINE, TY_CHAR) + + call rg_xstats (xc, REGIONS, Memc[regions], SZ_LINE) + if (rp < 1 || rp > MAX_NREGIONS || Memc[regions] == EOS) { + nregions = 0 + } else if (rg_xgrid (im, xc, rp, MAX_NREGIONS) > 0) { + nregions = rg_xstati (xc, NREGIONS) + } else if (rg_xgregions (im, xc, rp, MAX_NREGIONS) > 0) { + nregions = rg_xstati (xc, NREGIONS) + } else if (list != NULL) { + iferr { + if (fntgfnb (list, Memc[fname], SZ_FNAME) != EOF) { + fd = open (Memc[fname], READ_ONLY, TEXT_FILE) + nregions= rg_xrregions (fd, im, xc, rp, MAX_NREGIONS) + call close (fd) + } + } then + nregions = 0 + } else + nregions = 0 + + call sfree (sp) + + return (nregions) +end + + +# RG_XMKREGIONS -- Create a list of regions by marking image sections +# on the image display. + +int procedure rg_xmkregions (im, xc, rp, max_nregions, regions, maxch) + +pointer im #I pointer to the reference image +pointer xc #I pointer to the cross-correlation structure +int rp #I index of the current region +int max_nregions #I the maximum number of regions +char regions[ARB] #O the output regions string +int maxch #I maximum size of the output regions string + +int op, nregions, wcs, key +pointer sp, region, section, cmd +real xll, yll, xur, yur +int rg_xstati(), clgcur(), gstrcpy() +pointer rg_xstatp() + +begin + # Allocate working space. + call smark (sp) + call salloc (region, SZ_LINE, TY_CHAR) + call salloc (section, SZ_LINE, TY_CHAR) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Allocate the arrays to hold the regions information, + call rg_xrealloc (xc, max_nregions) + + # Initialize. + nregions = min (rp-1, rg_xstati (xc, NREGIONS)) + op = 1 + + # Mark the sections on the display. + while (nregions < max_nregions) { + + call printf ("Mark lower left corner of region %d [q to quit].\n") + call pargi (nregions + 1) + if (clgcur ("icommands", xll, yll, wcs, key, Memc[cmd], + SZ_LINE) == EOF) + break + if (key == 'q') + break + + call printf ("Mark upper right corner of region %d [q to quit].\n") + call pargi (nregions + 1) + if (clgcur ("icommands", xur, yur, wcs, key, Memc[cmd], + SZ_LINE) == EOF) + break + if (key == 'q') + break + + if (xll < 1.0 || xur > IM_LEN(im,1) || yll < 1.0 || yur > + IM_LEN(im,2)) + break + + Memi[rg_xstatp(xc,RC1)+nregions] = nint (xll) + Memi[rg_xstatp(xc,RC2)+nregions] = nint (xur) + Memi[rg_xstatp(xc,RL1)+nregions] = nint (yll) + Memi[rg_xstatp(xc,RL2)+nregions] = nint (yur) + Memr[rg_xstatp(xc,RZERO)+nregions] = INDEFR + Memr[rg_xstatp(xc,RXSLOPE)+nregions] = INDEFR + Memr[rg_xstatp(xc,RYSLOPE)+nregions] = INDEFR + Memr[rg_xstatp(xc,XSHIFTS)+nregions] = INDEFR + Memr[rg_xstatp(xc,YSHIFTS)+nregions] = INDEFR + nregions = nregions + 1 + + # Write the first 9 regions into the regions string. + call sprintf (Memc[cmd], SZ_LINE, "[%d:%d,%d:%d] ") + call pargi (nint (xll)) + call pargi (nint (xur)) + call pargi (nint (yll)) + call pargi (nint (yur)) + op = op + gstrcpy (Memc[cmd], regions[op], maxch - op + 1) + } + call printf ("\n") + + # Reallocate the correct amount of space. + call rg_xseti (xc, NREGIONS, nregions) + if (nregions > 0) + call rg_xrealloc (xc, nregions) + else + call rg_xrfree (xc) + + call sfree (sp) + + return (nregions) +end + + +# RG_XGRID - Decode the regions from a grid specification. + +int procedure rg_xgrid (im, xc, rp, max_nregions) + +pointer im #I pointer to the reference image +pointer xc #I pointer to the cross-correlation structure +int rp #I index of the current region +int max_nregions #I the maximum number of regions + +int i, istart, iend, j, jstart, jend, ncols, nlines, nxsample, nysample +int nxcols, nylines, nregions +pointer sp, region, section +int rg_xstati(), nscan(), strcmp() +pointer rg_xstatp() + +begin + # Allocate working space. + call smark (sp) + call salloc (region, SZ_LINE, TY_CHAR) + call salloc (section, SZ_LINE, TY_CHAR) + + # Allocate the arrays to hold the regions information, + call rg_xrealloc (xc, max_nregions) + + # Initialize. + call rg_xstats (xc, REGIONS, Memc[region], SZ_LINE) + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + nregions = min (rp - 1, rg_xstati (xc, NREGIONS)) + + # Decode the grid specification. + call sscan (Memc[region]) + call gargwrd (Memc[section], SZ_LINE) + call gargi (nxsample) + call gargi (nysample) + if ((nscan() != 3) || (strcmp (Memc[section], "grid") != 0)) { + call sfree (sp) + return (nregions) + } + + # Decode the regions. + if ((nxsample * nysample) > max_nregions) { + nxsample = nint (sqrt (real (max_nregions) * real (ncols) / + real (nlines))) + nysample = real (max_nregions) / real (nxsample) + } + nxcols = ncols / nxsample + nylines = nlines / nysample + jstart = 1 + (nlines - nysample * nylines) / 2 + jend = jstart + (nysample - 1) * nylines + do j = jstart, jend, nylines { + istart = 1 + (ncols - nxsample * nxcols) / 2 + iend = istart + (nxsample - 1) * nxcols + do i = istart, iend, nxcols { + Memi[rg_xstatp(xc,RC1)+nregions] = i + Memi[rg_xstatp(xc,RC2)+nregions] = i + nxcols - 1 + Memi[rg_xstatp(xc,RL1)+nregions] = j + Memi[rg_xstatp(xc,RL2)+nregions] = j + nylines - 1 + Memr[rg_xstatp(xc,RZERO)+nregions] = INDEFR + Memr[rg_xstatp(xc,RXSLOPE)+nregions] = INDEFR + Memr[rg_xstatp(xc,RYSLOPE)+nregions] = INDEFR + Memr[rg_xstatp(xc,XSHIFTS)+nregions] = INDEFR + Memr[rg_xstatp(xc,YSHIFTS)+nregions] = INDEFR + nregions = nregions + 1 + } + } + + call rg_xseti (xc, NREGIONS, nregions) + if (nregions > 0) + call rg_xrealloc (xc, nregions) + else + call rg_xrfree (xc) + call sfree (sp) + + return (nregions) +end + + +# RG_XRREGIONS -- Read and decode the regions from a file. + +int procedure rg_xrregions (fd, im, xc, rp, max_nregions) + +int fd #I regions file descriptor +pointer im #I pointer to the reference image +pointer xc #I pointer to the cross-correlation structure +int rp #I index of the current region +int max_nregions #I the maximum number of regions + +int ncols, nlines, nregions, x1, y1, x2, y2, step +pointer sp, line, section +int rg_xstati(), getline(), rg_xgsections() +pointer rg_xstatp() + +begin + # Allocate working space. + call smark (sp) + call salloc (line, SZ_LINE, TY_CHAR) + call salloc (section, SZ_LINE, TY_CHAR) + + # Allocate the arrays to hold the regions information, + call rg_xrealloc (xc, max_nregions) + + # Initialize. + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + nregions = min (rp - 1, rg_xstati (xc, NREGIONS)) + + # Decode the regions string. + while ((getline (fd, Memc[line]) != EOF) && nregions < max_nregions) { + call sscan (Memc[line]) + call gargwrd (Memc[section], SZ_LINE) + while ((Memc[section] != EOS) && (nregions < max_nregions)) { + if (rg_xgsections (Memc[section], x1, x2, step, y1, y2, step, + ncols, nlines) == OK) { + Memi[rg_xstatp(xc,RC1)+nregions] = x1 + Memi[rg_xstatp(xc,RC2)+nregions] = x2 + Memi[rg_xstatp(xc,RL1)+nregions] = y1 + Memi[rg_xstatp(xc,RL2)+nregions] = y2 + Memr[rg_xstatp(xc,RZERO)+nregions] = INDEFR + Memr[rg_xstatp(xc,RXSLOPE)+nregions] = INDEFR + Memr[rg_xstatp(xc,RYSLOPE)+nregions] = INDEFR + Memr[rg_xstatp(xc,XSHIFTS)+nregions] = INDEFR + Memr[rg_xstatp(xc,YSHIFTS)+nregions] = INDEFR + nregions = nregions + 1 + } + call gargwrd (Memc[section], SZ_LINE) + } + } + + # Reallocate the correct amount of space. + call rg_xseti (xc, NREGIONS, nregions) + if (nregions > 0) + call rg_xrealloc (xc, nregions) + else + call rg_xrfree (xc) + + call sfree (sp) + + return (nregions) +end + + +# RG_XGREGIONS -- Decode a list of regions from a string containing +# a list of sections. + +int procedure rg_xgregions (im, xc, rp, max_nregions) + +pointer im #I pointer to the reference image +pointer xc #I pointer to cross-correlation structure +int rp #I the index of the current region +int max_nregions #I the maximum number of regions + +int ncols, nlines, nregions, x1, x2, y1, y2, step +pointer sp, section, region +int rg_xstati(), rg_xgsections() +pointer rg_xstatp() + +begin + # Allocate working space. + call smark (sp) + call salloc (region, SZ_LINE, TY_CHAR) + call salloc (section, SZ_LINE, TY_CHAR) + + # Allocate the arrays to hold the regions information. + call rg_xrealloc (xc, max_nregions) + + # Initialize. + call rg_xstats (xc, REGIONS, Memc[region], SZ_LINE) + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + nregions = min (rp - 1, rg_xstati (xc, NREGIONS)) + + # Decode the sections + call sscan (Memc[region]) + call gargwrd (Memc[section], SZ_LINE) + while ((Memc[section] != EOS) && (nregions < max_nregions)) { + if (rg_xgsections (Memc[section], x1, x2, step, y1, y2, step, + ncols, nlines) == OK) { + Memi[rg_xstatp(xc,RC1)+nregions] = x1 + Memi[rg_xstatp(xc,RC2)+nregions] = x2 + Memi[rg_xstatp(xc,RL1)+nregions] = y1 + Memi[rg_xstatp(xc,RL2)+nregions] = y2 + Memr[rg_xstatp(xc,RZERO)+nregions] = INDEFR + Memr[rg_xstatp(xc,RXSLOPE)+nregions] = INDEFR + Memr[rg_xstatp(xc,RYSLOPE)+nregions] = INDEFR + Memr[rg_xstatp(xc,XSHIFTS)+nregions] = INDEFR + Memr[rg_xstatp(xc,YSHIFTS)+nregions] = INDEFR + nregions = nregions + 1 + } + call gargwrd (Memc[section], SZ_LINE) + } + + + # Reallocate the correct amount of space. + call rg_xseti (xc, NREGIONS, nregions) + if (nregions > 0) + call rg_xrealloc (xc, nregions) + else + call rg_xrfree (xc) + + call sfree (sp) + + return (nregions) +end + + +# RG_XGSECTIONS -- Decode an image section into column and line limits +# and a step size. Sections which describe the whole image are decoded into +# a block ncols * nlines long. + +int procedure rg_xgsections (section, x1, x2, xstep, y1, y2, ystep, ncols, + nlines) + +char section[ARB] #I the input section string +int x1, x2 #O the output column section limits +int xstep #O the output column step size +int y1, y2 #O the output line section limits +int ystep #O the output line step size +int ncols, nlines #I the maximum number of lines and columns + +int ip +int rg_xgdim() + +begin + ip = 1 + if (rg_xgdim (section, ip, x1, x2, xstep, ncols) == ERR) + return (ERR) + if (rg_xgdim (section, ip, y1, y2, ystep, nlines) == ERR) + return (ERR) + + return (OK) +end + + +# RG_XGDIM -- Decode a single subscript expression to produce the +# range of values for that subscript (X1:X2), and the sampling step size, STEP. +# Note that X1 may be less than, greater than, or equal to X2, and STEP may +# be a positive or negative nonzero integer. Various shorthand notations are +# permitted, as is embedded whitespace. + +int procedure rg_xgdim (section, ip, x1, x2, step, limit) + +char section[ARB] #I the input image section +int ip #I/O pointer to the position in section string +int x1 #O first limit of dimension +int x2 #O second limit of dimension +int step #O step size of dimension +int limit #I maximum size of dimension + +int temp +int ctoi() + +begin + x1 = 1 + x2 = limit + step = 1 + + while (IS_WHITE(section[ip])) + ip = ip + 1 + + if (section[ip] =='[') + ip = ip + 1 + + while (IS_WHITE(section[ip])) + ip = ip + 1 + + # Get X1, X2. + if (ctoi (section, ip, temp) > 0) { # [x1 + x1 = max (1, min (temp, limit)) + if (section[ip] == ':') { + ip = ip + 1 + if (ctoi (section, ip, temp) == 0) # [x1:x2 + return (ERR) + x2 = max (1, min (temp, limit)) + } else + x2 = x1 + + } else if (section[ip] == '-') { + x1 = limit + x2 = 1 + ip = ip + 1 + if (section[ip] == '*') + ip = ip + 1 + + } else if (section[ip] == '*') # [* + ip = ip + 1 + + while (IS_WHITE(section[ip])) + ip = ip + 1 + + # Get sample step size, if give. + if (section[ip] == ':') { # ..:step + ip = ip + 1 + if (ctoi (section, ip, step) == 0) + return (ERR) + else if (step == 0) + return (ERR) + } + + # Allow notation such as "-*:5", (or even "-:5") where the step + # is obviously supposed to be negative. + + if (x1 > x2 && step > 0) + step = -step + + while (IS_WHITE(section[ip])) + ip = ip + 1 + + if (section[ip] == ',') { + ip = ip + 1 + return (OK) + } else if (section[ip] == ']') + return (OK) + else + return (ERR) +end diff --git a/pkg/images/immatch/src/xregister/rgxshow.x b/pkg/images/immatch/src/xregister/rgxshow.x new file mode 100644 index 00000000..3a746d9c --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxshow.x @@ -0,0 +1,172 @@ +include "xregister.h" + +# RG_XSHOW -- Show the XREGISTER parameters. + +procedure rg_xshow (xc) + +pointer xc #I pointer to the main xregister structure + +begin + call rg_xnshow (xc) + call printf ("\n") + call rg_xbshow (xc) + call printf ("\n") + call rg_xxshow (xc) + call printf ("\n") + call rg_xpshow (xc) +end + + +# RG_XNSHOW -- Show the input/output data XREGISTER parameters. + +procedure rg_xnshow (xc) + +pointer xc #I pointer to the main xregister structure + +pointer sp, str +int rg_xstati() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Set the object characteristics. + call printf ("\nInput/output data\n") + call rg_xstats (xc, IMAGE, Memc[str], SZ_FNAME) + call printf (" %s: %s\n") + call pargstr (KY_IMAGE) + call pargstr (Memc[str]) + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_FNAME) + call printf (" %s: %s\n") + call pargstr (KY_REFIMAGE) + call pargstr (Memc[str]) + call rg_xstats (xc, REGIONS, Memc[str], SZ_FNAME) + call printf (" %s: %s\n") + call pargstr (KY_REGIONS) + call pargstr (Memc[str]) + call printf (" %s = %d %s = %d\n") + call pargstr (KY_XLAG) + call pargi (rg_xstati (xc, XLAG)) + call pargstr (KY_YLAG) + call pargi (rg_xstati (xc, YLAG)) + call printf (" %s = %d %s = %d\n") + call pargstr (KY_DXLAG) + call pargi (rg_xstati (xc, DXLAG)) + call pargstr (KY_DYLAG) + call pargi (rg_xstati (xc, DYLAG)) + call rg_xstats (xc, DATABASE, Memc[str], SZ_FNAME) + call printf (" %s: %s\n") + call pargstr (KY_DATABASE) + call pargstr (Memc[str]) + call rg_xstats (xc, RECORD, Memc[str], SZ_FNAME) + call printf (" %s: %s\n") + call pargstr (KY_RECORD) + call pargstr (Memc[str]) + call rg_xstats (xc, REFFILE, Memc[str], SZ_FNAME) + call printf (" %s: %s\n") + call pargstr (KY_REFFILE) + call pargstr (Memc[str]) + call rg_xstats (xc, OUTIMAGE, Memc[str], SZ_FNAME) + call printf (" %s: %s\n") + call pargstr (KY_OUTIMAGE) + call pargstr (Memc[str]) + + call sfree (sp) +end + + +# RG_XBSHOW -- Show the background fitting parameters. + +procedure rg_xbshow (xc) + +pointer xc #I pointer to the main xregister structure + +int back +pointer sp, str +int rg_xstati() +real rg_xstatr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + back = rg_xstati (xc, BACKGRD) + call printf ("Background fitting parameters:\n") + call rg_xstats (xc, BSTRING, Memc[str], SZ_LINE) + call printf (" %s: %s\n") + call pargstr (KY_BACKGROUND) + call pargstr (Memc[str]) + call printf (" %s = %d\n") + call pargstr (KY_BORDER) + call pargi (rg_xstati (xc, BORDER)) + call printf (" %s = %g %s = %g\n") + call pargstr (KY_LOREJECT) + call pargr (rg_xstatr (xc, LOREJECT)) + call pargstr (KY_HIREJECT) + call pargr (rg_xstatr (xc, HIREJECT)) + call printf (" %s = %g\n") + call pargstr (KY_APODIZE) + call pargr (rg_xstatr (xc, APODIZE)) + call rg_xstats (xc, FSTRING, Memc[str], SZ_LINE) + call printf (" %s: %s\n") + call pargstr (KY_FILTER) + call pargstr (Memc[str]) + + call sfree (sp) +end + + +# RG_XXSHOW -- Show the cross-correlation function parameters. + +procedure rg_xxshow (xc) + +pointer xc #I pointer to the main xregister structure + +pointer sp, str +int rg_xstati() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + call printf ("Cross correlation function:\n") + call rg_xstats (xc, CSTRING, Memc[str], SZ_LINE) + call printf (" %s: %s\n") + call pargstr (KY_CORRELATION) + call pargstr (Memc[str]) + call printf (" %s = %d %s = %d\n") + call pargstr (KY_XWINDOW) + call pargi (rg_xstati (xc, XWINDOW)) + call pargstr (KY_YWINDOW) + call pargi (rg_xstati (xc, YWINDOW)) + + call sfree (sp) +end + + +# RG_XPSHOW -- Show the peak centering parameters. + +procedure rg_xpshow (xc) + +pointer xc #I pointer to the main xregister structure + +pointer sp, str +int rg_xstati() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + call printf ("Peak centering parameters:\n") + call rg_xstats (xc, PSTRING, Memc[str], SZ_LINE) + call printf (" %s: %s\n") + call pargstr (KY_PEAKCENTER) + call pargstr (Memc[str]) + call printf (" %s = %d %s = %d\n") + call pargstr (KY_XCBOX) + call pargi (rg_xstati (xc, XCBOX)) + call pargstr (KY_YCBOX) + call pargi (rg_xstati (xc, YCBOX)) + + call sfree (sp) +end diff --git a/pkg/images/immatch/src/xregister/rgxtools.x b/pkg/images/immatch/src/xregister/rgxtools.x new file mode 100644 index 00000000..e1fb921e --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxtools.x @@ -0,0 +1,685 @@ +include "xregister.h" + +# RG_XINIT -- Initialize the cross-correlation code fitting structure. + +procedure rg_xinit (xc, cfunc) + +pointer xc #O pointer to the cross-correlation structure +int cfunc #I the input cross-correlation function + +begin + call malloc (xc, LEN_XCSTRUCT, TY_STRUCT) + + # Initialize the regions pointers. + XC_RC1(xc) = NULL + XC_RC2(xc) = NULL + XC_RL1(xc) = NULL + XC_RL2(xc) = NULL + XC_RZERO(xc) = NULL + XC_RXSLOPE(xc) = NULL + XC_RYSLOPE(xc) = NULL + XC_XSHIFTS(xc) = NULL + XC_YSHIFTS(xc) = NULL + XC_TXSHIFT(xc) = 0.0 + XC_TYSHIFT(xc) = 0.0 + XC_NREGIONS(xc) = 0 + XC_CREGION(xc) = 1 + + # Set up transformation parameters. + XC_NREFPTS(xc) = 0 + call malloc (XC_XREF(xc), MAX_NREF, TY_REAL) + call malloc (XC_YREF(xc), MAX_NREF, TY_REAL) + call malloc (XC_TRANSFORM(xc), MAX_NTRANSFORM, TY_REAL) + + # Initialize the region offsets + XC_IXLAG(xc) = DEF_IXLAG + XC_IYLAG(xc) = DEF_IYLAG + XC_XLAG(xc) = DEF_IXLAG + XC_YLAG(xc) = DEF_IYLAG + XC_DXLAG(xc) = DEF_DXLAG + XC_DYLAG(xc) = DEF_DYLAG + + # Define the background fitting parameters. + XC_BACKGRD(xc) = XC_BNONE + call strcpy ("none", XC_BSTRING(xc), SZ_FNAME) + XC_BVALUER(xc) = 0.0 + XC_BVALUE(xc) = 0.0 + XC_BORDER(xc) = DEF_BORDER + XC_LOREJECT(xc) = DEF_LOREJECT + XC_HIREJECT(xc) = DEF_HIREJECT + XC_APODIZE(xc) = 0.0 + XC_FILTER(xc) = XC_FNONE + call strcpy ("none", XC_FSTRING(xc), SZ_FNAME) + + # Get the correlation parameters. + XC_CFUNC(xc) = cfunc + switch (cfunc) { + case XC_DISCRETE: + call strcpy ("discrete", XC_CSTRING(xc), SZ_FNAME) + case XC_FOURIER: + call strcpy ("fourier", XC_CSTRING(xc), SZ_FNAME) + case XC_FILE: + call strcpy ("file", XC_CSTRING(xc), SZ_FNAME) + case XC_DIFFERENCE: + call strcpy ("difference", XC_CSTRING(xc), SZ_FNAME) + default: + call strcpy ("unknown", XC_CSTRING(xc), SZ_FNAME) + } + XC_XWINDOW(xc) = DEF_XWINDOW + XC_YWINDOW(xc) = DEF_YWINDOW + XC_XCOR(xc) = NULL + + # Define the peak fitting function. + XC_PFUNC(xc) = DEF_PFUNC + call sprintf (XC_PSTRING(xc), SZ_FNAME, "%s") + call pargstr ("centroid") + XC_XCBOX(xc) = DEF_XCBOX + XC_YCBOX(xc) = DEF_YCBOX + + # Initialize the strings. + XC_IMAGE(xc) = EOS + XC_REFIMAGE(xc) = EOS + XC_REGIONS(xc) = EOS + XC_DATABASE(xc) = EOS + XC_OUTIMAGE(xc) = EOS + XC_REFFILE(xc) = EOS + XC_RECORD(xc) = EOS + + # Initialize the buffers. + call rg_xrinit (xc) + +end + + +# RG_XRINIT -- Initialize the regions definition portion of the +# cross correlation code fitting structure. + +procedure rg_xrinit (xc) + +pointer xc #I pointer to crosscor structure + +begin + call rg_xrfree (xc) + + XC_NREGIONS(xc) = 0 + XC_CREGION(xc) = 1 + + call malloc (XC_RC1(xc), MAX_NREGIONS, TY_INT) + call malloc (XC_RC2(xc), MAX_NREGIONS, TY_INT) + call malloc (XC_RL1(xc), MAX_NREGIONS, TY_INT) + call malloc (XC_RL2(xc), MAX_NREGIONS, TY_INT) + call malloc (XC_RZERO(xc), MAX_NREGIONS, TY_REAL) + call malloc (XC_RXSLOPE(xc), MAX_NREGIONS, TY_REAL) + call malloc (XC_RYSLOPE(xc), MAX_NREGIONS, TY_REAL) + call malloc (XC_XSHIFTS(xc), MAX_NREGIONS, TY_REAL) + call malloc (XC_YSHIFTS(xc), MAX_NREGIONS, TY_REAL) + + call amovki (INDEFI, Memi[XC_RC1(xc)], MAX_NREGIONS) + call amovki (INDEFI, Memi[XC_RC2(xc)], MAX_NREGIONS) + call amovki (INDEFI, Memi[XC_RL1(xc)], MAX_NREGIONS) + call amovki (INDEFI, Memi[XC_RL2(xc)], MAX_NREGIONS) + call amovkr (INDEFR, Memr[XC_RZERO(xc)], MAX_NREGIONS) + call amovkr (INDEFR, Memr[XC_RXSLOPE(xc)], MAX_NREGIONS) + call amovkr (INDEFR, Memr[XC_RYSLOPE(xc)], MAX_NREGIONS) + call amovkr (INDEFR, Memr[XC_XSHIFTS(xc)], MAX_NREGIONS) + call amovkr (INDEFR, Memr[XC_YSHIFTS(xc)], MAX_NREGIONS) + + XC_TXSHIFT(xc) = 0.0 + XC_TYSHIFT(xc) = 0.0 +end + + +# RG_XCINDEFR -- Re-initialize the background and answers regions portion of +# the cross-correlation fitting structure + +procedure rg_xcindefr (xc, creg) + +pointer xc #I pointer to the cross-correlation structure +int creg #I the current region + +int nregions +int rg_xstati() + +begin + nregions = rg_xstati (xc, NREGIONS) + if (creg < 1 || creg > nregions) + return + + if (nregions > 0) { + Memr[XC_RZERO(xc)+creg-1] = INDEFR + Memr[XC_RXSLOPE(xc)+creg-1] = INDEFR + Memr[XC_RYSLOPE(xc)+creg-1] = INDEFR + Memr[XC_XSHIFTS(xc)+creg-1] = INDEFR + Memr[XC_YSHIFTS(xc)+creg-1] = INDEFR + } + + XC_TXSHIFT(xc) = 0.0 + XC_TYSHIFT(xc) = 0.0 +end + + +# RG_XINDEFR -- Re-initialize the background and answers regions portion of +# the cross-correlation fitting structure for all regions and reset the +# current region to 1. + +procedure rg_xindefr (xc) + +pointer xc #I pointer to the cross-correlation structure + +int nregions +int rg_xstati() + +begin + nregions = rg_xstati (xc, NREGIONS) + + if (nregions > 0) { + call amovkr (INDEFR, Memr[XC_RZERO(xc)], nregions) + call amovkr (INDEFR, Memr[XC_RXSLOPE(xc)], nregions) + call amovkr (INDEFR, Memr[XC_RYSLOPE(xc)], nregions) + call amovkr (INDEFR, Memr[XC_XSHIFTS(xc)], nregions) + call amovkr (INDEFR, Memr[XC_YSHIFTS(xc)], nregions) + } + + XC_CREGION(xc) = 1 + XC_TXSHIFT(xc) = 0.0 + XC_TYSHIFT(xc) = 0.0 +end + + +# RG_XREALLOC -- Reallocate the regions bufffers and initialize if necessary. + +procedure rg_xrealloc (xc, nregions) + +pointer xc #I pointer to crosscor structure +int nregions #I number of regions + +int nr +int rg_xstati() + +begin + nr = rg_xstati (xc, NREGIONS) + + call realloc (XC_RC1(xc), nregions, TY_INT) + call realloc (XC_RC2(xc), nregions, TY_INT) + call realloc (XC_RL1(xc), nregions, TY_INT) + call realloc (XC_RL2(xc), nregions, TY_INT) + call realloc (XC_RZERO(xc), nregions, TY_REAL) + call realloc (XC_RXSLOPE(xc), nregions, TY_REAL) + call realloc (XC_RYSLOPE(xc), nregions, TY_REAL) + call realloc (XC_XSHIFTS(xc), nregions, TY_REAL) + call realloc (XC_YSHIFTS(xc), nregions, TY_REAL) + + call amovki (INDEFI, Memi[XC_RC1(xc)+nr], nregions - nr) + call amovki (INDEFI, Memi[XC_RC2(xc)+nr], nregions - nr) + call amovki (INDEFI, Memi[XC_RL1(xc)+nr], nregions - nr) + call amovki (INDEFI, Memi[XC_RL2(xc)+nr], nregions - nr) + call amovkr (INDEFR, Memr[XC_RZERO(xc)+nr], nregions - nr) + call amovkr (INDEFR, Memr[XC_RXSLOPE(xc)+nr], nregions - nr) + call amovkr (INDEFR, Memr[XC_RYSLOPE(xc)+nr], nregions - nr) + call amovkr (INDEFR, Memr[XC_XSHIFTS(xc)+nr], nregions - nr) + call amovkr (INDEFR, Memr[XC_YSHIFTS(xc)+nr], nregions - nr) +end + + +# RG_XFREE -- Free the cross-correlation fitting structure. + +procedure rg_xfree (xc) + +pointer xc #I pointer to the cross-correlation structure + +begin + # Free the region descriptors. + call rg_xrfree (xc) + + # Free the transformation descriptors. + if (XC_XREF(xc) != NULL) + call mfree (XC_XREF(xc), TY_REAL) + if (XC_YREF(xc) != NULL) + call mfree (XC_YREF(xc), TY_REAL) + if (XC_TRANSFORM(xc) != NULL) + call mfree (XC_TRANSFORM(xc), TY_REAL) + + # Free the correlation function. + if (XC_XCOR(xc) != NULL) + call mfree (XC_XCOR(xc), TY_REAL) + + call mfree (xc, TY_STRUCT) +end + + +# RG_XRFREE -- Free the regions portion of the cross-correlation structure. + +procedure rg_xrfree (xc) + +pointer xc #I pointer to the cross-correlation structure + +begin + call rg_xseti (xc, NREGIONS, 0) + if (XC_RC1(xc) != NULL) + call mfree (XC_RC1(xc), TY_INT) + XC_RC1(xc) = NULL + if (XC_RC2(xc) != NULL) + call mfree (XC_RC2(xc), TY_INT) + XC_RC2(xc) = NULL + if (XC_RL1(xc) != NULL) + call mfree (XC_RL1(xc), TY_INT) + XC_RL1(xc) = NULL + if (XC_RL2(xc) != NULL) + call mfree (XC_RL2(xc), TY_INT) + XC_RL2(xc) = NULL + if (XC_RZERO(xc) != NULL) + call mfree (XC_RZERO(xc), TY_REAL) + XC_RZERO(xc) = NULL + if (XC_RXSLOPE(xc) != NULL) + call mfree (XC_RXSLOPE(xc), TY_REAL) + XC_RXSLOPE(xc) = NULL + if (XC_RYSLOPE(xc) != NULL) + call mfree (XC_RYSLOPE(xc), TY_REAL) + XC_RYSLOPE(xc) = NULL + if (XC_XSHIFTS(xc) != NULL) + call mfree (XC_XSHIFTS(xc), TY_REAL) + XC_XSHIFTS(xc) = NULL + if (XC_YSHIFTS(xc) != NULL) + call mfree (XC_YSHIFTS(xc), TY_REAL) + XC_YSHIFTS(xc) = NULL +end + + +# RG_XSTATI -- Fetch the value of a cross-correlation fitting structure +# integer parameter. + +int procedure rg_xstati (xc, param) + +pointer xc #I pointer to the cross-correlation fitting structure +int param #I parameter to be fetched + +begin + switch (param) { + case CFUNC: + return (XC_CFUNC(xc)) + case IXLAG: + return (XC_IXLAG(xc)) + case IYLAG: + return (XC_IYLAG(xc)) + case XLAG: + return (XC_XLAG(xc)) + case YLAG: + return (XC_YLAG(xc)) + case DXLAG: + return (XC_DXLAG(xc)) + case DYLAG: + return (XC_DYLAG(xc)) + case XWINDOW: + return (XC_XWINDOW(xc)) + case YWINDOW: + return (XC_YWINDOW(xc)) + case CREGION: + return (XC_CREGION(xc)) + case NREGIONS: + return (XC_NREGIONS(xc)) + case BACKGRD: + return (XC_BACKGRD(xc)) + case BORDER: + return (XC_BORDER(xc)) + case FILTER: + return (XC_FILTER(xc)) + case XCBOX: + return (XC_XCBOX(xc)) + case YCBOX: + return (XC_YCBOX(xc)) + case PFUNC: + return (XC_PFUNC(xc)) + case NREFPTS: + return (XC_NREFPTS(xc)) + default: + call error (0, "RG_XSTATI: Undefined integer parameter.") + } +end + + +# RG_XSTATP -- Fetch the value of a pointer parameter. + +pointer procedure rg_xstatp (xc, param) + +pointer xc #I pointer to the cross-correlation structure +int param #I parameter to be fetched + +begin + switch (param) { + case RC1: + return (XC_RC1(xc)) + case RC2: + return (XC_RC2(xc)) + case RL1: + return (XC_RL1(xc)) + case RL2: + return (XC_RL2(xc)) + case RZERO: + return (XC_RZERO(xc)) + case RXSLOPE: + return (XC_RXSLOPE(xc)) + case RYSLOPE: + return (XC_RYSLOPE(xc)) + case XSHIFTS: + return (XC_XSHIFTS(xc)) + case YSHIFTS: + return (XC_YSHIFTS(xc)) + case XCOR: + return (XC_XCOR(xc)) + case XREF: + return (XC_XREF(xc)) + case YREF: + return (XC_YREF(xc)) +# case CORAPODIZE: +# return (XC_CORAPODIZE(xc)) + case TRANSFORM: + return (XC_TRANSFORM(xc)) + default: + call error (0, "RG_XSTATP: Undefined pointer parameter.") + } +end + + +# RG_XSTATR -- Fetch the value of a real parameter. + +real procedure rg_xstatr (xc, param) + +pointer xc #I pointer to the cross-correlation structure +int param #I parameter to be fetched + +begin + switch (param) { + case BVALUER: + return (XC_BVALUER(xc)) + case BVALUE: + return (XC_BVALUE(xc)) + case LOREJECT: + return (XC_LOREJECT(xc)) + case HIREJECT: + return (XC_HIREJECT(xc)) + case APODIZE: + return (XC_APODIZE(xc)) + case TXSHIFT: + return (XC_TXSHIFT(xc)) + case TYSHIFT: + return (XC_TYSHIFT(xc)) + default: + call error (0, "RG_XSTATR: Undefined real parameter.") + } +end + + +# RG_XSTATS -- Fetch the value of a string parameter. + +procedure rg_xstats (xc, param, str, maxch) + +pointer xc #I pointer to the cross-correlation structure +int param #I parameter to be fetched +char str[ARB] #O output value of string parameter +int maxch #I maximum number of characters in output string + +begin + switch (param) { + case BSTRING: + call strcpy (XC_BSTRING(xc), str, maxch) + case FSTRING: + call strcpy (XC_FSTRING(xc), str, maxch) + case CSTRING: + call strcpy (XC_CSTRING(xc), str, maxch) + case PSTRING: + call strcpy (XC_PSTRING(xc), str, maxch) + case REFIMAGE: + call strcpy (XC_REFIMAGE(xc), str, maxch) + case IMAGE: + call strcpy (XC_IMAGE(xc), str, maxch) + case OUTIMAGE: + call strcpy (XC_OUTIMAGE(xc), str, maxch) + case REGIONS: + call strcpy (XC_REGIONS(xc), str, maxch) + case DATABASE: + call strcpy (XC_DATABASE(xc), str, maxch) + case RECORD: + call strcpy (XC_RECORD(xc), str, maxch) + case REFFILE: + call strcpy (XC_REFFILE(xc), str, maxch) + default: + call error (0, "RG_XSTATS: Undefined string parameter.") + } +end + + +# RG_XSETI -- Set the value of an integer parameter. + +procedure rg_xseti (xc, param, value) + +pointer xc #I pointer to the cross-correlation structure +int param #I parameter to be set +int value #O value of the integer parameter + +begin + switch (param) { + case CFUNC: + XC_CFUNC(xc) = value + switch (value) { + case XC_DISCRETE: + call strcpy ("discrete", XC_CSTRING(xc), SZ_FNAME) + case XC_FOURIER: + call strcpy ("fourier", XC_CSTRING(xc), SZ_FNAME) + case XC_FILE: + call strcpy ("file", XC_CSTRING(xc), SZ_FNAME) + case XC_DIFFERENCE: + call strcpy ("difference", XC_CSTRING(xc), SZ_FNAME) + default: + call strcpy ("unknown", XC_CSTRING(xc), SZ_FNAME) + } + case IXLAG: + XC_IXLAG(xc) = value + case IYLAG: + XC_IYLAG(xc) = value + case XLAG: + XC_XLAG(xc) = value + case YLAG: + XC_YLAG(xc) = value + case DXLAG: + XC_DXLAG(xc) = value + case DYLAG: + XC_DYLAG(xc) = value + case XWINDOW: + XC_XWINDOW(xc) = value + case YWINDOW: + XC_YWINDOW(xc) = value + case BACKGRD: + XC_BACKGRD(xc) = value + switch (value) { + case XC_BNONE: + call strcpy ("none", XC_BSTRING(xc), SZ_FNAME) + case XC_MEAN: + call strcpy ("mean", XC_BSTRING(xc), SZ_FNAME) + case XC_MEDIAN: + call strcpy ("median", XC_BSTRING(xc), SZ_FNAME) + case XC_SLOPE: + call strcpy ("plane", XC_BSTRING(xc), SZ_FNAME) + default: + call strcpy ("none", XC_BSTRING(xc), SZ_FNAME) + } + case BORDER: + XC_BORDER(xc) = value + case FILTER: + XC_FILTER(xc) = value + switch (value) { + case XC_FNONE: + call strcpy ("none", XC_FSTRING(xc), SZ_FNAME) + case XC_LAPLACE: + call strcpy ("laplace", XC_FSTRING(xc), SZ_FNAME) + default: + call strcpy ("none", XC_FSTRING(xc), SZ_FNAME) + } + case XCBOX: + XC_XCBOX(xc) = value + case YCBOX: + XC_YCBOX(xc) = value + case PFUNC: + XC_PFUNC(xc) = value + switch (value) { + case XC_PNONE: + call strcpy ("none", XC_PSTRING(xc), SZ_FNAME) + case XC_CENTROID: + call strcpy ("centroid", XC_PSTRING(xc), SZ_FNAME) + case XC_PARABOLA: + call strcpy ("parabolic", XC_PSTRING(xc), SZ_FNAME) + case XC_SAWTOOTH: + call strcpy ("sawtooth", XC_PSTRING(xc), SZ_FNAME) +# case XC_MARK: +# call strcpy ("mark", XC_PSTRING(xc), SZ_FNAME) + default: + ; + } + case NREFPTS: + XC_NREFPTS(xc) = value + case CREGION: + XC_CREGION(xc) = value + case NREGIONS: + XC_NREGIONS(xc) = value + default: + call error (0, "RG_XSETI: Undefined integer parameter.") + } +end + + +# RG_XSETP -- Set the value of a pointer parameter. + +procedure rg_xsetp (xc, param, value) + +pointer xc #I pointer to the cross-correlation structure +int param #I parameter to be set +pointer value #O value of the pointer parameter + +begin + switch (param) { + case RC1: + XC_RC1(xc) = value + case RC2: + XC_RC2(xc) = value + case RL1: + XC_RL1(xc) = value + case RL2: + XC_RL2(xc) = value + case RZERO: + XC_RZERO(xc) = value + case RXSLOPE: + XC_RXSLOPE(xc) = value + case RYSLOPE: + XC_RYSLOPE(xc) = value + case XSHIFTS: + XC_XSHIFTS(xc) = value + case YSHIFTS: + XC_YSHIFTS(xc) = value + case XCOR: + XC_XCOR(xc) = value + case XREF: + XC_XREF(xc) = value + case YREF: + XC_YREF(xc) = value + case TRANSFORM: + XC_TRANSFORM(xc) = value +# case CORAPODIZE: +# XC_CORAPODIZE(xc) = value + default: + call error (0, "RG_XSETP: Undefined pointer parameter.") + } +end + + +# RG_XSETR -- Set the value of a real parameter. + +procedure rg_xsetr (xc, param, value) + +pointer xc #I pointer to the cross-correlation structure +int param #I parameter to be set +real value #O value of real parameter + +begin + switch (param) { + case BVALUER: + XC_BVALUER(xc) = value + case BVALUE: + XC_BVALUE(xc) = value + case LOREJECT: + XC_LOREJECT(xc) = value + case HIREJECT: + XC_HIREJECT(xc) = value + case APODIZE: + XC_APODIZE(xc) = value + case TXSHIFT: + XC_TXSHIFT(xc) = value + case TYSHIFT: + XC_TYSHIFT(xc) = value + default: + call error (0, "RG_XSETR: Undefined real parameter.") + } +end + + +# RG_XSETS -- Set the value of a string parameter. + +procedure rg_xsets (xc, param, str) + +pointer xc #I pointer to the cross-correlation structure +int param #I parameter to be set +char str[ARB] #O value of string parameter + +int index +pointer sp, temp +int strdic(), fnldir() + +begin + call smark (sp) + call salloc (temp, SZ_FNAME, TY_CHAR) + + switch (param) { + case BSTRING: + index = strdic (str, str, SZ_LINE, XC_BTYPES) + if (index > 0) { + call strcpy (str, XC_BSTRING(xc), SZ_FNAME) + call rg_xseti (xc, BACKGRD, index) + } + case FSTRING: + index = strdic (str, str, SZ_LINE, XC_FTYPES) + if (index > 0) { + call strcpy (str, XC_FSTRING(xc), SZ_FNAME) + call rg_xseti (xc, FILTER, index) + } + case CSTRING: + index = strdic (str, str, SZ_LINE, XC_CTYPES) + if (index > 0) { + call strcpy (str, XC_CSTRING(xc), SZ_FNAME) + call rg_xseti (xc, CFUNC, index) + } + case PSTRING: + call strcpy (str, XC_PSTRING(xc), SZ_FNAME) + case REFIMAGE: + call imgcluster (str, Memc[temp], SZ_FNAME) + index = fnldir (Memc[temp], XC_REFIMAGE(xc), SZ_FNAME) + call strcpy (Memc[temp+index], XC_REFIMAGE(xc), SZ_FNAME) + case IMAGE: + call imgcluster (str, Memc[temp], SZ_FNAME) + index = fnldir (Memc[temp], XC_IMAGE(xc), SZ_FNAME) + call strcpy (Memc[temp+index], XC_IMAGE(xc), SZ_FNAME) + case OUTIMAGE: + call strcpy (str, XC_OUTIMAGE(xc), SZ_FNAME) + case REGIONS: + call strcpy (str, XC_REGIONS(xc), SZ_FNAME) + case DATABASE: + index = fnldir (str, XC_DATABASE(xc), SZ_FNAME) + call strcpy (str[index+1], XC_DATABASE(xc), SZ_FNAME) + case RECORD: + call strcpy (str, XC_RECORD(xc), SZ_FNAME) + case REFFILE: + index = fnldir (str, XC_REFFILE(xc), SZ_FNAME) + call strcpy (str[index+1], XC_REFFILE(xc), SZ_FNAME) + default: + call error (0, "RG_XSETS: Undefined string parameter.") + } + + call sfree (sp) +end diff --git a/pkg/images/immatch/src/xregister/rgxtransform.x b/pkg/images/immatch/src/xregister/rgxtransform.x new file mode 100644 index 00000000..63ee5f24 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxtransform.x @@ -0,0 +1,446 @@ +include +include +include "xregister.h" + +# RG_GXTRANSFORM -- Open the reference points file and the read the +# coordinates of the reference points in the reference image. Return +# the reference points file name and descriptor. + +int procedure rg_gxtransform (list, xc, reffile) + +int list #I list of reference points files +pointer xc #I pointer to the cross-correlation structure +char reffile[ARB] #O the output reference points file name + +int tdf +pointer sp, line, pxref, pyref +real x1, y1, x2, y2, x3, y3 +int fntgfnb(), open(), getline(), nscan() +pointer rg_xstatp() + +begin + # Get some working memory. + call smark (sp) + call salloc (line, SZ_LINE, TY_CHAR) + + # Get the points to the reference point lists. + pxref = rg_xstatp (xc, XREF) + pyref = rg_xstatp (xc, YREF) + call aclrr (Memr[rg_xstatp(xc, XREF)], MAX_NREF) + call aclrr (Memr[rg_xstatp(xc, YREF)], MAX_NREF) + + # Open the reference points file and read the coordinates. + while (fntgfnb (list, reffile, SZ_FNAME) != EOF) { + + iferr { + + # Open the reference file. + tdf = open (reffile, READ_ONLY, TEXT_FILE) + call aclrr (Memr[pxref], MAX_NREF) + call aclrr (Memr[pyref], MAX_NREF) + + # Read up to three valid reference points from the list. + while (getline (tdf, Memc[line]) != EOF) { + call sscan (Memc[line]) + call gargr (x1) + call gargr (y1) + call gargr (x2) + call gargr (y2) + call gargr (x3) + call gargr (y3) + if (nscan () >= 2) + break + } + + # Store the reference points. + if (nscan () == 2) { + Memr[pxref] = x1 + Memr[pyref] = y1 + call rg_xseti (xc, NREFPTS, 1) + } else if (nscan () == 4) { + Memr[pxref] = x1 + Memr[pyref] = y1 + Memr[pxref+1] = x2 + Memr[pyref+1] = y2 + call rg_xseti (xc, NREFPTS, 2) + } else if (nscan () == 6) { + Memr[pxref] = x1 + Memr[pyref] = y1 + Memr[pxref+1] = x2 + Memr[pyref+1] = y2 + Memr[pxref+2] = x3 + Memr[pyref+2] = y3 + call rg_xseti (xc, NREFPTS, 2) + } else + call rg_xseti (xc, NREFPTS, 0) + + } then { + call rg_xseti (xc, NREFPTS, 0) + } + } + + call sfree (sp) + + return (tdf) +end + + +# RG_ITRANSFORM -- Compute the transformation from the input image to the +# reference image interactively. + +procedure rg_itransform (xc, imr, im, id) + +pointer xc #I pointer to the cross-correlation stucture +pointer imr #I pointer to the reference image +pointer im #I pointer to the input image +pointer id #I pointer to the display device + +int nref, nstar, wcs, key +pointer sp, cmd, x, y, pxref, pyref, ptrans +real wx, wy +int clgcur() +pointer rg_xstatp() + +begin + # Allocate working space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (x, MAX_NREF, TY_REAL) + call salloc (y, MAX_NREF, TY_REAL) + call aclrr (Memr[x], MAX_NREF) + call aclrr (Memr[y], MAX_NREF) + + # Get the pointers. + pxref = rg_xstatp (xc, XREF) + pyref = rg_xstatp (xc, YREF) + ptrans = rg_xstatp (xc, TRANSFORM) + + # Mark up to three reference stars. + nref = 0 + call printf ("Mark reference star %d with the image cursor [q=quit]: ") + call pargi (nref + 1) + while ((nref < MAX_NREF) && clgcur ("icommands", wx, wy, wcs, key, + Memc[cmd], SZ_LINE) != EOF) { + if (key == 'q') { + call printf ("\n") + break + } + if (wx < 0.5 || wx > IM_LEN(imr,1) + 0.5) { + call printf ("\n") + next + } + if (wy < 0.5 || wy > IM_LEN(imr,2) + 0.5) { + call printf ("\n") + next + } + call printf ("%g %g\n") + call pargr (wx) + call pargr (wy) + Memr[pxref+nref] = wx + Memr[pyref+nref] = wy + nref = nref + 1 + call rg_xseti (xc, NREFPTS, nref) + if (nref >= MAX_NREF) + break + call printf ( + "Mark reference star %d with the image cursor [q=quit]: ") + call pargi (nref + 1) + } + + # Mark the corresponding input image stars. + if (nref > 0) { + + nstar = 0 + call printf ("Mark image star %d with the image cursor [q=quit]: ") + call pargi (nstar + 1) + while ((nstar < nref) && clgcur ("icommands", wx, wy, wcs, key, + Memc[cmd], SZ_LINE) != EOF) { + if (key == 'q') { + call printf ("\n") + break + } + if (wx < 0.5 || wx > IM_LEN(im,1) + 0.5) { + call printf ("\n") + next + } + if (wy < 0.5 || wy > IM_LEN(im,2) + 0.5) { + call printf ("\n") + next + } + call printf ("%g %g\n") + call pargr (wx) + call pargr (wy) + Memr[x+nstar] = wx + Memr[y+nstar] = wy + nstar = nstar + 1 + if (nstar >= MAX_NREF) + break + call printf ( + "Mark image star %d with the image cursor [q=quit]: ") + call pargi (nstar + 1) + } + + # Compute the transformation. + if (nstar > 0) { + switch (nstar) { + case 0: + call rg_xshift (Memr[pxref], Memr[pyref], Memr[pxref], + Memr[pyref], Memr[ptrans]) + case 1: + call rg_xshift (Memr[x], Memr[y], Memr[pxref], Memr[pyref], + Memr[ptrans]) + #case 2: + #call rg_xtwostar (Memr[x], Memr[y], Memr[pxref], + #Memr[pyref], Memr[ptrans]) + #case 3: + #call rg_xthreestar (Memr[x], Memr[y], Memr[pxref], + #Memr[pyref], Memr[ptrans]) + + default: + call rg_xshift (Memr[pxref], Memr[pyref], Memr[pxref], + Memr[pyref], Memr[ptrans]) + } + } + } + + call sfree (sp) +end + + +# RG_XTRANSFORM -- Compute the transformation from the input image to +# the reference image + +procedure rg_xtransform (tfd, xc) + +int tfd #I the reference points file descriptor +pointer xc #I the cross-correlation file descriptor + +int nref +pointer sp, line, x, y, pxref, pyref, ptrans +int getline(), rg_xstati(), nscan() +pointer rg_xstatp() + +begin + # Allocate working space. + call smark (sp) + call salloc (line, SZ_LINE, TY_CHAR) + call salloc (x, MAX_NREF, TY_REAL) + call salloc (y, MAX_NREF, TY_REAL) + call aclrr (Memr[x], MAX_NREF) + call aclrr (Memr[y], MAX_NREF) + + # Get the pointers to the reference image data. + nref = rg_xstati (xc, NREFPTS) + pxref = rg_xstatp (xc, XREF) + pyref = rg_xstatp (xc, YREF) + ptrans = rg_xstatp (xc, TRANSFORM) + + # Read the input image reference points. + while ((nref > 0) && getline (tfd, Memc[line]) != EOF) { + call sscan (Memc[line]) + call gargr (Memr[x]) + call gargr (Memr[y]) + call gargr (Memr[x+1]) + call gargr (Memr[y+1]) + call gargr (Memr[x+2]) + call gargr (Memr[y+2]) + if (nscan() >= 2 * nref) + break + } + + # Compute the transform. + if (nscan () < 2 * nref) { + call rg_xshift (Memr[pxref], Memr[pyref], Memr[pxref], Memr[pyref], + Memr[ptrans]) + } else { + switch (nref) { + case 0: + call rg_xshift (Memr[pxref], Memr[pyref], Memr[pxref], + Memr[pyref], Memr[ptrans]) + case 1: + call rg_xshift (Memr[x], Memr[y], Memr[pxref], Memr[pyref], + Memr[ptrans]) + case 2: + call rg_xtwostar (Memr[x], Memr[y], Memr[pxref], Memr[pyref], + Memr[ptrans]) + case 3: + call rg_xthreestar (Memr[x], Memr[y], Memr[pxref], Memr[pyref], + Memr[ptrans]) + } + } + + call sfree (sp) +end + + +# RG_ETRANSFORM -- Evaulate the current transform at a single point. + +procedure rg_etransform (xc, xin, yin, xout, yout) + +pointer xc #I pointer to the cross-correlation structure +real xin, yin #I the input x and y values +real xout, yout #O the output x and y values + +pointer ptrans +pointer rg_xstatp + +begin + ptrans = rg_xstatp (xc, TRANSFORM) + xout = Memr[ptrans] * xin + Memr[ptrans+1] * yin + Memr[ptrans+2] + yout = Memr[ptrans+3] * xin + Memr[ptrans+4] * yin + Memr[ptrans+5] +end + + +# RG_XSHIFT -- Compute the transformation coefficients required to define a +# simple shift using a single data point. + +procedure rg_xshift (xref, yref, xlist, ylist, coeff) + +real xref[ARB] #I x reference coordinates +real yref[ARB] #I y reference coordinates +real xlist[ARB] #I x input coordinates +real ylist[ARB] #I y input coordinates +real coeff[ARB] #O output coefficient array + +begin + # Compute the x transformation. + coeff[1] = 1.0 + coeff[2] = 0.0 + coeff[3] = xref[1] - xlist[1] + + # Compute the y transformation. + coeff[4] = 0.0 + coeff[5] = 1.0 + coeff[6] = yref[1] - ylist[1] +end + + +# RG_XTWOSTAR -- Compute the transformation coefficients required to +# define a simple shift, magnification which is the same in x and y, +# and rotation using two data points. + +procedure rg_xtwostar (xref, yref, xlist, ylist, coeff) + +real xref[ARB] #I x reference coordinates +real yref[ARB] #I y reference coordinates +real xlist[ARB] #I x input coordinates +real ylist[ARB] #I y input coordinates +real coeff[ARB] #O coefficient array + +real rot, mag, dxlis, dylis, dxref, dyref, cosrot, sinrot +real rg_xposangle() + +begin + # Compute the deltas. + dxlis = xlist[2] - xlist[1] + dylis = ylist[2] - ylist[1] + dxref = xref[2] - xref[1] + dyref = yref[2] - yref[1] + + # Compute the required rotation angle. + rot = rg_xposangle (dxref, dyref) - rg_xposangle (dxlis, dylis) + cosrot = cos (rot) + sinrot = sin (rot) + + # Compute the required magnification factor. + mag = dxlis ** 2 + dylis ** 2 + if (mag <= 0.0) + mag = 0.0 + else + mag = sqrt ((dxref ** 2 + dyref ** 2) / mag) + + # Compute the transformation coefficicents. + coeff[1] = mag * cosrot + coeff[2] = - mag * sinrot + coeff[3] = xref[1] - mag * cosrot * xlist[1] + mag * sinrot * ylist[1] + coeff[4] = mag * sinrot + coeff[5] = mag * cosrot + coeff[6] = yref[1] - mag * sinrot * xlist[1] - mag * cosrot * ylist[1] +end + + +# RG_THREESTAR -- Compute the transformation coefficients required to define +# x and y shifts, x and ymagnifications, a rotation and skew, and a possible +# axis flip using three tie points. + +procedure rg_xthreestar (xref, yref, xlist, ylist, coeff) + +real xref[ARB] #I x reference coordinates +real yref[ARB] #I y reference coordinates +real xlist[ARB] #I x input coordinates +real ylist[ARB] #I y input coordinates +real coeff[ARB] #O coefficient array + +real dx23, dx13, dx12, dy23, dy13, dy12, det +bool fp_equalr() + +begin + # Compute the deltas. + dx23 = xlist[2] - xlist[3] + dx13 = xlist[1] - xlist[3] + dx12 = xlist[1] - xlist[2] + dy23 = ylist[2] - ylist[3] + dy13 = ylist[1] - ylist[3] + dy12 = ylist[1] - ylist[2] + + # Compute the determinant. + det = xlist[1] * dy23 - xlist[2] * dy13 + xlist[3] * dy12 + if (fp_equalr (det, 0.0)) { + call rg_xtwostar (xref, yref, xlist, ylist, coeff) + return + } + + # Compute the x transformation. + coeff[1] = (xref[1] * dy23 - xref[2] * dy13 + xref[3] * dy12) / det + coeff[2] = (-xref[1] * dx23 + xref[2] * dx13 - xref[3] * dx12) / det + coeff[3] = (xref[1] * (xlist[2] * ylist[3] - xlist[3] * ylist[2]) + + xref[2] * (ylist[1] * xlist[3] - xlist[1] * ylist[3]) + + xref[3] * (xlist[1] * ylist[2] - ylist[1] * xlist[2])) / det + + # Compute the y transformation. + coeff[4] = (yref[1] * dy23 - yref[2] * dy13 + yref[3] * dy12) / det + coeff[5] = (-yref[1] * dx23 + yref[2] * dx13 - yref[3] * dx12) / det + coeff[6] = (yref[1] * (xlist[2] * ylist[3] - xlist[3] * ylist[2]) + + yref[2] * (ylist[1] * xlist[3] - xlist[1] * ylist[3]) + + yref[3] * (xlist[1] * ylist[2] - ylist[1] * xlist[2])) / det +end + + +# RG_XPOSANGLE -- Compute the position angle of a 2D vector. The angle is +# measured counter-clockwise from the positive x axis. + +real procedure rg_xposangle (x, y) + +real x #I x vector component +real y #I y vector component + +real theta +bool fp_equalr() + +begin + if (fp_equalr (y, 0.0)) { + if (x > 0.0) + theta = 0.0 + else if (x < 0.0) + theta = PI + else + theta = 0.0 + } else if (fp_equalr (x, 0.0)) { + if (y > 0.0) + theta = PI / 2.0 + else if (y < 0.0) + theta = 3.0 * PI / 2.0 + else + theta = 0.0 + } else if (x > 0.0 && y > 0.0) { # 1st quadrant + theta = atan (y / x) + } else if (x > 0.0 && y < 0.0) { # 4th quadrant + theta = 2.0 * PI + atan (y / x) + } else if (x < 0.0 && y > 0.0) { # 2nd quadrant + theta = PI + atan (y / x) + } else if (x < 0.0 && y < 0.0) { # 3rd quadrant + theta = PI + atan (y / x) + } + + return (theta) +end diff --git a/pkg/images/immatch/src/xregister/t_xregister.x b/pkg/images/immatch/src/xregister/t_xregister.x new file mode 100644 index 00000000..f9fc9b22 --- /dev/null +++ b/pkg/images/immatch/src/xregister/t_xregister.x @@ -0,0 +1,440 @@ +include +include +include +include +include "xregister.h" + +# T_XREGISTER -- Register a list of images using cross-correlation techniques. + +procedure t_xregister() + +pointer freglist # reference regions list +pointer database # the shifts database +int dformat # use the database format for the shifts file ? +int interactive # interactive mode ? +int verbose # verbose mode +pointer interpstr # interpolant type +int boundary # boundary extension type +real constant # constant for boundary extension + +int list1, listr, list2, reglist, reflist, reclist, tfd, stat, nregions +int c1, c2, l1, l2, ncols, nlines +pointer sp, image1, image2, imtemp, str, coords +pointer gd, id, imr, im1, im2, sdb, xc, mw +real shifts[2] +bool clgetb() +int imtopen(), imtlen(), imtgetim(), fntopnb(), clgwrd(), btoi() +int rg_xregions(), fntlenb(), rg_gxtransform(), rg_xstati() +int rg_xcorr(), rg_xicorr(), fntgfnb(), access(), open() +pointer gopen(), immap(), dtmap(), mw_openim() +real clgetr(), rg_xstatr() +errchk fntopnb(), gopen() + +begin + # Set STDOUT to flush on a newline character + call fseti (STDOUT, F_FLUSHNL, YES) + + # Allocate temporary working space. + call smark (sp) + + call salloc (freglist, SZ_LINE, TY_CHAR) + call salloc (image1, SZ_FNAME, TY_CHAR) + call salloc (image2, SZ_FNAME, TY_CHAR) + call salloc (imtemp, SZ_FNAME, TY_CHAR) + call salloc (database, SZ_FNAME, TY_CHAR) + call salloc (coords, SZ_FNAME, TY_CHAR) + call salloc (interpstr, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Get task parameters and open lists. + call clgstr ("input", Memc[str], SZ_LINE) + list1 = imtopen (Memc[str]) + call clgstr ("reference", Memc[str], SZ_LINE) + listr = imtopen (Memc[str]) + call clgstr ("regions", Memc[freglist], SZ_LINE) + call clgstr ("shifts", Memc[database], SZ_FNAME) + call clgstr ("output", Memc[str], SZ_LINE) + list2 = imtopen (Memc[str]) + call clgstr ("records", Memc[str], SZ_LINE) + if (Memc[str] == EOS) + reclist = NULL + else + reclist = fntopnb (Memc[str], NO) + call clgstr ("coords", Memc[coords], SZ_LINE) + + # Open the cross correlation fitting structure. + call rg_xgpars (xc) + + # Test the reference image list length. + if (rg_xstati (xc, CFUNC) != XC_FILE) { + 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 (Memc[coords] == EOS) + reflist = NULL + else { + reflist = fntopnb (Memc[coords], NO) + if (imtlen (listr) != fntlenb (reflist)) + call error (0, + "The number of reference point files and images is not the same.") + } + iferr { + reglist = fntopnb (Memc[freglist], NO) + } then { + reglist = NULL + } + call rg_xsets (xc, REGIONS, Memc[freglist]) + + } else { + call imtclose (listr) + listr = NULL + reflist = NULL + reglist = NULL + call rg_xsets (xc, REGIONS, "") + } + + # Close the output image list if it is empty. + if (imtlen (list2) == 0) { + call imtclose (list2) + list2 = NULL + } + + # Check that the output image list is the same size as the input + # image list. + if (list2 != NULL) { + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + if (list2 != NULL) + call imtclose (list2) + call error (0, + "The number of input and output images is not the same.") + } + } + + # Check that the record list is the same length as the input + # image list length. + if (reclist != NULL) { + if (fntlenb (reclist) != imtlen (list1)) + call error (0, + "Input image and record lists are not the same length.") + } + + + # Open the database file. + dformat = btoi (clgetb ("databasefmt")) + if (rg_xstati (xc, CFUNC) == XC_FILE) { + if (dformat == YES) + sdb = dtmap (Memc[database], READ_ONLY) + else + sdb = open (Memc[database], READ_ONLY, TEXT_FILE) + } else if (clgetb ("append")) { + if (dformat == YES) + sdb = dtmap (Memc[database], APPEND) + else + sdb = open (Memc[database], NEW_FILE, TEXT_FILE) + } else if (access (Memc[database], 0, 0) == YES) { + call error (0, "The shifts database file already exists") + } else { + if (dformat == YES) + sdb = dtmap (Memc[database], NEW_FILE) + else + sdb = open (Memc[database], NEW_FILE, TEXT_FILE) + } + call rg_xsets (xc, DATABASE, Memc[database]) + + # Get the boundary extension parameters for the image shift. + call clgstr ("interp_type", Memc[interpstr], SZ_FNAME) + boundary = clgwrd ("boundary_type", Memc[str], SZ_LINE, + "|constant|nearest|reflect|wrap|") + constant = clgetr ("constant") + + if (rg_xstati (xc, CFUNC) == XC_FILE) + interactive = NO + else + interactive = btoi (clgetb ("interactive")) + 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 { + if (rg_xstati (xc, PFUNC) == XC_MARK) + call rg_xseti (xc, PFUNC, XC_CENTROID) + gd = NULL + id = NULL + verbose = btoi (clgetb ("verbose")) + } + + # Initialize the reference image filter descriptors + imr = NULL + tfd = NULL + + # Initialize the overlap section. + c1 = INDEFI + c2 = INDEFI + l1 = INDEFI + l2 = INDEFI + ncols = INDEFI + nlines = INDEFI + + # Do each set of input, reference, and output images. + while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF)) { + + # Open the reference image, and associated regions and coordinates + # files if the correlation function is not file. + + if (rg_xstati (xc, CFUNC) != XC_FILE) { + if (imtgetim (listr, Memc[str], SZ_FNAME) != EOF) { + if (imr != NULL) + call imunmap (imr) + imr = immap (Memc[str], READ_ONLY, 0) + if (IM_NDIM(imr) > 2) + call error (0, "Reference images must be 1D or 2D") + call rg_xsets (xc, REFIMAGE, Memc[str]) + nregions = rg_xregions (reglist, imr, xc, 1) + if (nregions <= 0 && interactive == NO) + call error (0, "The regions list is empty.") + if (reflist != NULL) { + if (tfd != NULL) + call close (tfd) + tfd = rg_gxtransform (reflist, xc, Memc[str]) + call rg_xsets (xc, REFFILE, Memc[str]) + } + } + } else + call rg_xsets (xc, REFIMAGE, "reference") + + # Open the input image. + im1 = immap (Memc[image1], READ_ONLY, 0) + if (IM_NDIM(im1) > 2) { + call error (0, "Input images must be 1D or 2D") + } else if (imr != NULL) { + if (IM_NDIM(im1) != IM_NDIM(imr)) + call error (0, + "Input images must have same dimensionality as reference images") + } + 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))) + call rg_xsets (xc, IMAGE, Memc[image1]) + + # Open the output image if any. + if (list2 == NULL) { + im2 = NULL + Memc[image2] = EOS + } 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] = EOS + } + call rg_xsets (xc, OUTIMAGE, Memc[image2]) + + # Get the image record name for the shifts database. + if (reclist == NULL) + call strcpy (Memc[image1], Memc[str], SZ_FNAME) + else if (fntgfnb (reclist, Memc[str], SZ_FNAME) == EOF) + call strcpy (Memc[image1], Memc[str], SZ_FNAME) + call rg_xsets (xc, RECORD, Memc[str]) + + # Compute the initial coordinate shift. + if (tfd != NULL) + call rg_xtransform (tfd, xc) + + # Perform the cross correlation function. + if (interactive == YES) { + stat = rg_xicorr (imr, im1, im2, sdb, dformat, reglist, tfd, + xc, gd, id) + } else { + stat = rg_xcorr (imr, im1, sdb, dformat, xc) + if (verbose == YES) { + call rg_xstats (xc, REFIMAGE, Memc[str], SZ_LINE) + call printf ( + "Average shift from %s to %s is %g %g pixels\n") + call pargstr (Memc[image1]) + call pargstr (Memc[str]) + call pargr (rg_xstatr (xc, TXSHIFT)) + call pargr (rg_xstatr (xc, TYSHIFT)) + } + } + + # Compute the overlap region for the images. + call rg_overlap (im1, rg_xstatr (xc, TXSHIFT), + rg_xstatr (xc,TYSHIFT), c1, c2, l1, l2, ncols, nlines) + + # Shift the image and update the wcs. + if (im2 != NULL && stat == NO) { + if (verbose == YES) { + call printf ( + "\tShifting image %s to image %s ...\n") + call pargstr (Memc[image1]) + call pargstr (Memc[imtemp]) + } + + call rg_xshiftim (im1, im2, rg_xstatr (xc, TXSHIFT), + rg_xstatr (xc, TYSHIFT), Memc[interpstr], boundary, + constant) + mw = mw_openim (im1) + shifts[1] = rg_xstatr (xc, TXSHIFT) + shifts[2] = rg_xstatr (xc, TYSHIFT) + call mw_shift (mw, shifts, 03B) + call mw_saveim (mw, im2) + call mw_close (mw) + } + + # Close up the input and output images. + call imunmap (im1) + if (im2 != NULL) { + call imunmap (im2) + if (stat == YES) + call imdelete (Memc[image2]) + else + call xt_delimtemp (Memc[image2], Memc[imtemp]) + } + + if (stat == YES) + break + call rg_xindefr (xc) + } + + if (verbose == YES) + call rg_poverlap (c1, c2, l1, l2, ncols, nlines) + + call rg_xfree (xc) + + # Close up the lists. + if (imr != NULL) + call imunmap (imr) + call imtclose (list1) + if (listr != NULL) + call imtclose (listr) + if (reglist != NULL) + call fntclsb (reglist) + if (list2 != NULL) + call imtclose (list2) + if (tfd != NULL) + call close (tfd) + if (reflist != NULL) + call fntclsb (reflist) + if (reclist != NULL) + call fntclsb (reclist) + if (dformat == YES) + call dtunmap (sdb) + else + call close (sdb) + + # Close up the graphics and display devices. + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + + call sfree (sp) +end + + +# RG_OVERLAP -- Compute the overlap region of the list of images. + +procedure rg_overlap (im1, xshift, yshift, x1, x2, y1, y2, ncols, nlines) + +pointer im1 # pointer to the input image +real xshift # the computed x shift of the input image +real yshift # the computed y shift of the input image +int x1, x2 # the input/output column limits +int y1, y2 # the input/output line limits +int ncols, nlines # the input/output size limits + +int ixlo, ixhi, iylo, iyhi +real xlo, xhi, ylo, yhi + +begin + if (IS_INDEFR(xshift) || IS_INDEFR(yshift)) + return + + # Compute the limits of the shifted image. + xlo = 1.0 + xshift + xhi = IM_LEN(im1,1) + xshift + ylo = 1.0 + yshift + yhi = IM_LEN(im1,2) + yshift + + # Round up or down as appropriate. + ixlo = int (xlo) + if (xlo > ixlo) + ixlo = ixlo + 1 + ixhi = int (xhi) + if (xhi < ixhi) + ixhi = ixhi - 1 + iylo = int (ylo) + if (ylo > iylo) + iylo = iylo + 1 + iyhi = int (yhi) + if (yhi < iyhi) + iyhi = iyhi - 1 + + # Determine the new limits. + if (IS_INDEFI(x1)) + x1 = ixlo + else + x1 = max (ixlo, x1) + if (IS_INDEFI(x2)) + x2 = ixhi + else + x2 = min (ixhi, x2) + if (IS_INDEFI(y1)) + y1 = iylo + else + y1 = max (iylo, y1) + if (IS_INDEFI(y2)) + y2 = iyhi + else + y2 = min (iyhi, y2) + if (IS_INDEFI(ncols)) + ncols = IM_LEN(im1,1) + else + ncols = min (ncols, IM_LEN(im1,1)) + if (IS_INDEFI(nlines)) + nlines = IM_LEN(im1,2) + else + nlines = min (nlines, IM_LEN(im1,2)) +end + + +# RG_POVERLAP -- Procedure to print the overlap and/or vignetted region. + +procedure rg_poverlap (x1, x2, y1, y2, ncols, nlines) + +int x1, x2 # the input column limits +int y1, y2 # the input line limits +int ncols, nlines # the number of lines and columns + +int vx1, vx2, vy1, vy2 + +begin + vx1 = max (1, min (x1, ncols)) + vx2 = max (1, min (x2, ncols)) + vy1 = max (1, min (y1, nlines)) + vy2 = max (1, min (y2, nlines)) + + call printf ("Overlap region: [%d:%d,%d:%d]\n") + call pargi (x1) + call pargi (x2) + call pargi (y1) + call pargi (y2) + if (vx1 != x1 || vx2 != x2 || vy1 != y1 || vy2 != y2) { + call printf ("Vignetted overlap region: [%d:%d,%d:%d]\n") + call pargi (vx1) + call pargi (vx2) + call pargi (vy1) + call pargi (vy2) + } +end diff --git a/pkg/images/immatch/src/xregister/xregister.h b/pkg/images/immatch/src/xregister/xregister.h new file mode 100644 index 00000000..16c88b1e --- /dev/null +++ b/pkg/images/immatch/src/xregister/xregister.h @@ -0,0 +1,250 @@ +# Header file for XREGISTER + +# Define the cross correlation structure + +define LEN_XCSTRUCT (50 + 12 * SZ_FNAME + 12) + +define XC_RC1 Memi[$1] # pointers to 1st column of ref regions +define XC_RC2 Memi[$1+1] # pointers to 2nd column of ref regions +define XC_RL1 Memi[$1+2] # pointers to 1st line of ref regions +define XC_RL2 Memi[$1+3] # pointers to 2nd line of ref regions +define XC_RZERO Memi[$1+4] # pointers to zero pts of ref regions +define XC_RXSLOPE Memi[$1+5] # pointers to x slopes of ref regions +define XC_RYSLOPE Memi[$1+6] # pointers to y slopes of ref regions +define XC_XSHIFTS Memi[$1+7] # pointers to x shifts of ref regions +define XC_YSHIFTS Memi[$1+8] # pointers to y shifts of ref regions +define XC_NREGIONS Memi[$1+9] # total number of regions +define XC_CREGION Memi[$1+10] # the current region + +define XC_NREFPTS Memi[$1+11] # number of reference points +define XC_XREF Memi[$1+12] # pointer to x reference points +define XC_YREF Memi[$1+13] # pointer to y reference points +define XC_TRANSFORM Memi[$1+14] # pointer to the transform +define XC_IXLAG Memi[$1+15] # initial shift in x +define XC_IYLAG Memi[$1+16] # initial shift in y +define XC_XLAG Memi[$1+17] # current shift in x +define XC_YLAG Memi[$1+18] # current shift in y +define XC_DXLAG Memi[$1+19] # incremental shift in x +define XC_DYLAG Memi[$1+20] # incremental shift in y + +define XC_BACKGRD Memi[$1+21] # type of background subtraction +define XC_BORDER Memi[$1+22] # width of background border +define XC_BVALUER Memr[P2R($1+23)] # reference background value +define XC_BVALUE Memr[P2R($1+24)] # image bacground value +define XC_LOREJECT Memr[P2R($1+25)] # low side rejection +define XC_HIREJECT Memr[P2R($1+26)] # high side rejection +define XC_APODIZE Memr[P2R($1+27)] # fraction of apodized region +define XC_FILTER Memi[$1+28] # filter type + +define XC_CFUNC Memi[$1+30] # crosscor function +define XC_XWINDOW Memi[$1+31] # width of correlation window in x +define XC_YWINDOW Memi[$1+32] # width of correlation window in y +define XC_XCOR Memi[$1+33] # pointer to cross-correlation function + +define XC_PFUNC Memi[$1+34] # correlation peak fitting function +define XC_XCBOX Memi[$1+35] # x width of cor fitting box +define XC_YCBOX Memi[$1+36] # y width of cor fitting box + +define XC_TXSHIFT Memr[P2R($1+37)] # total x shift +define XC_TYSHIFT Memr[P2R($1+38)] # total y shift + +define XC_BSTRING Memc[P2C($1+50)] # background type +define XC_FSTRING Memc[P2C($1+50+SZ_FNAME+1)] # filter string +define XC_CSTRING Memc[P2C($1+50+2*SZ_FNAME+2)] # cross-correlation type +define XC_PSTRING Memc[P2C($1+50+3*SZ_FNAME+3)] # peak centering + +define XC_IMAGE Memc[P2C($1+50+4*SZ_FNAME+4)] # input image +define XC_REFIMAGE Memc[P2C($1+50+5*SZ_FNAME+5)] # reference image +define XC_REGIONS Memc[P2C($1+50+6*SZ_FNAME+6)] # regions list +define XC_DATABASE Memc[P2C($1+50+7*SZ_FNAME+7)] # shifts database +define XC_OUTIMAGE Memc[P2C($1+50+8*SZ_FNAME+8)] # output image +define XC_REFFILE Memc[P2C($1+50+9*SZ_FNAME+9)] # coordinates file +define XC_RECORD Memc[P2C($1+50+10*SZ_FNAME+10)] # record + +# Define the id strings + +define RC1 1 +define RC2 2 +define RL1 3 +define RL2 4 +define RZERO 5 +define RXSLOPE 6 +define RYSLOPE 7 +define XSHIFTS 8 +define YSHIFTS 9 +define NREGIONS 10 +define CREGION 11 + +define NREFPTS 12 +define XREF 13 +define YREF 14 +define TRANSFORM 15 +define IXLAG 16 +define IYLAG 17 +define XLAG 18 +define YLAG 19 +define DXLAG 20 +define DYLAG 21 + +define BACKGRD 22 +define BVALUER 23 +define BVALUE 24 +define BORDER 25 +define LOREJECT 26 +define HIREJECT 27 +define APODIZE 28 +define FILTER 29 + +define CFUNC 30 +define XWINDOW 31 +define YWINDOW 32 +define XCOR 33 + +define PFUNC 34 +define XCBOX 35 +define YCBOX 36 + +define TXSHIFT 37 +define TYSHIFT 38 + +define CSTRING 39 +define BSTRING 40 +define PSTRING 41 +define FSTRING 42 + +define IMAGE 43 +define REFIMAGE 44 +define REGIONS 45 +define OUTIMAGE 46 +define REFFILE 47 +define DATABASE 48 +define RECORD 49 + +# Define the default parameter values + +define DEF_IXLAG 0 +define DEF_IYLAG 0 +define DEF_DXLAG 0 +define DEF_DYLAG 0 +define DEF_XWINDOW 5 +define DEF_YWINDOW 5 + +define DEF_BACKGRD XC_BNONE +define DEF_BORDER INDEFI +define DEF_LOREJECT INDEFR +define DEF_HIREJECT INDEFR + +define DEF_XCBOX 5 +define DEF_YCBOX 5 +define DEF_PFUNC XC_CENTROID + +# Define the background fitting techniques + +define XC_BNONE 1 +define XC_MEAN 2 +define XC_MEDIAN 3 +define XC_SLOPE 4 + +define XC_BTYPES "|none|mean|median|plane|" + +# Define the filtering options + +define XC_FNONE 1 +define XC_LAPLACE 2 + +define XC_FTYPES "|none|laplace|" + +# Define the cross correlation techniques + +define XC_DISCRETE 1 +define XC_FOURIER 2 +define XC_DIFFERENCE 3 +define XC_FILE 4 + +define XC_CTYPES "|discrete|fourier|difference|file|" + +# Define the peak fitting functions + +define XC_PNONE 1 +define XC_CENTROID 2 +define XC_SAWTOOTH 3 +define XC_PARABOLA 4 +define XC_MARK 5 + +define XC_PTYPES "|none|centroid|sawtooth|parabola|mark|" + +# Miscellaneous + +define MAX_NREGIONS 100 +define MAX_NREF 3 +define MAX_NTRANSFORM 6 + +# Commands + +define XCMDS "|reference|input|regions|shifts|output|records|transform|\ +cregion|xlag|ylag|dxlag|dylag|background|border|loreject|hireject|apodize|\ +filter|correlation|xwindow|ywindow|function|xcbox|ycbox|show|mark|" + +define XSHOW "|data|background|correlation|center|" + +define XSHOW_DATA 1 +define XSHOW_BACKGROUND 2 +define XSHOW_CORRELATION 3 +define XSHOW_PEAKCENTER 4 + +define XCMD_REFIMAGE 1 +define XCMD_IMAGE 2 +define XCMD_REGIONS 3 +define XCMD_DATABASE 4 +define XCMD_OUTIMAGE 5 +define XCMD_RECORD 6 +define XCMD_REFFILE 7 +define XCMD_CREGION 8 +define XCMD_XLAG 9 +define XCMD_YLAG 10 +define XCMD_DXLAG 11 +define XCMD_DYLAG 12 +define XCMD_BACKGROUND 13 +define XCMD_BORDER 14 +define XCMD_LOREJECT 15 +define XCMD_HIREJECT 16 +define XCMD_APODIZE 17 +define XCMD_FILTER 18 +define XCMD_CORRELATION 19 +define XCMD_XWINDOW 20 +define XCMD_YWINDOW 21 +define XCMD_PEAKCENTER 22 +define XCMD_XCBOX 23 +define XCMD_YCBOX 24 +define XCMD_SHOW 25 +define XCMD_MARK 26 + +# Keywords + +define KY_REFIMAGE "reference" +define KY_IMAGE "input" +define KY_REGIONS "regions" +define KY_DATABASE "shifts" +define KY_OUTIMAGE "output" +define KY_RECORD "record" +define KY_REFFILE "coords" +define KY_NREGIONS "nregions" +define KY_CREGION "region" +define KY_XLAG "xlag" +define KY_YLAG "ylag" +define KY_DXLAG "dxlag" +define KY_DYLAG "dylag" +define KY_BACKGROUND "background" +define KY_BORDER "border" +define KY_LOREJECT "loreject" +define KY_HIREJECT "hireject" +define KY_APODIZE "apodize" +define KY_FILTER "filter" +define KY_CORRELATION "correlation" +define KY_XWINDOW "xwindow" +define KY_YWINDOW "ywindow" +define KY_PEAKCENTER "function" +define KY_XCBOX "xcbox" +define KY_YCBOX "ycbox" +define KY_TXSHIFT "xshift" +define KY_TYSHIFT "yshift" diff --git a/pkg/images/immatch/src/xregister/xregister.key b/pkg/images/immatch/src/xregister/xregister.key new file mode 100644 index 00000000..1956c88f --- /dev/null +++ b/pkg/images/immatch/src/xregister/xregister.key @@ -0,0 +1,47 @@ + Interactive Keystroke Commands + +? Print help +: Colon commands +t Define the offset between the reference and input images +c Draw a contour plot of the cross-correlation function +x Draw a column plot of the cross-correlation function +y Draw a line plot of the cross-correlation function +r Redraw the current plot +f Recompute the cross-correlation function +o Enter the image overlay plot submenu +w Update the task parameters +q Exit + + + Colon Commands + +:mark Mark regions on the display +:show Show current values of all the parameters + + + Show/set Parameters + +:reference [string] Show/set the current reference image name +:input [string] Show/set the current input image name +:regions [string] Show/set the regions to be cross-correlated +:shifts {string] Show/set the shifts database file name +:coords [string] Show/set the current coordinates file name +:output [string] Show/set the current output image name +:records [string] Show/set the current database record name +:xlag [value] Show/set the initial lag in x +:ylag [value] Show/set the initial lag in y +:dxlag [value] Show/set the incremental lag in x +:dylag [value] Show/set the incremental lag in y +:cregion [value] Show/set the current region +:background [string] Show/set the background fitting function +:border [value] Show/set border region for background fitting +: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 end points to apodize +:filter [string] Show/set the default spatial filter +:correlation [string] Show/set the cross-correlation function +:xwindow [value] Show/set width of cross-correlation window in x +:ywindow [value] Show/set width of cross-correlation window in y +:function [string] Show/set correlation peak centering function +:xcbox [value] Show/set the centering box width in x +:ycbox [value] Show/set the centering box width in y -- cgit