aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/immatch/src/psfmatch/rgpisfm.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/immatch/src/psfmatch/rgpisfm.x')
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpisfm.x556
1 files changed, 556 insertions, 0 deletions
diff --git a/pkg/images/immatch/src/psfmatch/rgpisfm.x b/pkg/images/immatch/src/psfmatch/rgpisfm.x
new file mode 100644
index 00000000..24df8fd7
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpisfm.x
@@ -0,0 +1,556 @@
+include <imhdr.h>
+include <ctype.h>
+include <gset.h>
+include "psfmatch.h"
+
+define HELPFILE "immatch$src/psfmatch/psfmatch.key"
+
+# Define the plot functions
+
+define PM_PPOWER 1
+define PM_PKERNEL 2
+
+# Define the plot types
+
+define PM_PCONTOUR 1
+define PM_PLINE 2
+define PM_PCOL 3
+
+# RG_PISFM -- Procedure to compute the shifts interactively.
+
+int procedure rg_pisfm (pm, imr, reglist, impsf, im1, imk, imp, im2, gd, id)
+
+pointer pm #I pointer to the psfmatch structure
+pointer imr #I/O pointer to the reference image/psf
+pointer reglist #I/O pointer to the regions list
+pointer impsf #I/O pointer to the input psf
+pointer im1 #I/O pointer to the input image
+pointer imp #I/O pointer to the fourier spectrum image
+pointer imk #I/O pointer to the kernel image
+pointer im2 #I/O pointer to the output image
+pointer gd #I graphics stream pointer
+pointer id #I display stream pointer
+
+int newref, newimage, newfourier, newfilter, plotfunc, plottype, wcs, key
+int newplot, ncolr, nliner, ip
+pointer sp, cmd
+real wx, wy
+int rg_pstati(), rg_psfm(), clgcur(), rg_pgqverify(), rg_pgtverify()
+int ctoi(), rg_pregions()
+pointer rg_pstatp()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ newref = YES
+ newimage = YES
+ newfourier = YES
+ newfilter = YES
+ ncolr = INDEFI
+ nliner = INDEFI
+ plotfunc = PM_PKERNEL
+ plottype = PM_PCONTOUR
+
+ # Compute the convolution kernel for the current image.
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONIMAGE && rg_pstati (pm,
+ NREGIONS) <= 0) {
+ call gclear (gd)
+ call gflush (gd)
+ call printf ("The objects list is empty\n")
+ } else {
+ if (rg_psfm (pm, imr, im1, impsf, imk, newref) == OK) {
+ call rg_pplot (gd, pm, ncolr, nliner, plotfunc, plottype)
+ newref = NO
+ newimage = NO
+ newfourier = NO
+ newfilter = NO
+ } else {
+ call gclear (gd)
+ call gflush (gd)
+ call rg_pstats (pm, IMAGE, Memc[cmd], SZ_FNAME)
+ call printf ("Error computing kernel for image %s\n")
+ call pargstr (Memc[cmd])
+ }
+ }
+ newplot = NO
+
+ # Loop over the cursor commands.
+ while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+ EOF) {
+
+ switch (key) {
+
+ # Print the help page.
+ case '?':
+ call gpagefile (gd, HELPFILE, "")
+
+ # Quit the task gracefully.
+ case 'q':
+ if (rg_pgqverify ("psfmatch", pm, imk, key) == YES) {
+ call sfree (sp)
+ return (rg_pgtverify (key))
+ }
+
+ # Process colon commands.
+ case ':':
+ for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1)
+ ;
+ switch (Memc[cmd+ip-1]) {
+
+ case 'x':
+ if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') {
+ call rg_pcolon (gd, pm, imr, reglist, impsf, im1, imk,
+ NULL, im2, Memc[cmd], newref, newimage,
+ newfourier, newfilter)
+ } else {
+ ip = ip + 1
+ if (ctoi (Memc[cmd], ip, ncolr) <= 0) {
+ switch (plotfunc) {
+ case PM_PPOWER:
+ ncolr = rg_pstati (pm, NXFFT) / 2 + 1
+ case PM_PKERNEL:
+ ncolr = rg_pstati (pm, KNX) / 2 + 1
+ default:
+ ncolr = rg_pstati (pm, KNX) / 2 + 1
+ }
+ }
+ plottype = PM_PCOL
+ newplot = YES
+ }
+
+ case 'y':
+ if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') {
+ call rg_pcolon (gd, pm, imr, reglist, impsf, im1, imk,
+ NULL, im2, Memc[cmd], newref, newimage,
+ newfourier, newfilter)
+ } else {
+ ip = ip + 1
+ if (ctoi (Memc[cmd], ip, nliner) <= 0) {
+ switch (plotfunc) {
+ case PM_PPOWER:
+ nliner = rg_pstati (pm, NYFFT) / 2 + 1
+ case PM_PKERNEL:
+ nliner = rg_pstati (pm, KNY) / 2 + 1
+ default:
+ nliner = rg_pstati (pm, KNY) / 2 + 1
+ }
+ }
+ plottype = PM_PLINE
+ newplot = YES
+ }
+
+
+ default:
+ call rg_pcolon (gd, pm, imr, reglist, impsf, im1, imk, NULL,
+ im2, Memc[cmd], newref, newimage, newfourier,
+ newfilter)
+ }
+
+ # Write the parameters to the parameter file.
+ case 'w':
+ call rg_pppars (pm)
+
+ # Recompute the convolution kernel function.
+ case 'f':
+
+ if (rg_pstati(pm,CONVOLUTION) == PM_CONIMAGE) {
+ if (newref == YES)
+ if (rg_pregions (reglist, imr, pm, 1, YES) > 0)
+ ;
+ else if (newimage == YES)
+ call rg_pindefr (pm)
+ }
+
+ if (rg_pstati (pm, NREGIONS) > 0 || rg_pstati (pm,
+ CONVOLUTION) != PM_CONIMAGE) {
+
+ if (newfourier == YES) {
+ call printf (
+ "\nRecomputing convolution kernel ...\n")
+ if (rg_psfm (pm, imr, im1, impsf, imk,
+ newref) == OK) {
+ ncolr = INDEFI
+ nliner = INDEFI
+ call rg_pplot (gd, pm, ncolr, nliner, plotfunc,
+ plottype)
+ newref = NO
+ newimage = NO
+ newfourier = NO
+ newfilter = NO
+ newplot = NO
+ } else
+ call printf (
+ "\nError computing new kernel ...\n")
+ }
+
+ if (newfilter == YES) {
+ if (Memr[rg_pstatp(pm,FFT)] != NULL) {
+ call rg_pfilter (pm)
+ ncolr = INDEFI
+ nliner = INDEFI
+ call rg_pplot (gd, pm, ncolr, nliner, plotfunc,
+ plottype)
+ newfilter = NO
+ newplot = NO
+ } else
+ call printf (
+ "The kernel fourier spectrum is undefined\n")
+ }
+
+ } else
+ call printf ("The objects list is empty\n")
+
+ # Draw a contour plot of the kernel.
+ case 'k':
+ if (plotfunc != PM_PKERNEL)
+ newplot = YES
+ if (plottype != PM_PCONTOUR)
+ newplot = YES
+ plotfunc = PM_PKERNEL
+ plottype = PM_PCONTOUR
+ ncolr = (1 + rg_pstati (pm, KNX)) / 2
+ nliner = (1 + rg_pstati (pm, KNY)) / 2
+
+ # Draw a contour plot of the fourier spectrum.
+ case 'p':
+ if (plotfunc != PM_PPOWER)
+ newplot = YES
+ if (plottype != PM_PCONTOUR)
+ newplot = YES
+ plotfunc = PM_PPOWER
+ plottype = PM_PCONTOUR
+ ncolr = (1 + rg_pstati (pm, NXFFT)) / 2
+ nliner = (1 + rg_pstati (pm, NYFFT)) / 2
+
+ # Plot a line of the current plot.
+ case 'x':
+ if (plottype != PM_PCOL)
+ newplot = YES
+ if (plottype == PM_PCONTOUR) {
+ ncolr = nint (wx)
+ nliner = nint (wy)
+ } else if (plottype == PM_PLINE) {
+ ncolr = nint (wx)
+ }
+ plottype = PM_PCOL
+
+ # Plot a line of the current plot.
+ case 'y':
+ if (plottype != PM_PLINE)
+ newplot = YES
+ if (plottype == PM_PCONTOUR) {
+ ncolr = nint (wx)
+ nliner = nint (wy)
+ } else if (plottype == PM_PCOL) {
+ ncolr = nint (wx)
+ }
+ plottype = PM_PLINE
+
+ # Redraw the current plot.
+ case 'r':
+ newplot = YES
+
+ # Do nothing gracefully.
+ default:
+ ;
+
+ }
+
+ if (newplot == YES) {
+ if (rg_pstati (pm, CONVOLUTION) == PM_CONIMAGE &&
+ rg_pstati (pm, NREGIONS) <= 0) {
+ call printf ("Warning: The objects list is empty\n")
+ } else if (newref == YES || newimage == YES ||
+ newfourier == YES || newfilter == YES) {
+ call printf (
+ "Warning: Convolution kernel should be refit\n")
+ } else if (rg_pstatp (pm, CONV) != NULL) {
+ call rg_pplot (gd, pm, ncolr, nliner, plotfunc, plottype)
+ newplot = NO
+ } else {
+ call printf (
+ "Warning: The convolution kernel is undefined\n")
+ }
+ }
+
+ }
+
+ call sfree (sp)
+end
+
+
+define QUERY "[Hit return to continue, n next image, q quit, w quit and update parameters]"
+
+# RG_PGQVERIFY -- Print a message in the status line asking the user if they
+# really want to quit, returning YES if they really want to quit, NO otherwise.
+
+int procedure rg_pgqverify (task, pm, imk, ch)
+
+char task[ARB] # task name
+pointer pm # pointer to psfmatch structure
+pointer imk # pointer to kernel image
+int ch # character keystroke command
+
+int wcs, stat
+pointer sp, cmd
+real wx, wy
+bool streq()
+int clgcur(), rg_pstati()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ # Print the status line query in reverse video and get the keystroke.
+ call printf (QUERY)
+ if (clgcur ("gcommands", wx, wy, wcs, ch, Memc[cmd], SZ_LINE) == EOF)
+ ;
+
+ # Process the command.
+ if (ch == 'q') {
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL)
+ call rg_pwrite (pm, imk, NULL)
+ stat = YES
+ } else if (ch == 'w') {
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL)
+ call rg_pwrite (pm, imk, NULL)
+ if (streq ("psfmatch", task))
+ call rg_pppars (pm)
+ stat = YES
+ } else if (ch == 'n') {
+ if (rg_pstati (pm, CONVOLUTION) != PM_CONKERNEL)
+ call rg_pwrite (pm, imk, NULL)
+ stat = YES
+ } else {
+ stat = NO
+ }
+
+ call sfree (sp)
+
+ return (stat)
+end
+
+
+# RG_PGTVERIFY -- Verify whether or not the user truly wishes to quit the
+# task.
+
+int procedure rg_pgtverify (ch)
+
+int ch #I the input keystroke command
+
+begin
+ if (ch == 'q') {
+ return (YES)
+ } else if (ch == 'w') {
+ return (YES)
+ } else if (ch == 'n') {
+ return (NO)
+ } else {
+ return (NO)
+ }
+end
+
+
+# RG_PPLOT -- Draw the default plot of the kernel fourier spectrum or the
+# kernel itself.
+
+procedure rg_pplot (gd, pm, col, line, plotfunc, plottype)
+
+pointer gd #I pointer to the graphics stream
+pointer pm #I pointer to the psfmatch structure
+int col #I column of cross-correlation function to plot
+int line #I line of cross-correlation function to plot
+int plotfunc #I the default plot function type
+int plottype #I the default plot type
+
+int nx, ny
+pointer sp, title, str, data
+int rg_pstati(), strlen()
+pointer rg_pstatp()
+
+begin
+ if (gd == NULL)
+ return
+
+ # Allocate working space.
+ call smark (sp)
+ call salloc (title, SZ_LINE, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Initialize the plot title and data.
+ switch (plotfunc) {
+ case PM_PPOWER:
+ call sprintf (Memc[title], SZ_LINE,
+ "Fourier Spectrum for Reference: %s Image: %s")
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ data = rg_pstatp (pm, ASFFT)
+ nx = rg_pstati (pm, NXFFT)
+ ny = rg_pstati (pm, NYFFT)
+ case PM_PKERNEL:
+ call sprintf (Memc[title], SZ_LINE,
+ "Convolution Kernel for Reference: %s Image: %s")
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ data = rg_pstatp (pm, CONV)
+ nx = rg_pstati (pm, KNX)
+ ny = rg_pstati (pm, KNY)
+ default:
+ call sprintf (Memc[title], SZ_LINE,
+ "Convolution Kernel for Reference: %s Image: %s")
+ call rg_pstats (pm, REFIMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ call rg_pstats (pm, IMAGE, Memc[str], SZ_FNAME)
+ call pargstr (Memc[str])
+ data = rg_pstatp (pm, CONV)
+ nx = rg_pstati (pm, KNX)
+ nx = rg_pstati (pm, KNY)
+ }
+ if (IS_INDEFI(col))
+ col = 1 + nx / 2
+ if (IS_INDEFI(line))
+ line = 1 + ny / 2
+
+ # Draw the plot.
+ if (ny == 1) {
+ switch (plotfunc) {
+ case PM_PPOWER:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nLine %d")
+ call pargi (1)
+ call rg_pcpline (gd, Memc[title], Memr[rg_pstatp(pm,ASFFT)],
+ nx, ny, 1)
+ case PM_PKERNEL:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nLine %d")
+ call pargi (1)
+ call rg_pcpline (gd, Memc[title], Memr[rg_pstatp(pm,CONV)],
+ nx, ny, 1)
+ default:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nLine %d")
+ call pargi (1)
+ call rg_pcpline (gd, Memc[title], Memr[rg_pstatp(pm,CONV)],
+ nx, ny, 1)
+ }
+ } else {
+ switch (plottype) {
+ case PM_PCONTOUR:
+ call rg_contour (gd, Memc[title], "", Memr[data], nx, ny)
+ case PM_PLINE:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nLine %d")
+ call pargi (line)
+ call rg_pcpline (gd, Memc[title], Memr[data], nx, ny, line)
+ case PM_PCOL:
+ call sprintf (Memc[title+strlen(Memc[title])], SZ_LINE,
+ "\nColumn %d")
+ call pargi (col)
+ call rg_pcpcol (gd, Memc[title], Memr[data], nx, ny, col)
+ default:
+ call rg_contour (gd, Memc[title], "", Memr[data], nx, ny)
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# RG_PCPLINE -- Plot a line of a 2D function.
+
+procedure rg_pcpline (gd, title, data, nx, ny, nline)
+
+pointer gd #I pointer to the graphics stream
+char title[ARB] #I title for the plot
+real data[nx,ARB] #I the input data array
+int nx, ny #I dimensions of the input data array
+int nline #I the line number
+
+int i
+pointer sp, str, x
+real ymin, ymax
+
+begin
+ # Return if no graphics stream.
+ if (gd == NULL)
+ return
+
+ # Check for valid line number.
+ if (nline < 1 || nline > ny)
+ return
+
+ # Allocate some working space.
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (x, nx, TY_REAL)
+
+ # Initialize the data.
+ do i = 1, nx
+ Memr[x+i-1] = i
+ call alimr (data[1,nline], nx, ymin, ymax)
+
+ # Set up the labels and the axes.
+ call gclear (gd)
+ call gswind (gd, 1.0, real (nx), ymin, ymax)
+ call glabax (gd, title, "X Lag", "X-Correlation Function")
+
+ # Plot the line profile.
+ call gseti (gd, G_PLTYPE, GL_SOLID)
+ call gpline (gd, Memr[x], data[1,nline], nx)
+ call gflush (gd)
+
+ call sfree (sp)
+end
+
+
+# RG_PCPCOL -- Plot a column of the cross-correlation function.
+
+procedure rg_pcpcol (gd, title, data, nx, ny, ncol)
+
+pointer gd #I pointer to the graphics stream
+char title[ARB] #I title of the column plot
+real data[nx,ARB] #I the input data array
+int nx, ny #I the dimensions of the input data array
+int ncol #I line number
+
+int i
+pointer sp, x, y
+real ymin, ymax
+
+begin
+ # Return if no graphics stream.
+ if (gd == NULL)
+ return
+
+ # Check for valid column number.
+ if (ncol < 1 || ncol > nx)
+ return
+
+ # Initialize.
+ call smark (sp)
+ call salloc (x, ny, TY_REAL)
+ call salloc (y, ny, TY_REAL)
+
+ # Get the data to be plotted.
+ do i = 1, ny {
+ Memr[x+i-1] = i
+ Memr[y+i-1] = data[ncol,i]
+ }
+ call alimr (Memr[y], ny, ymin, ymax)
+
+ # Set up the labels and the axes.
+ call gclear (gd)
+ call gswind (gd, 1.0, real (ny), ymin, ymax)
+ call glabax (gd, title, "Y Lag", "X-Correlation Function")
+
+ # Plot the profile.
+ call gseti (gd, G_PLTYPE, GL_SOLID)
+ call gpline (gd, Memr[x], Memr[y], ny)
+
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 13:47:47 -0400