diff options
Diffstat (limited to 'pkg/images/lib/geogmap.gx')
| -rw-r--r-- | pkg/images/lib/geogmap.gx | 494 |
1 files changed, 494 insertions, 0 deletions
diff --git a/pkg/images/lib/geogmap.gx b/pkg/images/lib/geogmap.gx new file mode 100644 index 00000000..e52a129e --- /dev/null +++ b/pkg/images/lib/geogmap.gx @@ -0,0 +1,494 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <math.h> +include <math/gsurfit.h> +include <gset.h> +include "geomap.h" +include "geogmap.h" + +define GHELPFILE "images$lib/geomap.key" +define CHELPFILE "images$lib/coomap.key" + +$for (rd) + +# GEO_MGFIT -- Fit the surface using interactive graphics. + +procedure geo_mgfit$t (gd, fit, sx1, sy1, sx2, sy2, xref, yref, xin, + yin, wts, npts, xerrmsg, yerrmsg, maxch) + +pointer gd #I graphics file descriptor +pointer fit #I pointer to the fit structure +pointer sx1 #I pointer to the linear x surface fit +pointer sy1 #I pointer to the linear y surface fit +pointer sx2 #I pointer to higher order x surface fit +pointer sy2 #I pointer to higher order y surface fit +PIXEL xref[npts] #I the x reference coordinates +PIXEL yref[npts] #I the y reference coordinates +PIXEL xin[npts] #I input x coordinates +PIXEL yin[npts] #I input y coordinates +PIXEL wts[npts] #I array of weights +int npts #I number of data points +char xerrmsg[ARB] #O the output x fit error message +char yerrmsg[ARB] #O the output x fit error message +int maxch #I the size of the error messages + +char errstr[SZ_LINE] +int newgraph, delete, wcs, key, errcode +pointer sp, w, gfit, xresid, yresid, cmd +pointer gt1, gt2, gt3, gt4, gt5 +real wx, wy +PIXEL xshift, yshift, xscale, yscale, thetax, thetay + +int clgcur(), errget() +pointer gt_init() + +errchk geo_fxy$t(), geo_mreject$t(), geo_ftheta$t() +errchk geo_fmagnify$t(), geo_flinear$t() + +begin + # Initialize gfit structure and working space. + call smark (sp) + call salloc (gfit, LEN_GEOGRAPH, TY_STRUCT) + call salloc (xresid, npts, TY_PIXEL) + call salloc (yresid, npts, TY_PIXEL) + call salloc (w, npts, TY_PIXEL) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Do initial fit. + iferr { + switch (GM_FIT(fit)) { + case GM_ROTATE: + call geo_ftheta$t (fit, sx1, sy1, xref, yref, xin, yin, wts, + Mem$t[xresid], Mem$t[yresid], npts, xerrmsg, maxch, + yerrmsg, maxch) + sx2 = NULL + sy2 = NULL + case GM_RSCALE: + call geo_fmagnify$t (fit, sx1, sy1, xref, yref, xin, yin, wts, + Mem$t[xresid], Mem$t[yresid], npts, xerrmsg, maxch, + yerrmsg, maxch) + sx2 = NULL + sy2 = NULL + case GM_RXYSCALE: + call geo_flinear$t (fit, sx1, sy1, xref, yref, xin, yin, wts, + Mem$t[xresid], Mem$t[yresid], npts, xerrmsg, maxch, + yerrmsg, maxch) + sx2 = NULL + sy2 = NULL + default: + call geo_fxy$t (fit, sx1, sx2, xref, yref, xin, wts, + Mem$t[xresid], npts, YES, xerrmsg, maxch) + call geo_fxy$t (fit, sy1, sy2, xref, yref, yin, wts, + Mem$t[yresid], npts, NO, yerrmsg, maxch) + } + if (GM_MAXITER(fit) <= 0 || IS_INDEFD(GM_REJECT(fit))) + GM_NREJECT(fit) = 0 + else + call geo_mreject$t (fit, sx1, sy1, sx2, sy2, xref, yref, xin, + yin, wts, Mem$t[xresid], Mem$t[yresid], npts, xerrmsg, + maxch, yerrmsg, maxch) + } then { + call sfree (sp) + if (GM_PROJECTION(fit) == GM_NONE) + call error (2, "Too few points for X and Y fits.") + else + call error (2, "Too few points for XI and ETA fits.") + } + + GG_NEWFUNCTION(gfit) = NO + GG_FITERROR(gfit) = NO + errcode = OK + + # Set up plotting defaults. + GG_PLOTTYPE(gfit) = FIT + GG_OVERPLOT(gfit) = NO + GG_CONSTXY(gfit) = YES + newgraph = NO + + # Allocate graphics tools. + gt1 = gt_init () + gt2 = gt_init () + gt3 = gt_init () + gt4 = gt_init () + gt5 = gt_init () + + # Set the plot title and x and y axis labels. + call geo_gtset (FIT, gt1, fit) + call geo_gtset (XXRESID, gt2, fit) + call geo_gtset (XYRESID, gt3, fit) + call geo_gtset (YXRESID, gt4, fit) + call geo_gtset (YYRESID, gt5, fit) + + # Make the first plot. + call gclear (gd) + call geo_label (FIT, gt1, fit) + call geo_1graph$t (gd, gt1, fit, gfit, xref, yref, xin, yin, wts, + npts) + if (GG_CONSTXY(gfit) == YES) + call geo_conxy$t (gd, fit, sx1, sy1, sx2, sy2) + call printf ("%s %s\n") + call pargstr (xerrmsg) + call pargstr (yerrmsg) + + # Read the cursor commands. + call amov$t (wts, Mem$t[w], npts) + while (clgcur ("cursor", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != EOF) { + + switch (key) { + + case 'q': + call amov$t (Mem$t[w], wts, npts) + break + + case '?': + if (GM_PROJECTION(fit) == GM_NONE) + call gpagefile (gd, GHELPFILE, "") + else + call gpagefile (gd, CHELPFILE, "") + + case ':': + call geo_colon (gd, fit, gfit, Memc[cmd], newgraph) + switch (GG_PLOTTYPE(gfit)) { + case FIT: + call gt_colon (Memc[cmd], gd, gt1, newgraph) + case XXRESID: + call gt_colon (Memc[cmd], gd, gt2, newgraph) + case XYRESID: + call gt_colon (Memc[cmd], gd, gt3, newgraph) + case YXRESID: + call gt_colon (Memc[cmd], gd, gt4, newgraph) + case YYRESID: + call gt_colon (Memc[cmd], gd, gt5, newgraph) + } + + case 'l': + if (GG_FITERROR(gfit) == NO) { + call geo_lcoeff$t (sx1, sy1, xshift, yshift, xscale, yscale, + thetax, thetay) + call printf ("xshift: %.2f yshift: %.2f ") + call parg$t (xshift) + call parg$t (yshift) + call printf ("xmag: %0.3g ymag: %0.3g ") + call parg$t (xscale) + call parg$t (yscale) + call printf ("xrot: %.2f yrot: %.2f\n") + call parg$t (thetax) + call parg$t (thetay) + } + + case 't': + if (GG_FITERROR(gfit) == NO && GG_PLOTTYPE(gfit) == FIT) + call geo_lxy$t (gd, fit, sx1, sy1, sx2, sy2, xref, yref, + xin, yin, npts, wx, wy) + + case 'c': + if (GG_CONSTXY(gfit) == YES) + GG_CONSTXY(gfit) = NO + else if (GG_CONSTXY(gfit) == NO) + GG_CONSTXY(gfit) = YES + + case 'd', 'u': + if (key == 'd') + delete = YES + else + delete = NO + + switch (GG_PLOTTYPE(gfit)) { + case FIT: + call geo_1delete$t (gd, xin, yin, Mem$t[w], wts, npts, wx, + wy, delete) + case XXRESID: + call geo_2delete$t (gd, xref, Mem$t[xresid], Mem$t[w], wts, + npts, wx, wy, delete) + case XYRESID: + call geo_2delete$t (gd, yref, Mem$t[xresid], Mem$t[w], wts, + npts, wx, wy, delete) + case YXRESID: + call geo_2delete$t (gd, xref, Mem$t[yresid], Mem$t[w], wts, + npts, wx, wy, delete) + case YYRESID: + call geo_2delete$t (gd, yref, Mem$t[yresid], Mem$t[w], wts, + npts, wx, wy, delete) + } + + GG_NEWFUNCTION(gfit) = YES + + case 'g': + if (GG_PLOTTYPE(gfit) != FIT) + newgraph = YES + GG_PLOTTYPE(gfit) = FIT + + case 'x': + if (GG_PLOTTYPE(gfit) != XXRESID) + newgraph = YES + GG_PLOTTYPE(gfit) = XXRESID + + case 'r': + if (GG_PLOTTYPE(gfit) != XYRESID) + newgraph = YES + GG_PLOTTYPE(gfit) = XYRESID + + case 'y': + if (GG_PLOTTYPE(gfit) != YXRESID) + newgraph = YES + GG_PLOTTYPE(gfit) = YXRESID + + case 's': + if (GG_PLOTTYPE(gfit) != YYRESID) + newgraph = YES + GG_PLOTTYPE(gfit) = YYRESID + + case 'f': + # do fit + if (GG_NEWFUNCTION(gfit) == YES) { + iferr { + switch (GM_FIT(fit)) { + case GM_ROTATE: + call geo_ftheta$t (fit, sx1, sy1, xref, yref, xin, + yin, Mem$t[w], Mem$t[xresid], Mem$t[yresid], + npts, xerrmsg, maxch, yerrmsg, maxch) + sx2 = NULL + sy2 = NULL + case GM_RSCALE: + call geo_fmagnify$t (fit, sx1, sy1, xref, yref, xin, + yin, Mem$t[w], Mem$t[xresid], Mem$t[yresid], + npts, xerrmsg, maxch, yerrmsg, maxch) + sx2 = NULL + sy2 = NULL + case GM_RXYSCALE: + call geo_flinear$t (fit, sx1, sy1, xref, yref, xin, + yin, Mem$t[w], Mem$t[xresid], Mem$t[yresid], + npts, xerrmsg, maxch, yerrmsg, maxch) + sx2 = NULL + sy2 = NULL + default: + call geo_fxy$t (fit, sx1, sx2, xref, yref, xin, + Mem$t[w], Mem$t[xresid], npts, YES, + xerrmsg, maxch) + call geo_fxy$t (fit, sy1, sy2, xref, yref, yin, + Mem$t[w], Mem$t[yresid], npts, NO, + yerrmsg, maxch) + } + if (GM_MAXITER(fit) <= 0 || IS_INDEFD(GM_REJECT(fit))) + GM_NREJECT(fit) = 0 + else + call geo_mreject$t (fit, sx1, sy1, sx2, sy2, xref, + yref, xin, yin, Mem$t[w], Mem$t[xresid], + Mem$t[yresid], npts, xerrmsg, maxch, + yerrmsg, maxch) + GG_NEWFUNCTION(gfit) = NO + GG_FITERROR(gfit) = NO + errcode = OK + } then { + errcode = errget (errstr, SZ_LINE) + call printf ("%s\n") + call pargstr (errstr) + GG_FITERROR(gfit) = YES + } + } + + # plot new graph + if (GG_FITERROR(gfit) == YES) + newgraph = NO + else + newgraph = YES + + case 'o': + GG_OVERPLOT(gfit) = YES + + default: + call printf ("\07") + + } + + if (newgraph == YES) { + switch (GG_PLOTTYPE(gfit)) { + case FIT: + call geo_label (FIT, gt1, fit) + call geo_1graph$t (gd, gt1, fit, gfit, xref, yref, xin, yin, + Mem$t[w], npts) + if (GG_CONSTXY(gfit) == YES) + call geo_conxy$t (gd, fit, sx1, sy1, sx2, sy2) + case XXRESID: + call geo_label (XXRESID, gt2, fit) + call geo_2graph$t (gd, gt2, fit, gfit, xref, Mem$t[xresid], + Mem$t[w], npts) + case XYRESID: + call geo_label (XYRESID, gt3, fit) + call geo_2graph$t (gd, gt3, fit, gfit, yref, Mem$t[xresid], + Mem$t[w], npts) + case YXRESID: + call geo_label (YXRESID, gt4, fit) + call geo_2graph$t (gd, gt4, fit, gfit, xref, Mem$t[yresid], + Mem$t[w], npts) + case YYRESID: + call geo_label (YYRESID, gt5, fit) + call geo_2graph$t (gd, gt5, fit, gfit, yref, Mem$t[yresid], + Mem$t[w], npts) + } + call printf ("%s %s\n") + call pargstr (xerrmsg) + call pargstr (yerrmsg) + newgraph = NO + } + } + + # Free space. + call gt_free (gt1) + call gt_free (gt2) + call gt_free (gt3) + call gt_free (gt4) + call gt_free (gt5) + call sfree (sp) + + # Call an error if appropriate. + if (errcode > 0) + call error (2, errstr) +end + +# GEO_LCOEFF -- Print the coefficents of the linear portion of the +# fit, xshift, yshift, xexpansion, yexpansion, x and y rotations. + +procedure geo_lcoeff$t (sx, sy, xshift, yshift, xscale, yscale, xrot, yrot) + +pointer sx #I pointer to the x surface fit +pointer sy #I pointer to the y surface fit +PIXEL xshift #O output x shift +PIXEL yshift #O output y shift +PIXEL xscale #O output x scale +PIXEL yscale #O output y scale +PIXEL xrot #O rotation of point on x axis +PIXEL yrot #O rotation of point on y axis + +int nxxcoeff, nxycoeff, nyxcoeff, nyycoeff +pointer sp, xcoeff, ycoeff +PIXEL xxrange, xyrange, xxmaxmin, xymaxmin +PIXEL yxrange, yyrange, yxmaxmin, yymaxmin +PIXEL a, b, c, d + +bool fp_equal$t() +$if (datatype == r) +int gsgeti() +real gsgetr() +$else +int dgsgeti() +double dgsgetd() +$endif + +begin + # Allocate working space. + call smark (sp) +$if (datatype == r) + call salloc (xcoeff, gsgeti (sx, GSNCOEFF), TY_PIXEL) + call salloc (ycoeff, gsgeti (sy, GSNCOEFF), TY_PIXEL) +$else + call salloc (xcoeff, dgsgeti (sx, GSNCOEFF), TY_PIXEL) + call salloc (ycoeff, dgsgeti (sy, GSNCOEFF), TY_PIXEL) +$endif + + # Get coefficients and numbers of coefficients. +$if (datatype == r) + call gscoeff (sx, Mem$t[xcoeff], nxxcoeff) + call gscoeff (sy, Mem$t[ycoeff], nyycoeff) + nxxcoeff = gsgeti (sx, GSNXCOEFF) + nxycoeff = gsgeti (sx, GSNYCOEFF) + nyxcoeff = gsgeti (sy, GSNXCOEFF) + nyycoeff = gsgeti (sy, GSNYCOEFF) +$else + call dgscoeff (sx, Mem$t[xcoeff], nxxcoeff) + call dgscoeff (sy, Mem$t[ycoeff], nyycoeff) + nxxcoeff = dgsgeti (sx, GSNXCOEFF) + nxycoeff = dgsgeti (sx, GSNYCOEFF) + nyxcoeff = dgsgeti (sy, GSNXCOEFF) + nyycoeff = dgsgeti (sy, GSNYCOEFF) +$endif + + # Get the data range. +$if (datatype == r) + if (gsgeti (sx, GSTYPE) != GS_POLYNOMIAL) { + xxrange = (gsgetr (sx, GSXMAX) - gsgetr (sx, GSXMIN)) / 2.0 + xxmaxmin = - (gsgetr (sx, GSXMAX) + gsgetr (sx, GSXMIN)) / 2.0 + xyrange = (gsgetr (sx, GSYMAX) - gsgetr (sx, GSYMIN)) / 2.0 + xymaxmin = - (gsgetr (sx, GSYMAX) + gsgetr (sx, GSYMIN)) / 2.0 +$else + if (dgsgeti (sx, GSTYPE) != GS_POLYNOMIAL) { + xxrange = (dgsgetd (sx, GSXMAX) - dgsgetd (sx, GSXMIN)) / 2.0d0 + xxmaxmin = - (dgsgetd (sx, GSXMAX) + dgsgetd (sx, GSXMIN)) / 2.0d0 + xyrange = (dgsgetd (sx, GSYMAX) - dgsgetd (sx, GSYMIN)) / 2.0d0 + xymaxmin = - (dgsgetd (sx, GSYMAX) + dgsgetd (sx, GSYMIN)) / 2.0d0 +$endif + } else { + xxrange = PIXEL(1.0) + xxmaxmin = PIXEL(0.0) + xyrange = PIXEL(1.0) + xymaxmin = PIXEL(0.0) + } + +$if (datatype == r) + if (gsgeti (sy, GSTYPE) != GS_POLYNOMIAL) { + yxrange = (gsgetr (sy, GSXMAX) - gsgetr (sy, GSXMIN)) / 2.0 + yxmaxmin = - (gsgetr (sy, GSXMAX) + gsgetr (sy, GSXMIN)) / 2.0 + yyrange = (gsgetr (sy, GSYMAX) - gsgetr (sy, GSYMIN)) / 2.0 + yymaxmin = - (gsgetr (sy, GSYMAX) + gsgetr (sy, GSYMIN)) / 2.0 +$else + if (dgsgeti (sy, GSTYPE) != GS_POLYNOMIAL) { + yxrange = (dgsgetd (sy, GSXMAX) - dgsgetd (sy, GSXMIN)) / 2.0d0 + yxmaxmin = - (dgsgetd (sy, GSXMAX) + dgsgetd (sy, GSXMIN)) / 2.0d0 + yyrange = (dgsgetd (sy, GSYMAX) - dgsgetd (sy, GSYMIN)) / 2.0d0 + yymaxmin = - (dgsgetd (sy, GSYMAX) + dgsgetd (sy, GSYMIN)) / 2.0d0 +$endif + } else { + yxrange = PIXEL(1.0) + yxmaxmin = PIXEL(0.0) + yyrange = PIXEL(1.0) + yymaxmin = PIXEL(0.0) + } + + # Get the shifts. + xshift = Mem$t[xcoeff] + Mem$t[xcoeff+1] * xxmaxmin / xxrange + + Mem$t[xcoeff+2] * xymaxmin / xyrange + yshift = Mem$t[ycoeff] + Mem$t[ycoeff+1] * yxmaxmin / yxrange + + Mem$t[ycoeff+2] * yymaxmin / yyrange + + # Get the rotation and scaling parameters and correct for normalization. + if (nxxcoeff > 1) + a = Mem$t[xcoeff+1] / xxrange + else + a = PIXEL(0.0) + if (nxycoeff > 1) + b = Mem$t[xcoeff+nxxcoeff] / xyrange + else + b = PIXEL(0.0) + if (nyxcoeff > 1) + c = Mem$t[ycoeff+1] / yxrange + else + c = PIXEL(0.0) + if (nyycoeff > 1) + d = Mem$t[ycoeff+nyxcoeff] / yyrange + else + d = PIXEL(0.0) + + # Get the magnification factors. + xscale = sqrt (a * a + c * c) + yscale = sqrt (b * b + d * d) + + # Get the x and y axes 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) + + call sfree (sp) +end + +$endfor |