include include include include # for SZ_TIME include include define HELPFILE "focgeom$rmarkerx.key" define LEN_RMSTRUCT 20 define RM_X1 Memr[$1] define RM_X2 Memr[$1+1] define RM_Y1 Memr[$1+2] define RM_Y2 Memr[$1+3] define RM_HBOXX Memi[$1+4] define RM_HBOXY Memi[$1+5] define RM_BOUNDARY MAX_REAL # value for out of bounds pixels define RM_BOUNDWIDTH 50 # width of allowed boundary extension # RMARKERX -- obtain approximate reseau locations interactively # # D. Giaretta STScI Mar1988 Initial code # Phil Hodge, 9-Dec-1993 Delete rm_markpos; include option to delete reseau; # allow indef positions; change parameter names. procedure t_rmarkerx() char inres[SZ_FNAME] # input reseau table with initial locations char outres[SZ_FNAME] # output reseau table for better locations char entry[SZ_R_PATTERN] # entry name for template pos char newentry[SZ_R_ENTNAME] # new entry name for better positions char image[SZ_FNAME] # input image name, if any char device[SZ_FNAME] # graphics device char logfile[SZ_FNAME] # output log file, if any #-- pointer sp pointer im pointer rpin, rpout, rentpin, rentpout pointer gd, gopen() pointer res pointer i_foc, x_foc_immap() char datetime[SZ_TIME] # date and time for log file long old_time, new_time, clktime() bool inplace int fd, open() int nowhite() int rs_readr(), rs_list() pointer rlist, rs_template() begin call smark (sp) call clgstr ("inres", inres, SZ_FNAME) if (nowhite (inres, inres, SZ_FNAME) < 1) call error (1, "input reseau table must be specified") call clgstr ("entry" , entry, SZ_R_PATTERN) call clgstr ("outres", outres, SZ_FNAME) if (nowhite (outres, outres, SZ_FNAME) < 1) call error (1, "output reseau table must be specified") call clgstr ("newentry", newentry, SZ_R_ENTNAME) call clgstr ("image", image, SZ_FNAME) call clgstr ("device", device, SZ_FNAME) gd = gopen (device, NEW_FILE, STDGRAPH) call clgstr ("logfile", logfile, SZ_FNAME) if (nowhite (logfile, logfile, SZ_FNAME) > 0) fd = open (logfile, APPEND, TEXT_FILE) else fd = NULL # initialise structure call rm_init (res) # Open the input and output reseau tables. call rm_open (inres, outres, rpin, rpout, inplace) call rs_alloc (rpout, TY_REAL, rentpout) call rs_alloc (rpin, TY_REAL, rentpin) # Get the input positions which the user may interactively modify. rlist = rs_template (rpin, entry) if (rs_readr (rpin, rentpin, EOS, rs_list(rlist)) == RES_F_NONEXTENTRY) call error (1, "unable to read entry positions") # If newentry was not specified, copy the entry name from inres. if (nowhite (newentry, newentry, SZ_R_ENTNAME) == 0) call strcpy (RES_ENT_ENTRY(rentpin), newentry, SZ_R_ENTNAME) # Copy reseau entry info. call strcpy (newentry, RES_ENT_ENTRY(rentpout), SZ_R_ENTRY) RES_ENT_NROWS(rentpout) = RES_ENT_NROWS(rentpin) RES_ENT_NCOLS(rentpout) = RES_ENT_NCOLS(rentpin) RES_ENT_SFIELD1(rentpout) = RES_ENT_SFIELD1(rentpin) RES_ENT_SFIELD2(rentpout) = RES_ENT_SFIELD2(rentpin) RES_ENT_SFIELD3(rentpout) = RES_ENT_SFIELD3(rentpin) RES_ENT_SFIELD4(rentpout) = RES_ENT_SFIELD4(rentpin) call strcpy (RES_ENT_TRACKING(rentpin), RES_ENT_TRACKING(rentpout), SZ_R_TRACKING) call rs_t_update (rentpout, "M", SZ_R_TRACKING) # append "M" # special treatment for FOC coords call salloc (i_foc, SZ_FOC_COORD, TY_STRUCT) FOC_COORD_INIT(i_foc) = false if (nowhite (image, image, SZ_FNAME) == 0) { im = NULL i_foc = NULL } else { im = x_foc_immap (image, i_foc, READ_ONLY, NULL) call imsetr( im, IM_BNDRYPIXVAL, RM_BOUNDARY) call imseti( im, IM_NBNDRYPIX, RM_BOUNDWIDTH) call imseti( im, IM_TYBNDRY, BT_CONSTANT) } # Log file and entry names. if (fd != NULL) { old_time = 0 new_time = clktime (old_time) call cnvtime (new_time, datetime, SZ_TIME) call fprintf (fd, "# %s\n") call pargstr (datetime) call fprintf (fd, "Input reseau table %s, entry %s\n") call pargstr (inres) call pargstr (RES_ENT_ENTRY(rentpin)) call fprintf (fd, "Output reseau table %s, entry %s\n") call pargstr (outres) call pargstr (RES_ENT_ENTRY(rentpout)) if (im != NULL) { call fprintf (fd, "Input image %s\n") call pargstr (image) } } # Interactively move the input reseau points, writing to output. call rm_inter (res, i_foc, RES_ENT_ENTRY(rentpin), gd, fd, RES_ENT_NCOLS(rentpin), RES_ENT_NROWS(rentpin), Memr[RES_XPT(rentpin)], Memr[RES_YPT(rentpin)], Memr[RES_XPT(rentpout)], Memr[RES_YPT(rentpout)]) call rs_write (rpout, rentpout, 0) if (im != NULL) call imunmap (im) # close the rest of the files if (fd != NULL) call close (fd) call gclose (gd) if (!inplace) call rs_close (rpout) call rs_close (rpin) call rm_free (res) call sfree (sp) end # RM_INIT -- initialise the rmarkerx structure procedure rm_init (res) pointer res # o: res structure #-- real clgetr() int clgeti() begin call malloc (res, LEN_RMSTRUCT, TY_STRUCT) RM_X1(res) = clgetr ("x1") RM_X2(res) = clgetr ("x2") RM_Y1(res) = clgetr ("y1") RM_Y2(res) = clgetr ("y2") RM_HBOXX(res) = min ((clgeti("xsearch") + 1) / 2, RM_BOUNDWIDTH) RM_HBOXY(res) = min ((clgeti("ysearch") + 1) / 2, RM_BOUNDWIDTH) end # RM_FREE -- free space from rm structure procedure rm_free (res) pointer res begin call mfree (res, TY_STRUCT) end procedure rm_inter (res, i_foc, entry, gd, fd, nx, ny, tx, ty, ax, ay) pointer res # i: rm structure pointer i_foc # i: input FOC image structure char entry[ARB] # i: input entry name pointer gd # i: GIO pointer int fd # i: fd for log file, or NULL int nx, ny # i: number of template marks in X and Y real tx[ARB] # i: x coords of template marks real ty[ARB] # i: y coords of template marks real ax[ARB] # o: better approx x real ay[ARB] # o: better approx y #-- int nres int i, j, key, wcs, nfit int junk, tipt[4], ipoint, ix, iy, minpos, xside, yside real wx, wy, xa[4], ya[4], sampoff, lineoff real oldx, oldy # for passing to rm_log pointer sp pointer strval pointer nax, nay, imgs2r(), buf, sax, say int nindex, scan(), nscan() int clgcur(), rm_gcur(), rm_mpoint() begin nres = nx * ny call smark (sp) call salloc (strval, SZ_LINE, TY_CHAR) call salloc (nax, nres, TY_REAL) call salloc (nay, nres, TY_REAL) call salloc (sax, nres, TY_REAL) call salloc (say, nres, TY_REAL) nfit = 0 if (i_foc == NULL) { sampoff = 0 lineoff = 0 } else { sampoff = FOC_SAMPBEG(i_foc) + FOC_SOFF1(i_foc) lineoff = FOC_LINEBEG(i_foc) + FOC_SOFF2(i_foc) } call amovr (tx, ax, nres) call amovr (ty, ay, nres) call amovr (tx, Memr[nax], nres) call amovr (ty, Memr[nay], nres) call amovr (tx, Memr[sax], nres) call amovr (ty, Memr[say], nres) # display image data positions on plot # deleted 12/9/93 by PEH # Draw the reseau positions. call sprintf (Memc[strval], SZ_LINE, "%s %d x %d") call pargstr (entry) call pargi (nx) call pargi (ny) call gclear (gd) call gswind (gd, RM_X1(res), RM_X2(res), RM_Y1(res), RM_Y2(res)) call glabax (gd, Memc[strval], "axis1", "axis2") # plot the template positions call gseti (gd, G_PMLTYPE, GL_SOLID) call gpmark (gd, ax, ay, nres, GM_PLUS, 2.0, 2.0) while (clgcur ("cursor", wx, wy, wcs, key, Memc[strval], SZ_LINE) != EOF) { switch ( key ) { case ':' : # process colon commands call rm_colon (gd, res, key, Memc[strval]) case 'q' : break case '?': call gpagefile (gd, HELPFILE, "") # Assign a position (e.g. replace an INDEF position). case 'a': call printf (" give column and row indices: ") call flush (STDOUT) if (scan() == EOF) { call printf (" no action taken\n") next } call gargi (i) call gargi (j) nindex = nscan() if (nindex < 1) { call printf (" no action taken\n") next } else if (nindex < 2) { j = 1 } ipoint = (j - 1) * nx + i oldx = ax[ipoint] oldy = ay[ipoint] call printf (" where should this be?\n") junk = rm_gcur ("cursor", ax[ipoint], ay[ipoint], wcs, key, Memc[strval], SZ_LINE) call gmark (gd, ax[ipoint], ay[ipoint], GM_BOX, 1.0, 1.0) call printf (" \n") Memr[nax+ipoint-1] = ax[ipoint] Memr[nay+ipoint-1] = ay[ipoint] call rm_log (fd, 'a', ipoint, nx, oldx, oldy, ax[ipoint], ay[ipoint]) # back to the original template case 'b': call amovr (tx, Memr[nax], nres) call amovr (ty, Memr[nay], nres) call rm_show (gd, ax, ay, Memr[nax], Memr[nay], Memr[sax], Memr[say], xa, ya, tipt, nfit, nres) call amovr (tx, Memr[sax], nres) call amovr (ty, Memr[say], nres) if (fd != NULL) call fprintf (fd, "Back to the original input positions.\n") # Delete the nearest reseau mark. case 'd': ipoint = rm_mpoint (wx, wy, ax, ay, nres) oldx = ax[ipoint] oldy = ay[ipoint] call gseti (gd, G_PMLTYPE, GL_CLEAR) call gmark (gd, ax[ipoint], ay[ipoint], GM_PLUS, 2.0, 2.0) call gseti (gd, G_PMLTYPE, GL_SOLID) call gmark (gd, ax[ipoint], ay[ipoint], GM_CROSS, 2.0, 2.0) ax[ipoint] = INDEFR ay[ipoint] = INDEFR Memr[nax+ipoint-1] = INDEFR Memr[nay+ipoint-1] = INDEFR i = mod (ipoint, nx) if (i < 1) i = nx j = (ipoint - i) / nx + 1 call printf (" [%d,%d] --> (%.1f,%.1f) deleted\n") call pargi (i) call pargi (j) call pargr (oldx) call pargr (oldy) call rm_log (fd, 'd', ipoint, nx, oldx, oldy, ax[ipoint], ay[ipoint]) # Lock on to minimum case 'l': if (i_foc == NULL) { call printf (" no image was specified\n") } else { xside = 2*RM_HBOXX(res) + 1 yside = 2*RM_HBOXY(res) + 1 call amovr (ax, Memr[nax], nres) call amovr (ay, Memr[nay], nres) do j = 1, nres { if (IS_INDEFR(Memr[nax+j-1]) || IS_INDEFR(Memr[nay+j-1])) next ix = int(Memr[nax+j-1] - sampoff) iy = int(Memr[nay+j-1] - lineoff) if ( ix < 1-RM_BOUNDWIDTH || ix > IM_LEN(FOC_IMPT(i_foc),1)+RM_BOUNDWIDTH-1 || iy < 1-RM_BOUNDWIDTH || iy > IM_LEN(FOC_IMPT(i_foc),1)+RM_BOUNDWIDTH-1) { next } buf = imgs2r (FOC_IMPT(i_foc), ix-RM_HBOXX(res), ix+RM_HBOXX(res), iy-RM_HBOXY(res), iy+RM_HBOXY(res)) minpos = xside*RM_HBOXY(res)+RM_HBOXX(res)+1 do i = 1, xside*yside { if (Memr[buf+i-1] < Memr[buf+minpos-1]) minpos = i } Memr[nax+j-1] = ix-RM_HBOXX(res)-1 + mod(minpos, xside) + sampoff Memr[nay+j-1] = iy-RM_HBOXY(res)-1 + (minpos+xside-1)/xside + lineoff } call rm_show (gd, ax, ay, Memr[nax], Memr[nay], Memr[sax], Memr[say], xa, ya, tipt, nfit, nres) } # specify 4 (Four) points to do fit case 'f' : do i = 1, 4 { if (i > 1) { call printf (" next reseau to mark for move\n") junk = clgcur("cursor", wx, wy, wcs, key, Memc[strval], SZ_LINE) } tipt[i] = rm_mpoint (wx, wy, ax, ay, nres) call gmark (gd, ax[tipt[i]], ay[tipt[i]], GM_PLUS + GM_CROSS, 2.0, 2.0) call printf (" where should this be?\n") junk = rm_gcur ("cursor", xa[i], ya[i], wcs, key, Memc[strval], SZ_LINE) call gmark (gd, xa[i], ya[i], GM_BOX, 1.0, 1.0) call rm_log (fd, 'f', tipt[i], nx, ax[tipt[i]], ay[tipt[i]], xa[i], ya[i]) } call printf (" \n") call rm_fit4 (ax, ay, xa, ya, tipt, nres, Memr[nax], Memr[nay]) nfit = 4 call rm_show (gd, ax, ay, Memr[nax], Memr[nay], Memr[sax], Memr[say], xa, ya, tipt, nfit, nres) # specify 1 (One) point, then shift case 'o' : tipt[1] = rm_mpoint (wx, wy, ax, ay, nres) call gmark (gd, ax[tipt[1]], ay[tipt[1]], GM_PLUS + GM_CROSS, 2.0, 2.0) call printf (" where should this be?\n") junk = rm_gcur ("cursor", xa[1], ya[1], wcs, key, Memc[strval], SZ_LINE) call gmark (gd, xa[1], ya[1], GM_BOX, 1.0, 1.0) call rm_log (fd, 'o', tipt[1], nx, ax[tipt[1]], ay[tipt[1]], xa[1], ya[1]) call printf (" \n") call rm_fit1 (ax, ay, xa, ya, tipt, nres, Memr[nax], Memr[nay]) nfit = 1 call rm_show (gd, ax, ay, Memr[nax], Memr[nay], Memr[sax], Memr[say], xa, ya, tipt, nfit, nres) # change a particular Point case 'p' : ipoint = rm_mpoint (wx, wy, ax, ay, nres) oldx = ax[ipoint] oldy = ay[ipoint] call gmark (gd, ax[ipoint], ay[ipoint], GM_PLUS + GM_CROSS, 2.0, 2.0) call printf (" where should this be?\n") junk = rm_gcur ("cursor", ax[ipoint], ay[ipoint], wcs, key, Memc[strval], SZ_LINE) call gmark (gd, ax[ipoint], ay[ipoint], GM_BOX, 1.0, 1.0) call printf (" \n") call amovr (ax, Memr[nax], nres) call amovr (ay, Memr[nay], nres) call rm_log (fd, 'p', ipoint, nx, oldx, oldy, ax[ipoint], ay[ipoint]) # Redraw the marked positions case 'r' : call sprintf (Memc[strval], SZ_LINE, "%s %d x %d") call pargstr (entry) call pargi (nx) call pargi (ny) call gclear (gd) call gswind (gd, RM_X1(res), RM_X2(res), RM_Y1(res), RM_Y2(res)) call glabax (gd, Memc[strval], "axis1", "axis2") call gseti (gd, G_PMLTYPE, GL_SOLID) call gpmark (gd, Memr[nax], Memr[nay], nres, GM_PLUS, 2.0, 2.0) # Save positions. case 's' : call rm_show (gd, ax, ay, Memr[nax], Memr[nay], Memr[sax], Memr[say], xa, ya, tipt, nfit, nres) call printf (" save current positions\n") if (fd != NULL) call fprintf (fd, "Save current positions.\n") # specify 2 (Two) points to do fit case 't' : do i = 1, 2 { if (i > 1) { call printf (" next reseau to mark for move\n") junk = clgcur ("cursor", wx, wy, wcs, key, Memc[strval], SZ_LINE) } tipt[i] = rm_mpoint (wx, wy, ax, ay, nres) call gmark (gd, ax[tipt[i]], ay[tipt[i]], GM_PLUS + GM_CROSS, 2.0, 2.0) call printf (" where should this be?\n") junk = rm_gcur ("cursor", xa[i], ya[i], wcs, key, Memc[strval], SZ_LINE) call gmark (gd, xa[i], ya[i], GM_BOX, 1.0, 1.0) call rm_log (fd, 't', tipt[i], nx, ax[tipt[i]], ay[tipt[i]], xa[i], ya[i]) } call printf (" \n") call rm_fit2 (ax, ay, xa, ya, tipt, nres, Memr[nax], Memr[nay]) nfit = 2 call rm_show (gd, ax, ay, Memr[nax], Memr[nay], Memr[sax], Memr[say], xa, ya, tipt, nfit, nres) # Undo changes since last call to rm_show. case 'u': call amovr (Memr[sax], Memr[nax], nres) call amovr (Memr[say], Memr[nay], nres) call rm_show (gd, ax, ay, Memr[nax], Memr[nay], Memr[sax], Memr[say], xa, ya, tipt, nfit, nres) call amovr (Memr[nax], Memr[sax], nres) call amovr (Memr[nay], Memr[say], nres) call printf (" restore to saved positions\n") if (fd != NULL) call fprintf (fd, "Restore to saved positions.\n") # Where are we? Print coordinates of reseau nearest cursor. case 'w': ipoint = rm_mpoint (wx, wy, ax, ay, nres) i = mod (ipoint, nx) if (i < 1) i = nx j = (ipoint - i) / nx + 1 call printf (" index: %d %d pixels: %.2f %.2f\n") call pargi (i) call pargi (j) call pargr (ax[ipoint]) call pargr (ay[ipoint]) default: if (key > ' ') call printf (" command `%c' not found\n") else call printf (" command `%xX' not found\n") call pargi (key) } } call sfree (sp) end # RM_SHOW -- show the new locations # Show where the current positions are, and update the ax,ay arrays. procedure rm_show (gd, ax, ay, nax, nay, sax, say, xa, ya, tipt, nfit, nres) pointer gd # i: graphics pointer real ax[ARB] # io: current x positions real ay[ARB] # io: current y positions real nax[ARB] # i: new x pos real nay[ARB] # i: new y pos real sax[ARB] # io: saved x pos real say[ARB] # io: saved y pos real xa[ARB] # i: markers x pos real ya[ARB] # i: markers y pos int tipt[ARB] # i: index to markers int nfit # i: number of marker pairs ( to be erased ) int nres # i: number of reseau #-- int i begin # remember where we were call amovr (ax, sax, nres) call amovr (ay, say, nres) call gseti (gd, G_PMLTYPE, GL_CLEAR) call gpmark (gd, ax, ay, nres, GM_PLUS, 2.0, 2.0) do i = 1, nfit { call gmark (gd, ax[tipt[i]], ay[tipt[i]], GM_PLUS + GM_CROSS, 2.0, 2.0) call gmark (gd, xa[i], ya[i], GM_BOX, 1.0, 1.0) } nfit = 0 call gseti (gd, G_PMLTYPE, GL_SOLID) call gpmark (gd, nax, nay, nres, GM_PLUS, 2.0, 2.0) call amovr (nax, ax, nres) call amovr (nay, ay, nres) end # RM_MPOINT -- choose positions closest to the cursor int procedure rm_mpoint (wx, wy, ax, ay, nres) real wx # i: x cursor coord real wy # i: y cursor coord real ax[ARB] # i: points to choose from, x coords real ay[ARB] # i: points to choose from, y coords int nres # i: number of positions #-- real dmin, dist int imin, i begin imin = nres dmin = INDEFR do i = 1, nres { if (!IS_INDEFR(ax[i]) && !IS_INDEFR(ay[i])) { dist = sqrt ((wx-ax[i])**2 + (wy-ay[i])**2) if ( dist < dmin ) { dmin = dist imin = i } } } return (imin) end # RM_FIT1 -- move reseaux by simple shift, based on info for one point procedure rm_fit1 (ax, ay, xa, ya, tipt, nres, newax, neway) real ax[ARB] # i: old x positions real ay[ARB] # i: old y positions real xa[4] # i: new selected x positions real ya[4] # i: new selected y positions int tipt[4] # i: index array for ax, ay int nres # i: number of reseau positions real newax[ARB] # o: new x positions real neway[ARB] # o: new y positions #-- int i real a, b begin a = xa[1] - ax[tipt[1]] b = ya[1] - ay[tipt[1]] do i = 1, nres { if (IS_INDEFR(ax[i]) || IS_INDEFR(ay[i])) { newax[i] = INDEFR neway[i] = INDEFR } else { newax[i] = ax[i] + a neway[i] = ay[i] + b } } end # RM_FIT2 -- move reseau positions based on info for two points # simple rotation, and scale procedure rm_fit2 (ax, ay, xa, ya, tipt, nres, newax, neway) real ax[ARB] # i: old x positions real ay[ARB] # i: old y positions real xa[4] # i: new selected x positions real ya[4] # i: new selected y positions int tipt[4] # i: index array for ax, ay int nres # i: number of reseau positions real newax[ARB] # o: new x positions real neway[ARB] # o: new y positions #-- int i real a, b real p1, p2 real x1, x2, y1, y2 begin x1 = ax[tipt[1]] x2 = ax[tipt[2]] y1 = ay[tipt[1]] y2 = ay[tipt[2]] p1 = xa[1] p2 = xa[2] a = (p1*y2-p2*y1)/(x1*y2-x2*y1) b = (p2*x1-p1*x2)/(x1*y2-x2*y1) do i = 1, nres { if (IS_INDEFR(ax[i]) || IS_INDEFR(ay[i])) newax[i] = INDEFR else newax[i] = a*ax[i] + b*ay[i] } p1 = ya[1] p2 = ya[2] a = (p1*y2-p2*y1)/(x1*y2-x2*y1) b = (p2*x1-p1*x2)/(x1*y2-x2*y1) do i = 1, nres { if (IS_INDEFR(ax[i]) || IS_INDEFR(ay[i])) neway[i] = INDEFR else neway[i] = a*ax[i] + b*ay[i] } end # RM_FIT4 -- move reseau based on info for 4 points # blinear fit these shifts and evaluate for remaining reseau procedure rm_fit4 (ax, ay, xa, ya, tipt, nres, newax, neway) real ax[ARB] # i: old x positions real ay[ARB] # i: old y positions real xa[4] # i: new selected x positions real ya[4] # i: new selected y positions int tipt[4] # i: index array for ax, ay int nres # i: number of reseau positions real newax[ARB] # o: new x positions real neway[ARB] # o: new y positions #-- int i real a, b, c, d real p1, p2, p3, p4 real x1, x2, x3, x4, y1, y2, y3, y4 real at, bt, ct, dt, mt, nt begin x1 = ax[tipt[1]] x2 = ax[tipt[2]] x3 = ax[tipt[3]] x4 = ax[tipt[4]] y1 = ay[tipt[1]] y2 = ay[tipt[2]] y3 = ay[tipt[3]] y4 = ay[tipt[4]] p1 = xa[1] p2 = xa[2] p3 = xa[3] p4 = xa[4] at = (x3-x1)*(y2-y1) - (x2-x1)*(y3-y1) bt = (x2*y2-x1*y1)*(x3-x1) - (x3*y3-x1*y1)*(x2-x1) ct = (x4-x1)*(y2-y1) - (x2-x1)*(y4-y1) dt = (x2*y2-x1*y1)*(x4-x1) - (x4*y4-x1*y1)*(x2-x1) mt = (x3-x1)*(p2-p1) - (x2-x1)*(p3-p1) nt = (x4-x1)*(p2-p1) - (x2-x1)*(p4-p1) c = (mt*dt - nt*bt)/(dt*at - bt*ct) d = (mt*ct - nt*at)/(bt*ct - dt*at) if (x2 != x1) b = ( (p2-p1) - c*(y2-y1) - d*(x2*y2 - x1*y1) )/ ( x2-x1) else if ( x3 != x1 ) b = ( (p3-p1) - c*(y3-y1) - d*(x3*y3 - x1*y1) )/ ( x2-x1) else if ( x4 != x1 ) b = ( (p4-p1) - c*(y4-y1) - d*(x4*y4 - x1*y1) )/ ( x2-x1) else { call printf (" 4 points must not be colinear\n") call amovr (ax, newax, nres) call amovr (ay, neway, nres) return } a = p1 - b*x1 - c*y1 - d*x1*y1 do i = 1, nres { if (IS_INDEFR(ax[i]) || IS_INDEFR(ay[i])) newax[i] = INDEFR else newax[i] = a + b*ax[i] + c*ay[i] + d*ax[i]*ay[i] } p1 = ya[1] p2 = ya[2] p3 = ya[3] p4 = ya[4] at = (x3-x1)*(y2-y1) - (x2-x1)*(y3-y1) bt = (x2*y2-x1*y1)*(x3-x1) - (x3*y3-x1*y1)*(x2-x1) ct = (x4-x1)*(y2-y1) - (x2-x1)*(y4-y1) dt = (x2*y2-x1*y1)*(x4-x1) - (x4*y4-x1*y1)*(x2-x1) mt = (x3-x1)*(p2-p1) - (x2-x1)*(p3-p1) nt = (x4-x1)*(p2-p1) - (x2-x1)*(p4-p1) c = (mt*dt - nt*bt)/(dt*at - bt*ct) d = (mt*ct - nt*at)/(bt*ct - dt*at) if (x2 != x1) b = ( (p2-p1) - c*(y2-y1) - d*(x2*y2 - x1*y1) )/ ( x2-x1) else if ( x3 != x1 ) b = ( (p3-p1) - c*(y3-y1) - d*(x3*y3 - x1*y1) )/ ( x2-x1) else if ( x4 != x1 ) b = ( (p4-p1) - c*(y4-y1) - d*(x4*y4 - x1*y1) )/ ( x2-x1) else { call printf (" 4 points must not be colinear\n") call amovr (ax, newax, nres) call amovr (ay, neway, nres) return } a = p1 - b*x1 - c*y1 - d*x1*y1 do i = 1, nres { if (IS_INDEFR(ax[i]) || IS_INDEFR(ay[i])) neway[i] = INDEFR else neway[i] = a + b*ax[i] + c*ay[i] + d*ax[i]*ay[i] } end # RM_COLON -- process colon commands. We are careful about saying we need to # reread the image - only done if the value is actually altered procedure rm_colon (gd, res, key, cmdstr) pointer gd # i: GIO pointer pointer res # i: RM pointer int key # i: key from clgcur char cmdstr[ARB] # i: Command string #-- char cmd[SZ_LINE] int ncmd, ival real rval int nscan(), strdic() string cmds "|show|x1|x2|y1|y2|xsearch|ysearch|" begin call sscan (cmdstr) call gargwrd (cmd, SZ_LINE) if (nscan() == 0) return ncmd = strdic (cmd, cmd, SZ_LINE, cmds) switch (ncmd) { case 1: call gdeactivate (gd, AW_CLEAR) call printf ("Current Parameters\n\n") call printf ("\tx1 \t= %g\n") call pargr (RM_X1(res)) call printf ("\tx2 \t= %g\n") call pargr (RM_X2(res)) call printf ("\ty1 \t= %g\n") call pargr (RM_Y1(res)) call printf ("\ty2 \t= %g\n") call pargr (RM_Y2(res)) call printf ("\txsearch \t= %d\n") call pargi (2*RM_HBOXX(res)+1) call printf ("\tysearch \t= %d\n") call pargi (2*RM_HBOXY(res)+1) call greactivate (gd, AW_PAUSE) case 2: call gargr (rval) if (nscan () == 1) { call printf ("x1 = %g\n") call pargr (RM_X1(res)) } else { RM_X1(res) = rval } case 3: call gargr (rval) if (nscan () == 1) { call printf ("x2 = %g\n") call pargr (RM_X2(res)) } else { RM_X2(res) = rval } case 4: call gargr (rval) if (nscan () == 1) { call printf ("y1 = %g\n") call pargr (RM_Y1(res)) } else { RM_Y1(res) = rval } case 5: call gargr (rval) if (nscan () == 1) { call printf ("y2 = %g\n") call pargr (RM_Y2(res)) } else { RM_Y2(res) = rval } case 6: call gargi (ival) if (nscan () == 1) { call printf ("xsearch = %d\n") call pargi (2*RM_HBOXX(res)+1) } else { if (ival >= 0) { RM_HBOXX(res) = min ((ival + 1) / 2, RM_BOUNDWIDTH) } } case 7: call gargi (ival) if (nscan () == 1) { call printf ("ysearch = %d\n") call pargi (2*RM_HBOXY(res)+1) } else { if (ival >= 0) { RM_HBOXY(res) = min ((ival + 1) / 2, RM_BOUNDWIDTH) } } } end # This routine opens the input and output reseau tables. procedure rm_open (inres, outres, irp, orp, inplace) char inres[ARB] # i: name of input reseau table char outres[ARB] # i: output reseau table pointer irp # o: pointer to struct for inres pointer orp # o: pointer to struct for outres; # this will be the same as irp if inplace is true bool inplace # o: true if outres is same as inres #-- pointer rs_open() int tbtacc() bool streq() begin if (outres[1] == EOS || streq (inres, outres)) { # Write the result into the input reseau table. inplace = true irp = rs_open (inres, READ_WRITE, NULL) orp = irp } else { # Write the result into a different output table. inplace = false irp = rs_open (inres, READ_ONLY, NULL) if (tbtacc (outres) == YES) # does it exist? orp = rs_open (outres, READ_WRITE, NULL) else orp = rs_open (outres, NEW_COPY, irp) } end # rm_gcur -- read the cursor position # This routine calls clgcur to get the cursor position. If the cursor key # is a colon, the X and Y coordinates are read from the string; otherwise, # the coordinates returned by clgcur are copied to the output parameters. # If the cursor key is a colon, either a comma or whitespace may be used # as a separator between the X and Y positions. If either X or Y is not # specified, EOF will be returned as the function value. int procedure rm_gcur (param, x, y, wcs, key, strval, maxch) char param[ARB] # i: parameter name for graphics cursor real x, y # o: cursor position int wcs # o: wcs to which coordinates belong int key # o: keystroke value of cursor event char strval[ARB] # o: string value, if any; used as scratch int maxch # i: size of strval #-- real tx, ty # temporary values of x, y int nitems int ip, ctor() int clgcur() begin nitems = clgcur (param, tx, ty, wcs, key, strval, maxch) if (nitems == EOF) return (EOF) if (key == ':') { ip = 1 if (ctor (strval, ip, x) < 1) return (EOF) if (strval[ip] == ',') ip = ip + 1 if (ctor (strval, ip, y) < 1) return (EOF) } else { x = tx y = ty } return (nitems) end # rm_log -- log a change to the logfile procedure rm_log (fd, key, ipoint, nx, oldx, oldy, newx, newy) int fd # i: fd for log file, or NULL if none int key # i: cursor key int ipoint # i: index in arrays of pixel coordinates int nx # i: number of columns in reseau grid real oldx, oldy # i: previous X and Y pixel coordinates real newx, newy # i: new X and Y pixel coordinates #-- int i, j # column and row indices of reseau mark begin if (fd == NULL) return # Convert ipoint to [i,j]. i = mod (ipoint, nx) if (i < 1) i = nx j = (ipoint - i) / nx + 1 if (key > ' ') call fprintf (fd, "'%c' [%d,%d] (%.2f,%.2f) --> (%.2f,%.2f)\n") else call fprintf (fd, "'%xX' [%d,%d] (%.2f,%.2f) --> (%.2f,%.2f)\n") call pargi (key) call pargi (i) call pargi (j) call pargr (oldx) call pargr (oldy) call pargr (newx) call pargr (newy) end