aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/immatch/src/psfmatch/rgpregions.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/immatch/src/psfmatch/rgpregions.x')
-rw-r--r--pkg/images/immatch/src/psfmatch/rgpregions.x464
1 files changed, 464 insertions, 0 deletions
diff --git a/pkg/images/immatch/src/psfmatch/rgpregions.x b/pkg/images/immatch/src/psfmatch/rgpregions.x
new file mode 100644
index 00000000..c04dcf97
--- /dev/null
+++ b/pkg/images/immatch/src/psfmatch/rgpregions.x
@@ -0,0 +1,464 @@
+include <fset.h>
+include <imhdr.h>
+include "psfmatch.h"
+
+# RG_PREGIONS -- Decoode the regions specification. If the sections
+# string is NULL then a default region dnx by dny pixels wide centered
+# on the reference image is used. Otherwise the section centers are
+# read from the regions string or from the objects list.
+
+int procedure rg_pregions (list, im, pm, rp, reread)
+
+int list #I pointer to regions file list
+pointer im #I pointer to the image
+pointer pm #I pointer to the psfmatch structure
+int rp #I region pointer
+int reread #I reread the current file
+
+char fname[SZ_FNAME]
+int nregions, fd
+int open(), rg_prregions(), rg_pgregions(), fntgfnb()
+int rg_pstati()
+data fname[1] /EOS/
+errchk open(), fntgfnb(), close()
+
+begin
+ if (rp < 1 || rp > MAX_NREGIONS) {
+ nregions = 0
+ } else if (rg_pgregions (im, pm, rp, MAX_NREGIONS) > 0) {
+ nregions = rg_pstati (pm, NREGIONS)
+ } else if (list != NULL) {
+ if (reread == NO) {
+ iferr {
+ if (fntgfnb (list, fname, SZ_FNAME) != EOF) {
+ fd = open (fname, READ_ONLY, TEXT_FILE)
+ nregions= rg_prregions (fd, im, pm, rp, MAX_NREGIONS)
+ call close (fd)
+ }
+ } then
+ nregions = 0
+ } else if (fname[1] != EOS) {
+ iferr {
+ fd = open (fname, READ_ONLY, TEXT_FILE)
+ nregions= rg_prregions (fd, im, pm, rp, MAX_NREGIONS)
+ call close (fd)
+ } then
+ nregions = 0
+ }
+ } else
+ nregions = 0
+
+ return (nregions)
+end
+
+
+# RG_PMKREGIONS -- Create a list of psf objects by selecting objects with
+# the image display cursor.
+
+int procedure rg_pmkregions (fd, im, pm, rp, max_nregions)
+
+int fd #I the output coordinates file descriptor
+pointer im #I pointer to the image
+pointer pm #I pointer to the psf matching structure
+int rp #I pointer to current region
+int max_nregions #I maximum number of regions
+
+int nregions, wcs, key, x1, x2, y1, y2
+pointer sp, region, cmd
+real x, y, xc, yc
+int clgcur(), rg_pstati()
+pointer rg_pstatp()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (region, SZ_FNAME, TY_CHAR)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ # Allocate the arrays to hold the regions information,
+ call rg_prealloc (pm, max_nregions)
+
+ nregions = min (rp-1, rg_pstati (pm, NREGIONS))
+ while (nregions < max_nregions) {
+
+ # Identify the object.
+ call printf ("Mark object %d [any key=mark,q=quit]:\n")
+ call pargi (nregions + 1)
+ if (clgcur ("icommands", x, y, wcs, key, Memc[cmd], SZ_LINE) == EOF)
+ break
+ if (key == 'q')
+ break
+
+ # Center the object.
+ if (rg_pstati (pm, CENTER) == YES) {
+ call rg_pcntr (im, x, y, max (rg_pstati(pm, PNX),
+ rg_pstati(pm, PNY)), xc, yc)
+ } else {
+ xc = x
+ yc = y
+ }
+
+ # Compute the data section.
+ x1 = xc - rg_pstati (pm, DNX) / 2
+ x2 = x1 + rg_pstati (pm, DNX) - 1
+ y1 = yc - rg_pstati (pm, DNY) / 2
+ y2 = y1 + rg_pstati (pm, DNY) - 1
+
+ # Make sure that the region is on the image.
+ if (x1 < 1 || x2 > IM_LEN(im,1) || y1 < 1 || y2 >
+ IM_LEN(im,2))
+ next
+
+ if (fd != NULL) {
+ call fprintf (fd, "%0.3f %0.3f\n")
+ call pargr (xc)
+ call pargr (yc)
+ }
+
+ Memi[rg_pstatp(pm,RC1)+nregions] = x1
+ Memi[rg_pstatp(pm,RC2)+nregions] = x2
+ Memi[rg_pstatp(pm,RL1)+nregions] = y1
+ Memi[rg_pstatp(pm,RL2)+nregions] = y2
+ Memr[rg_pstatp(pm,RZERO)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RXSLOPE)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RYSLOPE)+nregions] = INDEFR
+ nregions = nregions + 1
+
+ }
+
+ # Reallocate the correct amount of space.
+ call rg_pseti (pm, NREGIONS, nregions)
+ if (nregions > 0) {
+ call rg_prealloc (pm, nregions)
+ if (fd != NULL) {
+ call fstats (fd, F_FILENAME, Memc[region], SZ_FNAME)
+ call rg_psets (pm, PSFDATA, Memc[region])
+ } else
+ call rg_psets (pm, PSFDATA, "")
+ } else {
+ call rg_prfree (pm)
+ call rg_psets (pm, PSFDATA, "")
+ }
+
+ call sfree (sp)
+ return (nregions)
+end
+
+
+# RG_PRREGIONS -- Procedure to read the regions from a file.
+
+int procedure rg_prregions (fd, im, pm, rp, max_nregions)
+
+int fd #I regions file descriptor
+pointer im #I pointer to the image
+pointer pm #I pointer to psf matching structure
+int rp #I pointer to current region
+int max_nregions #I maximum number of regions
+
+int nregions, x1, y1, x2, y2
+pointer sp, line
+real x, y, xc, yc
+int rg_pstati(), getline()
+pointer rg_pstatp()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (line, SZ_LINE, TY_CHAR)
+
+ # Allocate the arrays to hold the regions information,
+ call rg_prealloc (pm, max_nregions)
+
+ # Decode the regions string.
+ nregions = min (rp - 1, rg_pstati (pm, NREGIONS))
+ while (getline (fd, Memc[line]) != EOF) {
+
+ if (nregions >= max_nregions)
+ break
+
+ call sscan (Memc[line])
+ call gargr (x)
+ call gargr (y)
+ if (rg_pstati (pm, CENTER) == YES) {
+ call rg_pcntr (im, x, y, max (rg_pstati(pm, PNX),
+ rg_pstati(pm, PNY)), xc, yc)
+ } else {
+ xc = x
+ yc = y
+ }
+
+ # Compute the data section.
+ x1 = xc - rg_pstati (pm, DNX) / 2
+ x2 = x1 + rg_pstati (pm, DNX) - 1
+ if (IM_NDIM(im) == 1) {
+ y1 = 1
+ y2 = 1
+ } else {
+ y1 = yc - rg_pstati (pm, DNY) / 2
+ y2 = y1 + rg_pstati (pm, DNY) - 1
+ }
+
+ # Make sure that the region is on the image.
+ if (x1 < 1 || x2 > IM_LEN(im,1) || y1 < 1 || y2 >
+ IM_LEN(im,2))
+ next
+
+ # Add the new region to the list.
+ Memi[rg_pstatp(pm,RC1)+nregions] = x1
+ Memi[rg_pstatp(pm,RC2)+nregions] = x2
+ Memi[rg_pstatp(pm,RL1)+nregions] = y1
+ Memi[rg_pstatp(pm,RL2)+nregions] = y2
+ Memr[rg_pstatp(pm,RZERO)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RXSLOPE)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RYSLOPE)+nregions] = INDEFR
+ nregions = nregions + 1
+ }
+
+ call rg_pseti (pm, NREGIONS, nregions)
+ if (nregions > 0)
+ call rg_prealloc (pm, nregions)
+ else
+ call rg_prfree (pm)
+
+ call sfree (sp)
+ return (nregions)
+end
+
+
+# RG_PGREGIONS -- Procedure to compute the column and line limits given
+# an x and y position and a default size.
+
+int procedure rg_pgregions (im, pm, rp, max_nregions)
+
+pointer im #I pointer to the image
+pointer pm #I pointer to psf matching structure
+int rp #I pointer to the current region
+int max_nregions #I maximum number of regions
+
+int ncols, nlines, nregions
+int x1, x2, y1, y2
+pointer sp, region
+real x, y, xc, yc
+int rg_pstati(), nscan()
+pointer rg_pstatp()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (region, SZ_LINE, TY_CHAR)
+
+ # Allocate the arrays to hold the regions information.
+ call rg_prealloc (pm, max_nregions)
+
+ # Get the constants.
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+
+ # Decode the center.
+ call rg_pstats (pm, PSFDATA, Memc[region], SZ_LINE)
+ nregions = min (rp - 1, rg_pstati (pm, NREGIONS))
+ call sscan (Memc[region])
+ call gargr (x)
+ call gargr (y)
+
+ # Compute the data region.
+ if (nscan() >= 2) {
+
+ # Compute a more accurate center.
+ if (rg_pstati (pm, CENTER) == YES) {
+ call rg_pcntr (im, x, y, max (rg_pstati(pm, PNX),
+ rg_pstati(pm, PNY)), xc, yc)
+ } else {
+ xc = x
+ yc = y
+ }
+
+ # Compute the data section.
+ x1 = xc - rg_pstati (pm, DNX) / 2
+ x2 = x1 + rg_pstati (pm, DNX) - 1
+ if (IM_NDIM(im) == 1) {
+ y1 = 1
+ y2 = 1
+ } else {
+ y1 = yc - rg_pstati (pm, DNY) / 2
+ y2 = y1 + rg_pstati (pm, DNY) - 1
+ }
+
+ # Make sure that the region is on the image.
+ if (x1 >= 1 && x2 <= IM_LEN(im,1) && y1 >= 1 &&
+ y2 <= IM_LEN(im,2)) {
+ Memi[rg_pstatp(pm,RC1)+nregions] = x1
+ Memi[rg_pstatp(pm,RC2)+nregions] = x2
+ Memi[rg_pstatp(pm,RL1)+nregions] = y1
+ Memi[rg_pstatp(pm,RL2)+nregions] = y2
+ Memr[rg_pstatp(pm,RZERO)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RXSLOPE)+nregions] = INDEFR
+ Memr[rg_pstatp(pm,RYSLOPE)+nregions] = INDEFR
+ nregions = nregions + 1
+ }
+ }
+
+
+ # Reallocate the correct amount of space.
+ call rg_pseti (pm, NREGIONS, nregions)
+ if (nregions > 0)
+ call rg_prealloc (pm, nregions)
+ else
+ call rg_prfree (pm)
+
+ call sfree (sp)
+
+ return (nregions)
+end
+
+
+# RG_PCNTR -- Compute star center using MPC algorithm.
+
+procedure rg_pcntr (im, xstart, ystart, boxsize, xcntr, ycntr)
+
+pointer im #I pointer to the input image
+real xstart, ystart #I initial position
+int boxsize #I width of the centering box
+real xcntr, ycntr #O computed center
+
+int x1, x2, y1, y2, half_box
+int ncols, nrows, nx, ny, try
+real xinit, yinit
+pointer bufptr, sp, x_vect, y_vect
+int imgs2r()
+
+begin
+ # Inialize.
+ half_box = (boxsize - 1) / 2
+ xinit = xstart
+ ncols = IM_LEN (im, 1)
+ if (IM_NDIM(im) == 1) {
+ yinit = 1
+ nrows = 1
+ } else {
+ yinit = ystart
+ nrows = IM_LEN (im, 2)
+ }
+ try = 0
+
+ # Iterate until pixel shifts are less than one.
+ repeat {
+
+ # Define region to extract.
+ x1 = max (xinit - half_box, 1.0) +0.5
+ x2 = min (xinit + half_box, real(ncols)) +0.5
+ y1 = max (yinit - half_box, 1.0) +0.5
+ y2 = min (yinit + half_box, real(nrows)) +0.5
+ nx = x2 - x1 + 1
+ ny = y2 - y1 + 1
+
+ # Extract region around center
+ bufptr = imgs2r (im, x1, x2, y1, y2)
+
+ # Compute the new center.
+ call smark (sp)
+ if (IM_NDIM(im) == 1) {
+ call salloc (x_vect, nx, TY_REAL)
+ call aclrr (Memr[x_vect], nx)
+ call rg_prowsum (Memr[bufptr], Memr[x_vect], nx, ny)
+ call rg_pcenter (Memr[x_vect], nx, xcntr)
+ ycntr = 1
+ } else {
+ call salloc (x_vect, nx, TY_REAL)
+ call salloc (y_vect, ny, TY_REAL)
+ call aclrr (Memr[x_vect], nx)
+ call aclrr (Memr[y_vect], ny)
+ call rg_prowsum (Memr[bufptr], Memr[x_vect], nx, ny)
+ call rg_pcolsum (Memr[bufptr], Memr[y_vect], nx, ny)
+ call rg_pcenter (Memr[x_vect], nx, xcntr)
+ call rg_pcenter (Memr[y_vect], ny, ycntr)
+ }
+ call sfree (sp)
+
+ # Check for INDEF centers.
+ if (IS_INDEFR(xcntr) || IS_INDEFR(ycntr)) {
+ xcntr = xinit
+ ycntr = yinit
+ break
+ }
+
+ # Add in offsets
+ xcntr = xcntr + x1
+ ycntr = ycntr + y1
+
+ try = try + 1
+ if (try == 1) {
+ if ((abs(xcntr-xinit) > 1.0) || (abs(ycntr-yinit) > 1.0)) {
+ xinit = xcntr
+ yinit = ycntr
+ }
+ } else
+ break
+ }
+end
+
+
+# RG_PROWSUM -- Sum all rows in a raster.
+
+procedure rg_prowsum (v, row, nx, ny)
+
+real v[nx,ny] #I the input subraster
+real row[ARB] #O the output row sum
+int nx, ny #I the dimensions of the subraster
+
+int i, j
+
+begin
+ do i = 1, ny
+ do j = 1, nx
+ row[j] = row[j] + v[j,i]
+end
+
+
+# RG_PCOLSUM -- Sum all columns in a raster.
+
+procedure rg_pcolsum (v, col, nx, ny)
+
+real v[nx,ny] #I the input subraster
+real col[ARB] #O the output column sum
+int nx, ny #I the dimensions of the subraster
+
+int i, j
+
+begin
+ do i = 1, ny
+ do j = 1, nx
+ col[j] = col[j] + v[i,j]
+end
+
+
+# RG_PCENTER -- Compute center of gravity of array.
+
+procedure rg_pcenter (v, nv, vc)
+
+real v[ARB] #I the input vector
+int nv #I the length of the vector
+real vc #O the output center
+
+int i
+real sum1, sum2, sigma, cont
+
+begin
+ # Compute first moment
+ sum1 = 0.0
+ sum2 = 0.0
+
+ call aavgr (v, nv, cont, sigma)
+
+ do i = 1, nv
+ if (v[i] > cont) {
+ sum1 = sum1 + (i-1) * (v[i] - cont)
+ sum2 = sum2 + (v[i] - cont)
+ }
+
+ # Determine center
+ if (sum2 == 0.0)
+ vc = INDEFR
+ else
+ vc = sum1 / sum2
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:30:27 -0400