diff options
Diffstat (limited to 'pkg/images/immatch/src/geometry')
-rw-r--r-- | pkg/images/immatch/src/geometry/geofunc.gx | 250 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/geofunc.x | 340 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/geotimtran.x | 543 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/geotran.h | 52 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/geotran.x | 1752 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/geoxytran.gx | 327 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/geoxytran.x | 446 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/mkpkg | 34 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/t_geomap.gx | 921 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/t_geomap.x | 1509 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/t_geotran.x | 880 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/t_geoxytran.x | 343 | ||||
-rw-r--r-- | pkg/images/immatch/src/geometry/trinvert.x | 163 |
13 files changed, 7560 insertions, 0 deletions
diff --git a/pkg/images/immatch/src/geometry/geofunc.gx b/pkg/images/immatch/src/geometry/geofunc.gx new file mode 100644 index 00000000..3b34a207 --- /dev/null +++ b/pkg/images/immatch/src/geometry/geofunc.gx @@ -0,0 +1,250 @@ +include <math.h> +include <math/gsurfit.h> + +$for (rd) + +# GEO_DROTMAG -- Adjust the coefficients of the fit using the database file. + +procedure geo_drotmag$t (dt, rec, sx1, sy1, xmag, ymag, xrot, yrot) + +pointer dt #I pointer to the text database file +int rec #I record number +pointer sx1, sy1 #I/O pointers to the x and y linear surfaces +PIXEL xmag, ymag #I/O the x and y magnification +PIXEL xrot, yrot #I/O the x and y axis rotation + +real dtgetr() + +begin + if (IS_$INDEF$T(xmag)) + xmag = PIXEL (dtgetr (dt, rec, "xmag")) + if (IS_$INDEF$T(ymag)) + ymag = PIXEL (dtgetr (dt, rec, "ymag")) + if (IS_$INDEF$T(xrot)) + xrot = DEGTORAD (PIXEL(dtgetr (dt, rec, "xrotation"))) + else + xrot = DEGTORAD(xrot) + if (IS_$INDEF$T(yrot)) + yrot = DEGTORAD (PIXEL (dtgetr (dt, rec, "yrotation"))) + else + yrot = DEGTORAD(yrot) + call geo_rotmag$t (sx1, sy1, xmag, ymag, xrot, yrot) +end + + +# GEO_DXYSHIFT -- Adjust the fitted xy shift using the database file. + +procedure geo_dxyshift$t (dt, rec, sx1, sy1, xout, yout, xref, yref, + xshift, yshift) + +pointer dt #I pointer to the text file database +int rec #I the database record +pointer sx1, sy1 #I/O pointers to the x and y linear surfaces +PIXEL xout, yout #I the input coordinate system origin +PIXEL xref, yref #I the reference coordinate system origin +PIXEL xshift, yshift #I the origin shift in input coordinates + +$if (datatype == r) +PIXEL gsgetr(), gseval() +$else +PIXEL dgsgetd(), dgseval() +$endif + +begin +$if (datatype == r) + if (IS_$INDEF$T(xref)) + xref = (gsgetr (sx1, GSXMIN) + gsgetr (sx1, GSXMAX)) / 2.0 + if (IS_$INDEF$T(yref)) + yref = (gsgetr (sy1, GSYMIN) + gsgetr (sy1, GSYMAX)) / 2.0 + + if (IS_$INDEF$T(xout)) + xout = gseval (sx1, xref, yref) + if (IS_$INDEF$T(yout)) + yout = gseval (sy1, xref, yref) + + if (IS_$INDEF$T(xshift)) + xshift = xout - gseval (sx1, xref, yref) + if (IS_$INDEF$T(yshift)) + yshift = yout - gseval (sy1, xref, yref) +$else + if (IS_$INDEF$T(xref)) + xref = (dgsgetd (sx1, GSXMIN) + dgsgetd (sx1, GSXMAX)) / 2.0d0 + if (IS_$INDEF$T(yref)) + yref = (dgsgetd (sy1, GSYMIN) + dgsgetd (sy1, GSYMAX)) / 2.0d0 + + if (IS_$INDEF$T(xout)) + xout = dgseval (sx1, xref, yref) + if (IS_$INDEF$T(yout)) + yout = dgseval (sy1, xref, yref) + + if (IS_$INDEF$T(xshift)) + xshift = xout - dgseval (sx1, xref, yref) + if (IS_$INDEF$T(yshift)) + yshift = yout - dgseval (sy1, xref, yref) +$endif + + call geo_xyshift$t (sx1, sy1, xshift, yshift) +end + + +# GEO_ROTMAG -- Edit the coefficients of the linear surface which determine +# magnification and rotation. + +procedure geo_rotmag$t (sx1, sy1, xscale, yscale, xrotation, yrotation) + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +PIXEL xscale, yscale #I the x and y scales +PIXEL xrotation,yrotation #I the x and y axis rotation angles in radians + +PIXEL cosx, sinx, cosy, siny, xrange, yrange +int ncoeff +pointer sp, xcoeff, ycoeff +$if (datatype == r) +real gsgetr() +int gsgeti() +$else +double dgsgetd() +int dgsgeti() +$endif + +begin + # Get the current solution. + call smark (sp) +$if (datatype == r) + ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE)) +$else + ncoeff = max (dgsgeti (sx1, GSNSAVE), dgsgeti (sy1, GSNSAVE)) +$endif + call salloc (xcoeff, ncoeff, TY_PIXEL) + call salloc (ycoeff, ncoeff, TY_PIXEL) +$if (datatype == r) + call gssave (sx1, Mem$t[xcoeff]) + call gssave (sy1, Mem$t[ycoeff]) +$else + call dgssave (sx1, Mem$t[xcoeff]) + call dgssave (sy1, Mem$t[ycoeff]) +$endif + + # Define the scaling parameters. + cosx = cos (xrotation) + sinx = sin (xrotation) + cosy = cos (yrotation) + siny = sin (yrotation) + + # Calculate coefficients. + Mem$t[xcoeff+GS_SAVECOEFF+1] = xscale * cosx + Mem$t[xcoeff+GS_SAVECOEFF+2] = yscale * siny + Mem$t[ycoeff+GS_SAVECOEFF+1] = -xscale * sinx + Mem$t[ycoeff+GS_SAVECOEFF+2] = yscale * cosy + + # Normalize coefficients for-non polynomial functions. +$if (datatype == r) + if (gsgeti (sx1, GSTYPE) != GS_POLYNOMIAL) { + xrange = gsget$t (sx1, GSXMAX) - gsget$t (sx1, GSXMIN) +$else + if (dgsgeti (sx1, GSTYPE) != GS_POLYNOMIAL) { + xrange = dgsget$t (sx1, GSXMAX) - dgsget$t (sx1, GSXMIN) +$endif + Mem$t[xcoeff+GS_SAVECOEFF+1] = Mem$t[xcoeff+GS_SAVECOEFF+1] * + xrange / 2.d0 + Mem$t[xcoeff+GS_SAVECOEFF+2] = Mem$t[xcoeff+GS_SAVECOEFF+2] * + yrange / 2.d0 + } +$if (datatype == r) + if (gsgeti (sy1, GSTYPE) != GS_POLYNOMIAL) { + yrange = gsget$t (sy1, GSYMAX) - gsget$t (sy1, GSYMIN) +$else + if (dgsgeti (sy1, GSTYPE) != GS_POLYNOMIAL) { + yrange = dgsget$t (sy1, GSYMAX) - dgsget$t (sy1, GSYMIN) +$endif + Mem$t[ycoeff+GS_SAVECOEFF+1] = Mem$t[ycoeff+GS_SAVECOEFF+1] * + xrange / 2.d0 + Mem$t[ycoeff+GS_SAVECOEFF+2] = Mem$t[ycoeff+GS_SAVECOEFF+2] * + yrange / 2.d0 + } + +$if (datatype == r) + # Free the original fit. + call gsfree (sx1) + call gsfree (sy1) + + # Restore the edited fit. + call gsrestore (sx1, Mem$t[xcoeff]) + call gsrestore (sy1, Mem$t[ycoeff]) +$else + # Free the original fit. + call dgsfree (sx1) + call dgsfree (sy1) + + # Restore the edited fit. + call dgsrestore (sx1, Mem$t[xcoeff]) + call dgsrestore (sy1, Mem$t[ycoeff]) +$endif + + call sfree (sp) +end + + +# GEO_XYSHIFT -- Shift the linear part of the fit in x and y. + +procedure geo_xyshift$t (sx1, sy1, xshift, yshift) + +pointer sx1, sy1 #I pointers to linear x and y surfaces +PIXEL xshift, yshift #I the input x and y shifts + +int ncoeff +pointer sp, xcoeff, ycoeff +$if (datatype == r) +int gsgeti() +$else +int dgsgeti() +$endif + +begin + call smark (sp) + + # Allocate working space. +$if (datatype == r) + ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE)) +$else + ncoeff = max (dgsgeti (sx1, GSNSAVE), dgsgeti (sy1, GSNSAVE)) +$endif + call salloc (xcoeff, ncoeff, TY_PIXEL) + call salloc (ycoeff, ncoeff, TY_PIXEL) + + # Get coefficients. +$if (datatype == r) + call gssave (sx1, Mem$t[xcoeff]) + call gssave (sy1, Mem$t[ycoeff]) +$else + call dgssave (sx1, Mem$t[xcoeff]) + call dgssave (sy1, Mem$t[ycoeff]) +$endif + + # Shift the coefficients. + Mem$t[xcoeff+GS_SAVECOEFF] = Mem$t[xcoeff+GS_SAVECOEFF] + xshift + Mem$t[ycoeff+GS_SAVECOEFF] = Mem$t[ycoeff+GS_SAVECOEFF] + yshift + +$if (datatype == r) + # Free original fit. + call gsfree (sx1) + call gsfree (sy1) + + # Restore fit. + call gsrestore (sx1, Mem$t[xcoeff]) + call gsrestore (sy1, Mem$t[ycoeff]) +$else + # Free original fit. + call dgsfree (sx1) + call dgsfree (sy1) + + # Restore fit. + call dgsrestore (sx1, Mem$t[xcoeff]) + call dgsrestore (sy1, Mem$t[ycoeff]) +$endif + + call sfree (sp) +end + + +$endfor diff --git a/pkg/images/immatch/src/geometry/geofunc.x b/pkg/images/immatch/src/geometry/geofunc.x new file mode 100644 index 00000000..c3be8fa5 --- /dev/null +++ b/pkg/images/immatch/src/geometry/geofunc.x @@ -0,0 +1,340 @@ +include <math.h> +include <math/gsurfit.h> + + + +# GEO_DROTMAG -- Adjust the coefficients of the fit using the database file. + +procedure geo_drotmagr (dt, rec, sx1, sy1, xmag, ymag, xrot, yrot) + +pointer dt #I pointer to the text database file +int rec #I record number +pointer sx1, sy1 #I/O pointers to the x and y linear surfaces +real xmag, ymag #I/O the x and y magnification +real xrot, yrot #I/O the x and y axis rotation + +real dtgetr() + +begin + if (IS_INDEFR(xmag)) + xmag = real (dtgetr (dt, rec, "xmag")) + if (IS_INDEFR(ymag)) + ymag = real (dtgetr (dt, rec, "ymag")) + if (IS_INDEFR(xrot)) + xrot = DEGTORAD (real(dtgetr (dt, rec, "xrotation"))) + else + xrot = DEGTORAD(xrot) + if (IS_INDEFR(yrot)) + yrot = DEGTORAD (real (dtgetr (dt, rec, "yrotation"))) + else + yrot = DEGTORAD(yrot) + call geo_rotmagr (sx1, sy1, xmag, ymag, xrot, yrot) +end + + +# GEO_DXYSHIFT -- Adjust the fitted xy shift using the database file. + +procedure geo_dxyshiftr (dt, rec, sx1, sy1, xout, yout, xref, yref, + xshift, yshift) + +pointer dt #I pointer to the text file database +int rec #I the database record +pointer sx1, sy1 #I/O pointers to the x and y linear surfaces +real xout, yout #I the input coordinate system origin +real xref, yref #I the reference coordinate system origin +real xshift, yshift #I the origin shift in input coordinates + +real gsgetr(), gseval() + +begin + if (IS_INDEFR(xref)) + xref = (gsgetr (sx1, GSXMIN) + gsgetr (sx1, GSXMAX)) / 2.0 + if (IS_INDEFR(yref)) + yref = (gsgetr (sy1, GSYMIN) + gsgetr (sy1, GSYMAX)) / 2.0 + + if (IS_INDEFR(xout)) + xout = gseval (sx1, xref, yref) + if (IS_INDEFR(yout)) + yout = gseval (sy1, xref, yref) + + if (IS_INDEFR(xshift)) + xshift = xout - gseval (sx1, xref, yref) + if (IS_INDEFR(yshift)) + yshift = yout - gseval (sy1, xref, yref) + + call geo_xyshiftr (sx1, sy1, xshift, yshift) +end + + +# GEO_ROTMAG -- Edit the coefficients of the linear surface which determine +# magnification and rotation. + +procedure geo_rotmagr (sx1, sy1, xscale, yscale, xrotation, yrotation) + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +real xscale, yscale #I the x and y scales +real xrotation,yrotation #I the x and y axis rotation angles in radians + +real cosx, sinx, cosy, siny, xrange, yrange +int ncoeff +pointer sp, xcoeff, ycoeff +real gsgetr() +int gsgeti() + +begin + # Get the current solution. + call smark (sp) + ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE)) + call salloc (xcoeff, ncoeff, TY_REAL) + call salloc (ycoeff, ncoeff, TY_REAL) + call gssave (sx1, Memr[xcoeff]) + call gssave (sy1, Memr[ycoeff]) + + # Define the scaling parameters. + cosx = cos (xrotation) + sinx = sin (xrotation) + cosy = cos (yrotation) + siny = sin (yrotation) + + # Calculate coefficients. + Memr[xcoeff+GS_SAVECOEFF+1] = xscale * cosx + Memr[xcoeff+GS_SAVECOEFF+2] = yscale * siny + Memr[ycoeff+GS_SAVECOEFF+1] = -xscale * sinx + Memr[ycoeff+GS_SAVECOEFF+2] = yscale * cosy + + # Normalize coefficients for-non polynomial functions. + if (gsgeti (sx1, GSTYPE) != GS_POLYNOMIAL) { + xrange = gsgetr (sx1, GSXMAX) - gsgetr (sx1, GSXMIN) + Memr[xcoeff+GS_SAVECOEFF+1] = Memr[xcoeff+GS_SAVECOEFF+1] * + xrange / 2.d0 + Memr[xcoeff+GS_SAVECOEFF+2] = Memr[xcoeff+GS_SAVECOEFF+2] * + yrange / 2.d0 + } + if (gsgeti (sy1, GSTYPE) != GS_POLYNOMIAL) { + yrange = gsgetr (sy1, GSYMAX) - gsgetr (sy1, GSYMIN) + Memr[ycoeff+GS_SAVECOEFF+1] = Memr[ycoeff+GS_SAVECOEFF+1] * + xrange / 2.d0 + Memr[ycoeff+GS_SAVECOEFF+2] = Memr[ycoeff+GS_SAVECOEFF+2] * + yrange / 2.d0 + } + + # Free the original fit. + call gsfree (sx1) + call gsfree (sy1) + + # Restore the edited fit. + call gsrestore (sx1, Memr[xcoeff]) + call gsrestore (sy1, Memr[ycoeff]) + + call sfree (sp) +end + + +# GEO_XYSHIFT -- Shift the linear part of the fit in x and y. + +procedure geo_xyshiftr (sx1, sy1, xshift, yshift) + +pointer sx1, sy1 #I pointers to linear x and y surfaces +real xshift, yshift #I the input x and y shifts + +int ncoeff +pointer sp, xcoeff, ycoeff +int gsgeti() + +begin + call smark (sp) + + # Allocate working space. + ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE)) + call salloc (xcoeff, ncoeff, TY_REAL) + call salloc (ycoeff, ncoeff, TY_REAL) + + # Get coefficients. + call gssave (sx1, Memr[xcoeff]) + call gssave (sy1, Memr[ycoeff]) + + # Shift the coefficients. + Memr[xcoeff+GS_SAVECOEFF] = Memr[xcoeff+GS_SAVECOEFF] + xshift + Memr[ycoeff+GS_SAVECOEFF] = Memr[ycoeff+GS_SAVECOEFF] + yshift + + # Free original fit. + call gsfree (sx1) + call gsfree (sy1) + + # Restore fit. + call gsrestore (sx1, Memr[xcoeff]) + call gsrestore (sy1, Memr[ycoeff]) + + call sfree (sp) +end + + + + +# GEO_DROTMAG -- Adjust the coefficients of the fit using the database file. + +procedure geo_drotmagd (dt, rec, sx1, sy1, xmag, ymag, xrot, yrot) + +pointer dt #I pointer to the text database file +int rec #I record number +pointer sx1, sy1 #I/O pointers to the x and y linear surfaces +double xmag, ymag #I/O the x and y magnification +double xrot, yrot #I/O the x and y axis rotation + +real dtgetr() + +begin + if (IS_INDEFD(xmag)) + xmag = double (dtgetr (dt, rec, "xmag")) + if (IS_INDEFD(ymag)) + ymag = double (dtgetr (dt, rec, "ymag")) + if (IS_INDEFD(xrot)) + xrot = DEGTORAD (double(dtgetr (dt, rec, "xrotation"))) + else + xrot = DEGTORAD(xrot) + if (IS_INDEFD(yrot)) + yrot = DEGTORAD (double (dtgetr (dt, rec, "yrotation"))) + else + yrot = DEGTORAD(yrot) + call geo_rotmagd (sx1, sy1, xmag, ymag, xrot, yrot) +end + + +# GEO_DXYSHIFT -- Adjust the fitted xy shift using the database file. + +procedure geo_dxyshiftd (dt, rec, sx1, sy1, xout, yout, xref, yref, + xshift, yshift) + +pointer dt #I pointer to the text file database +int rec #I the database record +pointer sx1, sy1 #I/O pointers to the x and y linear surfaces +double xout, yout #I the input coordinate system origin +double xref, yref #I the reference coordinate system origin +double xshift, yshift #I the origin shift in input coordinates + +double dgsgetd(), dgseval() + +begin + if (IS_INDEFD(xref)) + xref = (dgsgetd (sx1, GSXMIN) + dgsgetd (sx1, GSXMAX)) / 2.0d0 + if (IS_INDEFD(yref)) + yref = (dgsgetd (sy1, GSYMIN) + dgsgetd (sy1, GSYMAX)) / 2.0d0 + + if (IS_INDEFD(xout)) + xout = dgseval (sx1, xref, yref) + if (IS_INDEFD(yout)) + yout = dgseval (sy1, xref, yref) + + if (IS_INDEFD(xshift)) + xshift = xout - dgseval (sx1, xref, yref) + if (IS_INDEFD(yshift)) + yshift = yout - dgseval (sy1, xref, yref) + + call geo_xyshiftd (sx1, sy1, xshift, yshift) +end + + +# GEO_ROTMAG -- Edit the coefficients of the linear surface which determine +# magnification and rotation. + +procedure geo_rotmagd (sx1, sy1, xscale, yscale, xrotation, yrotation) + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +double xscale, yscale #I the x and y scales +double xrotation,yrotation #I the x and y axis rotation angles in radians + +double cosx, sinx, cosy, siny, xrange, yrange +int ncoeff +pointer sp, xcoeff, ycoeff +double dgsgetd() +int dgsgeti() + +begin + # Get the current solution. + call smark (sp) + ncoeff = max (dgsgeti (sx1, GSNSAVE), dgsgeti (sy1, GSNSAVE)) + call salloc (xcoeff, ncoeff, TY_DOUBLE) + call salloc (ycoeff, ncoeff, TY_DOUBLE) + call dgssave (sx1, Memd[xcoeff]) + call dgssave (sy1, Memd[ycoeff]) + + # Define the scaling parameters. + cosx = cos (xrotation) + sinx = sin (xrotation) + cosy = cos (yrotation) + siny = sin (yrotation) + + # Calculate coefficients. + Memd[xcoeff+GS_SAVECOEFF+1] = xscale * cosx + Memd[xcoeff+GS_SAVECOEFF+2] = yscale * siny + Memd[ycoeff+GS_SAVECOEFF+1] = -xscale * sinx + Memd[ycoeff+GS_SAVECOEFF+2] = yscale * cosy + + # Normalize coefficients for-non polynomial functions. + if (dgsgeti (sx1, GSTYPE) != GS_POLYNOMIAL) { + xrange = dgsgetd (sx1, GSXMAX) - dgsgetd (sx1, GSXMIN) + Memd[xcoeff+GS_SAVECOEFF+1] = Memd[xcoeff+GS_SAVECOEFF+1] * + xrange / 2.d0 + Memd[xcoeff+GS_SAVECOEFF+2] = Memd[xcoeff+GS_SAVECOEFF+2] * + yrange / 2.d0 + } + if (dgsgeti (sy1, GSTYPE) != GS_POLYNOMIAL) { + yrange = dgsgetd (sy1, GSYMAX) - dgsgetd (sy1, GSYMIN) + Memd[ycoeff+GS_SAVECOEFF+1] = Memd[ycoeff+GS_SAVECOEFF+1] * + xrange / 2.d0 + Memd[ycoeff+GS_SAVECOEFF+2] = Memd[ycoeff+GS_SAVECOEFF+2] * + yrange / 2.d0 + } + + # Free the original fit. + call dgsfree (sx1) + call dgsfree (sy1) + + # Restore the edited fit. + call dgsrestore (sx1, Memd[xcoeff]) + call dgsrestore (sy1, Memd[ycoeff]) + + call sfree (sp) +end + + +# GEO_XYSHIFT -- Shift the linear part of the fit in x and y. + +procedure geo_xyshiftd (sx1, sy1, xshift, yshift) + +pointer sx1, sy1 #I pointers to linear x and y surfaces +double xshift, yshift #I the input x and y shifts + +int ncoeff +pointer sp, xcoeff, ycoeff +int dgsgeti() + +begin + call smark (sp) + + # Allocate working space. + ncoeff = max (dgsgeti (sx1, GSNSAVE), dgsgeti (sy1, GSNSAVE)) + call salloc (xcoeff, ncoeff, TY_DOUBLE) + call salloc (ycoeff, ncoeff, TY_DOUBLE) + + # Get coefficients. + call dgssave (sx1, Memd[xcoeff]) + call dgssave (sy1, Memd[ycoeff]) + + # Shift the coefficients. + Memd[xcoeff+GS_SAVECOEFF] = Memd[xcoeff+GS_SAVECOEFF] + xshift + Memd[ycoeff+GS_SAVECOEFF] = Memd[ycoeff+GS_SAVECOEFF] + yshift + + # Free original fit. + call dgsfree (sx1) + call dgsfree (sy1) + + # Restore fit. + call dgsrestore (sx1, Memd[xcoeff]) + call dgsrestore (sy1, Memd[ycoeff]) + + call sfree (sp) +end + + + diff --git a/pkg/images/immatch/src/geometry/geotimtran.x b/pkg/images/immatch/src/geometry/geotimtran.x new file mode 100644 index 00000000..f84a794d --- /dev/null +++ b/pkg/images/immatch/src/geometry/geotimtran.x @@ -0,0 +1,543 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <imset.h> +include <mach.h> +include <math/gsurfit.h> +include <math/iminterp.h> +include "geotran.h" + +# GEO_IMTRAN -- Correct an entire image for geometric distortion using the +# transformed coordinates and image interpolation. + +procedure geo_imtran (input, output, geo, sx1, sy1, sx2, sy2) + +pointer input #I pointer to input image +pointer output #I pointer to output image +pointer geo #I pointer to geotran structure +pointer sx1, sy1 #I pointers to linear surface descriptors +pointer sx2, sy2 #I pointer to higher order surface descriptors + +int nincr +pointer sp, xref, yref, msi +real shift +real gsgetr() + +begin + # Initialize the interpolant and compute the out-of-bounds pixel + # margin required. + if (IM_NDIM(input) == 1) { + call asitype (GT_INTERPSTR(geo), GT_INTERPOLANT(geo), + GT_NSINC(geo), nincr, shift) + call asisinit (msi, GT_INTERPOLANT(geo), GT_NSINC(geo), + nincr, shift, 0.0) + } else { + call msitype (GT_INTERPSTR(geo), GT_INTERPOLANT(geo), + GT_NSINC(geo), nincr, shift) + call msisinit (msi, GT_INTERPOLANT(geo), GT_NSINC(geo), + nincr, nincr, shift, shift, 0.0) + } + call geo_margset (sx1, sy1, sx2, sy2, GT_XMIN(geo), GT_XMAX(geo), + GT_NCOLS(geo), GT_YMIN(geo), GT_YMAX(geo), GT_NLINES(geo), + GT_INTERPOLANT(geo), GT_NSINC(geo), GT_NXYMARGIN(geo)) + + # Allocate working space. + call smark (sp) + call salloc (xref, GT_NCOLS(geo), TY_REAL) + call salloc (yref, GT_NLINES(geo), TY_REAL) + + # Calculate the reference coordinates of the input image pixels. + call geo_ref (geo, Memr[xref], 1, GT_NCOLS(geo), GT_NCOLS(geo), + Memr[yref], 1, GT_NLINES(geo), GT_NLINES(geo), gsgetr (sx1, + GSXMIN), gsgetr (sx1, GSXMAX), gsgetr (sx1, GSYMIN), gsgetr (sx1, + GSYMAX), GT_ONE) + + # Configure the out-of-bounds pixel references for the input image. + call geo_imset (input, geo, sx1, sy1, sx2, sy2, Memr[xref], + GT_NCOLS(geo), Memr[yref], GT_NLINES(geo)) + + # Interpolate. + call geo_gsvector (input, output, geo, msi, Memr[xref], 1, + GT_NCOLS(geo), Memr[yref], 1, GT_NLINES(geo), sx1, sy1, sx2, sy2) + + # Clean up. + if (IM_NDIM(input) == 1) + call asifree (msi) + else + call msifree (msi) + call sfree (sp) +end + + +# GEO_SIMTRAN -- Correct an entire image for geometric distortion using +# nterpolated coordinate surfaces to speed up computation of the transformed +# coordinates and image interpolation. + +procedure geo_simtran (input, output, geo, sx1, sy1, sx2, sy2) + +pointer input #I pointer to input image +pointer output #I pointer to output image +pointer geo #I pointer to geotran structure +pointer sx1, sy1 #I pointer to linear surface descriptors +pointer sx2, sy2 #I pointer to higher order surface descriptors + +int nxsample, nysample, nincr +pointer sp, xsample, ysample, xinterp, yinterp +pointer xmsi, ymsi, jmsi, msi, xbuf, ybuf, jbuf +real shift +real gsgetr() + +begin + # Allocate working space and intialize the interpolant. + call smark (sp) + call salloc (xsample, GT_NCOLS(geo), TY_REAL) + call salloc (ysample, GT_NLINES(geo), TY_REAL) + call salloc (xinterp, GT_NCOLS(geo), TY_REAL) + call salloc (yinterp, GT_NLINES(geo), TY_REAL) + + # Set up sampling size. + if (GT_NCOLS(geo) == 1) + nxsample = 1 + else + nxsample = GT_NCOLS(geo) / GT_XSAMPLE(geo) + if (GT_NLINES(geo) == 1) + nysample = 1 + else + nysample = GT_NLINES(geo) / GT_YSAMPLE(geo) + + # Initialize interpolants. + if (IM_NDIM(input) == 1) { + call asiinit (xmsi, II_LINEAR) + call asiinit (ymsi, II_LINEAR) + call asitype (GT_INTERPSTR(geo), GT_INTERPOLANT(geo), + GT_NSINC(geo), nincr, shift) + call asisinit (msi, GT_INTERPOLANT(geo), GT_NSINC(geo), + nincr, shift, 0.0) + if (GT_FLUXCONSERVE(geo) == YES) + call asiinit (jmsi, II_LINEAR) + } else { + call msiinit (xmsi, II_BILINEAR) + call msiinit (ymsi, II_BILINEAR) + call msitype (GT_INTERPSTR(geo), GT_INTERPOLANT(geo), + GT_NSINC(geo), nincr, shift) + call msisinit (msi, GT_INTERPOLANT(geo), GT_NSINC(geo), + nincr, nincr, shift, shift, 0.0) + if (GT_FLUXCONSERVE(geo) == YES) + call msiinit (jmsi, II_BILINEAR) + } + call geo_margset (sx1, sy1, sx2, sy2, GT_XMIN(geo), GT_XMAX(geo), + GT_NCOLS(geo), GT_YMIN(geo), GT_YMAX(geo), GT_NLINES(geo), + GT_INTERPOLANT(geo), GT_NSINC(geo), GT_NXYMARGIN(geo)) + + # Setup input image boundary extension parameters. + call geo_ref (geo, Memr[xsample], 1, GT_NCOLS(geo), GT_NCOLS(geo), + Memr[ysample], 1, GT_NLINES(geo), GT_NLINES(geo), gsgetr (sx1, + GSXMIN), gsgetr (sx1, GSXMAX), gsgetr (sx1, GSYMIN), gsgetr (sx1, + GSYMAX), GT_ONE) + call geo_imset (input, geo, sx1, sy1, sx2, sy2, Memr[xsample], + GT_NCOLS(geo), Memr[ysample], GT_NLINES(geo)) + + # Calculate the sampled reference coordinates and the interpolated + # reference coordinates. + call geo_ref (geo, Memr[xsample], 1, nxsample, nxsample, Memr[ysample], + 1, nysample, nysample, gsgetr (sx1, GSXMIN), gsgetr (sx1, GSXMAX), + gsgetr (sx1, GSYMIN), gsgetr (sx1, GSYMAX), GT_ONE) + call geo_sample (geo, Memr[xinterp], 1, GT_NCOLS(geo), nxsample, + Memr[yinterp], 1, GT_NLINES(geo), nysample, GT_ONE) + + # Initialize the buffers + xbuf = NULL + ybuf = NULL + jbuf = NULL + + # Set up interpolants + call geo_xbuffer (sx1, sx2, xmsi, Memr[xsample], Memr[ysample], 1, + nxsample, 1, nysample, xbuf) + call geo_ybuffer (sy1, sy2, ymsi, Memr[xsample], Memr[ysample], 1, + nxsample, 1, nysample, ybuf) + if (GT_FLUXCONSERVE(geo) == YES && (sx2 != NULL || sy2 != NULL)) { + if (IM_NDIM(input) == 1) + call geo_jbuffer (sx1, NULL, sx2, NULL, jmsi, Memr[xsample], + Memr[ysample], 1, nxsample, 1, nysample, jbuf) + else + call geo_jbuffer (sx1, sy1, sx2, sy2, jmsi, Memr[xsample], + Memr[ysample], 1, nxsample, 1, nysample, jbuf) + } + + # Transform the image. + call geo_msivector (input, output, geo, xmsi, ymsi, jmsi, msi, + sx1, sy1, sx2, sy2, Memr[xinterp], 1, GT_NCOLS(geo), nxsample, + Memr[yinterp], 1, GT_NLINES(geo), nysample, 1, 1) + + # Free space. + if (IM_NDIM(input) == 1) { + call asifree (xmsi) + call asifree (ymsi) + call asifree (msi) + if (GT_FLUXCONSERVE(geo) == YES) + call asifree (jmsi) + } else { + call msifree (xmsi) + call msifree (ymsi) + call msifree (msi) + if (GT_FLUXCONSERVE(geo) == YES) + call msifree (jmsi) + } + call mfree (xbuf, TY_REAL) + call mfree (ybuf, TY_REAL) + if (jbuf != NULL) + call mfree (jbuf, TY_REAL) + call sfree (sp) +end + + +## GEO_IMSIVECTOR -- Evaluate the output image using interpolated surface +## coordinates. +# +#procedure geo_imsivector (in, out, geo, xmsi, ymsi, jmsi, msi, sx1, sy1, sx2, +# sy2, xref, yref, ncols, nlines) +# +#pointer in #I pointer to input image +#pointer out #I pointer to output image +#pointer geo #I pointer to geotran structure +#pointer xmsi, ymsi #I pointer to the interpolation xy surfaces +#pointer jmsi #I pointer to Jacobian surface +#pointer msi #I pointer to interpolation surface +#pointer sx1, sy1 #I linear surface descriptors +#pointer sx2, sy2 #I distortion surface pointers +#real xref[ARB] #I x reference coordinates +#real yref[ARB] #I y reference coordinates +#int ncols, nlines #I number of columns and rows +# +#int j +#pointer sp, x, y, xin, yin, xout, yout, inbuf, outbuf +#real factor +#pointer imgs1r(), imgs2r(), imps1r(), imps2r() +#real geo_jfactor() +# +#begin +# # Allocate working space. +# call smark (sp) +# call salloc (x, ncols, TY_REAL) +# call salloc (y, ncols, TY_REAL) +# call salloc (xin, ncols, TY_REAL) +# call salloc (yin, ncols, TY_REAL) +# call salloc (xout, ncols, TY_REAL) +# call salloc (yout, ncols, TY_REAL) +# +# # Fit the interpolant +# if (IM_NDIM(in) == 1) +# inbuf = imgs1r (in, 1, int (IM_LEN(in,1))) +# else +# inbuf = imgs2r (in, 1, int (IM_LEN(in,1)), 1, int (IM_LEN(in,2))) +# if (inbuf == EOF) +# call error (0, "Error reading image") +# if (IM_NDIM(in) == 1) +# call asifit (msi, Memr[inbuf], int (IM_LEN(in,1))) +# else +# call msifit (msi, Memr[inbuf], int (IM_LEN(in,1)), +# int (IM_LEN(in,2)), int (IM_LEN(in,1))) +# +# # Compute the output bufferr. +# do j = 1, nlines { +# +# # Compute coordinates. +# call amovkr (yref[j], Memr[y], ncols) +# if (IM_NDIM(in) == 1) { +# call asivector (xmsi, xref, Memr[xin], ncols) +# call asivector (ymsi, xref, Memr[yin], ncols) +# } else { +# call msivector (xmsi, xref, Memr[y], Memr[xin], ncols) +# call msivector (ymsi, xref, Memr[y], Memr[yin], ncols) +# } +# +# # Correct for out-of-bounds pixels. +# call geo_btran (in, geo, Memr[xin], Memr[yin], Memr[xout], +# Memr[yout], ncols) +# +# # Write to output image. +# if (IM_NDIM(in) == 1) +# outbuf = imps1r (out, 1, ncols) +# else +# outbuf = imps2r (out, 1, ncols, j, j) +# if (outbuf == EOF) +# call error (0, "Error writing output image") +# if (IM_NDIM(in) == 1) +# call asivector (msi, Memr[xout], Memr[outbuf], ncols) +# else +# call msivector (msi, Memr[xout], Memr[yout], Memr[outbuf], +# ncols) +# +# # Perform constant boundary extension. +# if (GT_BOUNDARY(geo) == BT_CONSTANT) +# call geo_bconstant (in, geo, Memr[xin], Memr[yin], +# Memr[outbuf], Memr[outbuf], ncols) +# +# # Preserve flux in image. +# if (GT_FLUXCONSERVE(geo) == YES) { +# factor = GT_XSCALE(geo) * GT_YSCALE(geo) +# if (GT_GEOMODE(geo) == GT_LINEAR || (sx2 == NULL && sy2 == +# NULL)) { +# if (IM_NDIM(in) == 1) +# call amulkr (Memr[outbuf], factor * geo_jfactor (sx1, +# NULL), Memr[outbuf], ncols) +# else +# call amulkr (Memr[outbuf], factor * geo_jfactor (sx1, +# sy1), Memr[outbuf], ncols) +# } else { +# if (IM_NDIM(in) == 1) +# call geo_msiflux (jmsi, xref, yref, Memr[outbuf], +# 1, ncols, 0, 1, 1) +# else +# call geo_msiflux (jmsi, xref, yref, Memr[outbuf], +# 1, ncols, j, 1, 1) +# call amulkr (Memr[outbuf], factor, Memr[outbuf], ncols) +# } +# } +# } +# +# call sfree (sp) +#end + + +## GEO_IGSVECTOR -- Evaluate the output image using fitted coordinates. +# +#procedure geo_igsvector (input, output, geo, msi, xref, yref, ncols, nlines, +# sx1, sy1, sx2, sy2) +# +#pointer input #I pointer to input image +#pointer output #I pointer to output image +#pointer geo #I pointer to geotran structure +#pointer msi #I pointer to interpolant +#real xref[ARB] #I x reference array +#real yref[ARB] #I y reference array +#int ncols, nlines #I number of columns and lines +#pointer sx1, sy1 #I pointer to linear surface +#pointer sx2, sy2 #I pointer to distortion surface +# +#int j +#pointer sp, y, xin, yin, xout, yout, temp, inbuf, outbuf +#real factor +#pointer imgs1r(), imgs2r(), imps1r(), imps2r() +#real geo_jfactor() +# +#begin +# # Allocate working space. +# call smark (sp) +# call salloc (y, ncols, TY_REAL) +# call salloc (xin, ncols, TY_REAL) +# call salloc (yin, ncols, TY_REAL) +# call salloc (xout, ncols, TY_REAL) +# call salloc (yout, ncols, TY_REAL) +# call salloc (temp, ncols, TY_REAL) +# +# # Fill image buffer. +# if (IM_NDIM(input) == 1) +# inbuf = imgs1r (input, 1, int (IM_LEN(input,1))) +# else +# inbuf = imgs2r (input, 1, int (IM_LEN(input,1)), 1, +# int (IM_LEN(input,2))) +# if (inbuf == EOF) +# call error (0, "Error reading image") +# +# # Fit the interpolant. +# if (IM_NDIM(input) == 1) +# call asifit (msi, Memr[inbuf], int (IM_LEN(input,1))) +# else +# call msifit (msi, Memr[inbuf], int (IM_LEN(input,1)), +# int (IM_LEN(input,2)), int (IM_LEN(input,1))) +# +# # Calculate the x and y input image coordinates. +# do j = 1, nlines { +# +# # Get output image buffer. +# if (IM_NDIM(input) == 1) +# outbuf = imps1r (output, 1, ncols) +# else +# outbuf = imps2r (output, 1, ncols, j, j) +# if (output == EOF) +# call error (0, "Error writing output image") +# +# # Fit x coords. +# call amovkr (yref[j], Memr[y], ncols) +# call gsvector (sx1, xref, Memr[y], Memr[xin], ncols) +# if (sx2 != NULL) { +# call gsvector (sx2, xref, Memr[y], Memr[temp], ncols) +# call aaddr (Memr[xin], Memr[temp], Memr[xin], ncols) +# } +# +# # Fit y coords. +# call gsvector (sy1, xref, Memr[y], Memr[yin], ncols) +# if (sy2 != NULL) { +# call gsvector (sy2, xref, Memr[y], Memr[temp], ncols) +# call aaddr (Memr[yin], Memr[temp], Memr[yin], ncols) +# } +# +# # Compute of of bounds pixels. +# call geo_btran (input, geo, Memr[xin], Memr[yin], Memr[xout], +# Memr[yout], ncols) +# +# # Interpolate in input image. +# if (IM_NDIM(input) == 1) +# call asivector (msi, Memr[xout], Memr[outbuf], ncols) +# else +# call msivector (msi, Memr[xout], Memr[yout], Memr[outbuf], +# ncols) +# +# # Correct for constant boundary extension. +# if (GT_BOUNDARY(geo) == BT_CONSTANT) +# call geo_bconstant (input, geo, Memr[xin], Memr[yin], +# Memr[outbuf], Memr[outbuf], ncols) +# +# # Preserve flux in image. +# if (GT_FLUXCONSERVE(geo) == YES) { +# factor = GT_XSCALE(geo) * GT_YSCALE(geo) +# if (GT_GEOMODE(geo) == GT_LINEAR || (sx2 == NULL && sy2 == +# NULL)) { +# if (IM_NDIM(input) == 1) +# call amulkr (Memr[outbuf], factor * geo_jfactor (sx1, +# NULL), Memr[outbuf], ncols) +# else +# call amulkr (Memr[outbuf], factor * geo_jfactor (sx1, +# sy1), Memr[outbuf], ncols) +# } else { +# if (IM_NDIM(input) == 1) +# call geo_gsflux (xref, yref, Memr[outbuf], 1, ncols, j, +# sx1, NULL, sx2, NULL) +# else +# call geo_gsflux (xref, yref, Memr[outbuf], 1, ncols, j, +# sx1, sy1, sx2, sy2) +# call amulkr (Memr[outbuf], factor, Memr[outbuf], ncols) +# } +# } +# } +# +# call sfree (sp) +#end + + +## GEO_BTRAN -- Map out-of-bounds pixel into the input image. +# +#procedure geo_btran (input, geo, xin, yin, xout, yout, ncols) +# +#pointer input #I pointer to the input image +#pointer geo #I pointer to geotran strcuture +#real xin[ARB] #I x input coords +#real yin[ARB] #I y input coords +#real xout[ARB] #O x output coords +#real yout[ARB] #O y output coords +#int ncols #I number of columns +# +#int i +#real xmax, ymax, xtemp, ytemp +# +#begin +# xmax = IM_LEN(input,1) +# if (IM_NDIM(input) == 1) +# ymax = 1.0 +# else +# ymax = IM_LEN(input,2) +# +# switch (GT_BOUNDARY(geo)) { +# case BT_CONSTANT, BT_NEAREST: +# do i = 1, ncols { +# if (xin[i] < 1.0) +# xout[i] = 1.0 +# else if (xin[i] > xmax) +# xout[i] = xmax +# else +# xout[i] = xin[i] +# if (yin[i] < 1.0) +# yout[i] = 1.0 +# else if (yin[i] > ymax) +# yout[i] = ymax +# else +# yout[i] = yin[i] +# } +# case BT_REFLECT: +# do i = 1, ncols { +# if (xin[i] < 1.0) +# xout[i] = 1.0 + (1.0 - xin[i]) +# else if (xin[i] > xmax) +# xout[i] = xmax - (xin[i] - xmax) +# else +# xout[i] = xin[i] +# if (yin[i] < 1.0) +# yout[i] = 1.0 + (1.0 - yin[i]) +# else if (yin[i] > ymax) +# yout[i] = ymax - (yin[i] - ymax) +# else +# yout[i] = yin[i] +# } +# case BT_WRAP: +# do i = 1, ncols { +# xtemp = xin[i] +# ytemp = yin[i] +# +# if (xtemp < 1.0) { +# while (xtemp < 1.0) +# xtemp = xtemp + xmax +# if (xtemp < 1.0) +# xtemp = xmax - xtemp +# else if (xtemp > xmax) +# xtemp = 2.0 + xmax - xtemp +# } else if (xtemp > xmax) { +# while (xtemp > xmax) +# xtemp = xtemp - xmax +# if (xtemp < 1.0) +# xtemp = xmax - xtemp +# else if (xtemp > xmax) +# xtemp = 2.0 + xmax - xtemp +# } +# xout[i] = xtemp +# +# if (ytemp < 1.0) { +# while (ytemp < 1.0) +# ytemp = ytemp + ymax +# if (ytemp < 1.0) +# ytemp = ymax - ytemp +# else if (ytemp > ymax) +# ytemp = 2.0 + ymax - ytemp +# } else if (ytemp > ymax) { +# while (ytemp > ymax) +# ytemp = ytemp - ymax +# if (ytemp < 1.0) +# ytemp = ymax - ytemp +# else if (ytemp > ymax) +# ytemp = 2.0 + ymax - ytemp +# } +# yout[i] = ytemp +# } +# } +#end + + +## GEO_BCONSTANT -- Map constant out-of-bounds pixels into the input image. +# +#procedure geo_bconstant (input, geo, xin, yin, inbuf, outbuf, ncols) +# +#pointer input #I pointer to the input image +#pointer geo #I pointer to geotran structure +#real xin[ARB] #I x input coords +#real yin[ARB] #I y input coords +#real inbuf[ARB] #I input buffer +#real outbuf[ARB] #O output buffer +#int ncols #I number of columns +# +#int i +#real xmax, ymax, constant +# +#begin +# xmax = IM_LEN(input,1) +# if (IM_NDIM(input) == 1) +# ymax = 1.0 +# else +# ymax = IM_LEN(input,2) +# constant = GT_CONSTANT(geo) +# do i = 1, ncols { +# if (xin[i] < 1.0 || xin[i] > xmax || yin[i] < 1.0 || yin[i] > ymax) +# outbuf[i] = constant +# else +# outbuf[i] = inbuf[i] +# } +#end diff --git a/pkg/images/immatch/src/geometry/geotran.h b/pkg/images/immatch/src/geometry/geotran.h new file mode 100644 index 00000000..d2fa6b55 --- /dev/null +++ b/pkg/images/immatch/src/geometry/geotran.h @@ -0,0 +1,52 @@ +# GEOTRAN Structure + +define LEN_GEOSTRUCT (30 + SZ_FNAME) + +# output picture formatting parameters + +define GT_NCOLS Memi[$1] # number of output columns +define GT_NLINES Memi[$1+1] # number of output lines +define GT_XMIN Memr[P2R($1+2)] # x minimum +define GT_XMAX Memr[P2R($1+3)] # x maximum +define GT_YMIN Memr[P2R($1+4)] # y minimun +define GT_YMAX Memr[P2R($1+5)] # y maximum +define GT_XSCALE Memr[P2R($1+6)] # x scale +define GT_YSCALE Memr[P2R($1+7)] # y scale + +# transformation parameters + +define GT_GEOMODE Memi[$1+8] # Type of transformation +define GT_XIN Memr[P2R($1+9)] # x input pixel +define GT_YIN Memr[P2R($1+10)] # y input pixel +define GT_XOUT Memr[P2R($1+11)] # x output pixel +define GT_YOUT Memr[P2R($1+12)] # y output pixel +define GT_XSHIFT Memr[P2R($1+13)] # x shift +define GT_YSHIFT Memr[P2R($1+14)] # y shift +define GT_XMAG Memr[P2R($1+15)] # input image x scale +define GT_YMAG Memr[P2R($1+16)] # input image y scale +define GT_XROTATION Memr[P2R($1+17)] # rotation angle +define GT_YROTATION Memr[P2R($1+18)] # scale angle + +# interpolation parameters +define GT_XSAMPLE Memr[P2R($1+19)] # x surface subsampling +define GT_YSAMPLE Memr[P2R($1+20)] # y surface subsampling +define GT_INTERPOLANT Memi[$1+21] # image interpolant +define GT_NSINC Memi[$1+22] # sinc width half-width +define GT_NXYMARGIN Memi[$1+23] # the interpolation margin +define GT_BOUNDARY Memi[$1+24] # boundary extension +define GT_CONSTANT Memr[P2R($1+25)] # constant boundary extension +define GT_FLUXCONSERVE Memi[$1+26] # conserve total flux +define GT_INTERPSTR Memc[P2C($1+27)] # interpolation string + +# GEOTRAN MODES + +define GT_NONE 1 # parameters defined by user +define GT_LINEAR 2 # use linear transformation +define GT_DISTORT 3 # distortion transformation only +define GT_GEOMETRIC 4 # use full transformation + +# GEOTRAN COORDINATE MODES + +define GT_ONE 1 +define GT_TWO 2 +define GT_FOUR 3 diff --git a/pkg/images/immatch/src/geometry/geotran.x b/pkg/images/immatch/src/geometry/geotran.x new file mode 100644 index 00000000..ee689d26 --- /dev/null +++ b/pkg/images/immatch/src/geometry/geotran.x @@ -0,0 +1,1752 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <imset.h> +include <mach.h> +include <math/gsurfit.h> +include <math/iminterp.h> +include "geotran.h" + +define NMARGIN 3 # number of boundary pixels +define NMARGIN_SPLINE3 16 # number of spline boundary pixels + +# GEO_TRAN -- Correct an image for geometric distortion block by block using +# fitted coordinates and image interpolation. + +procedure geo_tran (input, output, geo, sx1, sy1, sx2, sy2, nxblock, nyblock) + +pointer input #I pointer to input image +pointer output #I pointer to output image +pointer geo #I pointer to geotran structure +pointer sx1, sy1 #I pointers to linear surfaces +pointer sx2, sy2 #I pointers to higher order surfaces +int nxblock, nyblock #I working block size + +int l1, l2, c1, c2, nincr +pointer sp, xref, yref, msi +real shift +real gsgetr() + +begin + # Initialize the interpolant. + if (IM_NDIM(input) == 1) { + call asitype (GT_INTERPSTR(geo), GT_INTERPOLANT(geo), GT_NSINC(geo), + nincr, shift) + call asisinit (msi, GT_INTERPOLANT(geo), GT_NSINC(geo), nincr, + shift, 0.0) + } else { + call msitype (GT_INTERPSTR(geo), GT_INTERPOLANT(geo), + GT_NSINC(geo), nincr, shift) + call msisinit (msi, GT_INTERPOLANT(geo), GT_NSINC(geo), nincr, + nincr, shift, shift, 0.0) + } + call geo_margset (sx1, sy1, sx2, sy2, GT_XMIN(geo), GT_XMAX(geo), + GT_NCOLS(geo), GT_YMIN(geo), GT_YMAX(geo), GT_NLINES(geo), + GT_INTERPOLANT(geo), GT_NSINC(geo), GT_NXYMARGIN(geo)) + + # Allocate working space. + call smark (sp) + call salloc (xref, GT_NCOLS(geo), TY_REAL) + call salloc (yref, GT_NLINES(geo), TY_REAL) + + # Compute the reference coordinates corresponding to the center of + # the output image pixels. + call geo_ref (geo, Memr[xref], 1, GT_NCOLS(geo), GT_NCOLS(geo), + Memr[yref], 1, GT_NLINES(geo), GT_NLINES(geo), gsgetr (sx1, + GSXMIN), gsgetr (sx1, GSXMAX), gsgetr (sx1, GSYMIN), gsgetr (sx1, + GSYMAX), GT_ONE) + + # Configure the out-of-bounds pixel references for the input image. + call geo_imset (input, geo, sx1, sy1, sx2, sy2, Memr[xref], + GT_NCOLS(geo), Memr[yref], GT_NLINES(geo)) + + # Loop over the line blocks. + for (l1 = 1; l1 <= GT_NLINES(geo); l1 = l1 + nyblock) { + + # Set line limits in the output image. + l2 = min (l1 + nyblock - 1, GT_NLINES(geo)) + + # Loop over the column blocks + for (c1 = 1; c1 <= GT_NCOLS(geo); c1 = c1 + nxblock) { + + # Set column limits in the output image. + c2 = min (c1 + nxblock - 1, GT_NCOLS(geo)) + + # Interpolate + call geo_gsvector (input, output, geo, msi, Memr[xref], + c1, c2, Memr[yref], l1, l2, sx1, sy1, sx2, sy2) + } + } + + # Clean up. + if (IM_NDIM(input) == 1) + call asifree (msi) + else + call msifree (msi) + call sfree (sp) +end + + +# GEO_STRAN -- Correct an image for geometric distortion block by block using +# interpolated coordinates and image interpolation. + +procedure geo_stran (input, output, geo, sx1, sy1, sx2, sy2, nxblock, nyblock) + +pointer input #I pointer to input image +pointer output #I pointer to output image +pointer geo #I pointer to geotran structure +pointer sx1, sy1 #I pointers to linear surfaces +pointer sx2, sy2 #I pointers to higher order surfaces +int nxblock, nyblock #I working block size + +int nxsample, nysample, ncols, nlines, l1, l2, c1, c2 +int line1, line2, llast1, llast2, nincr +pointer sp, xsample, ysample, xinterp, yinterp +pointer xmsi, ymsi, jmsi, msi, xbuf, ybuf, jbuf +real shift +real gsgetr() + +begin + # Allocate working space and intialize the interpolant. + call smark (sp) + call salloc (xsample, GT_NCOLS(geo), TY_REAL) + call salloc (ysample, GT_NLINES(geo), TY_REAL) + call salloc (xinterp, GT_NCOLS(geo), TY_REAL) + call salloc (yinterp, GT_NLINES(geo), TY_REAL) + + # Compute the sample size. + if (GT_NCOLS(geo) == 1) + nxsample = 1 + else + nxsample = GT_NCOLS(geo) / GT_XSAMPLE(geo) + if (GT_NLINES(geo) == 1) + nysample = 1 + else + nysample = GT_NLINES(geo) / GT_YSAMPLE(geo) + + # Initialize interpolants. + if (IM_NDIM(input) == 1) { + call asiinit (xmsi, II_LINEAR) + call asiinit (ymsi, II_LINEAR) + call asitype (GT_INTERPSTR(geo), GT_INTERPOLANT(geo), + GT_NSINC(geo), nincr, shift) + call asisinit (msi, GT_INTERPOLANT(geo), GT_NSINC(geo), nincr, + shift, 0.0) + if (GT_FLUXCONSERVE(geo) == YES) + call asiinit (jmsi, II_LINEAR) + } else { + call msiinit (xmsi, II_BILINEAR) + call msiinit (ymsi, II_BILINEAR) + call msitype (GT_INTERPSTR(geo), GT_INTERPOLANT(geo), + GT_NSINC(geo), nincr, shift) + call msisinit (msi, GT_INTERPOLANT(geo), GT_NSINC(geo), nincr, + nincr, shift, shift, 0.0) + if (GT_FLUXCONSERVE(geo) == YES) + call msiinit (jmsi, II_BILINEAR) + } + call geo_margset (sx1, sy1, sx2, sy2, GT_XMIN(geo), GT_XMAX(geo), + GT_NCOLS(geo), GT_YMIN(geo), GT_YMAX(geo), GT_NLINES(geo), + GT_INTERPOLANT(geo), GT_NSINC(geo), GT_NXYMARGIN(geo)) + + # Setup input image boundary extension parameters. + call geo_ref (geo, Memr[xsample], 1, GT_NCOLS(geo), GT_NCOLS(geo), + Memr[ysample], 1, GT_NLINES(geo), GT_NLINES(geo), gsgetr (sx1, + GSXMIN), gsgetr (sx1, GSXMAX), gsgetr (sx1, GSYMIN), gsgetr (sx1, + GSYMAX), GT_ONE) + call geo_imset (input, geo, sx1, sy1, sx2, sy2, Memr[xsample], + GT_NCOLS(geo), Memr[ysample], GT_NLINES(geo)) + + # Calculate the sampled reference coordinates and the interpolated + # reference coordinates. + call geo_ref (geo, Memr[xsample], 1, nxsample, nxsample, Memr[ysample], + 1, nysample, nysample, gsgetr (sx1, GSXMIN), gsgetr (sx1, GSXMAX), + gsgetr (sx1, GSYMIN), gsgetr (sx1, GSYMAX), GT_ONE) + call geo_sample (geo, Memr[xinterp], 1, GT_NCOLS(geo), nxsample, + Memr[yinterp], 1, GT_NLINES(geo), nysample, GT_ONE) + + # Initialize the buffers. + xbuf = NULL + ybuf = NULL + jbuf = NULL + + # Loop over the line blocks. + for (l1 = 1; l1 <= GT_NLINES(geo); l1 = l1 + nyblock) { + + # Set line limits in the output image. + l2 = min (l1 + nyblock - 1, GT_NLINES(geo)) + nlines = l2 - l1 + 1 + + # Line1 and line2 are the coordinates in the interpolation surface + line1 = max (1, min (nysample - 1, int (Memr[yinterp+l1-1]))) + line2 = min (nysample, int (Memr[yinterp+l2-1] + 1.0)) + + if ((xbuf == NULL) || (ybuf == NULL) || (jbuf == NULL) || + (line1 < llast1) || (line2 > llast2)) { + call geo_xbuffer (sx1, sx2, xmsi, Memr[xsample], Memr[ysample], + 1, nxsample, line1, line2, xbuf) + call geo_ybuffer (sy1, sy2, ymsi, Memr[xsample], Memr[ysample], + 1, nxsample, line1, line2, ybuf) + if (GT_FLUXCONSERVE(geo) == YES) { + if (IM_NDIM(input) == 1) + call geo_jbuffer (sx1, NULL, sx2, NULL, jmsi, + Memr[xsample], Memr[ysample], 1, nxsample, + line1, line2, jbuf) + else + call geo_jbuffer (sx1, sy1, sx2, sy2, jmsi, + Memr[xsample], Memr[ysample], 1, nxsample, + line1, line2, jbuf) + } + llast1 = line1 + llast2 = line2 + } + + + # Loop over the column blocks. + for (c1 = 1; c1 <= GT_NCOLS(geo); c1 = c1 + nxblock) { + + # C1 and c2 are the column limits in the output image. + c2 = min (c1 + nxblock - 1, GT_NCOLS(geo)) + ncols = c2 - c1 + 1 + + # Calculate the coordinates of the output pixels in the input + # image. + call geo_msivector (input, output, geo, xmsi, ymsi, jmsi, msi, + sx1, sy1, sx2, sy2, Memr[xinterp], c1, c2, nxsample, + Memr[yinterp], l1, l2, nysample, 1, line1) + } + } + + # Free space. + if (IM_NDIM(input) == 1) { + call asifree (xmsi) + call asifree (ymsi) + call asifree (msi) + if (GT_FLUXCONSERVE(geo) == YES) + call asifree (jmsi) + } else { + call msifree (xmsi) + call msifree (ymsi) + call msifree (msi) + if (GT_FLUXCONSERVE(geo) == YES) + call msifree (jmsi) + } + call mfree (xbuf, TY_REAL) + call mfree (ybuf, TY_REAL) + if (GT_FLUXCONSERVE(geo) == YES) + call mfree (jbuf, TY_REAL) + call sfree (sp) +end + + +# GEO_REF -- Determine the x and y coordinates at which the coordinate +# surface will be subsampled. + +procedure geo_ref (geo, x, c1, c2, nx, y, l1, l2, ny, xmin, xmax, ymin, ymax, + cmode) + +pointer geo #I pointer to the geotran structure +real x[ARB] #O output x sample coordinates +int c1, c2, nx #I the column limits of the sampled array +real y[ARB] #O output y sample coordinates +int l1, l2, ny #I the line limits of the output coordinates +real xmin, xmax #I limits on x coordinates +real ymin, ymax #I limits on y coordinates +int cmode #I coordinate computation mode + +int i +real xtempmin, xtempmax, ytempmin, ytempmax, dx, dy + +begin + + switch (cmode) { + case GT_FOUR: + if (nx == 1) { + xtempmin = min (xmax, max (xmin, GT_XMIN(geo))) + xtempmax = min (xmax, max (xmin, GT_XMAX(geo))) + x[1] = xtempmin + x[2] = xtempmax + x[3] = xtempmax + x[4] = xtempmin + } else if (nx == GT_NCOLS(geo)) { + if (GT_XMIN(geo) > GT_XMAX(geo)) + dx = -GT_XSCALE(geo) + else + dx = GT_XSCALE(geo) + do i = c1, c2 { + xtempmin = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 1.5) * dx)) + xtempmax = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 0.5) * dx)) + x[4*(i-c1)+1] = xtempmin + x[4*(i-c1)+2] = xtempmax + x[4*(i-c1)+3] = xtempmax + x[4*(i-c1)+4] = xtempmin + } + } else { + if (GT_XMIN(geo) > GT_XMAX(geo)) + dx = -GT_XSCALE(geo) * (GT_NCOLS(geo) - 1.0) / (nx - 1.0) + else + dx = GT_XSCALE(geo) * (GT_NCOLS(geo) - 1.0) / (nx - 1.0) + do i = c1, c2 { + xtempmin = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 1.5) * dx)) + xtempmax = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 0.5) * dx)) + x[4*(i-c1)+1] = xtempmin + x[4*(i-c1)+2] = xtempmax + x[4*(i-c1)+3] = xtempmax + x[4*(i-c1)+4] = xtempmin + } + } + + case GT_TWO: + if (nx == 1) { + xtempmin = min (xmax, max (xmin, GT_XMIN(geo))) + xtempmax = min (xmax, max (xmin, GT_XMAX(geo))) + x[1] = xtempmin + x[2] = xtempmax + } else if (nx == GT_NCOLS(geo)) { + if (GT_XMIN(geo) > GT_XMAX(geo)) + dx = -GT_XSCALE(geo) + else + dx = GT_XSCALE(geo) + do i = c1, c2 { + xtempmin = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 1.5) * dx)) + xtempmax = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 0.5) * dx)) + x[2*(i-c1)+1] = xtempmin + x[2*(i-c1)+2] = xtempmax + } + } else { + if (GT_XMIN(geo) > GT_XMAX(geo)) + dx = -GT_XSCALE(geo) * (GT_NCOLS(geo) - 1.0) / (nx - 1.0) + else + dx = GT_XSCALE(geo) * (GT_NCOLS(geo) - 1.0) / (nx - 1.0) + do i = c1, c2 { + xtempmin = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 1.5) * dx)) + xtempmax = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 0.5) * dx)) + x[2*(i-c1)+1] = xtempmin + x[2*(i-c1)+2] = xtempmax + } + } + + case GT_ONE: + if (nx == 1) { + x[1] = min (xmax, max (xmin, + (GT_XMIN(geo) + GT_XMAX(geo)) / 2.0)) + } else if (nx == GT_NCOLS(geo)) { + if (GT_XMIN(geo) > GT_XMAX(geo)) + dx = -GT_XSCALE(geo) + else + dx = GT_XSCALE(geo) + do i = c1, c2 + x[i-c1+1] = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 1) * dx)) + } else { + if (GT_XMIN(geo) > GT_XMAX(geo)) + dx = -GT_XSCALE(geo) * (GT_NCOLS(geo) - 1.0) / (nx - 1.0) + else + dx = GT_XSCALE(geo) * (GT_NCOLS(geo) - 1.0) / (nx - 1.0) + do i = c1, c2 + x[i-c1+1] = min (xmax, max (xmin, GT_XMIN(geo) + + (i - 1) * dx)) + } + + } + + switch (cmode) { + case GT_FOUR: + if (ny == 1) { + ytempmin = min (ymax, max (ymin, GT_YMIN(geo))) + ytempmax = min (ymax, max (ymin, GT_YMAX(geo))) + y[1] = ytempmin + y[2] = ytempmin + y[3] = ytempmax + y[4] = ytempmax + } else if (ny == GT_NLINES(geo)) { + if (GT_YMIN(geo) > GT_YMAX(geo)) + dy = -GT_YSCALE(geo) + else + dy = GT_YSCALE(geo) + do i = l1, l2 { + ytempmin = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 1.5) * dy)) + ytempmax = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 0.5) * dy)) + y[4*(i-l1)+1] = ytempmin + y[4*(i-l1)+2] = ytempmin + y[4*(i-l1)+3] = ytempmax + y[4*(i-l1)+4] = ytempmax + } + } else { + if (GT_YMIN(geo) > GT_YMAX(geo)) + dy = -GT_YSCALE(geo) * (GT_NLINES(geo) - 1.0) / (ny - 1.0) + else + dy = GT_YSCALE(geo) * (GT_NLINES(geo) - 1.0) / (ny - 1.0) + do i = l1, l2 { + ytempmin = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 1.5) * dy)) + ytempmax = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 0.5) * dy)) + y[4*(i-l1)+1] = ytempmin + y[4*(i-l1)+2] = ytempmin + y[4*(i-l1)+3] = ytempmax + y[4*(i-l1)+4] = ytempmax + } + } + + case GT_TWO: + if (ny == 1) { + ytempmin = min (ymax, max (ymin, GT_YMIN(geo))) + ytempmax = min (ymax, max (ymin, GT_YMAX(geo))) + y[1] = ytempmin + y[2] = ytempmax + } else if (ny == GT_NLINES(geo)) { + if (GT_YMIN(geo) > GT_YMAX(geo)) + dy = -GT_YSCALE(geo) + else + dy = GT_YSCALE(geo) + do i = l1, l2 { + ytempmin = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 1.5) * dy)) + ytempmax = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 0.5) * dy)) + y[2*(i-l1)+1] = ytempmin + y[2*(i-l1)+2] = ytempmax + } + } else { + if (GT_YMIN(geo) > GT_YMAX(geo)) + dy = -GT_YSCALE(geo) * (GT_NLINES(geo) - 1.0) / (ny - 1.0) + else + dy = GT_YSCALE(geo) * (GT_NLINES(geo) - 1.0) / (ny - 1.0) + do i = l1, l2 { + ytempmin = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 1.5) * dy)) + ytempmax = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 0.5) * dy)) + y[2*(i-l1)+1] = ytempmin + y[2*(i-l1)+2] = ytempmax + } + } + case GT_ONE: + if (ny == 1) { + y[1] = min (ymax, max (ymin, + (GT_YMIN(geo) + GT_YMAX(geo)) / 2.0)) + } else if (ny == GT_NLINES(geo)) { + if (GT_YMIN(geo) > GT_YMAX(geo)) + dy = -GT_YSCALE(geo) + else + dy = GT_YSCALE(geo) + do i = l1, l2 + y[i-l1+1] = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 1) * dy)) + } else { + if (GT_YMIN(geo) > GT_YMAX(geo)) + dy = -GT_YSCALE(geo) * (GT_NLINES(geo) - 1.0) / (ny - 1.0) + else + dy = GT_YSCALE(geo) * (GT_NLINES(geo) - 1.0) / (ny - 1.0) + do i = l1, l2 + y[i-l1+1] = min (ymax, max (ymin, GT_YMIN(geo) + + (i - 1) * dy)) + } + + } +end + + +# GEO_SAMPLE -- Calculate the sampled reference points. + +procedure geo_sample (geo, xref, c1, c2, nxsample, yref, l1, l2, nysample, + cmode) + +pointer geo #I pointer to geotran structure +real xref[ARB] #O x reference values +int c1, c2, nxsample #I limits and number of sample points in x +real yref[ARB] #O y reference values +int l1, l2, nysample #I limits and number of sample points in y +int cmode #I coordinate computation mode + +int i +real xtempmin, xtempmax, ytempmin, ytempmax + +begin + switch (cmode) { + case GT_FOUR: + if (GT_NCOLS(geo) == 1) { + xref[1] = 0.5 + xref[2] = 1.5 + xref[3] = 1.5 + xref[4] = 0.5 + } else { + do i = c1, c2 { + xtempmin = min (real (nxsample), max (1., + real ((nxsample - 1) * (i - 0.5) + (GT_NCOLS(geo) - + nxsample)) / (GT_NCOLS(geo) - 1))) + xtempmax = min (real (nxsample), max (1., + real ((nxsample - 1) * (i + 0.5) + (GT_NCOLS(geo) - + nxsample)) / (GT_NCOLS(geo) - 1))) + xref[4*(i-c1)+1] = xtempmin + xref[4*(i-c1)+2] = xtempmax + xref[4*(i-c1)+3] = xtempmax + xref[4*(i-c1)+4] = xtempmin + } + + } + case GT_TWO: + if (GT_NCOLS(geo) == 1) { + xref[1] = 0.5 + xref[2] = 1.5 + } else { + do i = c1, c2 { + xtempmin = min (real (nxsample), max (1., + real ((nxsample - 1) * (i - 0.5) + (GT_NCOLS(geo) - + nxsample)) / (GT_NCOLS(geo) - 1))) + xtempmax = min (real (nxsample), max (1., + real ((nxsample - 1) * (i + 0.5) + (GT_NCOLS(geo) - + nxsample)) / (GT_NCOLS(geo) - 1))) + xref[2*(i-c1)+1] = xtempmin + xref[2*(i-c1)+2] = xtempmax + } + } + case GT_ONE: + if (GT_NCOLS(geo) == 1) + xref[1] = 1.0 + else { + do i = c1, c2 + xref[i-c1+1] = min (real (nxsample), max (1., + real ((nxsample - 1) * i + (GT_NCOLS(geo) - + nxsample)) / (GT_NCOLS(geo) - 1))) + } + } + + switch (cmode) { + case GT_FOUR: + if (GT_NLINES(geo) == 1) { + yref[1] = 0.5 + yref[2] = 0.5 + yref[3] = 1.5 + yref[4] = 1.5 + } else { + do i = l1, l2 { + ytempmin = min (real (nysample), max (1., + real ((nysample - 1) * (i - 0.5) + (GT_NLINES(geo) - + nysample)) / (GT_NLINES(geo) - 1))) + ytempmax = min (real (nysample), max (1., + real ((nysample - 1) * (i + 0.5) + (GT_NLINES(geo) - + nysample)) / (GT_NLINES(geo) - 1))) + yref[4*(i-l1)+1] = ytempmin + yref[4*(i-l1)+2] = ytempmin + yref[4*(i-l1)+3] = ytempmax + yref[4*(i-l1)+4] = ytempmax + } + } + case GT_TWO: + if (GT_NLINES(geo) == 1) { + yref[1] = 0.5 + yref[2] = 1.5 + } else { + do i = l1, l2 { + ytempmin = min (real (nysample), max (1., + real ((nysample - 1) * (i - 0.5) + (GT_NLINES(geo) - + nysample)) / (GT_NLINES(geo) - 1))) + ytempmax = min (real (nysample), max (1., + real ((nysample - 1) * (i + 0.5) + (GT_NLINES(geo) - + nysample)) / (GT_NLINES(geo) - 1))) + yref[2*(i-l1)+1] = ytempmin + yref[2*(i-l1)+2] = ytempmax + } + } + case GT_ONE: + if (GT_NLINES(geo) == 1) + yref[1] = 1.0 + else { + do i = l1, l2 + yref[i-l1+1] = min (real (nysample), max (1., + real ((nysample - 1) * i + (GT_NLINES(geo) - + nysample)) / (GT_NLINES(geo) - 1))) + } + } +end + + +# GEO_XBUFFER -- Compute the x interpolant and coordinates. + +procedure geo_xbuffer (s1, s2, msi, xsample, ysample, c1, c2, l1, l2, buf) + +pointer s1, s2 #I pointers to the x surface +pointer msi #I interpolant +real xsample[ARB] #I sampled x reference coordinates +real ysample[ARB] #I sampled y reference coordinates +int c1, c2 #I columns of interest in sampled image +int l1, l2 #I lines of interest in the sampled image +pointer buf #I pointer to output buffer + +int i, ncols, nlines, llast1, llast2, nclast, nllast +pointer sp, sf, y, z, buf1, buf2 + +begin + ncols = c2 - c1 + 1 + nlines = l2 - l1 + 1 + + # Combine surfaces. + if (s2 == NULL) + call gscopy (s1, sf) + else + call gsadd (s1, s2, sf) + + # Allocate working space. + call smark (sp) + call salloc (y, ncols, TY_REAL) + call salloc (z, ncols, TY_REAL) + + # If buffer undefined then allocate memory for the buffer. Reallocate + # the buffer if the number of lines or columns changes. + if (buf == NULL) { + call malloc (buf, ncols * nlines, TY_REAL) + llast1 = l1 - nlines + llast2 = l2 - nlines + } else if ((nlines != nllast) || (ncols != nclast)) { + call realloc (buf, ncols * nlines, TY_REAL) + llast1 = l1 - nlines + llast2 = l2 - nlines + } + + # Compute the coordinates. + if (l1 < llast1) { + do i = l2, l1, -1 { + if (i > llast1) + buf1 = buf + (i - llast1) * ncols + else { + buf1 = z + call amovkr (ysample[i], Memr[y], ncols) + call gsvector (sf, xsample[c1], Memr[y], Memr[buf1], ncols) + } + buf2 = buf + (i - l1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } else if (l2 > llast2) { + do i = l1, l2 { + if (i < llast2) + buf1 = buf + (i - llast1) * ncols + else { + buf1 = z + call amovkr (ysample[i], Memr[y], ncols) + call gsvector (sf, xsample[c1], Memr[y], Memr[buf1], ncols) + } + buf2 = buf + (i - l1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } + + llast1 = l1 + llast2 = l2 + nclast = ncols + nllast = nlines + + # Fit the interpolant. + if (nlines == 1) + call asifit (msi, Memr[buf], ncols) + else + call msifit (msi, Memr[buf], ncols, nlines, ncols) + + call gsfree (sf) + call sfree (sp) +end + + +# GEO_YBUFFER -- Compute the y interpolant and coordinates. + +procedure geo_ybuffer (s1, s2, msi, xsample, ysample, c1, c2, l1, l2, buf) + +pointer s1, s2 #I pointers to the y surface +pointer msi #I interpolant +real xsample[ARB] #I sampled x reference coordinates +real ysample[ARB] #I sampled y reference coordinates +int c1, c2 #I columns of interest in sampled image +int l1, l2 #I lines of interest in the sampled image +pointer buf #I pointer to output buffer + +int i, ncols, nlines, llast1, llast2, nclast, nllast +pointer sp, sf, y, z, buf1, buf2 + +begin + ncols = c2 - c1 + 1 + nlines = l2 - l1 + 1 + + # Combine surfaces. + if (s2 == NULL) + call gscopy (s1, sf) + else + call gsadd (s1, s2, sf) + + # Allocate working space. + call smark (sp) + call salloc (y, ncols, TY_REAL) + call salloc (z, ncols, TY_REAL) + + # If buffer undefined then allocate memory for the buffer. Reallocate + # the buffer if the number of lines or columns changes. + if (buf == NULL) { + call malloc (buf, ncols * nlines, TY_REAL) + llast1 = l1 - nlines + llast2 = l2 - nlines + } else if ((nlines != nllast) || (ncols != nclast)) { + call realloc (buf, ncols * nlines, TY_REAL) + llast1 = l1 - nlines + llast2 = l2 - nlines + } + + # Compute the coordinates. + if (l1 < llast1) { + do i = l2, l1, -1 { + if (i > llast1) + buf1 = buf + (i - llast1) * ncols + else { + buf1 = z + call amovkr (ysample[i], Memr[y], ncols) + call gsvector (sf, xsample[c1], Memr[y], Memr[buf1], ncols) + } + buf2 = buf + (i - l1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } else if (l2 > llast2) { + do i = l1, l2 { + if (i < llast2) + buf1 = buf + (i - llast1) * ncols + else { + buf1 = z + call amovkr (ysample[i], Memr[y], ncols) + call gsvector (sf, xsample[c1], Memr[y], Memr[buf1], ncols) + } + buf2 = buf + (i - l1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } + + llast1 = l1 + llast2 = l2 + nclast = ncols + nllast = nlines + + # Fit the interpolant. + if (nlines == 1) + call asifit (msi, Memr[buf], ncols) + else + call msifit (msi, Memr[buf], ncols, nlines, ncols) + + call gsfree (sf) + call sfree (sp) +end + + +# GEO_JBUFFER -- Fit the jacobian surface. + +procedure geo_jbuffer (sx1, sy1, sx2, sy2, jmsi, xsample, ysample, c1, c2, l1, + l2, jbuf) + +pointer sx1, sy1 #I pointers to the linear surfaces +pointer sx2, sy2 #I pointers to the distortion surfaces +pointer jmsi #I interpolant +real xsample[ARB] #I sampled x reference coordinates +real ysample[ARB] #I sampled y reference coordinates +int c1, c2 #I columns of interest in sampled image +int l1, l2 #I lines of interest in the sampled image +pointer jbuf #I pointer to output buffer + +int i, ncols, nlines, llast1, llast2, nclast, nllast +pointer sp, sx, sy, y, z, buf1, buf2 + +begin + ncols = c2 - c1 + 1 + nlines = l2 - l1 + 1 + + # Combine surfaces. + if (sx2 == NULL) + call gscopy (sx1, sx) + else + call gsadd (sx1, sx2, sx) + if (sy1 == NULL) + sy = NULL + else if (sy2 == NULL) + call gscopy (sy1, sy) + else + call gsadd (sy1, sy2, sy) + + call smark (sp) + call salloc (y, ncols, TY_REAL) + call salloc (z, ncols, TY_REAL) + + # If buffer undefined then allocate memory for the buffer. Reallocate + # the buffer if the number of lines or columns changes. + if (jbuf == NULL) { + call malloc (jbuf, ncols * nlines, TY_REAL) + llast1 = l1 - nlines + llast2 = l2 - nlines + } else if ((nlines != nllast) || (ncols != nclast)) { + call realloc (jbuf, ncols * nlines, TY_REAL) + llast1 = l1 - nlines + llast2 = l2 - nlines + } + + # Compute surface. + if (l1 < llast1) { + do i = l2, l1, -1 { + if (i > llast1) + buf1 = jbuf + (i - llast1) * ncols + else { + buf1 = z + call amovkr (ysample[i], Memr[y], ncols) + call geo_jgsvector (sx, sy, xsample[c1], Memr[y], + Memr[buf1], ncols) + } + buf2 = jbuf + (i - l1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } else if (l2 > llast2) { + do i = l1, l2 { + if (i < llast2) + buf1 = jbuf + (i - llast1) * ncols + else { + buf1 = z + call amovkr (ysample[i], Memr[y], ncols) + call geo_jgsvector (sx, sy, xsample[c1], Memr[y], + Memr[buf1], ncols) + } + buf2 = jbuf + (i - l1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } + + # Update buffer pointers. + llast1 = l1 + llast2 = l2 + nclast = ncols + nllast = nlines + + # Fit the interpolant. + if (nlines == 1) + call asifit (jmsi, Memr[jbuf], ncols) + else + call msifit (jmsi, Memr[jbuf], ncols, nlines, ncols) + + call gsfree (sx) + call gsfree (sy) + call sfree (sp) +end + + +# GEO_JGSVECTOR -- Procedure to compute the Jacobian of the transformation. + +procedure geo_jgsvector (sx, sy, x, y, out, ncols) + +pointer sx, sy #I surface descriptors +real x[ARB] #I x values +real y[ARB] #I y values +real out[ARB] #O output values +int ncols #I number of points + +pointer sp, der1, der2 + +begin + call smark (sp) + + if (sy == NULL) { + call gsder (sx, x, y, out, ncols, 1, 0) + } else { + call salloc (der1, ncols, TY_REAL) + call salloc (der2, ncols, TY_REAL) + call gsder (sx, x, y, Memr[der1], ncols, 1, 0) + call gsder (sy, x, y, Memr[der2], ncols, 0, 1) + call amulr (Memr[der1], Memr[der2], out, ncols) + call gsder (sx, x, y, Memr[der1], ncols, 0, 1) + call gsder (sy, x, y, Memr[der2], ncols, 1, 0) + call amulr (Memr[der1], Memr[der2], Memr[der1], ncols) + call asubr (out, Memr[der1], out, ncols) + } + + call sfree (sp) +end + + +# GEO_MSIVECTOR -- Procedure to interpolate the surface coordinates + +procedure geo_msivector (in, out, geo, xmsi, ymsi, jmsi, msi, sx1, sy1, sx2, + sy2, xref, c1, c2, nxsample, yref, l1, l2, nysample, x0, y0) + +pointer in #I pointer to input image +pointer out #I pointer to output image +pointer geo #I pointer to geotran structure +pointer xmsi, ymsi #I pointer to the interpolation cord surfaces +pointer jmsi #I pointer to Jacobian surface +pointer msi #I pointer to interpolation surface +pointer sx1, sy1 #I pointers to linear surfaces +pointer sx2, sy2 #I pointer to higher order surfaces +real xref[ARB] #I x reference coordinates +int c1, c2 #I column limits in output image +int nxsample #I the x sample size +real yref[ARB] #I y reference coordinates +int l1, l2 #I line limits in output image +int nysample #I the y sample size +int x0, y0 #I zero points of interpolation coordinates + +int j, ncols, nlines, ncols4, nlines4 +int imc1, imc2, iml1, iml2, nicols, nilines +pointer sp, txref, tyref, x, y, xin, yin, inbuf, outbuf +real xmin, xmax, ymin, ymax, factor +pointer imgs1r(), imgs2r(), imps1r(), imps2r() +real geo_jfactor() + +begin + ncols = c2 - c1 + 1 + nlines = l2 - l1 + 1 + + # Find min max of interpolation coords. + if (IM_NDIM(in) == 1) + call geo_iminmax (xref, yref, c1, c2, l1, l2, x0, 0, + xmsi, ymsi, xmin, xmax, ymin, ymax) + else + call geo_iminmax (xref, yref, c1, c2, l1, l2, x0, y0, + xmsi, ymsi, xmin, xmax, ymin, ymax) + + # Get the appropriate image section and fit the interpolant. + imc1 = int(xmin) - GT_NXYMARGIN(geo) + if (imc1 <= 0) + imc1 = imc1 - 1 + imc2 = nint (xmax) + GT_NXYMARGIN(geo) + 1 + nicols = imc2 - imc1 + 1 + if (IM_NDIM(in) == 1) { + ncols4 = 2 * ncols + nlines4 = 2 * nlines + iml1 = 1 + iml2 = 1 + nilines = 1 + inbuf = imgs1r (in, imc1, imc2) + if (inbuf == EOF) + call error (0, "Error reading image") + call asifit (msi, Memr[inbuf], nicols) + } else { + ncols4 = 4 * ncols + nlines4 = 4 * nlines + iml1 = int(ymin) - GT_NXYMARGIN(geo) + if (iml1 <= 0) + iml1 = iml1 - 1 + iml2 = nint (ymax) + GT_NXYMARGIN(geo) + 1 + nilines = iml2 - iml1 + 1 + inbuf = imgs2r (in, imc1, imc2, iml1, iml2) + if (inbuf == EOF) + call error (0, "Error reading image") + call msifit (msi, Memr[inbuf], nicols, nilines, nicols) + } + + # Allocate working space. + call smark (sp) + if (GT_INTERPOLANT(geo) == II_DRIZZLE || GT_INTERPOLANT(geo) == + II_BIDRIZZLE) { + call salloc (txref, ncols4, TY_REAL) + call salloc (tyref, nlines4, TY_REAL) + call salloc (x, ncols4, TY_REAL) + call salloc (y, ncols4, TY_REAL) + call salloc (xin, ncols4, TY_REAL) + call salloc (yin, ncols4, TY_REAL) + if (IM_NDIM(in) == 1) + call geo_sample (geo, Memr[txref], c1, c2, nxsample, + Memr[tyref], l1, l2, nysample, GT_TWO) + else + call geo_sample (geo, Memr[txref], c1, c2, nxsample, + Memr[tyref], l1, l2, nysample, GT_FOUR) + call aaddkr (Memr[txref], real (-x0 + 1), Memr[x], ncols4) + } else { + call salloc (x, ncols, TY_REAL) + call salloc (y, ncols, TY_REAL) + call salloc (xin, ncols, TY_REAL) + call salloc (yin, ncols, TY_REAL) + call aaddkr (xref[c1], real (-x0 + 1), Memr[x], ncols) + } + + # Compute the output buffer. + do j = l1, l2 { + + # Write the output image. + if (IM_NDIM(in) == 1) + outbuf = imps1r (out, c1, c2) + else + outbuf = imps2r (out, c1, c2, j, j) + if (outbuf == EOF) + call error (0, "Error writing output image") + + # Compute the interpolation coordinates. + if (GT_INTERPOLANT(geo) == II_DRIZZLE || GT_INTERPOLANT(geo) == + II_BIDRIZZLE) { + if (IM_NDIM(in) == 1) { + call asivector (xmsi, Memr[x], Memr[xin], ncols4) + call amovkr (1.0, Memr[yin], ncols4) + } else { + #call amovkr (yref[j] + real (-y0 + 1), Memr[y], ncols) + call geo_repeat (Memr[tyref+4*(j-l1)], 4, Memr[y], ncols) + call aaddkr (Memr[y], real(-y0 + 1), Memr[y], ncols4) + call msivector (xmsi, Memr[x], Memr[y], Memr[xin], ncols4) + call msivector (ymsi, Memr[x], Memr[y], Memr[yin], ncols4) + } + if (imc1 != 1) + call aaddkr (Memr[xin], real (-imc1 + 1), Memr[xin], ncols4) + if (iml1 != 1) + call aaddkr (Memr[yin], real (-iml1 + 1), Memr[yin], ncols4) + } else { + if (IM_NDIM(in) == 1) { + call asivector (xmsi, Memr[x], Memr[xin], ncols) + call amovkr (1.0, Memr[yin], ncols) + } else { + call amovkr (yref[j] + real (-y0 + 1), Memr[y], ncols) + call msivector (xmsi, Memr[x], Memr[y], Memr[xin], ncols) + call msivector (ymsi, Memr[x], Memr[y], Memr[yin], ncols) + } + if (imc1 != 1) + call aaddkr (Memr[xin], real (-imc1 + 1), Memr[xin], ncols) + if (iml1 != 1) + call aaddkr (Memr[yin], real (-iml1 + 1), Memr[yin], ncols) + } + + # Interpolate in the input image. + if (IM_NDIM(in) == 1) + call asivector (msi, Memr[xin], Memr[outbuf], ncols) + else + call msivector (msi, Memr[xin], Memr[yin], Memr[outbuf], ncols) + + # Preserve flux in image. + if (GT_FLUXCONSERVE(geo) == YES) { + factor = GT_XSCALE(geo) * GT_YSCALE(geo) + if (GT_GEOMODE(geo) == GT_LINEAR || (sx2 == NULL && sy2 == + NULL)) { + if (IM_NDIM(in) == 1) + call amulkr (Memr[outbuf], factor * geo_jfactor (sx1, + NULL), Memr[outbuf], ncols) + else + call amulkr (Memr[outbuf], factor * geo_jfactor (sx1, + sy1), Memr[outbuf], ncols) + } else { + if (IM_NDIM(in) == 1) + call geo_msiflux (jmsi, xref, yref, Memr[outbuf], + c1, c2, 0, x0, y0) + else + call geo_msiflux (jmsi, xref, yref, Memr[outbuf], + c1, c2, j, x0, y0) + call amulkr (Memr[outbuf], factor, Memr[outbuf], ncols) + } + } + } + + call sfree (sp) +end + + +# GEO_GSVECTOR -- Evaluate the output image pixels using fitted coordinate +# values and image interpolation. + +procedure geo_gsvector (input, output, geo, msi, xref, c1, c2, yref, l1, l2, + sx1, sy1, sx2, sy2) + +pointer input #I pointer to input image +pointer output #I pointer to output image +pointer geo #I pointer to geotran structure +pointer msi #I pointer to interpolant +real xref[ARB] #I x reference array +int c1, c2 #I columns of interest in output image +real yref[ARB] #I y reference array +int l1, l2 #I lines of interest in the output image +pointer sx1, sy1 #I linear surface descriptors +pointer sx2, sy2 #I distortion surface descriptors + +int j, ncols, nlines, ncols4, nlines4, nicols, nilines +int imc1, imc2, iml1, iml2 +pointer sp, txref, tyref, y, xin, yin, temp, inbuf, outbuf +real xmin, xmax, ymin, ymax, factor +pointer imgs1r(), imgs2r(), imps1r(), imps2r() +real gsgetr(), geo_jfactor() + +begin + # Compute the number of columns. + ncols = c2 - c1 + 1 + nlines = l2 - l1 + 1 + + # Compute the maximum and minimum coordinates. + call geo_minmax (xref, yref, c1, c2, l1, l2, sx1, sy1, sx2, sy2, + xmin, xmax, ymin, ymax) + + # Get the appropriate image section and fill the buffer. + imc1 = int(xmin) - GT_NXYMARGIN(geo) + if (imc1 <= 0) + imc1 = imc1 - 1 + imc2 = nint (xmax) + GT_NXYMARGIN(geo) + 1 + nicols = imc2 - imc1 + 1 + if (IM_NDIM(input) == 1) { + iml1 = 1 + iml2 = 1 + nilines = 1 + ncols4 = 2 * ncols + nlines4 = 2 * nlines + inbuf = imgs1r (input, imc1, imc2) + if (inbuf == EOF) + call error (0, "Error reading image") + call asifit (msi, Memr[inbuf], nicols) + } else { + iml1 = int(ymin) - GT_NXYMARGIN(geo) + if (iml1 <= 0) + iml1 = iml1 - 1 + iml2 = nint (ymax) + GT_NXYMARGIN(geo) + 1 + nilines = iml2 - iml1 + 1 + ncols4 = 4 * ncols + nlines4 = 4 * nlines + inbuf = imgs2r (input, imc1, imc2, iml1, iml2) + if (inbuf == EOF) + call error (0, "Error reading image") + call msifit (msi, Memr[inbuf], nicols, nilines, nicols) + } + + # Allocate working space. + call smark (sp) + if (GT_INTERPOLANT(geo) == II_DRIZZLE || GT_INTERPOLANT(geo) == + II_BIDRIZZLE) { + call salloc (txref, ncols4, TY_REAL) + call salloc (tyref, nlines4, TY_REAL) + call salloc (y, ncols4, TY_REAL) + call salloc (xin, ncols4, TY_REAL) + call salloc (yin, ncols4, TY_REAL) + call salloc (temp, ncols4, TY_REAL) + if (IM_NDIM(input) == 1) + call geo_ref (geo, Memr[txref], c1, c2, GT_NCOLS(geo), + Memr[tyref], l1, l2, GT_NLINES(geo), gsgetr (sx1, GSXMIN), + gsgetr (sx1, GSXMAX), gsgetr (sx1, GSYMIN), gsgetr (sx1, + GSYMAX), GT_TWO) + else + call geo_ref (geo, Memr[txref], c1, c2, GT_NCOLS(geo), + Memr[tyref], l1, l2, GT_NLINES(geo), gsgetr (sx1, GSXMIN), + gsgetr (sx1, GSXMAX), gsgetr (sx1, GSYMIN), gsgetr (sx1, + GSYMAX), GT_FOUR) + } else { + call salloc (y, ncols, TY_REAL) + call salloc (xin, ncols, TY_REAL) + call salloc (yin, ncols, TY_REAL) + call salloc (temp, ncols, TY_REAL) + } + + # Compute the pixels. + do j = l1, l2 { + + # Get output image buffer. + if (IM_NDIM(input) == 1) + outbuf = imps1r (output, c1, c2) + else + outbuf = imps2r (output, c1, c2, j, j) + if (output == EOF) + call error (0, "Error writing output image") + + # Compute the interpolation coordinates. + if (GT_INTERPOLANT(geo) == II_DRIZZLE || GT_INTERPOLANT(geo) == + II_BIDRIZZLE) { + + # Set the y coordinate. + if (IM_NDIM(input) == 1) + call geo_repeat (Memr[tyref+2*(j-l1)], 2, Memr[y], ncols) + else + call geo_repeat (Memr[tyref+4*(j-l1)], 4, Memr[y], ncols) + + # Fit x coords. + call gsvector (sx1, Memr[txref], Memr[y], Memr[xin], ncols4) + if (sx2 != NULL) { + call gsvector (sx2, Memr[txref], Memr[y], Memr[temp], + ncols4) + call aaddr (Memr[xin], Memr[temp], Memr[xin], ncols4) + } + if (imc1 != 1) + call aaddkr (Memr[xin], real (-imc1 + 1), Memr[xin], ncols4) + + # Fit y coords. + call gsvector (sy1, Memr[txref], Memr[y], Memr[yin], ncols4) + if (sy2 != NULL) { + call gsvector (sy2, Memr[txref], Memr[y], Memr[temp], + ncols4) + call aaddr (Memr[yin], Memr[temp], Memr[yin], ncols4) + } + if (iml1 != 1) + call aaddkr (Memr[yin], real (-iml1 + 1), Memr[yin], ncols4) + + } else { + + # Set the y coordinate. + call amovkr (yref[j], Memr[y], ncols) + + # Fit x coords. + call gsvector (sx1, xref[c1], Memr[y], Memr[xin], ncols) + if (sx2 != NULL) { + call gsvector (sx2, xref[c1], Memr[y], Memr[temp], ncols) + call aaddr (Memr[xin], Memr[temp], Memr[xin], ncols) + } + if (imc1 != 1) + call aaddkr (Memr[xin], real (-imc1 + 1), Memr[xin], ncols) + + # Fit y coords. + call gsvector (sy1, xref[c1], Memr[y], Memr[yin], ncols) + if (sy2 != NULL) { + call gsvector (sy2, xref[c1], Memr[y], Memr[temp], ncols) + call aaddr (Memr[yin], Memr[temp], Memr[yin], ncols) + } + if (iml1 != 1) + call aaddkr (Memr[yin], real (-iml1 + 1), Memr[yin], ncols) + } + + # Interpolate in input image. + if (IM_NDIM(input) == 1) + call asivector (msi, Memr[xin], Memr[outbuf], ncols) + else + call msivector (msi, Memr[xin], Memr[yin], Memr[outbuf], ncols) + + # Preserve flux in image. + if (GT_FLUXCONSERVE(geo) == YES) { + factor = GT_XSCALE(geo) * GT_YSCALE(geo) + if (GT_GEOMODE(geo) == GT_LINEAR || (sx2 == NULL && sy2 == + NULL)) { + if (IM_NDIM(input) == 1) + call amulkr (Memr[outbuf], factor * geo_jfactor (sx1, + NULL), Memr[outbuf], ncols) + else + call amulkr (Memr[outbuf], factor * geo_jfactor (sx1, + sy1), Memr[outbuf], ncols) + } else { + if (IM_NDIM(input) == 1) + call geo_gsflux (xref, yref, Memr[outbuf], c1, c2, j, + sx1, NULL, sx2, NULL) + else + call geo_gsflux (xref, yref, Memr[outbuf], c1, c2, j, + sx1, sy1, sx2, sy2) + call amulkr (Memr[outbuf], factor, Memr[outbuf], ncols) + } + } + } + + call sfree (sp) +end + + +# GEO_IMINMAX -- Find minimum and maximum interpolation coordinates. + +procedure geo_iminmax (xref, yref, c1, c2, l1, l2, x0, y0, xmsi, ymsi, xmin, + xmax, ymin, ymax) + +real xref[ARB] #I x reference coords +real yref[ARB] #I y reference coords +int c1, c2 #I columns limits +int l1, l2 #I line limits +int x0, y0 #I interpolation coord zero points +pointer xmsi, ymsi #I coord surfaces +real xmin, xmax #O output xmin and xmax +real ymin, ymax #O output ymin and ymax + +int j, ncols +pointer sp, x, y, xin, yin +real mintemp, maxtemp, x1, x2, y1, y2 +real asieval(), msieval() + +begin + call smark (sp) + ncols = c2 - c1 + 1 + call salloc (x, ncols, TY_REAL) + call salloc (y, ncols, TY_REAL) + call salloc (xin, ncols, TY_REAL) + call salloc (yin, ncols, TY_REAL) + + xmin = MAX_REAL + xmax = -MAX_REAL + ymin = MAX_REAL + ymax = -MAX_REAL + + # find the minimum and maximum + do j = l1, l2 { + + if (j == l1 || j == l2) { + + call aaddkr (xref[c1], real (-x0 + 1), Memr[x], ncols) + if (y0 <= 0) { + call asivector (xmsi, Memr[x], Memr[xin], ncols) + ymin = 1.0 + ymax = 1.0 + } else { + call amovkr (yref[j] + real (-y0 + 1), Memr[y], ncols) + call msivector (xmsi, Memr[x], Memr[y], Memr[xin], ncols) + call msivector (ymsi, Memr[x], Memr[y], Memr[yin], ncols) + call alimr (Memr[yin], ncols, mintemp, maxtemp) + ymin = min (ymin, mintemp) + ymax = max (ymax, maxtemp) + } + call alimr (Memr[xin], ncols, mintemp, maxtemp) + xmin = min (xmin, mintemp) + xmax = max (xmax, maxtemp) + } else { + if (y0 <= 0) { + x1 = asieval (xmsi, xref[c1] + real (-x0 + 1)) + x2 = asieval (xmsi, xref[c1+ncols-1] + real (-x0 + 1)) + ymin = 1.0 + ymax = 1.0 + } else { + x1 = msieval (xmsi, xref[c1] + real (-x0 + 1), + yref[j] + real (-y0 + 1)) + x2 = msieval (xmsi, xref[c1+ncols-1] + real (-x0 + 1), + yref[j] + real (-y0 + 1)) + y1 = msieval (ymsi, xref[c1] + real (-x0 + 1), + yref[j] + real (-y0 + 1)) + y2 = msieval (ymsi, xref[c1+ncols-1] + real (-x0 + 1), + yref[j] + real (-y0 + 1)) + ymin = min (ymin, y1, y2) + ymax = max (ymax, y1, y2) + } + xmin = min (xmin, x1, x2) + xmax = max (xmax, x1, x2) + + } + } + + call sfree (sp) + +end + + +# GEO_MINMAX -- Compute the minimum and maximum fitted coordinates. + +procedure geo_minmax (xref, yref, c1, c2, l1, l2, sx1, sy1, sx2, sy2, + xmin, xmax, ymin, ymax) + +real xref[ARB] #I x reference coords +real yref[ARB] #I y reference coords +int c1, c2 #I columns limits +int l1, l2 #I line limits +pointer sx1, sy1 #I linear surface descriptors +pointer sx2, sy2 #I distortion surface descriptors +real xmin, xmax #O output xmin and xmax +real ymin, ymax #O output ymin and ymax + +int j, ncols +pointer sp, y, xin, yin, temp +real x1, x2, y1, y2, mintemp, maxtemp +real gseval() + +begin + call smark (sp) + ncols = c2 - c1 + 1 + call salloc (y, ncols, TY_REAL) + call salloc (xin, ncols, TY_REAL) + call salloc (yin, ncols, TY_REAL) + call salloc (temp, ncols, TY_REAL) + + xmin = MAX_REAL + xmax = -MAX_REAL + ymin = MAX_REAL + ymax = -MAX_REAL + + # Find the maximum and minimum coordinates. + do j = l1, l2 { + + if (j == l1 || j == l2) { + + call amovkr (yref[j], Memr[y], ncols) + call gsvector (sx1, xref[c1], Memr[y], Memr[xin], ncols) + if (sx2 != NULL) { + call gsvector (sx2, xref[c1], Memr[y], Memr[temp], ncols) + call aaddr (Memr[xin], Memr[temp], Memr[xin], ncols) + } + call gsvector (sy1, xref[c1], Memr[y], Memr[yin], ncols) + if (sy2 != NULL) { + call gsvector (sy2, xref[c1], Memr[y], Memr[temp], ncols) + call aaddr (Memr[yin], Memr[temp], Memr[yin], ncols) + } + + call alimr (Memr[xin], ncols, mintemp, maxtemp) + xmin = min (xmin, mintemp) + xmax = max (xmax, maxtemp) + call alimr (Memr[yin], ncols, mintemp, maxtemp) + ymin = min (ymin, mintemp) + ymax = max (ymax, maxtemp) + + } else { + + x1 = gseval (sx1, xref[c1], yref[j]) + x2 = gseval (sx1, xref[c1+ncols-1], yref[j]) + if (sx2 != NULL) { + x1 = x1 + gseval (sx2, xref[c1], yref[j]) + x2 = x2 + gseval (sx2, xref[c1+ncols-1], yref[j]) + } + xmin = min (xmin, x1, x2) + xmax = max (xmax, x1, x2) + + y1 = gseval (sy1, xref[c1], yref[j]) + y2 = gseval (sy1, xref[c1+ncols-1], yref[j]) + if (sy2 != NULL) { + y1 = y1 + gseval (sy2, xref[c1], yref[j]) + y2 = y2 + gseval (sy2, xref[c1+ncols-1], yref[j]) + } + ymin = min (ymin, y1, y2) + ymax = max (ymax, y1, y2) + + } + } + + call sfree (sp) +end + + +# GEO_MARGSET -- Set up interpolation margin + +procedure geo_margset (sx1, sy1, sx2, sy2, xmin, xmax, ncols, ymin, ymax, + nlines, interpolant, nsinc, nxymargin) + +pointer sx1, sy1 #I linear surface descriptors +pointer sx2, sy2 #I distortion surface descriptors +real xmin, xmax #I the reference coordinate x limits +int ncols #I the number of output image columns +real ymin, ymax #I the reference coordinate y limits +int nlines #I the number of output image lines +int interpolant #I the interpolant type +int nsinc #I the sinc width +int nxymargin #O the interpolation margin + +int dist1, dist2, dist3, dist4, dist5, dist6 +pointer newsx, newsy +real x1, y1, x2, y2 +real gseval() + +begin + if (interpolant == II_SPLINE3 || interpolant == II_BISPLINE3) { + nxymargin = NMARGIN_SPLINE3 + } else if (interpolant == II_LSINC || interpolant == II_BILSINC) { + nxymargin = nsinc + } else if (interpolant == II_SINC || interpolant == II_BISINC) { + nxymargin = nsinc + } else if (interpolant == II_DRIZZLE || interpolant == II_BIDRIZZLE) { + if (sx2 == NULL) + call gscopy (sx1, newsx) + else + call gsadd (sx1, sx2, newsx) + if (sy2 == NULL) + call gscopy (sy1, newsy) + else + call gsadd (sy1, sy2, newsy) + x1 = gseval (newsx, xmin, ymin) + y1 = gseval (newsy, xmin, ymin) + x2 = gseval (newsx, xmax, ymin) + y2 = gseval (newsy, xmax, ymin) + dist1 = sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2) / ncols + x1 = gseval (newsx, xmax, ymax) + y1 = gseval (newsy, xmax, ymax) + dist2 = sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2) / nlines + x2 = gseval (newsx, xmin, ymax) + y2 = gseval (newsy, xmin, ymax) + dist3 = sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2) / ncols + x1 = gseval (newsx, xmin, ymin) + y1 = gseval (newsy, xmin, ymin) + dist4 = sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2) / nlines + x1 = gseval (newsx, xmin, (ymin + ymax) / 2.0) + y1 = gseval (newsy, xmin, (ymin + ymax) / 2.0) + x2 = gseval (newsx, xmax, (ymin + ymax) / 2.0) + y2 = gseval (newsy, xmax, (ymin + ymax) / 2.0) + dist5 = sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2) / ncols + x1 = gseval (newsx, (xmin + xmax) / 2.0, ymin) + y1 = gseval (newsy, (xmin + xmax) / 2.0, ymin) + x2 = gseval (newsx, (xmin + xmax) / 2.0, ymax) + y2 = gseval (newsy, (xmin + xmax) / 2.0, ymax) + dist6 = sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2) / nlines + nxymargin = max (NMARGIN, dist1, dist2, dist3, dist4, + dist5, dist6) + call gsfree (newsx) + call gsfree (newsy) + } else { + nxymargin = NMARGIN + } +end + + +# GEO_IMSET -- Set up input image boundary conditions. + +procedure geo_imset (im, geo, sx1, sy1, sx2, sy2, xref, nx, yref, ny) + +pointer im #I pointer to image +pointer geo #I pointer to geotran structure +pointer sx1, sy1 #I linear surface descriptors +pointer sx2, sy2 #I distortion surface descriptors +real xref[ARB] #I x reference coordinates +int nx #I number of x reference coordinates +real yref[ARB] #I y reference coordinates +int ny #I number of y reference coordinates + +int bndry, npts +pointer sp, x1, x2, y1, y2, xtemp, ytemp +real xn1, xn2, xn3, xn4, yn1, yn2, yn3, yn4, xmin, xmax, ymin, ymax +real gseval() + +begin + npts = max (nx, ny) + + xn1 = gseval (sx1, GT_XMIN(geo), GT_YMIN(geo)) + xn2 = gseval (sx1, GT_XMAX(geo), GT_YMIN(geo)) + xn3 = gseval (sx1, GT_XMAX(geo), GT_YMAX(geo)) + xn4 = gseval (sx1, GT_XMIN(geo), GT_YMAX(geo)) + + yn1 = gseval (sy1, GT_XMIN(geo), GT_YMIN(geo)) + yn2 = gseval (sy1, GT_XMAX(geo), GT_YMIN(geo)) + yn3 = gseval (sy1, GT_XMAX(geo), GT_YMAX(geo)) + yn4 = gseval (sy1, GT_XMIN(geo), GT_YMAX(geo)) + + xmin = min (xn1, xn2, xn3, xn4) + ymin = min (yn1, yn2, yn3, yn4) + xmax = max (xn1, xn2, xn3, xn4) + ymax = max (yn1, yn2, yn3, yn4) + + if (sx2 != NULL) { + call smark (sp) + call salloc (x1, npts, TY_REAL) + call salloc (x2, npts, TY_REAL) + call salloc (xtemp, npts, TY_REAL) + call salloc (ytemp, npts, TY_REAL) + + call amovkr (GT_YMIN(geo), Memr[ytemp], nx) + call gsvector (sx1, xref, Memr[ytemp], Memr[x1], nx) + call gsvector (sx2, xref, Memr[ytemp], Memr[x2], nx) + call aaddr (Memr[x1], Memr[x2], Memr[x1], nx) + call alimr (Memr[x1], nx, xn1, yn1) + + call amovkr (GT_XMAX(geo), Memr[xtemp], ny) + call gsvector (sx1, Memr[xtemp], yref, Memr[x1], ny) + call gsvector (sx2, Memr[xtemp], yref, Memr[x2], ny) + call aaddr (Memr[x1], Memr[x2], Memr[x1], ny) + call alimr (Memr[x1], ny, xn2, yn2) + + call amovkr (GT_YMAX(geo), Memr[ytemp], nx) + call gsvector (sx1, xref, Memr[ytemp], Memr[x1], nx) + call gsvector (sx2, xref, Memr[ytemp], Memr[x2], nx) + call aaddr (Memr[x1], Memr[x2], Memr[x1], nx) + call alimr (Memr[x1], nx, xn3, yn3) + + call amovkr (GT_XMIN(geo), Memr[xtemp], ny) + call gsvector (sx1, Memr[xtemp], yref, Memr[x1], ny) + call gsvector (sx2, Memr[xtemp], yref, Memr[x2], ny) + call aaddr (Memr[x1], Memr[x2], Memr[x1], ny) + call alimr (Memr[x1], ny, xn4, yn4) + + xmin = min (xn1, xn2, xn3, xn4) + xmax = max (yn1, yn2, yn3, yn4) + + call sfree (sp) + } + + if (sy2 != NULL) { + call smark (sp) + call salloc (y1, npts, TY_REAL) + call salloc (y2, npts, TY_REAL) + call salloc (xtemp, npts, TY_REAL) + call salloc (ytemp, npts, TY_REAL) + + call amovkr (GT_YMIN(geo), Memr[ytemp], nx) + call gsvector (sy1, xref, Memr[ytemp], Memr[y1], nx) + call gsvector (sy2, xref, Memr[ytemp], Memr[y2], nx) + call aaddr (Memr[y1], Memr[y2], Memr[y1], nx) + call alimr (Memr[y1], nx, xn1, yn1) + + call amovkr (GT_XMAX(geo), Memr[xtemp], ny) + call gsvector (sy1, Memr[xtemp], yref, Memr[y1], ny) + call gsvector (sy2, Memr[xtemp], yref, Memr[y2], ny) + call aaddr (Memr[y1], Memr[y2], Memr[y1], ny) + call alimr (Memr[y1], ny, xn2, yn2) + + call amovkr (GT_YMAX(geo), Memr[ytemp], nx) + call gsvector (sy1, xref, Memr[ytemp], Memr[y1], nx) + call gsvector (sy2, xref, Memr[ytemp], Memr[y2], nx) + call aaddr (Memr[y1], Memr[y2], Memr[y1], nx) + call alimr (Memr[y1], nx, xn3, yn3) + + call amovkr (GT_XMIN(geo), Memr[xtemp], ny) + call gsvector (sy1, Memr[xtemp], yref, Memr[y1], ny) + call gsvector (sy2, Memr[xtemp], yref, Memr[y2], ny) + call aaddr (Memr[y1], Memr[y2], Memr[y1], ny) + call alimr (Memr[y1], ny, xn4, yn4) + + ymin = min (xn1, xn2, xn3, xn4) + ymax = max (yn1, yn2, yn3, yn4) + + call sfree (sp) + } + + # Compute the out-of-bounds limit. + if (IM_NDIM(im) == 1) { + if (xmin < 1.0 || xmax > real (IM_LEN(im,1))) + bndry = max (1.0 - xmin, xmax - IM_LEN(im,1)) + 1 + else + bndry = 1 + } else { + if (xmin < 1.0 || ymin < 1.0 || xmax > real (IM_LEN(im,1)) || + ymax > real (IM_LEN(im,2))) + bndry = max (1.0 - xmin, 1.0 - ymin, xmax - IM_LEN(im,1), + ymax - IM_LEN(im,2)) + 1 + else + bndry = 1 + } + + call imseti (im, IM_NBNDRYPIX, bndry + GT_NXYMARGIN(geo) + 1) + call imseti (im, IM_TYBNDRY, GT_BOUNDARY(geo)) + call imsetr (im, IM_BNDRYPIXVAL, GT_CONSTANT(geo)) +end + + +# GEO_GSFLUX -- Preserve the image flux after a transformation. + +procedure geo_gsflux (xref, yref, buf, c1, c2, line, sx1, sy1, sx2, sy2) + +real xref[ARB] #I x reference coordinates +real yref[ARB] #I y reference coordinates +real buf[ARB] #O output image buffer +int c1, c2 #I column limits in the output image +int line #I line in the output image +pointer sx1, sy1 #I linear surface descriptors +pointer sx2, sy2 #I distortion surface descriptors + +int ncols +pointer sp, y, der1, der2, jacob, sx, sy + +begin + ncols = c2 - c1 + 1 + + # Get the reference coordinates. + call smark (sp) + call salloc (y, ncols, TY_REAL) + call salloc (jacob, ncols, TY_REAL) + + # Add the two surfaces together for efficiency. + if (sx2 != NULL) + call gsadd (sx1, sx2, sx) + else + call gscopy (sx1, sx) + if (sy1 == NULL) + sy = NULL + else if (sy2 != NULL) + call gsadd (sy1, sy2, sy) + else + call gscopy (sy1, sy) + + # Multiply the output buffer by the Jacobian. + call amovkr (yref[line], Memr[y], ncols) + if (sy == NULL) + call gsder (sx, xref[c1], Memr[y], Memr[jacob], ncols, 1, 0) + else { + call salloc (der1, ncols, TY_REAL) + call salloc (der2, ncols, TY_REAL) + call gsder (sx, xref[c1], Memr[y], Memr[der1], ncols, 1, 0) + call gsder (sy, xref[c1], Memr[y], Memr[der2], ncols, 0, 1) + call amulr (Memr[der1], Memr[der2], Memr[jacob], ncols) + call gsder (sx, xref[c1], Memr[y], Memr[der1], ncols, 0, 1) + call gsder (sy, xref[c1], Memr[y], Memr[der2], ncols, 1, 0) + call amulr (Memr[der1], Memr[der2], Memr[der1], ncols) + call asubr (Memr[jacob], Memr[der1], Memr[jacob], ncols) + } + call aabsr (Memr[jacob], Memr[jacob], ncols) + call amulr (buf, Memr[jacob], buf, ncols) + + # Clean up. + call gsfree (sx) + if (sy != NULL) + call gsfree (sy) + call sfree (sp) +end + + +# GEO_MSIFLUX -- Procedure to interpolate the surface coordinates + +procedure geo_msiflux (jmsi, xinterp, yinterp, outdata, c1, c2, line, x0, y0) + +pointer jmsi #I pointer to the jacobian interpolant +real xinterp[ARB] #I x reference coordinates +real yinterp[ARB] #I y reference coordinates +real outdata[ARB] #O output data +int c1, c2 #I column limits in output image +int line #I line to be flux corrected +int x0, y0 #I zero points of interpolation coordinates + +int ncols +pointer sp, x, y, jacob + +begin + # Allocate tempoaray space. + call smark (sp) + ncols = c2 - c1 + 1 + call salloc (x, ncols, TY_REAL) + call salloc (jacob, ncols, TY_REAL) + + # Calculate the x points. + if (x0 == 1) + call amovr (xinterp[c1], Memr[x], ncols) + else + call aaddkr (xinterp[c1], real (-x0 + 1), Memr[x], ncols) + + # Multiply the data by the Jacobian. + if (line == 0) { + call asivector (jmsi, Memr[x], Memr[jacob], ncols) + } else { + call salloc (y, ncols, TY_REAL) + call amovkr ((yinterp[line] + real (-y0 + 1)), Memr[y], ncols) + call msivector (jmsi, Memr[x], Memr[y], Memr[jacob], ncols) + } + call aabsr (Memr[jacob], Memr[jacob], ncols) + call amulr (outdata, Memr[jacob], outdata, ncols) + + call sfree (sp) +end + + +# GEO_JFACTOR -- Compute the Jacobian of a linear transformation. + +real procedure geo_jfactor (sx1, sy1) + +pointer sx1 #I pointer to x surface +pointer sy1 #I pointer to y surface + +real xval, yval, xx, xy, yx, yy +real gsgetr() + +begin + xval = (gsgetr (sx1, GSXMIN) + gsgetr (sx1, GSXMAX)) / 2.0 + if (sy1 == NULL) + yval = 1.0 + else + yval = (gsgetr (sy1, GSYMIN) + gsgetr (sy1, GSYMIN)) / 2.0 + + call gsder (sx1, xval, yval, xx, 1, 1, 0) + if (sy1 == NULL) { + xy = 0.0 + yy = 1.0 + yx = 0.0 + } else { + call gsder (sx1, xval, yval, xy, 1, 0, 1) + call gsder (sy1, xval, yval, yx, 1, 1, 0) + call gsder (sy1, xval, yval, yy, 1, 0, 1) + } + + return (abs (xx * yy - xy * yx)) +end + + +# GEO_REPEAT -- Copy a small repeated pattern into the output buffer. + +procedure geo_repeat (pat, npat, output, ntimes) + +real pat[ARB] #I the input pattern to be repeated +int npat #I the size of the pattern +real output[ARB] #O the output array +int ntimes #I the number of times the pattern is to be repeated + +int j, i, offset + +begin + do j = 1, ntimes { + offset = npat * j - npat + do i = 1, npat + output[offset+i] = pat[i] + } +end diff --git a/pkg/images/immatch/src/geometry/geoxytran.gx b/pkg/images/immatch/src/geometry/geoxytran.gx new file mode 100644 index 00000000..22d577f1 --- /dev/null +++ b/pkg/images/immatch/src/geometry/geoxytran.gx @@ -0,0 +1,327 @@ +include <ctype.h> +include <mach.h> +include <math.h> +include <math/gsurfit.h> + +define GEO_LINEAR 1 # Linear transformation only +define GEO_DISTORTION 2 # Distortion correction only +define GEO_GEOMETRIC 3 # Full transformation + +$for (rd) + +# GEO_LINIT -- Initialize the linear part of the transformation. + +$if (datatype == r) +procedure geo_linitr (sx1, sy1, sx2, sy2) +$else +procedure geo_linitd (sx1, sy1, sx2, sy2) +$endif + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +pointer sx2, sy2 #I/O pointer to the distortion x and y surfaces + +PIXEL xmag, ymag, xrot, yrot, xref, yref, xout, yout, xshift, yshift +$if (datatype == r) +real clgetr(), gseval() +$else +double clgetd(), dgseval() +$endif + +begin + # Initialize the surfaces. +$if (datatype == r) + call gsinit (sx1, GS_POLYNOMIAL, 2, 2, GS_XNONE, -MAX_REAL, MAX_REAL, + -MAX_REAL, MAX_REAL) + call gsinit (sy1, GS_POLYNOMIAL, 2, 2, GS_XNONE, -MAX_REAL, MAX_REAL, + -MAX_REAL, MAX_REAL) +$else + call dgsinit (sx1, GS_POLYNOMIAL, 2, 2, GS_XNONE, double (-MAX_REAL), + double (MAX_REAL), double (-MAX_REAL), double (MAX_REAL)) + call dgsinit (sy1, GS_POLYNOMIAL, 2, 2, GS_XNONE, double (-MAX_REAL), + double (MAX_REAL), double (-MAX_REAL), double (MAX_REAL)) +$endif + sx2 = NULL + sy2 = NULL + + # Get the magnification parameters. + xmag = clget$t ("xmag") + if (IS_$INDEF$T(xmag)) + xmag = PIXEL(1.0) + ymag = clget$t ("ymag") + if (IS_$INDEF$T(ymag)) + ymag = PIXEL(1.0) + + # Get the rotation parameters. + xrot = clget$t ("xrot") + if (IS_$INDEF$T(xrot)) + xrot = PIXEL(0.0) + xrot = -DEGTORAD(xrot) + yrot = clget$t ("yrot") + if (IS_$INDEF$T(yrot)) + yrot = PIXEL(0.0) + yrot = -DEGTORAD(yrot) + + # Set the magnification and rotation coefficients. + call geo_rotmag$t (sx1, sy1, xmag, ymag, xrot, yrot) + + # Compute the origin of the reference coordinates. + xref = clget$t ("xref") + if (IS_$INDEF$T(xref)) + xref = PIXEL(0.0) + yref = clget$t ("yref") + if (IS_$INDEF$T(yref)) + yref = PIXEL(0.0) + + # Compute the corresponding input coordinates. + xout = clget$t ("xout") + if (IS_$INDEF$T(xout)) +$if (datatype == r) + xout = gseval (sx1, xref, yref) +$else + xout = dgseval (sx1, xref, yref) +$endif + yout = clget$t ("yout") + if (IS_$INDEF$T(yout)) +$if (datatype == r) + yout = gseval (sy1, xref, yref) +$else + yout = dgseval (sy1, xref, yref) +$endif + + # Set the shifts. + xshift = clget$t ("xshift") + yshift = clget$t ("yshift") +$if (datatype == r) + if (IS_$INDEF$T(xshift)) + xshift = xout - gseval (sx1, xref, yref) + if (IS_$INDEF$T(yshift)) + yshift = yout - gseval (sy1, xref, yref) +$else + if (IS_$INDEF$T(xshift)) + xshift = xout - $tgseval (sx1, xref, yref) + if (IS_$INDEF$T(yshift)) + yshift = yout - $tgseval (sy1, xref, yref) +$endif + call geo_xyshift$t (sx1, sy1, xshift, yshift) +end + + +# GEO_SFREE -- Free the x and y surface fitting descriptors. + +$if (datatype == r) +procedure geo_sfreer (sx1, sy1, sx2, sy2) +$else +procedure geo_sfreed (sx1, sy1, sx2, sy2) +$endif + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +pointer sx2, sy2 #I/O pointer to the distortion x and y surfaces + +begin +$if (datatype == r) + call gsfree (sx1) + call gsfree (sy1) + if (sx2 != NULL) + call gsfree (sx2) + if (sy2 != NULL) + call gsfree (sy2) +$else + call dgsfree (sx1) + call dgsfree (sy1) + if (sx2 != NULL) + call dgsfree (sx2) + if (sy2 != NULL) + call dgsfree (sy2) +$endif +end + + +# GEO_SINIT -- Read the surface fits from the database file and make +# any requested changes. + +procedure geo_sinit$t (dt, record, geometry, sx1, sy1, sx2, sy2) + +pointer dt #I pointer to database file produced by geomap +char record[ARB] #I the name of the database record +int geometry #I the type of geometry to be computed +pointer sx1, sy1 #O pointers to the linear x and y surfaces +pointer sx2, sy2 #O pointers to the x and y distortion surfaces + +int i, rec, ncoeff, junk +PIXEL xmag, ymag, xrot, yrot, xref, yref, xout, yout, xshift, yshift +pointer newsx1, newsy1, xcoeff, ycoeff +int dtlocate(), dtscan(), dtgeti() +PIXEL clget$t() +$if (datatype == r) +errchk gsrestore +$else +errchk dgsrestore +$endif + +begin + # Locate record. + rec = dtlocate (dt, record) + + # Get linear part of fit. + ncoeff = dtgeti (dt, rec, "surface1") + call malloc (xcoeff, ncoeff, TY_PIXEL) + call malloc (ycoeff, ncoeff, TY_PIXEL) + do i = 1, ncoeff { + junk = dtscan (dt) + call garg$t (Mem$t[xcoeff+i-1]) + call garg$t (Mem$t[ycoeff+i-1]) + } + + # Restore linear part of fit. +$if (datatype == r) + call gsrestore (sx1, Mem$t[xcoeff]) + call gsrestore (sy1, Mem$t[ycoeff]) +$else + call dgsrestore (sx1, Mem$t[xcoeff]) + call dgsrestore (sy1, Mem$t[ycoeff]) +$endif + + # Get geometric transformation. + xmag = clget$t ("xmag") + ymag = clget$t ("ymag") + xrot = clget$t ("xrotation") + yrot = clget$t ("yrotation") + xout = clget$t ("xout") + yout = clget$t ("yout") + xref = clget$t ("xref") + yref = clget$t ("yref") + xshift = clget$t ("xshift") + yshift = clget$t ("yshift") + + # Get set to adjust linear part of the fit. +$if (datatype == r) + call gscopy (sx1, newsx1) + call gscopy (sy1, newsy1) +$else + call dgscopy (sx1, newsx1) + call dgscopy (sy1, newsy1) +$endif + + if (geometry == GEO_DISTORTION) + call geo_rotmag$t (newsx1, newsy1, PIXEL(1.0), PIXEL(1.0), + PIXEL(0.0), PIXEL(0.0)) + else if (! IS_$INDEF$T(xmag) || ! IS_$INDEF$T(ymag) || + ! IS_$INDEF$T(xrot) || ! IS_$INDEF$T(yrot)) + call geo_drotmag$t (dt, rec, newsx1, newsy1, xmag, ymag, + xrot, yrot) + call geo_dxyshift$t (dt, rec, newsx1, newsy1, xout, yout, xref, yref, + xshift, yshift) +$if (datatype == r) + call gssave (newsx1, Mem$t[xcoeff]) + call gssave (newsy1, Mem$t[ycoeff]) +$else + call dgssave (newsx1, Mem$t[xcoeff]) + call dgssave (newsy1, Mem$t[ycoeff]) +$endif + + # Get distortion part of fit. + ncoeff = dtgeti (dt, rec, "surface2") + if (ncoeff > 0 && (geometry == GEO_GEOMETRIC || + geometry == GEO_DISTORTION)) { + + call realloc (xcoeff, ncoeff, TY_PIXEL) + call realloc (ycoeff, ncoeff, TY_PIXEL) + do i = 1, ncoeff { + junk = dtscan (dt) + call garg$t (Mem$t[xcoeff+i-1]) + call garg$t (Mem$t[ycoeff+i-1]) + } + + # Restore distortion part of fit. +$if (datatype == r) + iferr { + call gsrestore (sx2, Mem$t[xcoeff]) + } then { + call mfree (sx2, TY_STRUCT) + sx2 = NULL + } + iferr { + call gsrestore (sy2, Mem$t[ycoeff]) + } then { + call mfree (sy2, TY_STRUCT) + sy2 = NULL + } +$else + iferr { + call dgsrestore (sx2, Mem$t[xcoeff]) + } then { + call mfree (sx2, TY_STRUCT) + sx2 = NULL + } + iferr { + call dgsrestore (sy2, Mem$t[ycoeff]) + } then { + call mfree (sy2, TY_STRUCT) + sy2 = NULL + } +$endif + + } else { + sx2 = NULL + sy2 = NULL + } + + # Redefine the linear surfaces. +$if (datatype == r) + call gsfree (sx1) + call gscopy (newsx1, sx1) + call gsfree (newsx1) + call gsfree (sy1) + call gscopy (newsy1, sy1) + call gsfree (newsy1) +$else + call dgsfree (sx1) + call dgscopy (newsx1, sx1) + call dgsfree (newsx1) + call dgsfree (sy1) + call dgscopy (newsy1, sy1) + call dgsfree (newsy1) +$endif + + # Cleanup. + call mfree (xcoeff, TY_PIXEL) + call mfree (ycoeff, TY_PIXEL) +end + + +# GEO_DO_TRANSFORM -- The linear transformation is performed in this procedure. +# First the coordinates are scaled, then rotated and translated. The +# transformed coordinates are returned. + +procedure geo_do_transform$t (x, y, xt, yt, sx1, sy1, sx2, sy2) + +PIXEL x, y # initial positions +PIXEL xt, yt # transformed positions +pointer sx1, sy1 # pointer to linear surfaces +pointer sx2, sy2 # pointer to distortion surfaces + +$if (datatype == r) +real gseval() +$else +double dgseval() +$endif + +begin +$if (datatype == r) + xt = gseval (sx1, x, y) + if (sx2 != NULL) + xt = xt + gseval (sx2, x, y) + yt = gseval (sy1, x, y) + if (sy2 != NULL) + yt = yt + gseval (sy2, x, y) +$else + xt = dgseval (sx1, x, y) + if (sx2 != NULL) + xt = xt + dgseval (sx2, x, y) + yt = dgseval (sy1, x, y) + if (sy2 != NULL) + yt = yt + dgseval (sy2, x, y) +$endif +end + +$endfor diff --git a/pkg/images/immatch/src/geometry/geoxytran.x b/pkg/images/immatch/src/geometry/geoxytran.x new file mode 100644 index 00000000..e8bb9f64 --- /dev/null +++ b/pkg/images/immatch/src/geometry/geoxytran.x @@ -0,0 +1,446 @@ +include <ctype.h> +include <mach.h> +include <math.h> +include <math/gsurfit.h> + +define GEO_LINEAR 1 # Linear transformation only +define GEO_DISTORTION 2 # Distortion correction only +define GEO_GEOMETRIC 3 # Full transformation + + + +# GEO_LINIT -- Initialize the linear part of the transformation. + +procedure geo_linitr (sx1, sy1, sx2, sy2) + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +pointer sx2, sy2 #I/O pointer to the distortion x and y surfaces + +real xmag, ymag, xrot, yrot, xref, yref, xout, yout, xshift, yshift +real clgetr(), gseval() + +begin + # Initialize the surfaces. + call gsinit (sx1, GS_POLYNOMIAL, 2, 2, GS_XNONE, -MAX_REAL, MAX_REAL, + -MAX_REAL, MAX_REAL) + call gsinit (sy1, GS_POLYNOMIAL, 2, 2, GS_XNONE, -MAX_REAL, MAX_REAL, + -MAX_REAL, MAX_REAL) + sx2 = NULL + sy2 = NULL + + # Get the magnification parameters. + xmag = clgetr ("xmag") + if (IS_INDEFR(xmag)) + xmag = real(1.0) + ymag = clgetr ("ymag") + if (IS_INDEFR(ymag)) + ymag = real(1.0) + + # Get the rotation parameters. + xrot = clgetr ("xrot") + if (IS_INDEFR(xrot)) + xrot = real(0.0) + xrot = -DEGTORAD(xrot) + yrot = clgetr ("yrot") + if (IS_INDEFR(yrot)) + yrot = real(0.0) + yrot = -DEGTORAD(yrot) + + # Set the magnification and rotation coefficients. + call geo_rotmagr (sx1, sy1, xmag, ymag, xrot, yrot) + + # Compute the origin of the reference coordinates. + xref = clgetr ("xref") + if (IS_INDEFR(xref)) + xref = real(0.0) + yref = clgetr ("yref") + if (IS_INDEFR(yref)) + yref = real(0.0) + + # Compute the corresponding input coordinates. + xout = clgetr ("xout") + if (IS_INDEFR(xout)) + xout = gseval (sx1, xref, yref) + yout = clgetr ("yout") + if (IS_INDEFR(yout)) + yout = gseval (sy1, xref, yref) + + # Set the shifts. + xshift = clgetr ("xshift") + yshift = clgetr ("yshift") + if (IS_INDEFR(xshift)) + xshift = xout - gseval (sx1, xref, yref) + if (IS_INDEFR(yshift)) + yshift = yout - gseval (sy1, xref, yref) + call geo_xyshiftr (sx1, sy1, xshift, yshift) +end + + +# GEO_SFREE -- Free the x and y surface fitting descriptors. + +procedure geo_sfreer (sx1, sy1, sx2, sy2) + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +pointer sx2, sy2 #I/O pointer to the distortion x and y surfaces + +begin + call gsfree (sx1) + call gsfree (sy1) + if (sx2 != NULL) + call gsfree (sx2) + if (sy2 != NULL) + call gsfree (sy2) +end + + +# GEO_SINIT -- Read the surface fits from the database file and make +# any requested changes. + +procedure geo_sinitr (dt, record, geometry, sx1, sy1, sx2, sy2) + +pointer dt #I pointer to database file produced by geomap +char record[ARB] #I the name of the database record +int geometry #I the type of geometry to be computed +pointer sx1, sy1 #O pointers to the linear x and y surfaces +pointer sx2, sy2 #O pointers to the x and y distortion surfaces + +int i, rec, ncoeff, junk +real xmag, ymag, xrot, yrot, xref, yref, xout, yout, xshift, yshift +pointer newsx1, newsy1, xcoeff, ycoeff +int dtlocate(), dtscan(), dtgeti() +real clgetr() +errchk gsrestore + +begin + # Locate record. + rec = dtlocate (dt, record) + + # Get linear part of fit. + ncoeff = dtgeti (dt, rec, "surface1") + call malloc (xcoeff, ncoeff, TY_REAL) + call malloc (ycoeff, ncoeff, TY_REAL) + do i = 1, ncoeff { + junk = dtscan (dt) + call gargr (Memr[xcoeff+i-1]) + call gargr (Memr[ycoeff+i-1]) + } + + # Restore linear part of fit. + call gsrestore (sx1, Memr[xcoeff]) + call gsrestore (sy1, Memr[ycoeff]) + + # Get geometric transformation. + xmag = clgetr ("xmag") + ymag = clgetr ("ymag") + xrot = clgetr ("xrotation") + yrot = clgetr ("yrotation") + xout = clgetr ("xout") + yout = clgetr ("yout") + xref = clgetr ("xref") + yref = clgetr ("yref") + xshift = clgetr ("xshift") + yshift = clgetr ("yshift") + + # Get set to adjust linear part of the fit. + call gscopy (sx1, newsx1) + call gscopy (sy1, newsy1) + + if (geometry == GEO_DISTORTION) + call geo_rotmagr (newsx1, newsy1, real(1.0), real(1.0), + real(0.0), real(0.0)) + else if (! IS_INDEFR(xmag) || ! IS_INDEFR(ymag) || + ! IS_INDEFR(xrot) || ! IS_INDEFR(yrot)) + call geo_drotmagr (dt, rec, newsx1, newsy1, xmag, ymag, + xrot, yrot) + call geo_dxyshiftr (dt, rec, newsx1, newsy1, xout, yout, xref, yref, + xshift, yshift) + call gssave (newsx1, Memr[xcoeff]) + call gssave (newsy1, Memr[ycoeff]) + + # Get distortion part of fit. + ncoeff = dtgeti (dt, rec, "surface2") + if (ncoeff > 0 && (geometry == GEO_GEOMETRIC || + geometry == GEO_DISTORTION)) { + + call realloc (xcoeff, ncoeff, TY_REAL) + call realloc (ycoeff, ncoeff, TY_REAL) + do i = 1, ncoeff { + junk = dtscan (dt) + call gargr (Memr[xcoeff+i-1]) + call gargr (Memr[ycoeff+i-1]) + } + + # Restore distortion part of fit. + iferr { + call gsrestore (sx2, Memr[xcoeff]) + } then { + call mfree (sx2, TY_STRUCT) + sx2 = NULL + } + iferr { + call gsrestore (sy2, Memr[ycoeff]) + } then { + call mfree (sy2, TY_STRUCT) + sy2 = NULL + } + + } else { + sx2 = NULL + sy2 = NULL + } + + # Redefine the linear surfaces. + call gsfree (sx1) + call gscopy (newsx1, sx1) + call gsfree (newsx1) + call gsfree (sy1) + call gscopy (newsy1, sy1) + call gsfree (newsy1) + + # Cleanup. + call mfree (xcoeff, TY_REAL) + call mfree (ycoeff, TY_REAL) +end + + +# GEO_DO_TRANSFORM -- The linear transformation is performed in this procedure. +# First the coordinates are scaled, then rotated and translated. The +# transformed coordinates are returned. + +procedure geo_do_transformr (x, y, xt, yt, sx1, sy1, sx2, sy2) + +real x, y # initial positions +real xt, yt # transformed positions +pointer sx1, sy1 # pointer to linear surfaces +pointer sx2, sy2 # pointer to distortion surfaces + +real gseval() + +begin + xt = gseval (sx1, x, y) + if (sx2 != NULL) + xt = xt + gseval (sx2, x, y) + yt = gseval (sy1, x, y) + if (sy2 != NULL) + yt = yt + gseval (sy2, x, y) +end + + + +# GEO_LINIT -- Initialize the linear part of the transformation. + +procedure geo_linitd (sx1, sy1, sx2, sy2) + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +pointer sx2, sy2 #I/O pointer to the distortion x and y surfaces + +double xmag, ymag, xrot, yrot, xref, yref, xout, yout, xshift, yshift +double clgetd(), dgseval() + +begin + # Initialize the surfaces. + call dgsinit (sx1, GS_POLYNOMIAL, 2, 2, GS_XNONE, double (-MAX_REAL), + double (MAX_REAL), double (-MAX_REAL), double (MAX_REAL)) + call dgsinit (sy1, GS_POLYNOMIAL, 2, 2, GS_XNONE, double (-MAX_REAL), + double (MAX_REAL), double (-MAX_REAL), double (MAX_REAL)) + sx2 = NULL + sy2 = NULL + + # Get the magnification parameters. + xmag = clgetd ("xmag") + if (IS_INDEFD(xmag)) + xmag = double(1.0) + ymag = clgetd ("ymag") + if (IS_INDEFD(ymag)) + ymag = double(1.0) + + # Get the rotation parameters. + xrot = clgetd ("xrot") + if (IS_INDEFD(xrot)) + xrot = double(0.0) + xrot = -DEGTORAD(xrot) + yrot = clgetd ("yrot") + if (IS_INDEFD(yrot)) + yrot = double(0.0) + yrot = -DEGTORAD(yrot) + + # Set the magnification and rotation coefficients. + call geo_rotmagd (sx1, sy1, xmag, ymag, xrot, yrot) + + # Compute the origin of the reference coordinates. + xref = clgetd ("xref") + if (IS_INDEFD(xref)) + xref = double(0.0) + yref = clgetd ("yref") + if (IS_INDEFD(yref)) + yref = double(0.0) + + # Compute the corresponding input coordinates. + xout = clgetd ("xout") + if (IS_INDEFD(xout)) + xout = dgseval (sx1, xref, yref) + yout = clgetd ("yout") + if (IS_INDEFD(yout)) + yout = dgseval (sy1, xref, yref) + + # Set the shifts. + xshift = clgetd ("xshift") + yshift = clgetd ("yshift") + if (IS_INDEFD(xshift)) + xshift = xout - dgseval (sx1, xref, yref) + if (IS_INDEFD(yshift)) + yshift = yout - dgseval (sy1, xref, yref) + call geo_xyshiftd (sx1, sy1, xshift, yshift) +end + + +# GEO_SFREE -- Free the x and y surface fitting descriptors. + +procedure geo_sfreed (sx1, sy1, sx2, sy2) + +pointer sx1, sy1 #I/O pointers to the linear x and y surfaces +pointer sx2, sy2 #I/O pointer to the distortion x and y surfaces + +begin + call dgsfree (sx1) + call dgsfree (sy1) + if (sx2 != NULL) + call dgsfree (sx2) + if (sy2 != NULL) + call dgsfree (sy2) +end + + +# GEO_SINIT -- Read the surface fits from the database file and make +# any requested changes. + +procedure geo_sinitd (dt, record, geometry, sx1, sy1, sx2, sy2) + +pointer dt #I pointer to database file produced by geomap +char record[ARB] #I the name of the database record +int geometry #I the type of geometry to be computed +pointer sx1, sy1 #O pointers to the linear x and y surfaces +pointer sx2, sy2 #O pointers to the x and y distortion surfaces + +int i, rec, ncoeff, junk +double xmag, ymag, xrot, yrot, xref, yref, xout, yout, xshift, yshift +pointer newsx1, newsy1, xcoeff, ycoeff +int dtlocate(), dtscan(), dtgeti() +double clgetd() +errchk dgsrestore + +begin + # Locate record. + rec = dtlocate (dt, record) + + # Get linear part of fit. + ncoeff = dtgeti (dt, rec, "surface1") + call malloc (xcoeff, ncoeff, TY_DOUBLE) + call malloc (ycoeff, ncoeff, TY_DOUBLE) + do i = 1, ncoeff { + junk = dtscan (dt) + call gargd (Memd[xcoeff+i-1]) + call gargd (Memd[ycoeff+i-1]) + } + + # Restore linear part of fit. + call dgsrestore (sx1, Memd[xcoeff]) + call dgsrestore (sy1, Memd[ycoeff]) + + # Get geometric transformation. + xmag = clgetd ("xmag") + ymag = clgetd ("ymag") + xrot = clgetd ("xrotation") + yrot = clgetd ("yrotation") + xout = clgetd ("xout") + yout = clgetd ("yout") + xref = clgetd ("xref") + yref = clgetd ("yref") + xshift = clgetd ("xshift") + yshift = clgetd ("yshift") + + # Get set to adjust linear part of the fit. + call dgscopy (sx1, newsx1) + call dgscopy (sy1, newsy1) + + if (geometry == GEO_DISTORTION) + call geo_rotmagd (newsx1, newsy1, double(1.0), double(1.0), + double(0.0), double(0.0)) + else if (! IS_INDEFD(xmag) || ! IS_INDEFD(ymag) || + ! IS_INDEFD(xrot) || ! IS_INDEFD(yrot)) + call geo_drotmagd (dt, rec, newsx1, newsy1, xmag, ymag, + xrot, yrot) + call geo_dxyshiftd (dt, rec, newsx1, newsy1, xout, yout, xref, yref, + xshift, yshift) + call dgssave (newsx1, Memd[xcoeff]) + call dgssave (newsy1, Memd[ycoeff]) + + # Get distortion part of fit. + ncoeff = dtgeti (dt, rec, "surface2") + if (ncoeff > 0 && (geometry == GEO_GEOMETRIC || + geometry == GEO_DISTORTION)) { + + call realloc (xcoeff, ncoeff, TY_DOUBLE) + call realloc (ycoeff, ncoeff, TY_DOUBLE) + do i = 1, ncoeff { + junk = dtscan (dt) + call gargd (Memd[xcoeff+i-1]) + call gargd (Memd[ycoeff+i-1]) + } + + # Restore distortion part of fit. + iferr { + call dgsrestore (sx2, Memd[xcoeff]) + } then { + call mfree (sx2, TY_STRUCT) + sx2 = NULL + } + iferr { + call dgsrestore (sy2, Memd[ycoeff]) + } then { + call mfree (sy2, TY_STRUCT) + sy2 = NULL + } + + } else { + sx2 = NULL + sy2 = NULL + } + + # Redefine the linear surfaces. + call dgsfree (sx1) + call dgscopy (newsx1, sx1) + call dgsfree (newsx1) + call dgsfree (sy1) + call dgscopy (newsy1, sy1) + call dgsfree (newsy1) + + # Cleanup. + call mfree (xcoeff, TY_DOUBLE) + call mfree (ycoeff, TY_DOUBLE) +end + + +# GEO_DO_TRANSFORM -- The linear transformation is performed in this procedure. +# First the coordinates are scaled, then rotated and translated. The +# transformed coordinates are returned. + +procedure geo_do_transformd (x, y, xt, yt, sx1, sy1, sx2, sy2) + +double x, y # initial positions +double xt, yt # transformed positions +pointer sx1, sy1 # pointer to linear surfaces +pointer sx2, sy2 # pointer to distortion surfaces + +double dgseval() + +begin + xt = dgseval (sx1, x, y) + if (sx2 != NULL) + xt = xt + dgseval (sx2, x, y) + yt = dgseval (sy1, x, y) + if (sy2 != NULL) + yt = yt + dgseval (sy2, x, y) +end + + diff --git a/pkg/images/immatch/src/geometry/mkpkg b/pkg/images/immatch/src/geometry/mkpkg new file mode 100644 index 00000000..e6e98b24 --- /dev/null +++ b/pkg/images/immatch/src/geometry/mkpkg @@ -0,0 +1,34 @@ +# Make the GEOMAP/GEOXYTRAN and CCMAP/CCSETWCS/CCTRAN tasks + +$checkout libpkg.a ../../../ +$update libpkg.a +$checkin libpkg.a ../../../ +$exit + +generic: + $set GEN = "$$generic -k" + + $ifolder (geofunc.x, geofunc.gx) + $(GEN) geofunc.gx -o geofunc.x $endif + $ifolder (t_geomap.x, t_geomap.gx) + $(GEN) t_geomap.gx -o t_geomap.x $endif + $ifolder (geoxytran.x,geoxytran.gx) + $(GEN) geoxytran.gx -o geoxytran.x $endif + ; + +libpkg.a: + $ifeq (USE_GENERIC, yes) $call generic $endif + + geofunc.x <math.h> <math/gsurfit.h> + geotimtran.x <imhdr.h> <imset.h> <mach.h> <math/gsurfit.h> \ + <math/iminterp.h> geotran.h + geotran.x <imhdr.h> <imset.h> <mach.h> <math/gsurfit.h> \ + <math/iminterp.h> geotran.h + geoxytran.x <mach.h> <ctype.h> <math.h> <math/gsurfit.h> + t_geomap.x <fset.h> <error.h> <mach.h> <math/gsurfit.h> \ + <math.h> "../../../lib/geomap.h" + t_geotran.x <imhdr.h> <mwset.h> <math.h> <math/gsurfit.h> \ + geotran.h + t_geoxytran.x <fset.h> <ctype.h> + trinvert.x + ; diff --git a/pkg/images/immatch/src/geometry/t_geomap.gx b/pkg/images/immatch/src/geometry/t_geomap.gx new file mode 100644 index 00000000..02d530e5 --- /dev/null +++ b/pkg/images/immatch/src/geometry/t_geomap.gx @@ -0,0 +1,921 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <fset.h> +include <error.h> +include <mach.h> +include <math.h> +include <math/gsurfit.h> +include "../../../lib/geomap.h" + +define GM_REAL 1 # computation type is real +define GM_DOUBLE 2 # computation type is double + +$for (r) + +# T_GEOMAP -- Procedure to calculate the transformation required to transform +# the coordinate system of a reference image to the coordinate system of +# an input image. The transformation is of the following form. +# +# xin = f (xref, yref) +# yin = g (xref, yref) + +procedure t_geomap () + +bool verbose, interactive +double xmin, xmax, ymin, ymax, reject +int geometry, function, calctype, nfiles, list, in, reclist, nrecords +int xxorder, xyorder, xxterms, yxorder, yyorder, yxterms, maxiter +int reslist, nresfiles, res +pointer sp, in_name, str, out, fit, gd, graphics +real rxmin, rxmax, rymin, rymax + +bool clgetb() +double clgetd() +int clgeti(), clgwrd(), clplen(), errget(), imtopenp(), imtlen() +int imtgetim() +pointer clpopnu(), clgfil(), dtmap(), gopen(), open() + +errchk geo_mapr(), geo_mapd() + +begin + # Get working space. + call smark (sp) + call salloc (in_name, SZ_FNAME, TY_CHAR) + call salloc (graphics, SZ_FNAME, TY_CHAR) + call salloc (str, max(SZ_LINE, SZ_FNAME), TY_CHAR) + + # Get input data file(s). + list = clpopnu ("input") + nfiles = clplen (list) + + # Open database output file. + call clgstr ("database", Memc[str], SZ_FNAME) + out = dtmap (Memc[str], APPEND) + + # Get minimum and maximum reference values. + xmin = clgetd ("xmin") + if (IS_INDEFD(xmin)) + rxmin = INDEFR + else + rxmin = xmin + xmax = clgetd ("xmax") + if (IS_INDEFD(xmax)) + rxmax = INDEFR + else + rxmax = xmax + ymin = clgetd ("ymin") + if (IS_INDEFD(ymin)) + rymin = INDEFR + else + rymin = ymin + ymax = clgetd ("ymax") + if (IS_INDEFD(ymax)) + rymax = INDEFR + else + rymax = ymax + + # Get the records list. + reclist = imtopenp ("transforms") + nrecords = imtlen (reclist) + if ((nrecords > 0) && (nrecords != nfiles)) { + call eprintf ( + "The number of records is not equal to the number of input files") + call clpcls (list) + call dtunmap (out) + call imtclose (reclist) + call sfree (sp) + return + } + + # Get the results file list. + reslist = clpopnu ("results") + nresfiles = clplen (reslist) + if (nresfiles > 1 && nresfiles != nfiles) { + call eprintf ("Error: there are too few results files\n") + call clpcls (list) + call dtunmap (out) + call imtclose (reclist) + call clpcls (reslist) + call sfree (sp) + return + } + + # Get the surface fitting parameters. + geometry = clgwrd ("fitgeometry", Memc[str], SZ_LINE, GM_GEOMETRIES) + function = clgwrd ("function", Memc[str], SZ_LINE, GM_FUNCS) + xxorder = clgeti ("xxorder") + xyorder = clgeti ("xyorder") + xxterms = clgwrd ("xxterms", Memc[str], SZ_LINE, GM_XFUNCS) - 1 + yxorder = clgeti ("yxorder") + yyorder = clgeti ("yyorder") + yxterms = clgwrd ("yxterms", Memc[str], SZ_LINE, GM_XFUNCS) - 1 + maxiter = clgeti ("maxiter") + reject = clgetd ("reject") + calctype = clgwrd ("calctype", Memc[str], SZ_LINE, ",real,double,") + + # Get the graphics parameters. + verbose = clgetb ("verbose") + interactive = clgetb ("interactive") + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + + # Flush standard output on newline. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Initialize the fit structure. + call geo_minit (fit, GM_NONE, geometry, function, xxorder, xyorder, + xxterms, yxorder, yyorder, yxterms, maxiter, reject) + + # Loop over the files. + while (clgfil (list, Memc[in_name], SZ_FNAME) != EOF) { + + # Open text file of coordinates. + in = open (Memc[in_name], READ_ONLY, TEXT_FILE) + + # Open the results files. + if (nresfiles <= 0) + res = NULL + else if (clgfil (reslist, Memc[str], SZ_FNAME) != EOF) + res = open (Memc[str], NEW_FILE, TEXT_FILE) + + # Set file name in structure. + if (nrecords > 0) { + if (imtgetim (reclist, GM_RECORD(fit), SZ_FNAME) != EOF) + ; + } else + call strcpy (Memc[in_name], GM_RECORD(fit), SZ_FNAME) + + if (verbose && res != STDOUT) { + call fstats (in, F_FILENAME, Memc[str], SZ_FNAME) + call printf ("\nCoordinate list: %s Transform: %s\n") + call pargstr (Memc[str]) + call pargstr (GM_RECORD(fit)) + if (res != NULL) + call fstats (res, F_FILENAME, Memc[str], SZ_FNAME) + else + call strcpy ("", Memc[str], SZ_FNAME) + call printf (" Results file: %s\n") + call pargstr (Memc[str]) + call flush (STDOUT) + } + if (res != NULL) { + call fstats (in, F_FILENAME, Memc[str], SZ_FNAME) + call fprintf (res, "\n# Coordinate list: %s Transform: %s\n") + call pargstr (Memc[str]) + call pargstr (GM_RECORD(fit)) + if (res != NULL) + call fstats (res, F_FILENAME, Memc[str], SZ_FNAME) + else + call strcpy ("", Memc[str], SZ_FNAME) + call fprintf (res, "# Results file: %s\n") + call pargstr (Memc[str]) + call flush (STDOUT) + } + + if (interactive) { + gd = gopen (Memc[graphics], NEW_FILE, STDGRAPH) + } else + gd = NULL + + iferr { + if (calctype == GM_REAL) + call geo_mapr (gd, in, out, res, fit, rxmin, rxmax, rymin, + rymax, verbose) + else + call geo_mapd (gd, in, out, res, fit, xmin, xmax, ymin, + ymax, verbose) + } then { + if (verbose && res != STDOUT) { + call printf ("Error fitting coordinate list: %s\n") + call pargstr (Memc[in_name]) + call flush (STDOUT) + if (errget (Memc[str], SZ_LINE) == 0) + ; + call printf ("\t%s\n") + call pargstr (Memc[str)) + } + if (res != NULL) { + call fprintf (res, "# Error fitting coordinate list: %s\n") + call pargstr (Memc[in_name]) + call flush (STDOUT) + if (errget (Memc[str], SZ_LINE) == 0) + ; + call fprintf (res, "# %s\n") + call pargstr (Memc[str)) + } + } + + call close (in) + if (nresfiles == nfiles) + call close ( res) + + if (gd != NULL) + call gclose (gd) + } + + # Close up. + call geo_free (fit) + if (nresfiles < nfiles) + call close ( res) + call dtunmap (out) + call imtclose (reclist) + call clpcls (list) + call sfree (sp) +end + +$endfor + +$for (rd) + +# GEO_MAP -- Procedure to calculate the coordinate transformations + +procedure geo_map$t (gd, in, out, res, fit, xmin, xmax, ymin, ymax, verbose) + +pointer gd #I the graphics stream +int in #I the input file descriptor +pointer out #I the output file descriptor +int res #I the results file descriptor +pointer fit #I pointer to fit parameters +PIXEL xmin, xmax #I max and min xref values +PIXEL ymin, ymax #I max and min yref values +bool verbose #I verbose mode + +int npts, ngood +pointer sp, str, xref, yref, xin, yin, wts, xfit, yfit, xerrmsg, yerrmsg +pointer sx1, sy1, sx2, sy2 +PIXEL mintemp, maxtemp + +PIXEL asum$t() +int geo_rdxy$t() +errchk geo_fit$t, geo_mgfit$t() + +begin + # Get working space. + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (xerrmsg, SZ_LINE, TY_CHAR) + call salloc (yerrmsg, SZ_LINE, TY_CHAR) + + # Initialize pointers. + xref = NULL + yref = NULL + xin = NULL + yin = NULL + wts = NULL + + # Read in data and check that data is in range. + npts = geo_rdxy$t (in, xref, yref, xin, yin, xmin, xmax, ymin, ymax) + if (npts <= 0) { + call fstats (in, F_FILENAME, Memc[str], SZ_FNAME) + call printf ("Coordinate list: %s has no data in range.\n") + call pargstr (Memc[str]) + call sfree (sp) + return + } + + # Compute the mean of the reference and input coordinates. + GM_XOREF(fit) = double (asum$t (Mem$t[xref], npts) / npts) + GM_YOREF(fit) = double (asum$t (Mem$t[yref], npts) / npts) + GM_XOIN(fit) = double (asum$t (Mem$t[xin], npts) / npts) + GM_YOIN(fit) = double (asum$t (Mem$t[yin], npts) / npts) + + # Set the reference point for the projections to INDEF. + GM_XREFPT(fit) = INDEFD + GM_YREFPT(fit) = INDEFD + + # Compute the weights. + call malloc (xfit, npts, TY_PIXEL) + call malloc (yfit, npts, TY_PIXEL) + call malloc (wts, npts, TY_PIXEL) + call amovk$t (PIXEL(1.), Mem$t[wts], npts) + + # Determine the x max and min. + if (IS_$INDEF$T(xmin) || IS_$INDEF$T(xmax)) { + call alim$t (Mem$t[xref], npts, mintemp, maxtemp) + if (! IS_$INDEF$T(xmin)) + GM_XMIN(fit) = double (xmin) + else + GM_XMIN(fit) = double (mintemp) + if (! IS_$INDEF$T(xmax)) + GM_XMAX(fit) = double (xmax) + else + GM_XMAX(fit) = double (maxtemp) + } else { + GM_XMIN(fit) = double (xmin) + GM_XMAX(fit) = double (xmax) + } + + # Determine the y max and min. + if (IS_$INDEF$T(ymin) || IS_$INDEF$T(ymax)) { + call alim$t (Mem$t[yref], npts, mintemp, maxtemp) + if (! IS_$INDEF$T(ymin)) + GM_YMIN(fit) = double (ymin) + else + GM_YMIN(fit) = double (mintemp) + if (! IS_$INDEF$T(ymax)) + GM_YMAX(fit) = double (ymax) + else + GM_YMAX(fit) = double (maxtemp) + } else { + GM_YMIN(fit) = double (ymin) + GM_YMAX(fit) = double (ymax) + } + + # Initalize surface pointers. + sx1 = NULL + sy1 = NULL + sx2 = NULL + sy2 = NULL + + # Fit the data. + if (gd != NULL) { + iferr { + call geo_mgfit$t (gd, fit, sx1, sy1, sx2, sy2, Mem$t[xref], + Mem$t[yref], Mem$t[xin], Mem$t[yin], Mem$t[wts], npts, + Memc[xerrmsg], Memc[yerrmsg], SZ_LINE) + } then { + call gdeactivate (gd, 0) + call mfree (xfit, TY_PIXEL) + call mfree (yfit, TY_PIXEL) + call mfree (wts, TY_PIXEL) + call geo_mmfree$t (sx1, sy1, sx2, sy2) + call sfree (sp) + call error (3, "Too few points for X or Y fits.") + } + call gdeactivate (gd, 0) + if (verbose && res != STDOUT) { + call printf ("Coordinate mapping status\n") + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "# Coordinate mapping status\n") + } + } else { + if (verbose && res != STDOUT) { + call printf ("Coordinate mapping status\n ") + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "# Coordinate mapping status\n# ") + } + iferr { + call geo_fit$t (fit, sx1, sy1, sx2, sy2, Mem$t[xref], + Mem$t[yref], Mem$t[xin], Mem$t[yin], Mem$t[wts], npts, + Memc[xerrmsg], Memc[yerrmsg], SZ_LINE) + } then { + call mfree (xfit, TY_PIXEL) + call mfree (yfit, TY_PIXEL) + call mfree (wts, TY_PIXEL) + call geo_mmfree$t (sx1, sy1, sx2, sy2) + call sfree (sp) + call error (3, "Too few points for X or Y fits.") + } + if (verbose && res != STDOUT) { + call printf ("%s %s\n") + call pargstr (Memc[xerrmsg]) + call pargstr (Memc[yerrmsg]) + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "%s %s\n") + call pargstr (Memc[xerrmsg]) + call pargstr (Memc[yerrmsg]) + call flush (STDOUT) + } + } + ngood = GM_NPTS(fit) - GM_NWTS0(fit) + if (verbose && res != STDOUT) { + call printf (" Xin and Yin fit rms: %0.7g %0.7g\n") + if (ngood <= 1) { + call pargd (0.0d0) + call pargd (0.0d0) + } else { + call pargd (sqrt (GM_XRMS(fit) / (ngood - 1))) + call pargd (sqrt (GM_YRMS(fit) / (ngood - 1))) + } + call geo_show$t (STDOUT, fit, sx1, sy1, NO) + } + if (res != NULL) { + call fprintf (res, "# Xin and Yin fit rms: %0.7g %0.7g\n") + if (ngood <= 1) { + call pargd (0.0) + call pargd (0.0) + } else { + call pargd (sqrt (GM_XRMS(fit) / (ngood - 1))) + call pargd (sqrt (GM_YRMS(fit) / (ngood - 1))) + } + call geo_show$t (res, fit, sx1, sy1, YES) + } + + # Compute and print the fitted x and y values. + if (res != NULL) { + call geo_eval$t (sx1, sy1, sx2, sy2, Mem$t[xref], Mem$t[yref], + Mem$t[xfit], Mem$t[yfit], npts) + call geo_plist$t (res, fit, Mem$t[xref], Mem$t[yref], Mem$t[xin], + Mem$t[yin], Mem$t[xfit], Mem$t[yfit], Mem$t[wts], npts) + } + + # Free the data + if (xref != NULL) + call mfree (xref, TY_PIXEL) + if (yref != NULL) + call mfree (yref, TY_PIXEL) + if (xin != NULL) + call mfree (xin, TY_PIXEL) + if (yin != NULL) + call mfree (yin, TY_PIXEL) + if (xfit != NULL) + call mfree (xfit, TY_PIXEL) + if (yfit != NULL) + call mfree (yfit, TY_PIXEL) + if (wts != NULL) + call mfree (wts, TY_PIXEL) + + # Output the data. + call geo_mout$t (fit, out, sx1, sy1, sx2, sy2) + + # Free the space and close files. + call geo_mmfree$t (sx1, sy1, sx2, sy2) + call sfree (sp) +end + + +define GEO_DEFBUFSIZE 1000 # default data buffer sizes + +# GEO_RDXY -- Read in the data points. + +int procedure geo_rdxy$t (fd, xref, yref, xin, yin, xmin, xmax, ymin, ymax) + +int fd # the input file descriptor +pointer xref # the x reference coordinates +pointer yref # the y reference coordinates +pointer xin # the x coordinates +pointer yin # the y coordinates +PIXEL xmin, xmax # the range of the x coordinates +PIXEL ymin, ymax # the range of the y coordinates + +int npts, bufsize +int fscan(), nscan() + +begin + bufsize = GEO_DEFBUFSIZE + call malloc (xref, bufsize, TY_PIXEL) + call malloc (yref, bufsize, TY_PIXEL) + call malloc (xin, bufsize, TY_PIXEL) + call malloc (yin, bufsize, TY_PIXEL) + + npts = 0 + while (fscan (fd) != EOF) { + + # Decode the data. + call garg$t (Mem$t[xref+npts]) + call garg$t (Mem$t[yref+npts]) + call garg$t (Mem$t[xin+npts]) + call garg$t (Mem$t[yin+npts]) + if (nscan() < 4) + next + + # Check the data limits. + if (! IS_$INDEF$T(xmin)) { + if (Mem$t[xref+npts] < xmin) + next + } + if (! IS_$INDEF$T(xmax)) { + if (Mem$t[xref+npts] > xmax) + next + } + if (! IS_$INDEF$T(ymin)) { + if (Mem$t[yref+npts] < ymin) + next + } + if (! IS_$INDEF$T(ymax)) { + if (Mem$t[yref+npts] > ymax) + next + } + + npts = npts + 1 + if (npts >= bufsize) { + bufsize = bufsize + GEO_DEFBUFSIZE + call realloc (xref, bufsize, TY_PIXEL) + call realloc (yref, bufsize, TY_PIXEL) + call realloc (xin, bufsize, TY_PIXEL) + call realloc (yin, bufsize, TY_PIXEL) + } + } + + if (npts <= 0) { + call mfree (xref, TY_PIXEL) + call mfree (yref, TY_PIXEL) + call mfree (xin, TY_PIXEL) + call mfree (yin, TY_PIXEL) + xref = NULL + yref = NULL + xin = NULL + yin = NULL + } else if (npts < bufsize) { + call realloc (xref, npts, TY_PIXEL) + call realloc (yref, npts, TY_PIXEL) + call realloc (xin, npts, TY_PIXEL) + call realloc (yin, npts, TY_PIXEL) + } + + return (npts) +end + + +# GEO_EVAL -- Evalute the fit. + +procedure geo_eval$t (sx1, sy1, sx2, sy2, xref, yref, xi, eta, npts) + +pointer sx1, sy1 #I pointer to linear surfaces +pointer sx2, sy2 #I pointer to higher order surfaces +PIXEL xref[ARB] #I the x reference coordinates +PIXEL yref[ARB] #I the y reference coordinates +PIXEL xi[ARB] #O the fitted xi coordinates +PIXEL eta[ARB] #O the fitted eta coordinates +int npts #I the number of points + +pointer sp, temp + +begin + call smark (sp) + call salloc (temp, npts, TY_PIXEL) + +$if (datatype == r) + call gsvector (sx1, xref, yref, xi, npts) +$else + call dgsvector (sx1, xref, yref, xi, npts) +$endif + if (sx2 != NULL) { +$if (datatype == r) + call gsvector (sx2, xref, yref, Mem$t[temp], npts) +$else + call dgsvector (sx2, xref, yref, Mem$t[temp], npts) +$endif + call aadd$t (Mem$t[temp], xi, xi, npts) + } +$if (datatype == r) + call gsvector (sy1, xref, yref, eta, npts) +$else + call dgsvector (sy1, xref, yref, eta, npts) +$endif + if (sy2 != NULL) { +$if (datatype == r) + call gsvector (sy2, xref, yref, Mem$t[temp], npts) +$else + call dgsvector (sy2, xref, yref, Mem$t[temp], npts) +$endif + + call aadd$t (Mem$t[temp], eta, eta, npts) + } + + call sfree (sp) +end + + +# GEO_MOUT -- Write the output database file. + +procedure geo_mout$t (fit, out, sx1, sy1, sx2, sy2) + +pointer fit #I pointer to fitting structure +int out #I pointer to database file +pointer sx1, sy1 #I pointer to linear surfaces +pointer sx2, sy2 #I pointer to distortion surfaces + +int i, npts, ncoeff +pointer sp, str, xcoeff, ycoeff +PIXEL xrms, yrms, xshift, yshift, xscale, yscale, xrot, yrot +$if (datatype == r) +int gsgeti() +$else +int dgsgeti() +$endif +int rg_wrdstr() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Compute the x and y fit rms. + #npts = max (0, GM_NPTS(fit) - GM_NREJECT(fit) - GM_NWTS0(fit)) + npts = max (0, GM_NPTS(fit) - GM_NWTS0(fit)) + xrms = max (0.0d0, GM_XRMS(fit)) + yrms = max (0.0d0, GM_YRMS(fit)) + if (npts > 1) { + xrms = sqrt (xrms / (npts - 1)) + yrms = sqrt (yrms / (npts - 1)) + } else { + xrms = 0.0d0 + yrms = 0.0d0 + } + + # Print title. + call dtptime (out) + call dtput (out, "begin\t%s\n") + call pargstr (GM_RECORD(fit)) + + # Print the x and y mean values. + call dtput (out, "\txrefmean\t%g\n") + call pargd (GM_XOREF(fit)) + call dtput (out, "\tyrefmean\t%g\n") + call pargd (GM_YOREF(fit)) + call dtput (out, "\txmean\t\t%g\n") + call pargd (GM_XOIN(fit)) + call dtput (out, "\tymean\t\t%g\n") + call pargd (GM_YOIN(fit)) + + # Print some of the fitting parameters. + if (rg_wrdstr (GM_FIT(fit), Memc[str], SZ_FNAME, GM_GEOMETRIES) <= 0) + call strcpy ("general", Memc[str], SZ_FNAME) + call dtput (out, "\tgeometry\t%s\n") + call pargstr (Memc[str]) + if (rg_wrdstr (GM_FUNCTION(fit), Memc[str], SZ_FNAME, GM_FUNCS) <= 0) + call strcpy ("polynomial", Memc[str], SZ_FNAME) + call dtput (out, "\tfunction\t%s\n") + call pargstr (Memc[str]) + + # Output the geometric parameters. + call geo_lcoeff$t (sx1, sy1, xshift, yshift, xscale, yscale, xrot, yrot) + call dtput (out, "\txshift\t\t%g\n") + call parg$t (xshift) + call dtput (out, "\tyshift\t\t%g\n") + call parg$t (yshift) + call dtput (out, "\txmag\t\t%g\n") + call parg$t (xscale) + call dtput (out, "\tymag\t\t%g\n") + call parg$t (yscale) + call dtput (out, "\txrotation\t%g\n") + call parg$t (xrot) + call dtput (out, "\tyrotation\t%g\n") + call parg$t (yrot) + + # Out the rms values. + call dtput (out, "\txrms\t\t%g\n") + call parg$t (PIXEL(xrms)) + call dtput (out, "\tyrms\t\t%g\n") + call parg$t (PIXEL(yrms)) + + # Allocate memory for linear coefficients. +$if (datatype == r) + ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE)) +$else + ncoeff = max (dgsgeti (sx1, GSNSAVE), dgsgeti (sy1, GSNSAVE)) +$endif + call calloc (xcoeff, ncoeff, TY_PIXEL) + call calloc (ycoeff, ncoeff, TY_PIXEL) + + # Output the linear coefficients. +$if (datatype == r) + call gssave (sx1, Mem$t[xcoeff]) + call gssave (sy1, Mem$t[ycoeff]) +$else + call dgssave (sx1, Mem$t[xcoeff]) + call dgssave (sy1, Mem$t[ycoeff]) +$endif + call dtput (out, "\tsurface1\t%d\n") + call pargi (ncoeff) + do i = 1, ncoeff { + call dtput (out, "\t\t\t%g\t%g\n") + call parg$t (Mem$t[xcoeff+i-1]) + call parg$t (Mem$t[ycoeff+i-1]) + } + + call mfree (xcoeff, TY_PIXEL) + call mfree (ycoeff, TY_PIXEL) + + # Allocate memory for higer order coefficients. + if (sx2 == NULL) + ncoeff = 0 + else +$if (datatype == r) + ncoeff = gsgeti (sx2, GSNSAVE) +$else + ncoeff = dgsgeti (sx2, GSNSAVE) +$endif + if (sy2 == NULL) + ncoeff = max (0, ncoeff) + else +$if (datatype == r) + ncoeff = max (gsgeti (sy2, GSNSAVE), ncoeff) +$else + ncoeff = max (dgsgeti (sy2, GSNSAVE), ncoeff) +$endif + call calloc (xcoeff, ncoeff, TY_PIXEL) + call calloc (ycoeff, ncoeff, TY_PIXEL) + + # Save the coefficients. +$if (datatype == r) + call gssave (sx2, Mem$t[xcoeff]) + call gssave (sy2, Mem$t[ycoeff]) +$else + call dgssave (sx2, Mem$t[xcoeff]) + call dgssave (sy2, Mem$t[ycoeff]) +$endif + + # Output the coefficients. + call dtput (out, "\tsurface2\t%d\n") + call pargi (ncoeff) + do i = 1, ncoeff { + call dtput (out, "\t\t\t%g\t%g\n") + call parg$t (Mem$t[xcoeff+i-1]) + call parg$t (Mem$t[ycoeff+i-1]) + } + + # Cleanup. + call mfree (xcoeff, TY_PIXEL) + call mfree (ycoeff, TY_PIXEL) + call sfree (sp) +end + + +# GEO_PLIST -- Print the input, output, and fitted data and the residuals. + +procedure geo_plist$t (fd, fit, xref, yref, xin, yin, xfit, yfit, wts, npts) + +int fd #I the results file descriptor +pointer fit #I pointer to the fit structure +PIXEL xref[ARB] #I the input x coordinates +PIXEL yref[ARB] #I the input y coordinates +PIXEL xin[ARB] #I the input ra / longitude coordinates +PIXEL yin[ARB] #I the input dec / latitude coordinates +PIXEL xfit[ARB] #I the fitted ra / longitude coordinates +PIXEL yfit[ARB] #I the fitted dec / latitude coordinates +PIXEL wts[ARB] #I the weights array +int npts #I the number of data points + +int i, index +pointer sp, fmtstr, twts + +begin + # Allocate working space. + call smark (sp) + call salloc (fmtstr, SZ_LINE, TY_CHAR) + call salloc (twts, npts, TY_PIXEL) + + # Compute the weights. + call amov$t (wts, Mem$t[twts], npts) + do i = 1, GM_NREJECT(fit) { + index = Memi[GM_REJ(fit)+i-1] + if (wts[index] > PIXEL(0.0)) + Mem$t[twts+index-1] = PIXEL(0.0) + } + + # Print banner. + call fprintf (fd, "\n# Input Coordinate Listing\n") + call fprintf (fd, "# Column 1: X (reference) \n") + call fprintf (fd, "# Column 2: Y (reference)\n") + call fprintf (fd, "# Column 3: X (input)\n") + call fprintf (fd, "# Column 4: Y (input)\n") + call fprintf (fd, "# Column 5: X (fit)\n") + call fprintf (fd, "# Column 6: Y (fit)\n") + call fprintf (fd, "# Column 7: X (residual)\n") + call fprintf (fd, "# Column 8: Y (residual)\n\n") + + # Create the format string. + call sprintf (Memc[fmtstr], SZ_LINE, "%s %s %s %s %s %s %s %s\n") +$if (datatype == r) + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") +$else + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") +$endif + + # Print the data. + do i = 1, npts { + call fprintf (fd, Memc[fmtstr]) + call parg$t (xref[i]) + call parg$t (yref[i]) + call parg$t (xin[i]) + call parg$t (yin[i]) + if (Mem$t[twts+i-1] > 0.0d0) { + call parg$t (xfit[i]) + call parg$t (yfit[i]) + call parg$t (xin[i] - xfit[i]) + call parg$t (yin[i] - yfit[i]) + } else { + call parg$t (INDEF) + call parg$t (INDEF) + call parg$t (INDEF) + call parg$t (INDEF) + } + + } + + call fprintf (fd, "\n") + + call sfree (sp) + +end + +# GEO_SHOW -- Print the coordinate mapping parameters. + +procedure geo_show$t (fd, fit, sx1, sy1, comment) + +int fd #I the output file descriptor +pointer fit #I pointer to the fit structure +pointer sx1, sy1 #I pointer to linear surfaces +int comment #I comment the output ? + +PIXEL xshift, yshift, a, b, c, d +PIXEL xscale, yscale, xrot, yrot +pointer sp, str +bool fp_equal$t() + +begin + # Allocate temporary space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Compute the geometric parameters. + call geo_gcoeff$t (sx1, sy1, xshift, yshift, a, b, c, d) + + if (comment == NO) { + call fprintf (fd, "Coordinate mapping parameters\n") + } else { + call fprintf (fd, "# Coordinate mapping parameters\n") + } + + if (comment == NO) { + call fprintf (fd, + " Mean Xref and Yref: %0.7g %0.7g\n") + call pargd (GM_XOREF(fit)) + call pargd (GM_YOREF(fit)) + call fprintf (fd, + " Mean Xin and Yin: %0.7g %0.7g\n") + call pargd (GM_XOIN(fit)) + call pargd (GM_YOIN(fit)) + call fprintf (fd, + " X and Y shift: %0.7g %0.7g (xin yin)\n") + call parg$t (xshift) + call parg$t (yshift) + } else { + call fprintf (fd, + "# Mean Xref and Yref: %0.7g %0.7g\n") + call pargd (GM_XOREF(fit)) + call pargd (GM_YOREF(fit)) + call fprintf (fd, + "# Mean Xin and Yin: %0.7g %g0.7\n") + call pargd (GM_XOIN(fit)) + call pargd (GM_YOIN(fit)) + call fprintf (fd, + "# X and Y shift: %0.7g %0.7g (xin yin)\n") + call parg$t (xshift) + call parg$t (yshift) + } + + # Output the scale factors. + xscale = sqrt (a * a + c * c) + yscale = sqrt (b * b + d * d) + if (comment == NO) { + call fprintf (fd, + " X and Y scale: %0.7g %0.7g (xin / xref yin / yref)\n") + call parg$t (xscale) + call parg$t (yscale) + } else { + call fprintf (fd, + "# X and Y scale: %0.7g %0.7g (xin / xref yin / yref)\n") + call parg$t (xscale) + call parg$t (yscale) + } + + # Output the rotation factors. + if (fp_equal$t (a, PIXEL(0.0)) && fp_equal$t (c, PIXEL(0.0))) + xrot = PIXEL(0.0) + else + xrot = RADTODEG (atan2 (-c, a)) + if (xrot < PIXEL(0.0)) + xrot = xrot + PIXEL(360.0) + if (fp_equal$t (b, PIXEL(0.0)) && fp_equal$t (d, PIXEL(0.0))) + yrot = PIXEL(0.0) + else + yrot = RADTODEG (atan2 (b, d)) + if (yrot < PIXEL(0.0)) + yrot = yrot + PIXEL(360.0) + if (comment == NO) { + call fprintf (fd, + " X and Y axis rotation: %0.5f %0.5f (degrees degrees)\n") + call parg$t (xrot) + call parg$t (yrot) + } else { + call fprintf (fd, + "# X and Y axis rotation: %0.5f %0.5f (degrees degrees)\n") + call parg$t (xrot) + call parg$t (yrot) + } + + call sfree (sp) +end + +$endfor diff --git a/pkg/images/immatch/src/geometry/t_geomap.x b/pkg/images/immatch/src/geometry/t_geomap.x new file mode 100644 index 00000000..6f1c20f0 --- /dev/null +++ b/pkg/images/immatch/src/geometry/t_geomap.x @@ -0,0 +1,1509 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <fset.h> +include <error.h> +include <mach.h> +include <math.h> +include <math/gsurfit.h> +include "../../../lib/geomap.h" + +define GM_REAL 1 # computation type is real +define GM_DOUBLE 2 # computation type is double + + + +# T_GEOMAP -- Procedure to calculate the transformation required to transform +# the coordinate system of a reference image to the coordinate system of +# an input image. The transformation is of the following form. +# +# xin = f (xref, yref) +# yin = g (xref, yref) + +procedure t_geomap () + +bool verbose, interactive +double xmin, xmax, ymin, ymax, reject +int geometry, function, calctype, nfiles, list, in, reclist, nrecords +int xxorder, xyorder, xxterms, yxorder, yyorder, yxterms, maxiter +int reslist, nresfiles, res +pointer sp, in_name, str, out, fit, gd, graphics +real rxmin, rxmax, rymin, rymax + +bool clgetb() +double clgetd() +int clgeti(), clgwrd(), clplen(), errget(), imtopenp(), imtlen() +int imtgetim() +pointer clpopnu(), clgfil(), dtmap(), gopen(), open() + +errchk geo_mapr(), geo_mapd() + +begin + # Get working space. + call smark (sp) + call salloc (in_name, SZ_FNAME, TY_CHAR) + call salloc (graphics, SZ_FNAME, TY_CHAR) + call salloc (str, max(SZ_LINE, SZ_FNAME), TY_CHAR) + + # Get input data file(s). + list = clpopnu ("input") + nfiles = clplen (list) + + # Open database output file. + call clgstr ("database", Memc[str], SZ_FNAME) + out = dtmap (Memc[str], APPEND) + + # Get minimum and maximum reference values. + xmin = clgetd ("xmin") + if (IS_INDEFD(xmin)) + rxmin = INDEFR + else + rxmin = xmin + xmax = clgetd ("xmax") + if (IS_INDEFD(xmax)) + rxmax = INDEFR + else + rxmax = xmax + ymin = clgetd ("ymin") + if (IS_INDEFD(ymin)) + rymin = INDEFR + else + rymin = ymin + ymax = clgetd ("ymax") + if (IS_INDEFD(ymax)) + rymax = INDEFR + else + rymax = ymax + + # Get the records list. + reclist = imtopenp ("transforms") + nrecords = imtlen (reclist) + if ((nrecords > 0) && (nrecords != nfiles)) { + call eprintf ( + "The number of records is not equal to the number of input files") + call clpcls (list) + call dtunmap (out) + call imtclose (reclist) + call sfree (sp) + return + } + + # Get the results file list. + reslist = clpopnu ("results") + nresfiles = clplen (reslist) + if (nresfiles > 1 && nresfiles != nfiles) { + call eprintf ("Error: there are too few results files\n") + call clpcls (list) + call dtunmap (out) + call imtclose (reclist) + call clpcls (reslist) + call sfree (sp) + return + } + + # Get the surface fitting parameters. + geometry = clgwrd ("fitgeometry", Memc[str], SZ_LINE, GM_GEOMETRIES) + function = clgwrd ("function", Memc[str], SZ_LINE, GM_FUNCS) + xxorder = clgeti ("xxorder") + xyorder = clgeti ("xyorder") + xxterms = clgwrd ("xxterms", Memc[str], SZ_LINE, GM_XFUNCS) - 1 + yxorder = clgeti ("yxorder") + yyorder = clgeti ("yyorder") + yxterms = clgwrd ("yxterms", Memc[str], SZ_LINE, GM_XFUNCS) - 1 + maxiter = clgeti ("maxiter") + reject = clgetd ("reject") + calctype = clgwrd ("calctype", Memc[str], SZ_LINE, ",real,double,") + + # Get the graphics parameters. + verbose = clgetb ("verbose") + interactive = clgetb ("interactive") + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + + # Flush standard output on newline. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Initialize the fit structure. + call geo_minit (fit, GM_NONE, geometry, function, xxorder, xyorder, + xxterms, yxorder, yyorder, yxterms, maxiter, reject) + + # Loop over the files. + while (clgfil (list, Memc[in_name], SZ_FNAME) != EOF) { + + # Open text file of coordinates. + in = open (Memc[in_name], READ_ONLY, TEXT_FILE) + + # Open the results files. + if (nresfiles <= 0) + res = NULL + else if (clgfil (reslist, Memc[str], SZ_FNAME) != EOF) + res = open (Memc[str], NEW_FILE, TEXT_FILE) + + # Set file name in structure. + if (nrecords > 0) { + if (imtgetim (reclist, GM_RECORD(fit), SZ_FNAME) != EOF) + ; + } else + call strcpy (Memc[in_name], GM_RECORD(fit), SZ_FNAME) + + if (verbose && res != STDOUT) { + call fstats (in, F_FILENAME, Memc[str], SZ_FNAME) + call printf ("\nCoordinate list: %s Transform: %s\n") + call pargstr (Memc[str]) + call pargstr (GM_RECORD(fit)) + if (res != NULL) + call fstats (res, F_FILENAME, Memc[str], SZ_FNAME) + else + call strcpy ("", Memc[str], SZ_FNAME) + call printf (" Results file: %s\n") + call pargstr (Memc[str]) + call flush (STDOUT) + } + if (res != NULL) { + call fstats (in, F_FILENAME, Memc[str], SZ_FNAME) + call fprintf (res, "\n# Coordinate list: %s Transform: %s\n") + call pargstr (Memc[str]) + call pargstr (GM_RECORD(fit)) + if (res != NULL) + call fstats (res, F_FILENAME, Memc[str], SZ_FNAME) + else + call strcpy ("", Memc[str], SZ_FNAME) + call fprintf (res, "# Results file: %s\n") + call pargstr (Memc[str]) + call flush (STDOUT) + } + + if (interactive) { + gd = gopen (Memc[graphics], NEW_FILE, STDGRAPH) + } else + gd = NULL + + iferr { + if (calctype == GM_REAL) + call geo_mapr (gd, in, out, res, fit, rxmin, rxmax, rymin, + rymax, verbose) + else + call geo_mapd (gd, in, out, res, fit, xmin, xmax, ymin, + ymax, verbose) + } then { + if (verbose && res != STDOUT) { + call printf ("Error fitting coordinate list: %s\n") + call pargstr (Memc[in_name]) + call flush (STDOUT) + if (errget (Memc[str], SZ_LINE) == 0) + ; + call printf ("\t%s\n") + call pargstr (Memc[str)) + } + if (res != NULL) { + call fprintf (res, "# Error fitting coordinate list: %s\n") + call pargstr (Memc[in_name]) + call flush (STDOUT) + if (errget (Memc[str], SZ_LINE) == 0) + ; + call fprintf (res, "# %s\n") + call pargstr (Memc[str)) + } + } + + call close (in) + if (nresfiles == nfiles) + call close ( res) + + if (gd != NULL) + call gclose (gd) + } + + # Close up. + call geo_free (fit) + if (nresfiles < nfiles) + call close ( res) + call dtunmap (out) + call imtclose (reclist) + call clpcls (list) + call sfree (sp) +end + + + + + +# GEO_MAP -- Procedure to calculate the coordinate transformations + +procedure geo_mapr (gd, in, out, res, fit, xmin, xmax, ymin, ymax, verbose) + +pointer gd #I the graphics stream +int in #I the input file descriptor +pointer out #I the output file descriptor +int res #I the results file descriptor +pointer fit #I pointer to fit parameters +real xmin, xmax #I max and min xref values +real ymin, ymax #I max and min yref values +bool verbose #I verbose mode + +int npts, ngood +pointer sp, str, xref, yref, xin, yin, wts, xfit, yfit, xerrmsg, yerrmsg +pointer sx1, sy1, sx2, sy2 +real mintemp, maxtemp + +real asumr() +int geo_rdxyr() +errchk geo_fitr, geo_mgfitr() + +begin + # Get working space. + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (xerrmsg, SZ_LINE, TY_CHAR) + call salloc (yerrmsg, SZ_LINE, TY_CHAR) + + # Initialize pointers. + xref = NULL + yref = NULL + xin = NULL + yin = NULL + wts = NULL + + # Read in data and check that data is in range. + npts = geo_rdxyr (in, xref, yref, xin, yin, xmin, xmax, ymin, ymax) + if (npts <= 0) { + call fstats (in, F_FILENAME, Memc[str], SZ_FNAME) + call printf ("Coordinate list: %s has no data in range.\n") + call pargstr (Memc[str]) + call sfree (sp) + return + } + + # Compute the mean of the reference and input coordinates. + GM_XOREF(fit) = double (asumr (Memr[xref], npts) / npts) + GM_YOREF(fit) = double (asumr (Memr[yref], npts) / npts) + GM_XOIN(fit) = double (asumr (Memr[xin], npts) / npts) + GM_YOIN(fit) = double (asumr (Memr[yin], npts) / npts) + + # Set the reference point for the projections to INDEF. + GM_XREFPT(fit) = INDEFD + GM_YREFPT(fit) = INDEFD + + # Compute the weights. + call malloc (xfit, npts, TY_REAL) + call malloc (yfit, npts, TY_REAL) + call malloc (wts, npts, TY_REAL) + call amovkr (real(1.), Memr[wts], npts) + + # Determine the x max and min. + if (IS_INDEFR(xmin) || IS_INDEFR(xmax)) { + call alimr (Memr[xref], npts, mintemp, maxtemp) + if (! IS_INDEFR(xmin)) + GM_XMIN(fit) = double (xmin) + else + GM_XMIN(fit) = double (mintemp) + if (! IS_INDEFR(xmax)) + GM_XMAX(fit) = double (xmax) + else + GM_XMAX(fit) = double (maxtemp) + } else { + GM_XMIN(fit) = double (xmin) + GM_XMAX(fit) = double (xmax) + } + + # Determine the y max and min. + if (IS_INDEFR(ymin) || IS_INDEFR(ymax)) { + call alimr (Memr[yref], npts, mintemp, maxtemp) + if (! IS_INDEFR(ymin)) + GM_YMIN(fit) = double (ymin) + else + GM_YMIN(fit) = double (mintemp) + if (! IS_INDEFR(ymax)) + GM_YMAX(fit) = double (ymax) + else + GM_YMAX(fit) = double (maxtemp) + } else { + GM_YMIN(fit) = double (ymin) + GM_YMAX(fit) = double (ymax) + } + + # Initalize surface pointers. + sx1 = NULL + sy1 = NULL + sx2 = NULL + sy2 = NULL + + # Fit the data. + if (gd != NULL) { + iferr { + call geo_mgfitr (gd, fit, sx1, sy1, sx2, sy2, Memr[xref], + Memr[yref], Memr[xin], Memr[yin], Memr[wts], npts, + Memc[xerrmsg], Memc[yerrmsg], SZ_LINE) + } then { + call gdeactivate (gd, 0) + call mfree (xfit, TY_REAL) + call mfree (yfit, TY_REAL) + call mfree (wts, TY_REAL) + call geo_mmfreer (sx1, sy1, sx2, sy2) + call sfree (sp) + call error (3, "Too few points for X or Y fits.") + } + call gdeactivate (gd, 0) + if (verbose && res != STDOUT) { + call printf ("Coordinate mapping status\n") + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "# Coordinate mapping status\n") + } + } else { + if (verbose && res != STDOUT) { + call printf ("Coordinate mapping status\n ") + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "# Coordinate mapping status\n# ") + } + iferr { + call geo_fitr (fit, sx1, sy1, sx2, sy2, Memr[xref], + Memr[yref], Memr[xin], Memr[yin], Memr[wts], npts, + Memc[xerrmsg], Memc[yerrmsg], SZ_LINE) + } then { + call mfree (xfit, TY_REAL) + call mfree (yfit, TY_REAL) + call mfree (wts, TY_REAL) + call geo_mmfreer (sx1, sy1, sx2, sy2) + call sfree (sp) + call error (3, "Too few points for X or Y fits.") + } + if (verbose && res != STDOUT) { + call printf ("%s %s\n") + call pargstr (Memc[xerrmsg]) + call pargstr (Memc[yerrmsg]) + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "%s %s\n") + call pargstr (Memc[xerrmsg]) + call pargstr (Memc[yerrmsg]) + call flush (STDOUT) + } + } + ngood = GM_NPTS(fit) - GM_NWTS0(fit) + if (verbose && res != STDOUT) { + call printf (" Xin and Yin fit rms: %0.7g %0.7g\n") + if (ngood <= 1) { + call pargd (0.0d0) + call pargd (0.0d0) + } else { + call pargd (sqrt (GM_XRMS(fit) / (ngood - 1))) + call pargd (sqrt (GM_YRMS(fit) / (ngood - 1))) + } + call geo_showr (STDOUT, fit, sx1, sy1, NO) + } + if (res != NULL) { + call fprintf (res, "# Xin and Yin fit rms: %0.7g %0.7g\n") + if (ngood <= 1) { + call pargd (0.0) + call pargd (0.0) + } else { + call pargd (sqrt (GM_XRMS(fit) / (ngood - 1))) + call pargd (sqrt (GM_YRMS(fit) / (ngood - 1))) + } + call geo_showr (res, fit, sx1, sy1, YES) + } + + # Compute and print the fitted x and y values. + if (res != NULL) { + call geo_evalr (sx1, sy1, sx2, sy2, Memr[xref], Memr[yref], + Memr[xfit], Memr[yfit], npts) + call geo_plistr (res, fit, Memr[xref], Memr[yref], Memr[xin], + Memr[yin], Memr[xfit], Memr[yfit], Memr[wts], npts) + } + + # Free the data + if (xref != NULL) + call mfree (xref, TY_REAL) + if (yref != NULL) + call mfree (yref, TY_REAL) + if (xin != NULL) + call mfree (xin, TY_REAL) + if (yin != NULL) + call mfree (yin, TY_REAL) + if (xfit != NULL) + call mfree (xfit, TY_REAL) + if (yfit != NULL) + call mfree (yfit, TY_REAL) + if (wts != NULL) + call mfree (wts, TY_REAL) + + # Output the data. + call geo_moutr (fit, out, sx1, sy1, sx2, sy2) + + # Free the space and close files. + call geo_mmfreer (sx1, sy1, sx2, sy2) + call sfree (sp) +end + + +define GEO_DEFBUFSIZE 1000 # default data buffer sizes + +# GEO_RDXY -- Read in the data points. + +int procedure geo_rdxyr (fd, xref, yref, xin, yin, xmin, xmax, ymin, ymax) + +int fd # the input file descriptor +pointer xref # the x reference coordinates +pointer yref # the y reference coordinates +pointer xin # the x coordinates +pointer yin # the y coordinates +real xmin, xmax # the range of the x coordinates +real ymin, ymax # the range of the y coordinates + +int npts, bufsize +int fscan(), nscan() + +begin + bufsize = GEO_DEFBUFSIZE + call malloc (xref, bufsize, TY_REAL) + call malloc (yref, bufsize, TY_REAL) + call malloc (xin, bufsize, TY_REAL) + call malloc (yin, bufsize, TY_REAL) + + npts = 0 + while (fscan (fd) != EOF) { + + # Decode the data. + call gargr (Memr[xref+npts]) + call gargr (Memr[yref+npts]) + call gargr (Memr[xin+npts]) + call gargr (Memr[yin+npts]) + if (nscan() < 4) + next + + # Check the data limits. + if (! IS_INDEFR(xmin)) { + if (Memr[xref+npts] < xmin) + next + } + if (! IS_INDEFR(xmax)) { + if (Memr[xref+npts] > xmax) + next + } + if (! IS_INDEFR(ymin)) { + if (Memr[yref+npts] < ymin) + next + } + if (! IS_INDEFR(ymax)) { + if (Memr[yref+npts] > ymax) + next + } + + npts = npts + 1 + if (npts >= bufsize) { + bufsize = bufsize + GEO_DEFBUFSIZE + call realloc (xref, bufsize, TY_REAL) + call realloc (yref, bufsize, TY_REAL) + call realloc (xin, bufsize, TY_REAL) + call realloc (yin, bufsize, TY_REAL) + } + } + + if (npts <= 0) { + call mfree (xref, TY_REAL) + call mfree (yref, TY_REAL) + call mfree (xin, TY_REAL) + call mfree (yin, TY_REAL) + xref = NULL + yref = NULL + xin = NULL + yin = NULL + } else if (npts < bufsize) { + call realloc (xref, npts, TY_REAL) + call realloc (yref, npts, TY_REAL) + call realloc (xin, npts, TY_REAL) + call realloc (yin, npts, TY_REAL) + } + + return (npts) +end + + +# GEO_EVAL -- Evalute the fit. + +procedure geo_evalr (sx1, sy1, sx2, sy2, xref, yref, xi, eta, npts) + +pointer sx1, sy1 #I pointer to linear surfaces +pointer sx2, sy2 #I pointer to higher order surfaces +real xref[ARB] #I the x reference coordinates +real yref[ARB] #I the y reference coordinates +real xi[ARB] #O the fitted xi coordinates +real eta[ARB] #O the fitted eta coordinates +int npts #I the number of points + +pointer sp, temp + +begin + call smark (sp) + call salloc (temp, npts, TY_REAL) + + call gsvector (sx1, xref, yref, xi, npts) + if (sx2 != NULL) { + call gsvector (sx2, xref, yref, Memr[temp], npts) + call aaddr (Memr[temp], xi, xi, npts) + } + call gsvector (sy1, xref, yref, eta, npts) + if (sy2 != NULL) { + call gsvector (sy2, xref, yref, Memr[temp], npts) + + call aaddr (Memr[temp], eta, eta, npts) + } + + call sfree (sp) +end + + +# GEO_MOUT -- Write the output database file. + +procedure geo_moutr (fit, out, sx1, sy1, sx2, sy2) + +pointer fit #I pointer to fitting structure +int out #I pointer to database file +pointer sx1, sy1 #I pointer to linear surfaces +pointer sx2, sy2 #I pointer to distortion surfaces + +int i, npts, ncoeff +pointer sp, str, xcoeff, ycoeff +real xrms, yrms, xshift, yshift, xscale, yscale, xrot, yrot +int gsgeti() +int rg_wrdstr() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Compute the x and y fit rms. + #npts = max (0, GM_NPTS(fit) - GM_NREJECT(fit) - GM_NWTS0(fit)) + npts = max (0, GM_NPTS(fit) - GM_NWTS0(fit)) + xrms = max (0.0d0, GM_XRMS(fit)) + yrms = max (0.0d0, GM_YRMS(fit)) + if (npts > 1) { + xrms = sqrt (xrms / (npts - 1)) + yrms = sqrt (yrms / (npts - 1)) + } else { + xrms = 0.0d0 + yrms = 0.0d0 + } + + # Print title. + call dtptime (out) + call dtput (out, "begin\t%s\n") + call pargstr (GM_RECORD(fit)) + + # Print the x and y mean values. + call dtput (out, "\txrefmean\t%g\n") + call pargd (GM_XOREF(fit)) + call dtput (out, "\tyrefmean\t%g\n") + call pargd (GM_YOREF(fit)) + call dtput (out, "\txmean\t\t%g\n") + call pargd (GM_XOIN(fit)) + call dtput (out, "\tymean\t\t%g\n") + call pargd (GM_YOIN(fit)) + + # Print some of the fitting parameters. + if (rg_wrdstr (GM_FIT(fit), Memc[str], SZ_FNAME, GM_GEOMETRIES) <= 0) + call strcpy ("general", Memc[str], SZ_FNAME) + call dtput (out, "\tgeometry\t%s\n") + call pargstr (Memc[str]) + if (rg_wrdstr (GM_FUNCTION(fit), Memc[str], SZ_FNAME, GM_FUNCS) <= 0) + call strcpy ("polynomial", Memc[str], SZ_FNAME) + call dtput (out, "\tfunction\t%s\n") + call pargstr (Memc[str]) + + # Output the geometric parameters. + call geo_lcoeffr (sx1, sy1, xshift, yshift, xscale, yscale, xrot, yrot) + call dtput (out, "\txshift\t\t%g\n") + call pargr (xshift) + call dtput (out, "\tyshift\t\t%g\n") + call pargr (yshift) + call dtput (out, "\txmag\t\t%g\n") + call pargr (xscale) + call dtput (out, "\tymag\t\t%g\n") + call pargr (yscale) + call dtput (out, "\txrotation\t%g\n") + call pargr (xrot) + call dtput (out, "\tyrotation\t%g\n") + call pargr (yrot) + + # Out the rms values. + call dtput (out, "\txrms\t\t%g\n") + call pargr (real(xrms)) + call dtput (out, "\tyrms\t\t%g\n") + call pargr (real(yrms)) + + # Allocate memory for linear coefficients. + ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE)) + call calloc (xcoeff, ncoeff, TY_REAL) + call calloc (ycoeff, ncoeff, TY_REAL) + + # Output the linear coefficients. + call gssave (sx1, Memr[xcoeff]) + call gssave (sy1, Memr[ycoeff]) + call dtput (out, "\tsurface1\t%d\n") + call pargi (ncoeff) + do i = 1, ncoeff { + call dtput (out, "\t\t\t%g\t%g\n") + call pargr (Memr[xcoeff+i-1]) + call pargr (Memr[ycoeff+i-1]) + } + + call mfree (xcoeff, TY_REAL) + call mfree (ycoeff, TY_REAL) + + # Allocate memory for higer order coefficients. + if (sx2 == NULL) + ncoeff = 0 + else + ncoeff = gsgeti (sx2, GSNSAVE) + if (sy2 == NULL) + ncoeff = max (0, ncoeff) + else + ncoeff = max (gsgeti (sy2, GSNSAVE), ncoeff) + call calloc (xcoeff, ncoeff, TY_REAL) + call calloc (ycoeff, ncoeff, TY_REAL) + + # Save the coefficients. + call gssave (sx2, Memr[xcoeff]) + call gssave (sy2, Memr[ycoeff]) + + # Output the coefficients. + call dtput (out, "\tsurface2\t%d\n") + call pargi (ncoeff) + do i = 1, ncoeff { + call dtput (out, "\t\t\t%g\t%g\n") + call pargr (Memr[xcoeff+i-1]) + call pargr (Memr[ycoeff+i-1]) + } + + # Cleanup. + call mfree (xcoeff, TY_REAL) + call mfree (ycoeff, TY_REAL) + call sfree (sp) +end + + +# GEO_PLIST -- Print the input, output, and fitted data and the residuals. + +procedure geo_plistr (fd, fit, xref, yref, xin, yin, xfit, yfit, wts, npts) + +int fd #I the results file descriptor +pointer fit #I pointer to the fit structure +real xref[ARB] #I the input x coordinates +real yref[ARB] #I the input y coordinates +real xin[ARB] #I the input ra / longitude coordinates +real yin[ARB] #I the input dec / latitude coordinates +real xfit[ARB] #I the fitted ra / longitude coordinates +real yfit[ARB] #I the fitted dec / latitude coordinates +real wts[ARB] #I the weights array +int npts #I the number of data points + +int i, index +pointer sp, fmtstr, twts + +begin + # Allocate working space. + call smark (sp) + call salloc (fmtstr, SZ_LINE, TY_CHAR) + call salloc (twts, npts, TY_REAL) + + # Compute the weights. + call amovr (wts, Memr[twts], npts) + do i = 1, GM_NREJECT(fit) { + index = Memi[GM_REJ(fit)+i-1] + if (wts[index] > real(0.0)) + Memr[twts+index-1] = real(0.0) + } + + # Print banner. + call fprintf (fd, "\n# Input Coordinate Listing\n") + call fprintf (fd, "# Column 1: X (reference) \n") + call fprintf (fd, "# Column 2: Y (reference)\n") + call fprintf (fd, "# Column 3: X (input)\n") + call fprintf (fd, "# Column 4: Y (input)\n") + call fprintf (fd, "# Column 5: X (fit)\n") + call fprintf (fd, "# Column 6: Y (fit)\n") + call fprintf (fd, "# Column 7: X (residual)\n") + call fprintf (fd, "# Column 8: Y (residual)\n\n") + + # Create the format string. + call sprintf (Memc[fmtstr], SZ_LINE, "%s %s %s %s %s %s %s %s\n") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + call pargstr ("%9.7g") + + # Print the data. + do i = 1, npts { + call fprintf (fd, Memc[fmtstr]) + call pargr (xref[i]) + call pargr (yref[i]) + call pargr (xin[i]) + call pargr (yin[i]) + if (Memr[twts+i-1] > 0.0d0) { + call pargr (xfit[i]) + call pargr (yfit[i]) + call pargr (xin[i] - xfit[i]) + call pargr (yin[i] - yfit[i]) + } else { + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + call pargr (INDEFR) + } + + } + + call fprintf (fd, "\n") + + call sfree (sp) + +end + +# GEO_SHOW -- Print the coordinate mapping parameters. + +procedure geo_showr (fd, fit, sx1, sy1, comment) + +int fd #I the output file descriptor +pointer fit #I pointer to the fit structure +pointer sx1, sy1 #I pointer to linear surfaces +int comment #I comment the output ? + +real xshift, yshift, a, b, c, d +real xscale, yscale, xrot, yrot +pointer sp, str +bool fp_equalr() + +begin + # Allocate temporary space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Compute the geometric parameters. + call geo_gcoeffr (sx1, sy1, xshift, yshift, a, b, c, d) + + if (comment == NO) { + call fprintf (fd, "Coordinate mapping parameters\n") + } else { + call fprintf (fd, "# Coordinate mapping parameters\n") + } + + if (comment == NO) { + call fprintf (fd, + " Mean Xref and Yref: %0.7g %0.7g\n") + call pargd (GM_XOREF(fit)) + call pargd (GM_YOREF(fit)) + call fprintf (fd, + " Mean Xin and Yin: %0.7g %0.7g\n") + call pargd (GM_XOIN(fit)) + call pargd (GM_YOIN(fit)) + call fprintf (fd, + " X and Y shift: %0.7g %0.7g (xin yin)\n") + call pargr (xshift) + call pargr (yshift) + } else { + call fprintf (fd, + "# Mean Xref and Yref: %0.7g %0.7g\n") + call pargd (GM_XOREF(fit)) + call pargd (GM_YOREF(fit)) + call fprintf (fd, + "# Mean Xin and Yin: %0.7g %g0.7\n") + call pargd (GM_XOIN(fit)) + call pargd (GM_YOIN(fit)) + call fprintf (fd, + "# X and Y shift: %0.7g %0.7g (xin yin)\n") + call pargr (xshift) + call pargr (yshift) + } + + # Output the scale factors. + xscale = sqrt (a * a + c * c) + yscale = sqrt (b * b + d * d) + if (comment == NO) { + call fprintf (fd, + " X and Y scale: %0.7g %0.7g (xin / xref yin / yref)\n") + call pargr (xscale) + call pargr (yscale) + } else { + call fprintf (fd, + "# X and Y scale: %0.7g %0.7g (xin / xref yin / yref)\n") + call pargr (xscale) + call pargr (yscale) + } + + # Output the rotation factors. + if (fp_equalr (a, real(0.0)) && fp_equalr (c, real(0.0))) + xrot = real(0.0) + else + xrot = RADTODEG (atan2 (-c, a)) + if (xrot < real(0.0)) + xrot = xrot + real(360.0) + if (fp_equalr (b, real(0.0)) && fp_equalr (d, real(0.0))) + yrot = real(0.0) + else + yrot = RADTODEG (atan2 (b, d)) + if (yrot < real(0.0)) + yrot = yrot + real(360.0) + if (comment == NO) { + call fprintf (fd, + " X and Y axis rotation: %0.5f %0.5f (degrees degrees)\n") + call pargr (xrot) + call pargr (yrot) + } else { + call fprintf (fd, + "# X and Y axis rotation: %0.5f %0.5f (degrees degrees)\n") + call pargr (xrot) + call pargr (yrot) + } + + call sfree (sp) +end + + + +# GEO_MAP -- Procedure to calculate the coordinate transformations + +procedure geo_mapd (gd, in, out, res, fit, xmin, xmax, ymin, ymax, verbose) + +pointer gd #I the graphics stream +int in #I the input file descriptor +pointer out #I the output file descriptor +int res #I the results file descriptor +pointer fit #I pointer to fit parameters +double xmin, xmax #I max and min xref values +double ymin, ymax #I max and min yref values +bool verbose #I verbose mode + +int npts, ngood +pointer sp, str, xref, yref, xin, yin, wts, xfit, yfit, xerrmsg, yerrmsg +pointer sx1, sy1, sx2, sy2 +double mintemp, maxtemp + +double asumd() +int geo_rdxyd() +errchk geo_fitd, geo_mgfitd() + +begin + # Get working space. + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (xerrmsg, SZ_LINE, TY_CHAR) + call salloc (yerrmsg, SZ_LINE, TY_CHAR) + + # Initialize pointers. + xref = NULL + yref = NULL + xin = NULL + yin = NULL + wts = NULL + + # Read in data and check that data is in range. + npts = geo_rdxyd (in, xref, yref, xin, yin, xmin, xmax, ymin, ymax) + if (npts <= 0) { + call fstats (in, F_FILENAME, Memc[str], SZ_FNAME) + call printf ("Coordinate list: %s has no data in range.\n") + call pargstr (Memc[str]) + call sfree (sp) + return + } + + # Compute the mean of the reference and input coordinates. + GM_XOREF(fit) = double (asumd (Memd[xref], npts) / npts) + GM_YOREF(fit) = double (asumd (Memd[yref], npts) / npts) + GM_XOIN(fit) = double (asumd (Memd[xin], npts) / npts) + GM_YOIN(fit) = double (asumd (Memd[yin], npts) / npts) + + # Set the reference point for the projections to INDEF. + GM_XREFPT(fit) = INDEFD + GM_YREFPT(fit) = INDEFD + + # Compute the weights. + call malloc (xfit, npts, TY_DOUBLE) + call malloc (yfit, npts, TY_DOUBLE) + call malloc (wts, npts, TY_DOUBLE) + call amovkd (double(1.), Memd[wts], npts) + + # Determine the x max and min. + if (IS_INDEFD(xmin) || IS_INDEFD(xmax)) { + call alimd (Memd[xref], npts, mintemp, maxtemp) + if (! IS_INDEFD(xmin)) + GM_XMIN(fit) = double (xmin) + else + GM_XMIN(fit) = double (mintemp) + if (! IS_INDEFD(xmax)) + GM_XMAX(fit) = double (xmax) + else + GM_XMAX(fit) = double (maxtemp) + } else { + GM_XMIN(fit) = double (xmin) + GM_XMAX(fit) = double (xmax) + } + + # Determine the y max and min. + if (IS_INDEFD(ymin) || IS_INDEFD(ymax)) { + call alimd (Memd[yref], npts, mintemp, maxtemp) + if (! IS_INDEFD(ymin)) + GM_YMIN(fit) = double (ymin) + else + GM_YMIN(fit) = double (mintemp) + if (! IS_INDEFD(ymax)) + GM_YMAX(fit) = double (ymax) + else + GM_YMAX(fit) = double (maxtemp) + } else { + GM_YMIN(fit) = double (ymin) + GM_YMAX(fit) = double (ymax) + } + + # Initalize surface pointers. + sx1 = NULL + sy1 = NULL + sx2 = NULL + sy2 = NULL + + # Fit the data. + if (gd != NULL) { + iferr { + call geo_mgfitd (gd, fit, sx1, sy1, sx2, sy2, Memd[xref], + Memd[yref], Memd[xin], Memd[yin], Memd[wts], npts, + Memc[xerrmsg], Memc[yerrmsg], SZ_LINE) + } then { + call gdeactivate (gd, 0) + call mfree (xfit, TY_DOUBLE) + call mfree (yfit, TY_DOUBLE) + call mfree (wts, TY_DOUBLE) + call geo_mmfreed (sx1, sy1, sx2, sy2) + call sfree (sp) + call error (3, "Too few points for X or Y fits.") + } + call gdeactivate (gd, 0) + if (verbose && res != STDOUT) { + call printf ("Coordinate mapping status\n") + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "# Coordinate mapping status\n") + } + } else { + if (verbose && res != STDOUT) { + call printf ("Coordinate mapping status\n ") + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "# Coordinate mapping status\n# ") + } + iferr { + call geo_fitd (fit, sx1, sy1, sx2, sy2, Memd[xref], + Memd[yref], Memd[xin], Memd[yin], Memd[wts], npts, + Memc[xerrmsg], Memc[yerrmsg], SZ_LINE) + } then { + call mfree (xfit, TY_DOUBLE) + call mfree (yfit, TY_DOUBLE) + call mfree (wts, TY_DOUBLE) + call geo_mmfreed (sx1, sy1, sx2, sy2) + call sfree (sp) + call error (3, "Too few points for X or Y fits.") + } + if (verbose && res != STDOUT) { + call printf ("%s %s\n") + call pargstr (Memc[xerrmsg]) + call pargstr (Memc[yerrmsg]) + call flush (STDOUT) + } + if (res != NULL) { + call fprintf (res, "%s %s\n") + call pargstr (Memc[xerrmsg]) + call pargstr (Memc[yerrmsg]) + call flush (STDOUT) + } + } + ngood = GM_NPTS(fit) - GM_NWTS0(fit) + if (verbose && res != STDOUT) { + call printf (" Xin and Yin fit rms: %0.7g %0.7g\n") + if (ngood <= 1) { + call pargd (0.0d0) + call pargd (0.0d0) + } else { + call pargd (sqrt (GM_XRMS(fit) / (ngood - 1))) + call pargd (sqrt (GM_YRMS(fit) / (ngood - 1))) + } + call geo_showd (STDOUT, fit, sx1, sy1, NO) + } + if (res != NULL) { + call fprintf (res, "# Xin and Yin fit rms: %0.7g %0.7g\n") + if (ngood <= 1) { + call pargd (0.0) + call pargd (0.0) + } else { + call pargd (sqrt (GM_XRMS(fit) / (ngood - 1))) + call pargd (sqrt (GM_YRMS(fit) / (ngood - 1))) + } + call geo_showd (res, fit, sx1, sy1, YES) + } + + # Compute and print the fitted x and y values. + if (res != NULL) { + call geo_evald (sx1, sy1, sx2, sy2, Memd[xref], Memd[yref], + Memd[xfit], Memd[yfit], npts) + call geo_plistd (res, fit, Memd[xref], Memd[yref], Memd[xin], + Memd[yin], Memd[xfit], Memd[yfit], Memd[wts], npts) + } + + # Free the data + if (xref != NULL) + call mfree (xref, TY_DOUBLE) + if (yref != NULL) + call mfree (yref, TY_DOUBLE) + if (xin != NULL) + call mfree (xin, TY_DOUBLE) + if (yin != NULL) + call mfree (yin, TY_DOUBLE) + if (xfit != NULL) + call mfree (xfit, TY_DOUBLE) + if (yfit != NULL) + call mfree (yfit, TY_DOUBLE) + if (wts != NULL) + call mfree (wts, TY_DOUBLE) + + # Output the data. + call geo_moutd (fit, out, sx1, sy1, sx2, sy2) + + # Free the space and close files. + call geo_mmfreed (sx1, sy1, sx2, sy2) + call sfree (sp) +end + + +define GEO_DEFBUFSIZE 1000 # default data buffer sizes + +# GEO_RDXY -- Read in the data points. + +int procedure geo_rdxyd (fd, xref, yref, xin, yin, xmin, xmax, ymin, ymax) + +int fd # the input file descriptor +pointer xref # the x reference coordinates +pointer yref # the y reference coordinates +pointer xin # the x coordinates +pointer yin # the y coordinates +double xmin, xmax # the range of the x coordinates +double ymin, ymax # the range of the y coordinates + +int npts, bufsize +int fscan(), nscan() + +begin + bufsize = GEO_DEFBUFSIZE + call malloc (xref, bufsize, TY_DOUBLE) + call malloc (yref, bufsize, TY_DOUBLE) + call malloc (xin, bufsize, TY_DOUBLE) + call malloc (yin, bufsize, TY_DOUBLE) + + npts = 0 + while (fscan (fd) != EOF) { + + # Decode the data. + call gargd (Memd[xref+npts]) + call gargd (Memd[yref+npts]) + call gargd (Memd[xin+npts]) + call gargd (Memd[yin+npts]) + if (nscan() < 4) + next + + # Check the data limits. + if (! IS_INDEFD(xmin)) { + if (Memd[xref+npts] < xmin) + next + } + if (! IS_INDEFD(xmax)) { + if (Memd[xref+npts] > xmax) + next + } + if (! IS_INDEFD(ymin)) { + if (Memd[yref+npts] < ymin) + next + } + if (! IS_INDEFD(ymax)) { + if (Memd[yref+npts] > ymax) + next + } + + npts = npts + 1 + if (npts >= bufsize) { + bufsize = bufsize + GEO_DEFBUFSIZE + call realloc (xref, bufsize, TY_DOUBLE) + call realloc (yref, bufsize, TY_DOUBLE) + call realloc (xin, bufsize, TY_DOUBLE) + call realloc (yin, bufsize, TY_DOUBLE) + } + } + + if (npts <= 0) { + call mfree (xref, TY_DOUBLE) + call mfree (yref, TY_DOUBLE) + call mfree (xin, TY_DOUBLE) + call mfree (yin, TY_DOUBLE) + xref = NULL + yref = NULL + xin = NULL + yin = NULL + } else if (npts < bufsize) { + call realloc (xref, npts, TY_DOUBLE) + call realloc (yref, npts, TY_DOUBLE) + call realloc (xin, npts, TY_DOUBLE) + call realloc (yin, npts, TY_DOUBLE) + } + + return (npts) +end + + +# GEO_EVAL -- Evalute the fit. + +procedure geo_evald (sx1, sy1, sx2, sy2, xref, yref, xi, eta, npts) + +pointer sx1, sy1 #I pointer to linear surfaces +pointer sx2, sy2 #I pointer to higher order surfaces +double xref[ARB] #I the x reference coordinates +double yref[ARB] #I the y reference coordinates +double xi[ARB] #O the fitted xi coordinates +double eta[ARB] #O the fitted eta coordinates +int npts #I the number of points + +pointer sp, temp + +begin + call smark (sp) + call salloc (temp, npts, TY_DOUBLE) + + call dgsvector (sx1, xref, yref, xi, npts) + if (sx2 != NULL) { + call dgsvector (sx2, xref, yref, Memd[temp], npts) + call aaddd (Memd[temp], xi, xi, npts) + } + call dgsvector (sy1, xref, yref, eta, npts) + if (sy2 != NULL) { + call dgsvector (sy2, xref, yref, Memd[temp], npts) + + call aaddd (Memd[temp], eta, eta, npts) + } + + call sfree (sp) +end + + +# GEO_MOUT -- Write the output database file. + +procedure geo_moutd (fit, out, sx1, sy1, sx2, sy2) + +pointer fit #I pointer to fitting structure +int out #I pointer to database file +pointer sx1, sy1 #I pointer to linear surfaces +pointer sx2, sy2 #I pointer to distortion surfaces + +int i, npts, ncoeff +pointer sp, str, xcoeff, ycoeff +double xrms, yrms, xshift, yshift, xscale, yscale, xrot, yrot +int dgsgeti() +int rg_wrdstr() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Compute the x and y fit rms. + #npts = max (0, GM_NPTS(fit) - GM_NREJECT(fit) - GM_NWTS0(fit)) + npts = max (0, GM_NPTS(fit) - GM_NWTS0(fit)) + xrms = max (0.0d0, GM_XRMS(fit)) + yrms = max (0.0d0, GM_YRMS(fit)) + if (npts > 1) { + xrms = sqrt (xrms / (npts - 1)) + yrms = sqrt (yrms / (npts - 1)) + } else { + xrms = 0.0d0 + yrms = 0.0d0 + } + + # Print title. + call dtptime (out) + call dtput (out, "begin\t%s\n") + call pargstr (GM_RECORD(fit)) + + # Print the x and y mean values. + call dtput (out, "\txrefmean\t%g\n") + call pargd (GM_XOREF(fit)) + call dtput (out, "\tyrefmean\t%g\n") + call pargd (GM_YOREF(fit)) + call dtput (out, "\txmean\t\t%g\n") + call pargd (GM_XOIN(fit)) + call dtput (out, "\tymean\t\t%g\n") + call pargd (GM_YOIN(fit)) + + # Print some of the fitting parameters. + if (rg_wrdstr (GM_FIT(fit), Memc[str], SZ_FNAME, GM_GEOMETRIES) <= 0) + call strcpy ("general", Memc[str], SZ_FNAME) + call dtput (out, "\tgeometry\t%s\n") + call pargstr (Memc[str]) + if (rg_wrdstr (GM_FUNCTION(fit), Memc[str], SZ_FNAME, GM_FUNCS) <= 0) + call strcpy ("polynomial", Memc[str], SZ_FNAME) + call dtput (out, "\tfunction\t%s\n") + call pargstr (Memc[str]) + + # Output the geometric parameters. + call geo_lcoeffd (sx1, sy1, xshift, yshift, xscale, yscale, xrot, yrot) + call dtput (out, "\txshift\t\t%g\n") + call pargd (xshift) + call dtput (out, "\tyshift\t\t%g\n") + call pargd (yshift) + call dtput (out, "\txmag\t\t%g\n") + call pargd (xscale) + call dtput (out, "\tymag\t\t%g\n") + call pargd (yscale) + call dtput (out, "\txrotation\t%g\n") + call pargd (xrot) + call dtput (out, "\tyrotation\t%g\n") + call pargd (yrot) + + # Out the rms values. + call dtput (out, "\txrms\t\t%g\n") + call pargd (double(xrms)) + call dtput (out, "\tyrms\t\t%g\n") + call pargd (double(yrms)) + + # Allocate memory for linear coefficients. + ncoeff = max (dgsgeti (sx1, GSNSAVE), dgsgeti (sy1, GSNSAVE)) + call calloc (xcoeff, ncoeff, TY_DOUBLE) + call calloc (ycoeff, ncoeff, TY_DOUBLE) + + # Output the linear coefficients. + call dgssave (sx1, Memd[xcoeff]) + call dgssave (sy1, Memd[ycoeff]) + call dtput (out, "\tsurface1\t%d\n") + call pargi (ncoeff) + do i = 1, ncoeff { + call dtput (out, "\t\t\t%g\t%g\n") + call pargd (Memd[xcoeff+i-1]) + call pargd (Memd[ycoeff+i-1]) + } + + call mfree (xcoeff, TY_DOUBLE) + call mfree (ycoeff, TY_DOUBLE) + + # Allocate memory for higer order coefficients. + if (sx2 == NULL) + ncoeff = 0 + else + ncoeff = dgsgeti (sx2, GSNSAVE) + if (sy2 == NULL) + ncoeff = max (0, ncoeff) + else + ncoeff = max (dgsgeti (sy2, GSNSAVE), ncoeff) + call calloc (xcoeff, ncoeff, TY_DOUBLE) + call calloc (ycoeff, ncoeff, TY_DOUBLE) + + # Save the coefficients. + call dgssave (sx2, Memd[xcoeff]) + call dgssave (sy2, Memd[ycoeff]) + + # Output the coefficients. + call dtput (out, "\tsurface2\t%d\n") + call pargi (ncoeff) + do i = 1, ncoeff { + call dtput (out, "\t\t\t%g\t%g\n") + call pargd (Memd[xcoeff+i-1]) + call pargd (Memd[ycoeff+i-1]) + } + + # Cleanup. + call mfree (xcoeff, TY_DOUBLE) + call mfree (ycoeff, TY_DOUBLE) + call sfree (sp) +end + + +# GEO_PLIST -- Print the input, output, and fitted data and the residuals. + +procedure geo_plistd (fd, fit, xref, yref, xin, yin, xfit, yfit, wts, npts) + +int fd #I the results file descriptor +pointer fit #I pointer to the fit structure +double xref[ARB] #I the input x coordinates +double yref[ARB] #I the input y coordinates +double xin[ARB] #I the input ra / longitude coordinates +double yin[ARB] #I the input dec / latitude coordinates +double xfit[ARB] #I the fitted ra / longitude coordinates +double yfit[ARB] #I the fitted dec / latitude coordinates +double wts[ARB] #I the weights array +int npts #I the number of data points + +int i, index +pointer sp, fmtstr, twts + +begin + # Allocate working space. + call smark (sp) + call salloc (fmtstr, SZ_LINE, TY_CHAR) + call salloc (twts, npts, TY_DOUBLE) + + # Compute the weights. + call amovd (wts, Memd[twts], npts) + do i = 1, GM_NREJECT(fit) { + index = Memi[GM_REJ(fit)+i-1] + if (wts[index] > double(0.0)) + Memd[twts+index-1] = double(0.0) + } + + # Print banner. + call fprintf (fd, "\n# Input Coordinate Listing\n") + call fprintf (fd, "# Column 1: X (reference) \n") + call fprintf (fd, "# Column 2: Y (reference)\n") + call fprintf (fd, "# Column 3: X (input)\n") + call fprintf (fd, "# Column 4: Y (input)\n") + call fprintf (fd, "# Column 5: X (fit)\n") + call fprintf (fd, "# Column 6: Y (fit)\n") + call fprintf (fd, "# Column 7: X (residual)\n") + call fprintf (fd, "# Column 8: Y (residual)\n\n") + + # Create the format string. + call sprintf (Memc[fmtstr], SZ_LINE, "%s %s %s %s %s %s %s %s\n") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + call pargstr ("%16.14g") + + # Print the data. + do i = 1, npts { + call fprintf (fd, Memc[fmtstr]) + call pargd (xref[i]) + call pargd (yref[i]) + call pargd (xin[i]) + call pargd (yin[i]) + if (Memd[twts+i-1] > 0.0d0) { + call pargd (xfit[i]) + call pargd (yfit[i]) + call pargd (xin[i] - xfit[i]) + call pargd (yin[i] - yfit[i]) + } else { + call pargd (INDEFD) + call pargd (INDEFD) + call pargd (INDEFD) + call pargd (INDEFD) + } + + } + + call fprintf (fd, "\n") + + call sfree (sp) + +end + +# GEO_SHOW -- Print the coordinate mapping parameters. + +procedure geo_showd (fd, fit, sx1, sy1, comment) + +int fd #I the output file descriptor +pointer fit #I pointer to the fit structure +pointer sx1, sy1 #I pointer to linear surfaces +int comment #I comment the output ? + +double xshift, yshift, a, b, c, d +double xscale, yscale, xrot, yrot +pointer sp, str +bool fp_equald() + +begin + # Allocate temporary space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Compute the geometric parameters. + call geo_gcoeffd (sx1, sy1, xshift, yshift, a, b, c, d) + + if (comment == NO) { + call fprintf (fd, "Coordinate mapping parameters\n") + } else { + call fprintf (fd, "# Coordinate mapping parameters\n") + } + + if (comment == NO) { + call fprintf (fd, + " Mean Xref and Yref: %0.7g %0.7g\n") + call pargd (GM_XOREF(fit)) + call pargd (GM_YOREF(fit)) + call fprintf (fd, + " Mean Xin and Yin: %0.7g %0.7g\n") + call pargd (GM_XOIN(fit)) + call pargd (GM_YOIN(fit)) + call fprintf (fd, + " X and Y shift: %0.7g %0.7g (xin yin)\n") + call pargd (xshift) + call pargd (yshift) + } else { + call fprintf (fd, + "# Mean Xref and Yref: %0.7g %0.7g\n") + call pargd (GM_XOREF(fit)) + call pargd (GM_YOREF(fit)) + call fprintf (fd, + "# Mean Xin and Yin: %0.7g %g0.7\n") + call pargd (GM_XOIN(fit)) + call pargd (GM_YOIN(fit)) + call fprintf (fd, + "# X and Y shift: %0.7g %0.7g (xin yin)\n") + call pargd (xshift) + call pargd (yshift) + } + + # Output the scale factors. + xscale = sqrt (a * a + c * c) + yscale = sqrt (b * b + d * d) + if (comment == NO) { + call fprintf (fd, + " X and Y scale: %0.7g %0.7g (xin / xref yin / yref)\n") + call pargd (xscale) + call pargd (yscale) + } else { + call fprintf (fd, + "# X and Y scale: %0.7g %0.7g (xin / xref yin / yref)\n") + call pargd (xscale) + call pargd (yscale) + } + + # Output the rotation factors. + if (fp_equald (a, double(0.0)) && fp_equald (c, double(0.0))) + xrot = double(0.0) + else + xrot = RADTODEG (atan2 (-c, a)) + if (xrot < double(0.0)) + xrot = xrot + double(360.0) + if (fp_equald (b, double(0.0)) && fp_equald (d, double(0.0))) + yrot = double(0.0) + else + yrot = RADTODEG (atan2 (b, d)) + if (yrot < double(0.0)) + yrot = yrot + double(360.0) + if (comment == NO) { + call fprintf (fd, + " X and Y axis rotation: %0.5f %0.5f (degrees degrees)\n") + call pargd (xrot) + call pargd (yrot) + } else { + call fprintf (fd, + "# X and Y axis rotation: %0.5f %0.5f (degrees degrees)\n") + call pargd (xrot) + call pargd (yrot) + } + + call sfree (sp) +end + + diff --git a/pkg/images/immatch/src/geometry/t_geotran.x b/pkg/images/immatch/src/geometry/t_geotran.x new file mode 100644 index 00000000..5e5cd2e3 --- /dev/null +++ b/pkg/images/immatch/src/geometry/t_geotran.x @@ -0,0 +1,880 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <mwset.h> +include <math.h> +include <math/gsurfit.h> +include "geotran.h" + +# T_GEOTRAN -- Geometrically transform a list of images either linearly or +# using a transformation computed by the GEOMAP task. + +procedure t_geotran () + +int ncols, nlines # output picture size +real xmin, xmax, ymin, ymax # minimum and maximum ref values +real xscale, yscale # output picture scale +real xin, yin # input picture origin +real xshift, yshift # x and y shifts +real xout, yout # output picture origin +real xmag, ymag # input picture scale +real xrotation, yrotation # rotation angle +int nxblock, nyblock # block size of image to be used + +bool verbose +int list1, list2, tflist, ndim, nc, nl, mode +pointer sp, imtlist1, imtlist2, database, transform, record +pointer image1, image2, imtemp, imroot, section, str +pointer geo, sx1, sy1, sx2, sy2, in, out, mw +real xs, ys, txshift, tyshift, txmag, tymag, txrot, tyrot +double oltv[2], nltv[2], oltm[2,2], nltm[2,2] + +bool clgetb(), envgetb(), streq() +int imtopen(), imtlen(), clgeti(), imtgetim(), clgwrd(), btoi() +pointer immap(), mw_openim() +real clgetr() +errchk immap() + +begin + # Set up the geotran structure. + call smark (sp) + call salloc (imtlist1, SZ_LINE, TY_CHAR) + call salloc (imtlist2, SZ_LINE, TY_CHAR) + call salloc (database, SZ_FNAME, TY_CHAR) + call salloc (transform, SZ_FNAME, TY_CHAR) + call salloc (record, SZ_FNAME, 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 (imroot, SZ_FNAME, TY_CHAR) + call salloc (section, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (geo, LEN_GEOSTRUCT, TY_STRUCT) + + # Get the input and output lists and database file. + call clgstr ("input", Memc[imtlist1], SZ_FNAME) + call clgstr ("output", Memc[imtlist2], SZ_FNAME) + call clgstr ("database", Memc[database], SZ_FNAME) + if (Memc[database] != EOS) { + call clgstr ("transforms", Memc[transform], SZ_FNAME) + tflist = imtopen (Memc[transform]) + GT_GEOMODE(geo) = clgwrd ("geometry", Memc[str], SZ_LINE, + ",junk,linear,distortion,geometric,") + } else { + tflist = NULL + GT_GEOMODE(geo) = GT_NONE + } + + # Get the output picture format parameters. + xmin = clgetr ("xmin") + xmax = clgetr ("xmax") + ymin = clgetr ("ymin") + ymax = clgetr ("ymax") + xscale = clgetr ("xscale") + yscale = clgetr ("yscale") + ncols= clgeti ("ncols") + nlines = clgeti ("nlines") + + # Get the geometric transformation parameters. + xin = clgetr ("xin") + yin = clgetr ("yin") + xshift = clgetr ("xshift") + yshift = clgetr ("yshift") + xout = clgetr ("xout") + yout = clgetr ("yout") + xmag = clgetr ("xmag") + ymag = clgetr ("ymag") + xrotation = clgetr ("xrotation") + yrotation = clgetr ("yrotation") + + # Get the interpolation parameters. + call clgstr ("interpolant", GT_INTERPSTR(geo), SZ_FNAME) + #GT_INTERPOLANT(geo) = clgwrd ("interpolant", Memc[str], SZ_LINE, + #",nearest,linear,poly3,poly5,spline3,") + GT_BOUNDARY(geo) = clgwrd ("boundary", Memc[str], SZ_LINE, + ",constant,nearest,reflect,wrap,") + GT_CONSTANT(geo) = clgetr ("constant") + GT_XSAMPLE(geo) = clgetr ("xsample") + GT_YSAMPLE(geo) = clgetr ("ysample") + GT_FLUXCONSERVE(geo) = btoi (clgetb("fluxconserve")) + + nxblock = clgeti ("nxblock") + nyblock = clgeti ("nyblock") + verbose = clgetb ("verbose") + + # Open the lists of images and check the scale lengths. + list1 = imtopen (Memc[imtlist1]) + list2 = imtopen (Memc[imtlist2]) + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + if (tflist != NULL) + call imtclose (tflist) + call error (0, "Input and output lists not the same length.") + } + + # Check the transform list. + if (tflist != NULL) { + if (imtlen (tflist) > 1 && imtlen (tflist) != imtlen (list1)) { + call imtclose (list1) + call imtclose (list2) + call imtclose (tflist) + call error (0, "Transform and input lists not the same length.") + } + } + + # Loop over the images. + if (verbose) { + call printf ("\n") + } + while (imtgetim (list1, Memc[image1], SZ_FNAME) != EOF && + imtgetim (list2, Memc[image2], SZ_FNAME) != EOF) { + + # Print messages. + if (verbose) { + call printf ("Transforming image %s to image %s\n") + call pargstr (Memc[image1]) + call pargstr (Memc[image2]) + call flush (STDOUT) + } + + # Open the images. + in = immap (Memc[image1], READ_ONLY, 0) + call imgimage (Memc[image1], Memc[str], SZ_FNAME) + call imgimage (Memc[image2], Memc[imroot], SZ_FNAME) + call imgsection (Memc[image2], Memc[section], SZ_FNAME) + if (streq (Memc[str], Memc[imroot])) { + call strcpy (Memc[imroot], Memc[imtemp], SZ_FNAME) + call mktemp ("tmp", Memc[image2], SZ_FNAME) + } else + call strcpy (Memc[image2], Memc[imtemp], SZ_FNAME) + ifnoerr (out = immap (Memc[image2], READ_WRITE, 0)) { + mode = READ_WRITE + nc = IM_LEN(out,1) + nl = IM_LEN(out,2) + xs = INDEF + ys = INDEF + } else if (Memc[section] != EOS) { + mode = NEW_IMAGE + out = immap (Memc[imroot], NEW_IMAGE, 0) + IM_NDIM(out) = IM_NDIM(in) + if (IS_INDEFI(ncols)) + IM_LEN(out,1) = IM_LEN(in,1) + else + IM_LEN(out,1) = ncols + if (IS_INDEFI(nlines)) + IM_LEN(out,2) = IM_LEN(in,2) + else + IM_LEN(out,2) = nlines + IM_PIXTYPE(out) = IM_PIXTYPE(in) + call geo_imzero (out, GT_CONSTANT(geo)) + call imunmap (out) + out = immap (Memc[image2], READ_WRITE, 0) + nc = IM_LEN(out,1) + nl = IM_LEN(out,2) + xs = INDEF + ys = INDEF + } else { + mode = NEW_COPY + out = immap (Memc[image2], NEW_COPY, in) + nc = ncols + nl = nlines + xs = xscale + ys = yscale + } + + # Set the geometry parameters. + call geo_set (geo, xmin, xmax, ymin, ymax, xs, ys, nc, nl, xin, + yin, xshift, yshift, xout, yout, xmag, ymag, xrotation, + yrotation) + + # Get the coordinate surfaces. + if (GT_GEOMODE(geo) == GT_NONE) { + call geo_format (in, out, geo, sx1, sy1, sx2, sy2) + if (verbose) { + call geo_lcoeffr (sx1, sy1, txshift, tyshift, txmag, + tymag, txrot, tyrot) + call printf (" xshift: %.2f yshift: %.2f ") + call pargr (txshift) + call pargr (tyshift) + call printf ("xmag: %.2f ymag: %.2f ") + call pargr (txmag) + call pargr (tymag) + call printf ("xrot: %.2f yrot: %.2f\n") + call pargr (txrot) + call pargr (tyrot) + call flush (STDOUT) + } + } else { + if (imtgetim (tflist, Memc[str], SZ_FNAME) != EOF) + call strcpy (Memc[str], Memc[record], SZ_FNAME) + call geo_dformat (in, out, geo, Memc[database], Memc[record], + sx1, sy1, sx2, sy2) + if (verbose) { + call printf (" Using transform %s in database %s\n") + call pargstr (Memc[record]) + call pargstr (Memc[database]) + call flush (STDOUT) + } + } + + # Transform the image. + if (IM_LEN(out,1) <= nxblock && IM_LEN(out,2) <= nyblock) { + if (GT_XSAMPLE(geo) > 1.0 || GT_YSAMPLE(geo) > 1.0) + call geo_simtran (in, out, geo, sx1, sy1, sx2, sy2) + else + call geo_imtran (in, out, geo, sx1, sy1, sx2, sy2) + } else { + if (GT_XSAMPLE(geo) > 1.0 || GT_YSAMPLE(geo) > 1.0) { + if (IM_NDIM(out) == 1) + call geo_stran (in, out, geo, sx1, sy1, sx2, sy2, + nxblock, 1) + else + call geo_stran (in, out, geo, sx1, sy1, sx2, sy2, + nxblock, nyblock) + } else { + if (IM_NDIM(out) == 1) + call geo_tran (in, out, geo, sx1, sy1, sx2, sy2, + nxblock, 1) + else + call geo_tran (in, out, geo, sx1, sy1, sx2, sy2, + nxblock, nyblock) + } + } + + # Update the linear part of the wcs. + if (!envgetb ("nomwcs") && mode == NEW_COPY) { + ndim = IM_NDIM(in) + mw = mw_openim (in) + call geo_gwcs (geo, sx1, sy1, oltm, oltv) + call mw_invertd (oltm, nltm, ndim) + call mw_vmuld (nltm, oltv, nltv, ndim) + call anegd (nltv, nltv, ndim) + call geo_swcs (mw, nltm, nltv, ndim) + call mw_saveim (mw, out) + call mw_close (mw) + } + + # Free the surfaces. + call gsfree (sx1) + call gsfree (sy1) + call gsfree (sx2) + call gsfree (sy2) + + # Close the images. + call imunmap (in) + call imunmap (out) + + call xt_delimtemp (Memc[image2], Memc[imtemp]) + } + + # Clean up. + call sfree (sp) + if (tflist != NULL) + call imtclose (tflist) + call imtclose (list1) + call imtclose (list2) +end + + +# GEO_IMZERO -- Create a dummy output image filled with the constant boundary +# extension value. + +procedure geo_imzero (im, constant) + +pointer im #I pointer to the input image +real constant #I the constant value to insert in the imagw + +int npix +pointer sp, v, buf +int impnls(), impnll(), impnlr(), impnld(), impnlx() + +begin + # Setup start vector for sequential reads and writes. + call smark (sp) + call salloc (v, IM_MAXDIM, TY_LONG) + call amovkl (long(1), Meml[v], IM_MAXDIM) + + # Initialize the image. + npix = IM_LEN(im, 1) + switch (IM_PIXTYPE(im)) { + case TY_SHORT: + while (impnls (im, buf, Meml[v]) != EOF) + call amovks (short (constant), Mems[buf], npix) + case TY_USHORT, TY_INT, TY_LONG: + while (impnll (im, buf, Meml[v]) != EOF) + call amovkl (long (constant), Meml[buf], npix) + case TY_REAL: + while (impnlr (im, buf, Meml[v]) != EOF) + call amovkr (constant, Memr[buf], npix) + case TY_DOUBLE: + while (impnld (im, buf, Meml[v]) != EOF) + call amovkd (double (constant), Memd[buf], npix) + case TY_COMPLEX: + while (impnlx (im, buf, Meml[v]) != EOF) + call amovkx (complex (constant, 0.0), Memx[buf], npix) + default: + call error (1, "Unknown pixel datatype") + } + + call sfree (sp) +end + + +# GEO_SET -- Set the image dependent task parameters individually for each +# image. + +procedure geo_set (geo, xmin, xmax, ymin, ymax, xscale, yscale, ncols, nlines, + xin, yin, xshift, yshift, xout, yout, xmag, ymag, xrotation, yrotation) + +pointer geo #I pointer to geotran structure +real xmin, xmax #I minimum and maximum reference values +real ymin, ymax #I minimum and maximum reference values +real xscale, yscale #I output picture scale +int ncols, nlines #I output picture size +real xin, yin #I input picture pixel coordinates +real xshift, yshift #I shift of origin +real xout, yout #I corresponding output picture coords +real xmag, ymag #I input picture scale +real xrotation, yrotation #I scale angle + +begin + # Set the output picture format parameters. + GT_XMIN(geo) = xmin + GT_XMAX(geo) = xmax + GT_YMIN(geo) = ymin + GT_YMAX(geo) = ymax + GT_XSCALE(geo) = xscale + GT_YSCALE(geo) = yscale + GT_NCOLS(geo) = ncols + GT_NLINES(geo) = nlines + + # Set the transformation parameters. + GT_XIN(geo) = xin + GT_YIN(geo) = yin + GT_XSHIFT(geo) = xshift + GT_YSHIFT(geo) = yshift + GT_XOUT(geo) = xout + GT_YOUT(geo) = yout + GT_XMAG(geo) = xmag + GT_YMAG(geo) = ymag + GT_XROTATION(geo) = xrotation + GT_YROTATION(geo) = yrotation +end + + +# GEO_FORMAT -- Format the output picture when there is no database file. + +procedure geo_format (in, out, geo, sx1, sy1, sx2, sy2) + +pointer in #I pointer to the input image +pointer out #I pointer to the ouput image +pointer geo #I pointer to the geotran structure +pointer sx1, sy1 #O pointer to linear surfaces +pointer sx2, sy2 #O pointer to distortion surfaces + +real xmax, ymax + +begin + # Get the scale transformation parameters. + if (IS_INDEFR(GT_XMAG(geo))) + GT_XMAG(geo) = 1. + if (IM_NDIM(in) == 1) + GT_YMAG(geo) = 1. + else if (IS_INDEFR(GT_YMAG(geo))) + GT_YMAG(geo) = 1. + + # Get the rotate transformation parameters. + if (IM_NDIM(in) == 1) + GT_XROTATION(geo) = DEGTORAD(0.) + else if (IS_INDEFR(GT_XROTATION(geo))) + GT_XROTATION(geo) = DEGTORAD(0.) + else + GT_XROTATION(geo) = DEGTORAD(GT_XROTATION(geo)) + if (IM_NDIM(in) == 1) + GT_YROTATION(geo) = DEGTORAD(0.) + else if (IS_INDEFR(GT_YROTATION(geo))) + GT_YROTATION(geo) = DEGTORAD(0.) + else + GT_YROTATION(geo) = DEGTORAD(GT_YROTATION(geo)) + + # Automatically compute the maximum extent of the image. + if (GT_XMAX(geo) <= 0.0 || GT_YMAX(geo) <= 0.0) { + + # Compute the size of the output image. + xmax = abs (cos(GT_XROTATION(geo)) * IM_LEN(in,1) / + GT_XMAG(geo)) + abs(sin(GT_YROTATION(geo)) * IM_LEN(in,2) / + GT_YMAG(geo)) + ymax = abs (sin(GT_XROTATION(geo)) * IM_LEN(in, 1) / + GT_XMAG(geo)) + abs (cos(GT_YROTATION(geo)) * IM_LEN(in,2) / + GT_YMAG(geo)) + } + + # Set up the x reference coordinate limits. + if (IS_INDEF(GT_XMIN(geo))) + GT_XMIN(geo) = 1. + else + GT_XMIN(geo) = max (1.0, GT_XMIN(geo)) + if (IS_INDEF(GT_XMAX(geo))) + GT_XMAX(geo) = IM_LEN(in,1) + else if (GT_XMAX(geo) <= 0.0) + #GT_XMAX(geo) = int (xmax + 1.0) + GT_XMAX(geo) = xmax + + # Set up the y reference coordinate limits. + if (IS_INDEF(GT_YMIN(geo))) + GT_YMIN(geo) = 1. + else + GT_YMIN(geo) = max (1.0, GT_YMIN(geo)) + if (IS_INDEF(GT_YMAX(geo))) + GT_YMAX(geo) = IM_LEN(in, 2) + else if (GT_YMAX(geo) <= 0.0) + #GT_YMAX(geo) = int (ymax + 1.0) + GT_YMAX(geo) = ymax + + # Set the number of columns and rows. + if (IS_INDEFI(GT_NCOLS(geo))) + GT_NCOLS(geo) = IM_LEN(in, 1) + if (IM_NDIM(in) == 1) + GT_NLINES(geo) = 1 + else if (IS_INDEFI(GT_NLINES(geo))) + GT_NLINES(geo) = IM_LEN(in, 2) + + # Set scale, overiding number of columns and rows if necessary. + if (IS_INDEFR(GT_XSCALE(geo))) + GT_XSCALE(geo) = (GT_XMAX(geo) - GT_XMIN(geo)) / (GT_NCOLS(geo) - 1) + else + GT_NCOLS(geo) = (GT_XMAX(geo) - GT_XMIN(geo)) / GT_XSCALE(geo) + 1 + if (IM_NDIM(in) == 1) + GT_YSCALE(geo) = 1.0 + else if (IS_INDEFR(GT_YSCALE(geo))) + GT_YSCALE(geo) = (GT_YMAX(geo) - GT_YMIN(geo)) / + (GT_NLINES(geo) - 1) + else + GT_NLINES(geo) = (GT_YMAX(geo) - GT_YMIN(geo)) / GT_YSCALE(geo) + 1 + IM_LEN(out, 1) = GT_NCOLS(geo) + IM_LEN(out, 2) = GT_NLINES(geo) + + # Set up the surfaces, distortion surfaces are NULL. + if (IM_NDIM(in) == 1) { + call gsinit (sx1, GS_POLYNOMIAL, 2, 2, GS_XNONE, GT_XMIN(geo), + GT_XMAX(geo), 0.5, 1.5) + call gsinit (sy1, GS_POLYNOMIAL, 2, 2, GS_XNONE, GT_XMIN(geo), + GT_XMAX(geo), 0.5, 1.5) + } else { + call gsinit (sx1, GS_POLYNOMIAL, 2, 2, GS_XNONE, GT_XMIN(geo), + GT_XMAX(geo), GT_YMIN(geo), GT_YMAX(geo)) + call gsinit (sy1, GS_POLYNOMIAL, 2, 2, GS_XNONE, GT_XMIN(geo), + GT_XMAX(geo), GT_YMIN(geo), GT_YMAX(geo)) + } + sx2 = NULL + sy2 = NULL + + # Adjust rotation, x and y scale, scale angle, and flip. + call geo_rotmagr (sx1, sy1, GT_XMAG(geo), GT_YMAG(geo), + GT_XROTATION(geo), GT_YROTATION(geo)) + + # Adjust the shift. + call geo_shift (in, out, geo, sx1, sy1) +end + + +# GEO_DFORMAT -- Get the coordinate transformation from a database file. + +procedure geo_dformat (in, out, geo, database, transform, sx1, sy1, sx2, sy2) + +pointer in, out #I pointers to input and output images +pointer geo #I pointer to geotran structure +char database[ARB] #I name of database file +char transform[ARB] #I name of transform +pointer sx1, sy1 #O pointer to linear part of surface fit +pointer sx2, sy2 #O pointer to higher order surface + +int i, dt, rec, ncoeff, junk +pointer xcoeff, ycoeff, newsx1, newsy1 +int dtmap(), dtlocate(), dtgeti(), dtscan() +errchk gsrestore + +begin + # Map the database and locate the transformation record. + dt = dtmap (database, READ_ONLY) + rec = dtlocate (dt, transform) + + # Get the linear part of the fit. + ncoeff = dtgeti (dt, rec, "surface1") + call malloc (xcoeff, ncoeff, TY_REAL) + call malloc (ycoeff, ncoeff, TY_REAL) + do i = 1, ncoeff { + junk = dtscan (dt) + call gargr (Memr[xcoeff+i-1]) + call gargr (Memr[ycoeff+i-1]) + } + call gsrestore (sx1, Memr[xcoeff]) + call gsrestore (sy1, Memr[ycoeff]) + + # Set the output image format parameters. + call geo_dout (in, out, geo, sx1, sy1) + + # Adjust the linear part of the fit. + call gscopy (sx1, newsx1) + call gscopy (sy1, newsy1) + if (GT_GEOMODE(geo) == GT_DISTORT) + call geo_rotmagr (newsx1, newsy1, 1.0, 1.0, 0.0, 0.0) + else if (! IS_INDEFR(GT_XMAG(geo)) || ! IS_INDEFR(GT_YMAG(geo)) || + ! IS_INDEFR(GT_XROTATION(geo)) || ! IS_INDEFR(GT_YROTATION(geo))) + call geo_dcoeff (geo, dt, rec, newsx1, newsy1) + call geo_dshift (in, out, dt, rec, geo, newsx1, newsy1) + + # Get the higher order part of the fit. + ncoeff = dtgeti (dt, rec, "surface2") + if (ncoeff > 0 && (GT_GEOMODE(geo) == GT_GEOMETRIC || GT_GEOMODE(geo) == + GT_DISTORT)) { + + # Get the distortion coefficients. + call realloc (xcoeff, ncoeff, TY_REAL) + call realloc (ycoeff, ncoeff, TY_REAL) + do i = 1, ncoeff { + junk = dtscan(dt) + call gargr (Memr[xcoeff+i-1]) + call gargr (Memr[ycoeff+i-1]) + } + iferr { + call gsrestore (sx2, Memr[xcoeff]) + } then { + call mfree (sx2, TY_STRUCT) + sx2 = NULL + } + iferr { + call gsrestore (sy2, Memr[ycoeff]) + } then { + call mfree (sy2, TY_STRUCT) + sy2 = NULL + } + + } else { + + sx2 = NULL + sy2 = NULL + } + + # Redefine the surfaces. + call gsfree (sx1) + call gscopy (newsx1, sx1) + call gsfree (newsx1) + call gsfree (sy1) + call gscopy (newsy1, sy1) + call gsfree (newsy1) + + # Cleanup. + call mfree (xcoeff, TY_REAL) + call mfree (ycoeff, TY_REAL) + call dtunmap (dt) +end + + +# GEO_DOUT -- Set the output image format using information in the database +# file. + +procedure geo_dout (in, out, geo, sx1, sy1) + +pointer in, out #I pointers to input and output image +pointer geo #I pointer to geotran sturcture +pointer sx1, sy1 #I pointers to linear surface descriptors + +real gsgetr () + +begin + # Set the reference coordinate limits. + if (IS_INDEFR(GT_XMIN(geo))) + GT_XMIN(geo) = gsgetr (sx1, GSXMIN) + if (IS_INDEFR(GT_XMAX(geo))) + GT_XMAX(geo) = gsgetr (sx1, GSXMAX) + if (IS_INDEFR(GT_YMIN(geo))) + GT_YMIN(geo) = gsgetr (sy1, GSYMIN) + if (IS_INDEFR(GT_YMAX(geo))) + GT_YMAX(geo) = gsgetr (sy1, GSYMAX) + + # Set the number of lines and columns. + if (IS_INDEFI(GT_NCOLS(geo))) + GT_NCOLS(geo) = IM_LEN(in, 1) + if (IM_NDIM(in) == 1) + GT_NLINES(geo) = 1 + else if (IS_INDEFI(GT_NLINES(geo))) + GT_NLINES(geo) = IM_LEN(in, 2) + + # Set scale, overiding the number of columns and rows if necessary. + if (IS_INDEFR(GT_XSCALE(geo))) + GT_XSCALE(geo) = (GT_XMAX(geo) - GT_XMIN(geo)) / (GT_NCOLS(geo) - 1) + else + GT_NCOLS(geo) = abs ((GT_XMAX(geo) - GT_XMIN(geo)) / + GT_XSCALE(geo)) + 1 + if (IM_NDIM(in) == 1) + GT_YSCALE(geo) = 1.0 + else if (IS_INDEFR(GT_YSCALE(geo))) + GT_YSCALE(geo) = (GT_YMAX(geo) - GT_YMIN(geo)) / + (GT_NLINES(geo) - 1) + else + GT_NLINES(geo) = abs ((GT_YMAX(geo) - GT_YMIN(geo)) / + GT_YSCALE(geo)) + 1 + + # Set the output image size. + IM_LEN(out,1) = GT_NCOLS(geo) + IM_LEN(out,2) = GT_NLINES(geo) +end + + +# GEO_DSHIFT -- Adjust the shifts using information in the database file. + +procedure geo_dshift (in, out, dt, rec, geo, sx1, sy1) + +pointer in, out #I pointer to input and output images +pointer dt #I pointer to database +int rec #I pointer to database record +pointer geo #I pointer to geotran structure +pointer sx1, sy1 #U pointers to linear surfaces + +real gseval() + +begin + # Define the output origin. + if (IS_INDEFR(GT_XOUT(geo))) + GT_XOUT(geo) = (GT_XMAX(geo) + GT_XMIN(geo)) / 2.0 + if (IS_INDEFR(GT_YOUT(geo))) + GT_YOUT(geo) = (GT_YMAX(geo) + GT_YMIN(geo)) / 2.0 + + # Define the input image origin. + if (IS_INDEFR(GT_XIN(geo))) + GT_XIN(geo) = gseval (sx1, GT_XOUT(geo), GT_YOUT(geo)) + if (IS_INDEFR(GT_YIN(geo))) + GT_YIN(geo) = gseval (sy1, GT_XOUT(geo), GT_YOUT(geo)) + + # Define the shifts. + if (IS_INDEFR(GT_XSHIFT(geo))) + GT_XSHIFT(geo) = GT_XIN(geo) - gseval (sx1, GT_XOUT(geo), + GT_YOUT(geo)) + if (IS_INDEFR(GT_YSHIFT(geo))) + GT_YSHIFT(geo) = GT_YIN(geo) - gseval (sy1, GT_XOUT(geo), + GT_YOUT(geo)) + + # Correct the coefficients. + call geo_xyshiftr (sx1, sy1, GT_XSHIFT(geo), GT_YSHIFT(geo)) +end + + +# GEO_SHIFT -- Compute the shift. + +procedure geo_shift (in, out, geo, sx1, sy1) + +pointer in, out #I pointer to input and output images +pointer geo #I pointer to geotran structure +pointer sx1, sy1 #I pointers to linear surfaces + +real gseval() + +begin + # Determine the output origin. + if (IS_INDEFR(GT_XOUT(geo))) + GT_XOUT(geo) = (GT_XMAX(geo) + GT_XMIN(geo)) / 2.0 + if (IS_INDEFR(GT_YOUT(geo))) + GT_YOUT(geo) = (GT_YMAX(geo) + GT_YMIN(geo)) / 2.0 + + # Determine the input origin. + if (IS_INDEFR(GT_XIN(geo))) + GT_XIN(geo) = (real (IM_LEN (in, 1)) + 1.) / 2. + if (IS_INDEFR(GT_YIN(geo))) + GT_YIN(geo) = (real (IM_LEN (in, 2)) + 1.) / 2. + + # Determine the final x and y shifts. + if (! IS_INDEFR(GT_XSHIFT(geo))) + GT_XOUT(geo) = GT_XIN(geo) + GT_XSHIFT(geo) + if (! IS_INDEFR(GT_YSHIFT(geo))) + GT_YOUT(geo) = GT_YIN(geo) + GT_YSHIFT(geo) + GT_XSHIFT(geo) = GT_XIN(geo) - gseval (sx1, GT_XOUT(geo), + GT_YOUT(geo)) + GT_YSHIFT(geo) = GT_YIN(geo) - gseval (sy1, GT_XOUT(geo), + GT_YOUT(geo)) + + # Alter coefficients. + call geo_xyshiftr (sx1, sy1, GT_XSHIFT(geo), GT_YSHIFT(geo)) +end + + +# GEO_DCOEFF -- Alter the linear componets of the surface fit after the fact. + +procedure geo_dcoeff (geo, dt, rec, sx1, sy1) + +pointer geo #I pointer to geotran structure +pointer dt #I pointer to database record +int rec #I database record +pointer sx1, sy1 #U pointers to the linear surface + +real dtgetr() +errchk dtgetr() + +begin + # Get the transformation parameters. + if (IS_INDEFR(GT_XMAG(geo))) { + iferr (GT_XMAG(geo) = dtgetr (dt, rec, "xmag")) + GT_XMAG(geo) = dtgetr (dt, rec, "xscale") + } + if (IS_INDEFR(GT_YMAG(geo))) { + iferr (GT_YMAG(geo) = dtgetr (dt, rec, "ymag")) + GT_YMAG(geo) = dtgetr (dt, rec, "yscale") + } + if (IS_INDEFR(GT_XROTATION(geo))) + GT_XROTATION(geo) = DEGTORAD(dtgetr (dt, rec, "xrotation")) + else + GT_XROTATION(geo) = DEGTORAD(GT_XROTATION(geo)) + if (IS_INDEFR(GT_YROTATION(geo))) + GT_YROTATION(geo) = DEGTORAD(dtgetr (dt, rec, "yrotation")) + else + GT_YROTATION(geo) = DEGTORAD(GT_YROTATION(geo)) + + call geo_rotmagr (sx1, sy1, GT_XMAG(geo), GT_YMAG(geo), + GT_XROTATION(geo), GT_YROTATION(geo)) +end + + +# GEO_GWCS -- Compute the ltm and ltv vectors using the GEOTRAN coordinate +# surfaces. + +procedure geo_gwcs (geo, sx1, sy1, ltm, ltv) + +pointer geo # pointer to the geotran structure +pointer sx1 # pointer to the linear x coordinate surface +pointer sy1 # pointer to the linear y coordinate surface +double ltm[2,2] # rotation matrix +double ltv[2] # shift vector + +double xscale, yscale, xmin, ymin +int ncoeff +pointer sp, xcoeff, ycoeff +real xrange, yrange +int gsgeti() +real gsgetr() + +begin + # Allocate space for the coefficients. + call smark (sp) + ncoeff = max (gsgeti (sx1, GSNSAVE), gsgeti (sy1, GSNSAVE)) + call salloc (xcoeff, ncoeff, TY_REAL) + call salloc (ycoeff, ncoeff, TY_REAL) + + # Fetch the coefficients. + call gssave (sx1, Memr[xcoeff]) + call gssave (sy1, Memr[ycoeff]) + + # Denormalize the coefficients for non-polynomial functions. + xrange = gsgetr (sx1, GSXMAX) - gsgetr (sx1, GSXMIN) + yrange = gsgetr (sy1, GSYMAX) - gsgetr (sy1, GSYMIN) + if (gsgeti (sx1, GSTYPE) != GS_POLYNOMIAL) { + Memr[xcoeff+GS_SAVECOEFF+1] = Memr[xcoeff+GS_SAVECOEFF+1] * 2. / + xrange + Memr[xcoeff+GS_SAVECOEFF+2] = Memr[xcoeff+GS_SAVECOEFF+2] * 2. / + yrange + } + if (gsgeti (sy1, GSTYPE) != GS_POLYNOMIAL) { + Memr[ycoeff+GS_SAVECOEFF+1] = Memr[ycoeff+GS_SAVECOEFF+1] * 2. / + xrange + Memr[ycoeff+GS_SAVECOEFF+2] = Memr[ycoeff+GS_SAVECOEFF+2] * 2. / + yrange + } + + # Set the shift vector. + ltv[1] = Memr[xcoeff+GS_SAVECOEFF] + ltv[2] = Memr[ycoeff+GS_SAVECOEFF] + + # Set the rotation vector. + ltm[1,1] = Memr[xcoeff+GS_SAVECOEFF+1] + ltm[2,1] = Memr[xcoeff+GS_SAVECOEFF+2] + ltm[1,2] = Memr[ycoeff+GS_SAVECOEFF+1] + ltm[2,2] = Memr[ycoeff+GS_SAVECOEFF+2] + + # Get the sign of the scale vector which is always +ve. + xmin = GT_XMIN(geo) + ymin = GT_YMIN(geo) + if (GT_XMIN(geo) > GT_XMAX(geo)) + xscale = -GT_XSCALE(geo) + else + xscale = GT_XSCALE(geo) + if (GT_YMIN(geo) > GT_YMAX(geo)) + yscale = -GT_YSCALE(geo) + else + yscale = GT_YSCALE(geo) + + # Correct for reference units that are not in pixels. + ltv[1] = ltv[1] + ltm[1,1] * xmin + ltm[2,1] * ymin - ltm[1,1] * + xscale - ltm[2,1] * yscale + ltv[2] = ltv[2] + ltm[1,2] * xmin + ltm[2,2] * ymin - ltm[1,2] * + xscale - ltm[2,2] * yscale + ltm[1,1] = ltm[1,1] * xscale + ltm[2,1] = ltm[2,1] * yscale + ltm[1,2] = ltm[1,2] * xscale + ltm[2,2] = ltm[2,2] * yscale + + call sfree (sp) +end + + +define LTM Memd[ltm+(($2)-1)*pdim+($1)-1] + +# GEO_SWCS -- Update the wcs and write it to the image header. + +procedure geo_swcs (mw, gltm, gltv, ldim) + +pointer mw # the mwcs descriptor +double gltm[ldim,ldim] # the input cd matrix from geotran +double gltv[ldim] # the input shift vector from geotran +int ldim # number of logical dimensions + +int axes[IM_MAXDIM], naxes, pdim, nelem, axmap, ax1, ax2 +pointer sp, ltm, ltv_1, ltv_2 +int mw_stati() + +begin + # Convert axis bitflags to the axis lists. + if (ldim == 1) { + call mw_gaxlist (mw, 01B, axes, naxes) + if (naxes < 1) + return + } else { + call mw_gaxlist (mw, 03B, axes, naxes) + if (naxes < 2) + return + } + + # Initialize the parameters. + pdim = mw_stati (mw, MW_NDIM) + nelem = pdim * pdim + axmap = mw_stati (mw, MW_USEAXMAP) + call mw_seti (mw, MW_USEAXMAP, NO) + + # Allocate working space. + call smark (sp) + call salloc (ltm, nelem, TY_DOUBLE) + call salloc (ltv_1, pdim, TY_DOUBLE) + call salloc (ltv_2, pdim, TY_DOUBLE) + + # Initialize the vectors and matrices. + call mw_mkidmd (Memd[ltm], pdim) + call aclrd (Memd[ltv_1], pdim) + call aclrd (Memd[ltv_2], pdim) + + # Enter the linear operation. + ax1 = axes[1] + Memd[ltv_2+ax1-1] = gltv[1] + LTM(ax1,ax1) = gltm[1,1] + if (ldim == 2) { + ax2 = axes[2] + Memd[ltv_2+ax2-1] = gltv[2] + LTM(ax2,ax1) = gltm[2,1] + LTM(ax1,ax2) = gltm[1,2] + LTM(ax2,ax2) = gltm[2,2] + } + + # Perform the translation. + call mw_translated (mw, Memd[ltv_1], Memd[ltm], Memd[ltv_2], pdim) + + call sfree (sp) + call mw_seti (mw, MW_USEAXMAP, axmap) +end diff --git a/pkg/images/immatch/src/geometry/t_geoxytran.x b/pkg/images/immatch/src/geometry/t_geoxytran.x new file mode 100644 index 00000000..c99b9a0c --- /dev/null +++ b/pkg/images/immatch/src/geometry/t_geoxytran.x @@ -0,0 +1,343 @@ +include <fset.h> +include <ctype.h> +include <math/gsurfit.h> + +define MAX_FIELDS 100 # Maximum number of fields in list +define TABSIZE 8 # Spacing of tab stops + +# Define the permitted computation types +define GEO_REAL 1 # Computation type is real +define GEO_DOUBLE 2 # Computation type is double + +# T_GEOXYTRAN -- Transform a list of x and y coordinates using the geometric +# transformation operations computed by the GEOMAP task. + +procedure t_geoxytran() + +int inlist, outlist, reclist, calctype, geometry, dir, xcolumn, ycolumn +int min_sigdigits, infd, outfd +pointer sp, in_fname, out_fname, record, xformat, yformat, str, dt +pointer sx1, sy1, sx2, sy2 +int clgwrd(), clgeti(), open() +bool streq() +int fntopnb(), fntlenb(), fntgfnb(), imtopenp(), imtlen(), imtgetim() +pointer dtmap() + +begin + # Allocate memory for transformation parameters structure + call smark (sp) + call salloc (in_fname, SZ_FNAME, TY_CHAR) + call salloc (out_fname, SZ_FNAME, TY_CHAR) + call salloc (record, SZ_FNAME, TY_CHAR) + call salloc (xformat, SZ_FNAME, TY_CHAR) + call salloc (yformat, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Open the input and output file lists. + call clgstr ("input", Memc[str], SZ_FNAME) + if (Memc[str] == EOS) + call strcpy ("STDIN", Memc[str], SZ_FNAME) + inlist = fntopnb(Memc[str], NO) + call clgstr ("output", Memc[str], SZ_FNAME) + if (Memc[str] == EOS) + call strcpy ("STDOUT", Memc[str], SZ_FNAME) + outlist = fntopnb (Memc[str], NO) + call clgstr ("database", Memc[str], SZ_FNAME) + if (Memc[str] != EOS) { + dt = dtmap (Memc[str], READ_ONLY) + reclist = imtopenp ("transforms") + } else { + dt = NULL + reclist = NULL + } + + # Test the input and out file and record lists for validity. + if (fntlenb(inlist) <= 0) + call error (0, "The input file list is empty") + if (fntlenb(outlist) <= 0) + call error (0, "The output file list is empty") + if (fntlenb(outlist) > 1 && fntlenb(outlist) != fntlenb(inlist)) + call error (0, + "Input and output file lists are not the same length") + if (dt != NULL && reclist != NULL) { + if (imtlen (reclist) > 1 && imtlen (reclist) != fntlenb (inlist)) + call error (0, + "Input file and record lists are not the same length.") + } + + # Get geometry and transformation direction. + geometry = clgwrd ("geometry", Memc[str], SZ_LINE, + ",linear,distortion,geometric,") + dir = clgwrd ("direction", Memc[str], SZ_LINE, + ",forward,backward,") + + # Get field numbers from cl + if (dir == 1) + calctype = clgwrd ("calctype", Memc[str], SZ_LINE, + ",real,double,") + else + calctype = GEO_DOUBLE + xcolumn = clgeti ("xcolumn") + ycolumn = clgeti ("ycolumn") + call clgstr ("xformat", Memc[xformat], SZ_FNAME) + call clgstr ("yformat", Memc[yformat], SZ_FNAME) + min_sigdigits = clgeti ("min_sigdigits") + + # Get the output file name. + if (fntgfnb (outlist, Memc[out_fname], SZ_FNAME) == EOF) + call strcpy ("STDOUT", Memc[out_fname], SZ_FNAME) + outfd = open (Memc[out_fname], NEW_FILE, TEXT_FILE) + if (streq (Memc[out_fname], "STDOUT") || outfd == STDOUT) + call fseti (outfd, F_FLUSHNL, YES) + + # Get the record name. + if (reclist == NULL) + Memc[record] = EOS + else if (imtgetim (reclist, Memc[record], SZ_FNAME) == EOF) + Memc[record] = EOS + + # Call procedure to get parameters and fill structure. + sx1 = NULL; sy1 = NULL; sx2 = NULL; sy2 = NULL + call geo_init_transform (dt, Memc[record], calctype, geometry, + sx1, sy1, sx2, sy2) + + # While input list is not depleted, open file and transform list. + while (fntgfnb (inlist, Memc[in_fname], SZ_FNAME) != EOF) { + + infd = open (Memc[in_fname], READ_ONLY, TEXT_FILE) + + # Transform the coordinates. + call geo_transform_file (infd, outfd, xcolumn, ycolumn, dir, + calctype, Memc[xformat], Memc[yformat], min_sigdigits, + sx1, sy1, sx2, sy2) + + # Do not get a new output file name if there is not output + # file list or if only one output file was specified. + # Otherwise fetch the new name. + if (fntlenb(outlist) > 1) { + call close (outfd) + if (fntgfnb (outlist, Memc[out_fname], SZ_FNAME) != EOF) + outfd = open (Memc[out_fname], NEW_FILE, TEXT_FILE) + if (streq (Memc[out_fname], "STDOUT") || outfd == STDOUT) + call fseti (outfd, F_FLUSHNL, YES) + } + + call close (infd) + + # Do not reset the transformation if there is no record list + # or only one record is specified. Otherwise fetch the next + # record name. + if (reclist != NULL && imtlen (reclist) > 1) { + if (imtgetim (reclist, Memc[record], SZ_FNAME) != EOF) { + call geo_free_transform (calctype, sx1, sy1, sx2, sy2) + call geo_init_transform (dt, Memc[record], calctype, + geometry, sx1, sy1, sx2, sy2) + } + } + } + + # Free the surface descriptors. + call geo_free_transform (calctype, sx1, sy1, sx2, sy2) + + # Close up file and record templates. + if (dt != NULL) + call dtunmap (dt) + call close (outfd) + call fntclsb (inlist) + call fntclsb (outlist) + if (reclist != NULL) + call imtclose (reclist) + call sfree (sp) +end + + +# GEO_INIT_TRANSFORM -- gets parameter values relevant to the +# transformation from the cl. List entries will be transformed +# in procedure rg_transform. + +procedure geo_init_transform (dt, record, calctype, geometry, sx1, sy1, + sx2, sy2) + +pointer dt #I pointer to database file produced by geomap +char record[ARB] #I the name of the database record +int calctype #I the computation data type +int geometry #I the type of geometry to be computed +pointer sx1, sy1 #O pointers to the linear x and y surfaces +pointer sx2, sy2 #O pointers to the x and y distortion surfaces + +begin + if (dt == NULL) { + + if (calctype == GEO_REAL) + call geo_linitr (sx1, sy1, sx2, sy2) + else + call geo_linitd (sx1, sy1, sx2, sy2) + + } else { + + if (calctype == GEO_REAL) + call geo_sinitr (dt, record, geometry, sx1, sy1, + sx2, sy2) + else + call geo_sinitd (dt, record, geometry, sx1, sy1, + sx2, sy2) + } +end + + +# GEO_FREE_TRANSFORM -- Free the previously defined transformation + +procedure geo_free_transform (calctype, sx1, sy1, sx2, sy2) + +int calctype #I the computation data type +pointer sx1, sy1 #O pointers to the linear x and y surfaces +pointer sx2, sy2 #O pointers to the x and y distortion surfaces + +begin + if (calctype == GEO_REAL) + call geo_sfreer (sx1, sy1, sx2, sy2) + else + call geo_sfreed (sx1, sy1, sx2, sy2) +end + + +# GEO_TRANSFORM_FILE -- This procedure is called once for each file +# in the input list. For each line in the input file that isn't +# blank or comment, the line is transformed. Blank and comment +# lines are output unaltered. + +procedure geo_transform_file (infd, outfd, xfield, yfield, dir, calctype, + xformat, yformat, min_sigdigits, sx1, sy1, sx2, sy2) + +int infd #I the input file descriptor +int outfd #I the output file descriptor +int xfield #I the x column number +int yfield #I the y column number +int dir #I transform direction +int calctype #I the computation type +char xformat[ARB] #I output format of the x coordinate +char yformat[ARB] #I output format of the y coordinate +int min_sigdigits #I the minimum number of digits to be output +pointer sx1, sy1 #I pointers to the linear x and y surfaces +pointer sx2, sy2 #I pointers to the x and y distortion surfaces + +double xd, yd, xtd, ytd +int max_fields, nline, nfields, nchars, nsdig_x, nsdig_y, offset +real xr, yr, xtr, ytr +pointer sp, inbuf, linebuf, field_pos, outbuf, ip +int getline(), li_get_numr(), li_get_numd() + +int nsx, nsy +double der[8], xmin, xmax, ymin, ymax, tol +pointer sx[2], sy[2] +double dgsgetd() + +#double x, y, xt, yt + +begin + call smark (sp) + call salloc (inbuf, SZ_LINE, TY_CHAR) + call salloc (linebuf, SZ_LINE, TY_CHAR) + call salloc (field_pos, MAX_FIELDS, TY_INT) + call salloc (outbuf, SZ_LINE, TY_CHAR) + + max_fields = MAX_FIELDS + + # Initialize for backward transform. + if (dir == 2) { + sx[1] = sx1; sy[1] = sy1; sx[2] = sx2; sy[2] = sy2 + nsx = 2; nsy = 2 + if (sx2 == NULL) + nsx = 1 + if (sy2 == NULL) + nsy = 1 + xmin = dgsgetd (sx1, GSXMIN) + xmax = dgsgetd (sx1, GSXMAX) + ymin = dgsgetd (sx1, GSYMIN) + ymax = dgsgetd (sx1, GSYMAX) + tol = abs (xmax - xmin) / 1E10 + xd = (xmin + xmax) / 2 + yd = (ymin + ymax) / 2 + call tr_init (sx, nsx, sy, nsy, xd, yd, der) + } + + for (nline=1; getline (infd, Memc[inbuf]) != EOF; nline = nline + 1) { + for (ip=inbuf; IS_WHITE(Memc[ip]); ip=ip+1) + ; + if (Memc[ip] == '#') { + # Pass comment lines on to the output unchanged. + call putline (outfd, Memc[inbuf]) + next + } else if (Memc[ip] == '\n' || Memc[ip] == EOS) { + # Blank lines too. + call putline (outfd, Memc[inbuf]) + next + } + + # Expand tabs into blanks, determine field offsets. + call strdetab (Memc[inbuf], Memc[linebuf], SZ_LINE, TABSIZE) + call li_find_fields (Memc[linebuf], Memi[field_pos], max_fields, + nfields) + + if (xfield > nfields || yfield > nfields) { + call fstats (infd, F_FILENAME, Memc[outbuf], SZ_LINE) + call eprintf ("Not enough fields in file %s line %d\n") + call pargstr (Memc[outbuf]) + call pargi (nline) + call putline (outfd, Memc[linebuf]) + next + } + + offset = Memi[field_pos+xfield-1] + if (calctype == GEO_REAL) + nchars = li_get_numr (Memc[linebuf+offset-1], xr, nsdig_x) + else + nchars = li_get_numd (Memc[linebuf+offset-1], xd, nsdig_x) + if (nchars == 0) { + call fstats (infd, F_FILENAME, Memc[outbuf], SZ_LINE) + call eprintf ("Bad x value in file '%s' at line %d:\n") + call pargstr (Memc[outbuf]) + call pargi (nline) + call putline (outfd, Memc[linebuf]) + next + } + + offset = Memi[field_pos+yfield-1] + if (calctype == GEO_REAL) + nchars = li_get_numr (Memc[linebuf+offset-1], yr, nsdig_y) + else + nchars = li_get_numd (Memc[linebuf+offset-1], yd, nsdig_y) + if (nchars == 0) { + call fstats (infd, F_FILENAME, Memc[outbuf], SZ_LINE) + call eprintf ("Bad y value in file '%s' at line %d:\n") + call pargstr (Memc[outbuf]) + call pargi (nline) + call putline (outfd, Memc[linebuf]) + next + } + + if (calctype == GEO_REAL) { + call geo_do_transformr (xr, yr, xtr, ytr, + sx1, sy1, sx2, sy2) + call li_pack_liner (Memc[linebuf], Memc[outbuf], SZ_LINE, + Memi[field_pos], nfields, xfield, yfield, xtr, ytr, + xformat, yformat, nsdig_x, nsdig_y, min_sigdigits) + + } else { + if (dir == 1) + call geo_do_transformd (xd, yd, xtd, ytd, + sx1, sy1, sx2, sy2) + else + call tr_invert (sx, nsx, sy, nsy, xd, yd, xtd, ytd, + der, xmin, xmax, ymin, ymax, tol) + call li_pack_lined (Memc[linebuf], Memc[outbuf], SZ_LINE, + Memi[field_pos], nfields, xfield, yfield, xtd, ytd, + xformat, yformat, nsdig_x, nsdig_y, min_sigdigits) + } + + call putline (outfd, Memc[outbuf]) + } + + call sfree (sp) +end + diff --git a/pkg/images/immatch/src/geometry/trinvert.x b/pkg/images/immatch/src/geometry/trinvert.x new file mode 100644 index 00000000..5f75cdc2 --- /dev/null +++ b/pkg/images/immatch/src/geometry/trinvert.x @@ -0,0 +1,163 @@ +# The code here is taken from t_transform.x in the longslit package. The +# changes are to use a sum instead of an average when multiple surfaces +# are given and not to use the xgs interface. Also the convergence +# tolerance is user specified since in this application the units might +# not be pixels. + + +define MAX_ITERATE 20 +define ERROR 0.05 +define FUDGE 0.5 + +# TR_INVERT -- Given user coordinate surfaces U(X,Y) and V(X,Y) +# (if none use one-to-one mapping and if more than one sum) +# corresponding to a given U and V and also the various partial +# derivatives. This is done using a gradient following interative +# method based on evaluating the partial derivative at each point +# and solving the linear Taylor expansions simultaneously. The last +# point sampled is used as the starting point. Thus, if the +# input U and V progress smoothly then the number of iterations +# can be small. The output is returned in x and y and in the derivative array +# DER. A point outside of the surfaces is returned as the nearest +# point at the edge of the surfaces in the DER array. + +procedure tr_invert (usf, nusf, vsf, nvsf, u, v, x, y, der, + xmin, xmax, ymin, ymax, tol) + +pointer usf[ARB], vsf[ARB] # User coordinate surfaces U(X,Y) and V(X,Y) +int nusf, nvsf # Number of surfaces for each coordinate +double u, v # Input U and V to determine X and Y +double x, y # Output X and Y +double der[8] # Last result as input, new result as output + # 1=X, 2=Y, 3=U, 4=DUDX, 5=DUDY, 6=V, + # 7=DVDX, 8=DVDY +double xmin, xmax, ymin, ymax # Limits of coordinate surfaces. +double tol # Tolerance + +int i, j, nedge +double fudge, du, dv, dx, dy, tmp[3] + +begin + # Use the last result as the starting point for the next position. + # If this is near the desired value then the interation will converge + # quickly. Allow a iteration to go off the surface twice. + # Quit when DX and DY are within tol. + + nedge = 0 + do i = 1, MAX_ITERATE { + du = u - der[3] + dv = v - der[6] + dx = (der[8] * du - der[5] * dv) / + (der[8] * der[4] - der[5] * der[7]) + dy = (dv - der[7] * dx) / der[8] + fudge = 1 - FUDGE / i + x = der[1] + fudge * dx + y = der[2] + fudge * dy + der[1] = max (xmin, min (xmax, x)) + der[2] = max (ymin, min (ymax, y)) + if ((abs (dx) < tol) && (abs (dy) < tol)) + break + + if (nusf == 0) + der[3] = der[1] + else if (nusf == 1) { + call dgsder (usf[1], der[1], der[2], der[3], 1, 0, 0) + call dgsder (usf[1], der[1], der[2], der[4], 1, 1, 0) + call dgsder (usf[1], der[1], der[2], der[5], 1, 0, 1) + } else { + call dgsder (usf[1], der[1], der[2], der[3], 1, 0, 0) + call dgsder (usf[1], der[1], der[2], der[4], 1, 1, 0) + call dgsder (usf[1], der[1], der[2], der[5], 1, 0, 1) + do j = 2, nusf { + call dgsder (usf[j], der[1], der[2], tmp[1], 1, 0, 0) + call dgsder (usf[j], der[1], der[2], tmp[2], 1, 1, 0) + call dgsder (usf[j], der[1], der[2], tmp[3], 1, 0, 1) + der[3] = der[3] + tmp[1] + der[4] = der[4] + tmp[2] + der[5] = der[5] + tmp[3] + } + } + + if (nvsf == 0) + der[6] = der[2] + else if (nvsf == 1) { + call dgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0) + call dgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0) + call dgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1) + } else { + call dgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0) + call dgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0) + call dgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1) + do j = 2, nvsf { + call dgsder (vsf[j], der[1], der[2], tmp[1], 1, 0, 0) + call dgsder (vsf[j], der[1], der[2], tmp[2], 1, 1, 0) + call dgsder (vsf[j], der[1], der[2], tmp[3], 1, 0, 1) + der[6] = der[6] + tmp[1] + der[7] = der[7] + tmp[2] + der[8] = der[8] + tmp[3] + } + } + } +end + + +# TR_INIT -- Since the inversion iteration always begins from the last +# point we need to initialize before the first call to TR_INVERT. + +procedure tr_init (usf, nusf, vsf, nvsf, x, y, der) + +pointer usf[ARB], vsf[ARB] # User coordinate surfaces +int nusf, nvsf # Number of surfaces for each coordinate +double x, y # Starting X and Y +double der[8] # Inversion data + +int j +double tmp[3] + +begin + der[1] = x + der[2] = y + if (nusf == 0) { + der[3] = der[1] + der[4] = 1. + der[5] = 0. + } else if (nusf == 1) { + call dgsder (usf[1], der[1], der[2], der[3], 1, 0, 0) + call dgsder (usf[1], der[1], der[2], der[4], 1, 1, 0) + call dgsder (usf[1], der[1], der[2], der[5], 1, 0, 1) + } else { + call dgsder (usf[1], der[1], der[2], der[3], 1, 0, 0) + call dgsder (usf[1], der[1], der[2], der[4], 1, 1, 0) + call dgsder (usf[1], der[1], der[2], der[5], 1, 0, 1) + do j = 2, nusf { + call dgsder (usf[j], der[1], der[2], tmp[1], 1, 0, 0) + call dgsder (usf[j], der[1], der[2], tmp[2], 1, 1, 0) + call dgsder (usf[j], der[1], der[2], tmp[3], 1, 0, 1) + der[3] = der[3] + tmp[1] + der[4] = der[4] + tmp[2] + der[5] = der[5] + tmp[3] + } + } + + if (nvsf == 0) { + der[6] = der[2] + der[7] = 0. + der[8] = 1. + } else if (nvsf == 1) { + call dgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0) + call dgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0) + call dgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1) + } else { + call dgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0) + call dgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0) + call dgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1) + do j = 2, nvsf { + call dgsder (vsf[j], der[1], der[2], tmp[1], 1, 0, 0) + call dgsder (vsf[j], der[1], der[2], tmp[2], 1, 1, 0) + call dgsder (vsf[j], der[1], der[2], tmp[3], 1, 0, 1) + der[6] = der[6] + tmp[1] + der[7] = der[7] + tmp[2] + der[8] = der[8] + tmp[3] + } + } +end |