aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/dtoi/hdicfit/hdicgfit.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/imred/dtoi/hdicfit/hdicgfit.x')
-rw-r--r--noao/imred/dtoi/hdicfit/hdicgfit.x402
1 files changed, 402 insertions, 0 deletions
diff --git a/noao/imred/dtoi/hdicfit/hdicgfit.x b/noao/imred/dtoi/hdicfit/hdicgfit.x
new file mode 100644
index 00000000..b50ed03c
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgfit.x
@@ -0,0 +1,402 @@
+include <error.h>
+include <pkg/gtools.h>
+include <mach.h>
+include <gset.h>
+include "hdicfit.h"
+
+define HELP "noao$lib/scr/hdicgfit.key"
+define PROMPT "hdicfit options"
+define EB_WTS 10
+define EB_SDEV 11
+
+# ICG_FIT -- Interactive curve fitting with graphics. This is the main
+# entry point for the interactive graphics part of the icfit package.
+
+procedure icg_fitd (ic, gp, cursor, gt, cv, density, loge, owts, osdev, npts)
+
+pointer ic # ICFIT pointer
+pointer gp # GIO pointer
+char cursor[ARB] # GIO cursor input
+pointer gt # GTOOLS pointer
+pointer cv # CURFIT pointer
+double density[npts] # Original density, not fog subtracted
+double loge[npts] # Original log exposure values
+double owts[npts] # Original weights array
+double osdev[npts] # Original standard deviation array
+int npts # Number of points
+
+real wx, wy, wwy, oldx, oldy
+int wcs, key, nptso
+char cmd[SZ_LINE]
+
+int i, newgraph, axes[2], linetype
+double x1, newwt
+pointer userwts, x, y, wts, den, whydel, sdev, ebw
+
+int clgcur(), stridxs(), scan(), nscan()
+int icg_nearestd(), gstati()
+double dcveval()
+errchk ic_fitd, malloc
+define redraw_ 91
+
+begin
+ # Allocate memory for the fit and a copy of the weights.
+ # The weights are copied because they are changed when points are
+ # deleted. The input x is the untransformed density, and is used
+ # to generate other types of transforms. Points can be added to
+ # the sample, so the y array and weights array can change as well.
+ # The original number of points is also remembered.
+
+ call malloc (userwts, npts, TY_DOUBLE)
+ call malloc (x, npts, TY_DOUBLE)
+ call malloc (y, npts, TY_DOUBLE)
+ call malloc (wts, npts, TY_DOUBLE)
+ call malloc (den, npts, TY_DOUBLE)
+ call malloc (whydel, npts, TY_INT)
+ call malloc (sdev, npts, TY_DOUBLE)
+ call malloc (ebw, npts, TY_DOUBLE)
+
+ call amovd (owts, Memd[userwts], npts)
+ call amovd (owts, Memd[wts], npts)
+ call amovd (loge, Memd[y], npts)
+ call amovd (density, Memd[den], npts)
+ call amovki (NDELETE, Memi[whydel], npts)
+ call amovd (osdev, Memd[sdev], npts)
+ nptso = npts
+
+ # Initialize
+ IC_OVERPLOT(ic) = NO
+ IC_NEWX(ic) = YES
+ IC_NEWY(ic) = YES
+ IC_NEWWTS(ic) = YES
+ IC_NEWFUNCTION(ic) = YES
+ IC_NEWTRANSFORM(ic) = YES
+ IC_UPDATE(ic) = YES
+ IC_EBARS(ic) = EB_SDEV
+
+ # Read cursor commands.
+
+ key = 'f'
+ axes[1] = IC_AXES(ic, IC_GKEY(ic), 1)
+ axes[2] = IC_AXES(ic, IC_GKEY(ic), 2)
+
+ repeat {
+ switch (key) {
+ case '?': # Print help text.
+ call gpagefile (gp, HELP, PROMPT)
+
+ case 'q': # Terminate cursor loop
+ break
+
+ case ':': # List or set parameters
+ if (stridxs ("/", cmd) == 1)
+ call gt_colon (cmd, gp, gt, newgraph)
+ else
+ call icg_colond (ic, cmd, gp, gt, cv, Memd[x],
+ Memd[y], Memd[wts], npts)
+
+ if (IC_RESET(ic) == YES) {
+ npts = nptso
+ call amovd (owts, Memd[userwts], npts)
+ call amovd (owts, Memd[wts], npts)
+ call amovd (loge, Memd[y], npts)
+ call amovd (density, Memd[den], npts)
+ call amovki (NDELETE, Memi[whydel], npts)
+ call amovd (osdev, Memd[sdev], npts)
+ call hdic_init (density, npts, 0.0)
+ }
+
+ # See if user wants to quit without updating
+ if (IC_UPDATE(ic) == NO) {
+ call mfree (x, TY_DOUBLE)
+ call mfree (y, TY_DOUBLE)
+ call mfree (wts, TY_DOUBLE)
+ call mfree (userwts, TY_DOUBLE)
+ call mfree (den, TY_DOUBLE)
+ call mfree (sdev, TY_DOUBLE)
+ return
+ }
+
+ case 'a': # Add data points to the sample. This is only possible
+ # from an HD curve plot.
+
+ if ((IC_AXES (ic, IC_GKEY(ic), 1) == 'y') &&
+ (IC_AXES (ic, IC_GKEY(ic), 2) == 'u')) {
+
+ # Query for weight after plotting current location
+ # call gt_plot (gp, gt, wx, wy, 1)
+ call gmark (gp, wx, wy, GM_CIRCLE, 2.0, 2.0)
+ newwt = 1.0D0
+ call printf ("Enter weight of new point (%g): ")
+ call pargd (newwt)
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call pargd (x1)
+ if (nscan() == 1) {
+ if (!IS_INDEFD (x1)) {
+ newwt = x1
+ }
+ }
+ }
+
+ } else {
+ call eprintf ("Points can be added only from an HD Curve\n")
+ next
+ }
+
+ # Add fog into "density above fog" value read from cursor
+ wwy = wy + IC_FOG (ic)
+ if (wwy < 0.0) {
+ call eprintf (
+ "New density (%g) is below fog and will not be added\n")
+ call pargr (wwy)
+ next
+ }
+
+ # Add point into sample
+ call eprintf ("New Point: density above fog = %.4f, log ")
+ call pargr (wwy)
+ call eprintf ("exposure = %.4f, weight = %.4f\n")
+ call pargr (wx)
+ call pargd (newwt)
+
+ call hdic_addpoint (ic, wwy, wx, newwt, den, y, wts, userwts,
+ x, whydel, sdev, npts)
+
+ call realloc (ebw, npts, TY_DOUBLE)
+
+ call hdic_transform (ic, Memd[den], Memd[userwts], Memd[x],
+ Memd[wts], Memi[whydel], npts)
+
+ case 'c': # Print the positions of data points.
+ i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts,
+ wx, wy)
+
+ if (i != 0) {
+ call printf ("den= %7.4g x= %7.4g exp= %7.4g fit= %7.4g")
+ call pargd (Memd[den+i-1])
+ call pargd (Memd[x+i-1])
+ call pargd (Memd[y+i-1])
+ call pargd (dcveval (cv, Memd[x+i-1]))
+ }
+
+ case 'd': # Delete data points.
+ call icg_deleted (ic, gp, gt, cv, Memd[x], Memd[y], Memd[wts],
+ npts, wx, wy)
+
+ case 'f': # Fit the function and reset the flags.
+ iferr {
+ # Copy new transformed vector, if necessary
+ if (IC_NEWTRANSFORM(ic) == YES || IC_NEWFOG(ic) == YES)
+ call hdic_transform (ic, Memd[den], Memd[userwts],
+ Memd[x], Memd[wts], Memi[whydel], npts)
+
+ call ic_fitd (ic, cv, Memd[x], Memd[y], Memd[wts], npts,
+ IC_NEWX(ic), IC_NEWY(ic), IC_NEWWTS(ic),
+ IC_NEWFUNCTION(ic))
+
+ IC_NEWX(ic) = NO
+ IC_NEWY(ic) = NO
+ IC_NEWWTS(ic) = NO
+ IC_NEWFUNCTION(ic) = NO
+ IC_NEWTRANSFORM(ic) = NO
+ IC_FITERROR(ic) = NO
+ IC_NEWFOG(ic) = NO
+ newgraph = YES
+ } then {
+ IC_FITERROR(ic) = YES
+ call erract (EA_WARN)
+ newgraph = NO
+ }
+
+ case 'g': # Set graph axes types.
+ call printf ("Graph key to be defined: ")
+ call flush (STDOUT)
+ if (scan() == EOF)
+ goto redraw_
+ call gargc (cmd[1])
+
+ switch (cmd[1]) {
+ case '\n':
+ case 'h', 'i', 'j', 'k', 'l':
+ switch (cmd[1]) {
+ case 'h':
+ key = 1
+ case 'i':
+ key = 2
+ case 'j':
+ key = 3
+ case 'k':
+ key = 4
+ case 'l':
+ key = 5
+ }
+
+ call printf ("Set graph axes types (%c, %c): ")
+ call pargc (IC_AXES(ic, key, 1))
+ call pargc (IC_AXES(ic, key, 2))
+ call flush (STDOUT)
+ if (scan() == EOF)
+ goto redraw_
+ call gargc (cmd[1])
+
+ switch (cmd[1]) {
+ case '\n':
+ default:
+ call gargc (cmd[2])
+ call gargc (cmd[2])
+ if (cmd[2] != '\n') {
+ IC_AXES(ic, key, 1) = cmd[1]
+ IC_AXES(ic, key, 2) = cmd[2]
+ }
+ }
+ default:
+ call printf ("Not a graph key")
+ }
+
+ case 'h':
+ if (IC_GKEY(ic) != 1) {
+ IC_GKEY(ic) = 1
+ newgraph = YES
+ }
+
+ case 'i':
+ if (IC_GKEY(ic) != 2) {
+ IC_GKEY(ic) = 2
+ newgraph = YES
+ }
+
+ case 'j':
+ if (IC_GKEY(ic) != 3) {
+ IC_GKEY(ic) = 3
+ newgraph = YES
+ }
+
+ case 'k':
+ if (IC_GKEY(ic) != 4) {
+ IC_GKEY(ic) = 4
+ newgraph = YES
+ }
+
+ case 'l':
+ if (IC_GKEY(ic) != 5) {
+ IC_GKEY(ic) = 5
+ newgraph = YES
+ }
+
+ case 'o': # Set overplot flag
+ IC_OVERPLOT(ic) = YES
+
+ case 'r': # Redraw the graph
+ newgraph = YES
+
+ case 'u': # Undelete data points.
+ call icg_undeleted (ic, gp, gt, cv, Memd[x], Memd[y],
+ Memd[wts], Memd[userwts], npts, wx, wy)
+
+ case 'w': # Window graph
+ call gt_window (gt, gp, cursor, newgraph)
+
+ case 'x': # Reset the value of the x point.
+ i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx,
+ wy)
+
+ if (i != 0) {
+ call printf ("Enter new x (%g): ")
+ call pargd (Memd[x+i-1])
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (x1)
+ if (nscan() == 1) {
+ if (!IS_INDEF (x1)) {
+ oldx = Memd[x+i-1]
+ oldy = Memd[y+i-1]
+ Memd[x+i-1] = x1
+ call hd_redraw (gp, oldx, oldy, x1, oldy)
+ IC_NEWX(ic) = YES
+ }
+ }
+ }
+ }
+
+ case 'y': # Reset the value of the y point.
+ i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx,
+ wy)
+
+ if (i != 0) {
+ call printf ("Enter new y (%g): ")
+ call pargd (Memd[y+i-1])
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (x1)
+ if (nscan() == 1) {
+ if (!IS_INDEF (x1)) {
+ oldx = Memd[x+i-1]
+ oldy = Memd[y+i-1]
+ Memd[y+i-1] = x1
+ call hd_redraw (gp, oldx, oldy, oldx, x1)
+ IC_NEWY(ic) = YES
+ }
+ }
+ }
+ }
+
+ case 'z': # Reset the weight value of the nearest point
+ i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx,
+ wy)
+
+ if (i != 0) {
+ call printf ("Enter new weight (%g): ")
+ call pargd (Memd[wts+i-1])
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (x1)
+ if (nscan() == 1) {
+ if (!IS_INDEF (x1)) {
+ Memd[wts+i-1] = x1
+ IC_NEWWTS(ic) = YES
+ }
+ }
+ }
+ }
+
+ default: # Let the user decide on any other keys.
+ call icg_user (ic, gp, gt, cv, wx, wy, wcs, key, cmd)
+ }
+
+ # Redraw the graph if necessary.
+redraw_ if (newgraph == YES) {
+ if (IC_AXES(ic, IC_GKEY(ic), 1) != axes[1]) {
+ axes[1] = IC_AXES(ic, IC_GKEY(ic), 1)
+ call gt_setr (gt, GTXMIN, INDEF)
+ call gt_setr (gt, GTXMAX, INDEF)
+ }
+ if (IC_AXES(ic, IC_GKEY(ic), 2) != axes[2]) {
+ axes[2] = IC_AXES(ic, IC_GKEY(ic), 2)
+ call gt_setr (gt, GTYMIN, INDEF)
+ call gt_setr (gt, GTYMAX, INDEF)
+ }
+
+ # Overplot with a different line type
+ if (IC_OVERPLOT(ic) == YES)
+ linetype = min ((gstati (gp, G_PLTYPE) + 1), 4)
+ else
+ linetype = GL_SOLID
+ call gseti (gp, G_PLTYPE, linetype)
+
+ call hdic_ebw (ic, Memd[den], Memd[x], Memd[sdev], Memd[ebw],
+ npts)
+
+ call icg_graphd (ic, gp, gt, cv, Memd[x], Memd[y], Memd[wts],
+ Memd[ebw], npts)
+
+ newgraph = NO
+ }
+ } until (clgcur (cursor, wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+
+ call mfree (x, TY_DOUBLE)
+ call mfree (y, TY_DOUBLE)
+ call mfree (wts, TY_DOUBLE)
+ call mfree (userwts, TY_DOUBLE)
+ call mfree (den, TY_DOUBLE)
+end
generated by cgit (git 2.34.1) at 2025-09-20 12:08:16 -0400