# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. include include include include include include "../lib/skywcs/skywcs.h" include "wcspix.h" # Image class data. define LEN_IMGDATA 15 define IMG_WP Memi[$1 ] # wcspix back-pointer define IMG_IM Memi[$1+1] # image pointer define IMG_BPM Memi[$1+2] # bad pixel mask pointer define IMG_MW Memi[$1+3] # image wcs pointer define IMG_CO Memi[$1+4] # skywcs transform pointer define IMG_CTW Memi[$1+5] # mwcs log->world transform ptr define IMG_CTP Memi[$1+6] # mwcs log->phys transform ptr define IMG_CTA Memi[$1+7] # mwcs log->amplifier transform define IMG_CTD Memi[$1+8] # mwcs log->detector transform define IMG_ROT Memr[$1+9] # rotation angle define IMG_SCALE Memr[$1+10] # plate scale define IMG_LINEAR Memi[$1+11] # linear coords define IMG_DEBUG FALSE # IMG_INIT -- Initialize the object structure. procedure img_init (cp, wp) pointer cp #i cache pointer pointer wp #i WCSPIX structure pointer img # data pointer begin if (IMG_DEBUG) call printf ("img_init: \n") # Allocate the image data structure if not previously allocated. if (C_DATA(cp) == NULL) { iferr (call calloc (C_DATA(cp), LEN_IMGDATA, TY_STRUCT)) return } img = C_DATA(cp) IMG_WP(img) = wp IMG_IM(img) = NULL IMG_MW(img) = NULL IMG_CO(img) = NULL IMG_CTW(img) = NULL IMG_CTP(img) = NULL IMG_ROT(img) = 0.0 IMG_SCALE(img) = 0.0 IMG_LINEAR(img) = YES end # IMG_CACHE -- Cache an image in the object cache. procedure img_cache (cp, objid, regid, ref) pointer cp #i cache pointer int objid #i object id int regid #i region id char ref[ARB] #i object reference pointer img, im, wp int stat char alert[SZ_LINE] pointer immap(), ds_pmmap(), mw_sctran() pointer img_amp_wcs(), img_det_wcs() int imaccf(), sk_decim() errchk immap, ds_pmmap(), mw_sctran, sk_decim begin if (IMG_DEBUG) call printf ("img_cache: \n") # Now map the image and WCS. img = C_DATA(cp) wp = IMG_WP(img) iferr (IMG_IM(img) = immap (ref, READ_ONLY, 0)) { # Send alert to the GUI. call sprintf (alert, SZ_FNAME, "Unable to cache\n%s") call pargstr (ref) call xim_alert (alert, "", "") return } IMG_CO(img) = NULL IMG_CTW(img) = NULL IMG_CTP(img) = NULL iferr { stat = sk_decim (IMG_IM(img), "world", IMG_MW(img), IMG_CO(img)) if (stat == ERR || IMG_MW(img) == NULL) IMG_LINEAR(img) = YES if (IMG_MW(img) != NULL) { IMG_CTW(img) = mw_sctran (IMG_MW(img), "logical", "world", 03B) IMG_CTP(img) = mw_sctran (IMG_MW(img), "logical", "physical", 03B) # Get the amplifier transformation values if present. im = IMG_IM(img) if (imaccf(im,"ATM1_1") == YES && imaccf(im,"ATM2_2") == YES && imaccf(im,"ATV1") == YES && imaccf(im,"ATV2") == YES) IMG_CTA(img) = img_amp_wcs (im, IMG_MW(img)) if (imaccf(im,"DTM1_1") == YES && imaccf(im,"DTM2_2") == YES && imaccf(im,"DTV1") == YES && imaccf(im,"DTV2") == YES) IMG_CTD(img) = img_det_wcs (im, IMG_MW(img)) } } then { # Send alert to the GUI. call sprintf (alert, SZ_FNAME, "Unable to decode image WCS\n%s") call pargstr (ref) call xim_alert (alert, "", "") IMG_LINEAR(img) = YES } # See if we can find a bad pixel mask. if (WP_BPM(wp) == YES) { iferr (IMG_BPM(img) = ds_pmmap ("BPM", IMG_IM(img))) IMG_BPM(img) = NULL } C_OBJID(cp) = objid C_REGID(cp) = regid C_NREF(cp) = C_NREF(cp) + 1 call strcpy (ref, C_REF(cp), 128) end # IMG_UNCACHE -- Uncache an image in the object cache. procedure img_uncache (cp, id) pointer cp #i cache pointer int id #i image id pointer img begin if (IMG_DEBUG) call printf ("img_uncache: \n") C_OBJID(cp) = NULL C_NREF(cp) = 0 call strcpy ("", C_REF(cp), SZ_FNAME) img = C_DATA(cp) if (IMG_MW(img) != NULL) call mw_close (IMG_MW(img)) if (IMG_BPM(img) != NULL) call imunmap (IMG_BPM(img)) if (IMG_IM(img) != NULL) call imunmap (IMG_IM(img)) IMG_IM(img) = NULL IMG_BPM(img) = NULL IMG_MW(img) = NULL IMG_CTW(img) = NULL IMG_CTP(img) = NULL IMG_CO(img) = NULL IMG_ROT(img) = 0.0 IMG_SCALE(img) = 0.0 IMG_LINEAR(img) = NO call mfree (C_DATA(cp), TY_STRUCT) C_DATA(cp) = NULL end # IMG_WCSTRAN -- Translate object source (x,y) coordinates to the # desired output WCSs. Message is returned as something like: # # set value { # { object } { region } # { pixval [] } # { bpm } # { coord [ ] } # { coord [ ] } # : # } procedure img_wcstran (cp, id, x, y) pointer cp #i cache pointer int id #i image id real x, y #i source coords pointer img, im, wp, co double dx, dy, wx, wy, pixval real rx, ry int i, bpm # Use static storage to avoid allocation overhead. char buf[SZ_LINE] char msg[SZ_LINE], wcs[LEN_WCSNAME], xc[LEN_WCSNAME], yc[LEN_WCSNAME] char xunits[LEN_WCSNAME], yunits[LEN_WCSNAME] double sk_statd() begin if (IMG_DEBUG) call printf ("img_wcstran: \n") img = C_DATA(cp) # initialize co = IMG_CO(img) wp = IMG_WP(img) im = IMG_IM(img) # Get the translation to the image section. dx = (double(x) - sk_statd(co,S_VXOFF)) / sk_statd(co,S_VXSTEP) dy = (double(y) - sk_statd(co,S_VYOFF)) / sk_statd(co,S_VYSTEP) rx = dx ry = dy # Read the pixel data. call img_get_data (cp, id, rx, ry, pixval, bpm) # Begin formatting the message. call aclrc (msg, SZ_LINE) call sprintf (msg, SZ_LINE, "wcstran { object %d } { region %d } ") call pargi (C_OBJID(cp)) call pargi (C_REGID(cp)) call sprintf (buf, SZ_LINE, "{ pixval %9.9g } { bpm %d }\n") call pargd (pixval) call pargi (bpm) call strcat (buf, msg, SZ_LINE) # Now loop over the requested systems and generate a coordinate # for each. for (i=1; i <= MAX_WCSLINES; i=i+1) { # Get the coordinate value. call img_get_coord (img, dx, dy, SYSTEMS(wp,i), WCSNAME(wp,i), wx, wy) # Get the system name, labels, and formats strings for the WCS. call img_coord_labels (cp, i, wcs, xunits, yunits) # Format the values as requested. call img_coord_fmt (cp, i, wx, wy, xc, yc) # Format the coord buffer and append it to the message. call sprintf (buf, SZ_LINE, "{coord {%9s} {%12s} {%12s} {%4s} {%4s}}\n") call pargstr (wcs) call pargstr (xc) call pargstr (yc) call pargstr (xunits) call pargstr (yunits) call strcat (buf, msg, SZ_LINE) } # Now send the completed message. call wcspix_message (msg); end # IMG_WCSLIST -- List the WCSs available for the given image. procedure img_wcslist (cp, id) pointer cp #i cache pointer int id #i image id pointer img, im, mw char msg[SZ_LINE] begin if (IMG_DEBUG) call printf ("img_wcslist: \n") img = C_DATA(cp) # initialize mw = IMG_MW(img) im = IMG_IM(img) call strcpy ("wcslist {None Logical World Physical line ", msg, SZ_LINE) # See if we can do amplifier/detector coords by checking for ATM/ATV # and DTM/DTV keywords. if (IMG_CTA(img) != NULL) call strcat (" Amplifier ", msg, SZ_LINE) if (IMG_CTD(img) != NULL) call strcat (" Detector ", msg, SZ_LINE) if (IMG_CTA(img) != NULL || IMG_CTD(img) != NULL) call strcat (" CCD ", msg, SZ_LINE) call strcat (" line ", msg, SZ_LINE) # If we have a MWCS pointer list the sky projections. if (mw != NULL) call strcat (SKYPROJ, msg, SZ_LINE) # Close the message. call strcat ("}", msg, SZ_LINE) call wcspix_message (msg) end # IMG_GET_DATA -- Get data from the image. procedure img_get_data (cp, id, x, y, pixval, bpm_pix) pointer cp #i cache pointer int id #i image id real x, y #i source coords double pixval #o central pixel value int bpm_pix #o bad pixel mask value pointer img, wp, im, bpm, pix int nl, nc, ix, iy int size, x1, x2, y1, y2 pointer imgs2r(), imgs2i() begin if (IMG_DEBUG) call printf ("img_get_data: \n") img = C_DATA(cp) wp = IMG_WP(img) im = IMG_IM(img) bpm = IMG_BPM(img) nc = IM_LEN(im,1) nl = IM_LEN(im,2) size = WP_PTABSZ(wp) # Sanity check on the cursor image position. if (x < 0.0 || y < 0.0 || x > nc || y > nl) return # Bounds checking. Rather than deal with out of bounds pixels we'll # adjust the center pixel so we get the same size raster up to each # boundary. ix = int (x + 0.5) ; iy = int (y + 0.5) ix = max (size/2+1, ix) ; iy = max (size/2+1, iy) ix = min (ix, (nc-(size/2)-1)) ; iy = min (iy, (nl-(size/2)-1)) # Compute the box offset given the center and size. x1 = ix - size / 2 + 0.5 x2 = ix + size / 2 + 0.5 y1 = iy - size / 2 + 0.5 y2 = iy + size / 2 + 0.5 # Get the image pixels x1 = max (1, x1) x2 = min (nc, x2) y1 = max (1, y1) y2 = min (nl, y2) pix = imgs2r (im, int(x1), int(x2), int(y1), int(y2)) if (bpm != NULL && WP_BPM(wp) == YES) bpm_pix = Memi[imgs2i (bpm, ix, ix, iy, iy)] else bpm_pix = 0 # Compute the image pixel associated with the requested coords. pixval = Memr[pix + ((size/2)*size) + (size/2)] * 1.0d0 # Send the pixel table. if (WP_PTABSZ(wp) > 1) call img_send_pixtab (Memr[pix], WP_PTABSZ(wp), x1, x2, y1, y2) end # IMG_OBJINFO -- Get header information from the image. procedure img_objinfo (cp, id, template) pointer cp #i cache pointer int id #i image id char template[ARB] #i keyword template pointer im, img define WCS_TEMPLATE "WCSDIM,CTYPE*,CRPIX*,CRVAL*,CD*,CROTA2,LTV*,LTM*,WSV*,WAT*,RA*,DEC*,EQUINOX,EPOCH,MJD*,DATE-OBS" begin if (IMG_DEBUG) call printf ("img_objinfo: \n") # Send the full header (or keyword filtered header), only the WCS # keywords, and a plain-text explanation of the WCS. img = C_DATA(cp) im = IMG_IM(img) call img_send_header (im, "imghdr", template) call img_send_header (im, "wcshdr", WCS_TEMPLATE) call img_send_wcsinfo (im, cp) call img_send_compass (im, cp) end #============================================================================== # IMG_SEND_HEADER -- Send an image header to the named GUI object. Keywords # are filtered according to a specified template procedure img_send_header (im, object, template) pointer im #i image descriptor char object[ARB] #i object for the message char template[ARB] #i keyword template pointer sp, hdr, lbuf, line, field, keyw, dict pointer ip, lp, list int nlines, in, out, i, hdr_size bool keyw_filter int stropen(), getline(), stridx(), imgnfn(), strdic() pointer imofnlu() bool streq() errchk stropen, getline, putci, putline, imgnfn, imofnlu, strdic define USER_AREA Memc[($1+IMU-1)*SZ_STRUCT + 1] define SZ_KEYW 8 begin hdr_size = (LEN_IMDES + IM_LENHDRMEM(im) - IMU) * SZ_STRUCT - 1 hdr_size = hdr_size + SZ_LINE call smark (sp) call salloc (hdr, hdr_size, TY_CHAR) call salloc (dict, hdr_size, TY_CHAR) call salloc (field, SZ_LINE, TY_CHAR) call salloc (lbuf, SZ_LINE, TY_CHAR) call salloc (line, SZ_LINE, TY_CHAR) call salloc (keyw, SZ_KEYW, TY_CHAR) in = stropen (USER_AREA(im), hdr_size, READ_ONLY) out = stropen (Memc[hdr], hdr_size, WRITE_ONLY) call fprintf (out, "%s {") call pargstr (object) # Build up a dictionary of header keywords based on the template. keyw_filter = (!streq (template, "*")) if (keyw_filter) { list = imofnlu (im, template) call strcpy ("|", Memc[dict], hdr_size) while (imgnfn (list, Memc[field], SZ_FNAME) != EOF) { call strcat (Memc[field], Memc[dict], hdr_size) call strcat ("|", Memc[dict], hdr_size) } call imcfnl (list) } # Copy header records to the output, stripping any trailing # whitespace and clipping at the right margin. We also filter # against the keyword dictionary found above. nlines = 0 while (getline (in, Memc[lbuf]) != EOF) { call aclrc (Memc[line], SZ_LINE) # Escape any brackets passed to the Tcl. ip = lbuf lp = line while (Memc[ip] != EOS && Memc[ip] != '\n') { if (stridx (Memc[ip], "[{") > 0) { Memc[lp] = '\\' lp = lp + 1 } Memc[lp] = Memc[ip] ip = ip + 1 lp = lp + 1 } Memc[lp] = '\n' Memc[lp+1] = EOS # See whether the line matches a keyword we want to output. if (keyw_filter) { for (i=0; i < SZ_KEYW && !IS_WHITE(Memc[line+i]); i=i+1) Memc[keyw+i] = Memc[line+i] Memc[keyw+i] = '\0' # If not in the dictionary skip to the next line. if (strdic (Memc[keyw], Memc[keyw], SZ_KEYW, Memc[dict]) == 0) next } call putci (out, ' ') call putline (out, Memc[line]) # Send the header in small chunks so we don't overflow the # message buffer. nlines = nlines + 1 if (mod(nlines,10) == 0) { call fprintf (out, "}") call close (out) call wcspix_message (Memc[hdr]); call aclrc (Memc[hdr], hdr_size) out = stropen (Memc[hdr], hdr_size, WRITE_ONLY) call fprintf (out, "%s {") call pargstr (object) } } call fprintf (out, "}") call close (in) call close (out) # Send the final message. call wcspix_message (Memc[hdr]) # Pad a few lines for the GUI call sprintf (Memc[hdr], SZ_LINE, "%d { \n\n\n }") call pargstr (object) call wcspix_message (Memc[hdr]) call sfree (sp) end # IMG_SEND_COMPASS -- Send information about the image WCS in a plain-english # string. procedure img_send_compass (im, cp) pointer im #i image descriptor pointer cp #i cache element pointer pointer sp, buf, img, co double cx, cy, cx1, cy1, dx, dy, x1, y1 double cosa, sina, angle int i, j, comp_x, comp_y long axis[IM_MAXDIM], lv[IM_MAXDIM], pv1[IM_MAXDIM], pv2[IM_MAXDIM] begin call smark (sp) call salloc (buf, SZ_LINE, TY_CHAR) call aclrc (Memc[buf], SZ_LINE) # Get the data pointer. img = C_DATA(cp) co = IMG_CO(img) # Get world coords at the image corners. if (IMG_CTW(img) != NULL) { if (IMG_ROT(img) > 0.0) angle = -IMG_ROT(img) else angle = IMG_ROT(img) + 360.0 cosa = cos (DEGTORAD(angle)) sina = sin (DEGTORAD(angle)) # Image center position cx = IM_LEN(im,1) / 2.0d0 cy = IM_LEN(im,2) / 2.0d0 call mw_c2trand (IMG_CTW(img), cx, cy, cx1, cy1) # Extend a unit vector up from the center assuming it's North # and rotate it by the wcs angle. dx = cx + ( 10.0 * sina) dy = cy + ( 10.0 * cosa) call mw_c2trand (IMG_CTW(img), dx, dy, x1, y1) # Check new point Y value relative to the center position. if (y1 >= cy1) comp_y = 1 # North is up else comp_y = -1 # North is down # Extend a unit vector left from the center assuming it's East # and rotate it by the wcs angle. dx = cx + (-10.0 * cosa) dy = cy + ( 10.0 * sina) call mw_c2trand (IMG_CTW(img), dx, dy, x1, y1) # Check new point X value relative to the center position. if (x1 >= cx1) comp_x = 1 # East is left and we have a WCS else comp_x = -1 # East is right } else { # Determine the logical to physical mapping by evaluating two # points and determining the axis reduction if any. pv1 will be # the offset and pv2-pv1 will be the scale. lv[1] = 0; lv[2] = 0; call imaplv (im, lv, pv1, 2) lv[1] = 1; lv[2] = 1; call imaplv (im, lv, pv2, 2) i = 1 axis[1] = 1; axis[2] = 2 do j = 1, IM_MAXDIM { if (pv1[j] != pv2[j]) { axis[i] = j i = i + 1 } } comp_x = - (pv2[axis[1]] - pv1[axis[1]]) comp_y = (pv2[axis[2]] - pv1[axis[2]]) } call sprintf (Memc[buf], SZ_LINE, "compass %d %g %d %d %s\0") call pargi (C_OBJID(cp)) call pargr (IMG_ROT(img)) call pargi (comp_x) call pargi (comp_y) if (IMG_MW(img) != NULL) call pargstr ("E N") else call pargstr ("X Y") call wcspix_message (Memc[buf]) call sfree (sp) end # IMG_SEND_WCSINFO -- Send information about the image WCS in a plain-english # string. procedure img_send_wcsinfo (im, cp) pointer im #i image descriptor pointer cp #i cache element pointer pointer sp, co, img, mw pointer buf, proj, radecstr int fd, radecsys, ctype, wtype, ndim double crpix1, crpix2, crval1, crval2, cval1, cval2 double xscale, yscale, xrot, yrot double r[IM_MAXDIM], w[IM_MAXDIM], cd[IM_MAXDIM,IM_MAXDIM], int idxstr(), sk_stati(), stropen(), mw_stati() double sk_statd(), sl_epj(), sl_epb() bool fp_equald() errchk stropen begin call smark (sp) call salloc (buf, SZ_LINE, TY_CHAR) call salloc (proj, SZ_FNAME, TY_CHAR) call salloc (radecstr, SZ_FNAME, TY_CHAR) # Open a string on a file. fd = stropen (Memc[buf], SZ_LINE, WRITE_ONLY) # Get the data pointer. img = C_DATA(cp) # Get the coordinate transform descriptor. co = IMG_CO(img) radecsys = sk_stati (co, S_RADECSYS) ctype = sk_stati (co, S_CTYPE) wtype = sk_stati (co, S_WTYPE) mw = IMG_MW(img) if (mw != NULL) { # Now get the mwcs Rterm (CRPIXi), Wterm (CRVALi), and CD matrix. ndim = mw_stati (mw, MW_NPHYSDIM) call wcs_gfterm (mw, r, w, cd, ndim) crpix1 = r[1] crpix2 = r[2] crval1 = w[1] crval2 = w[2] xscale = sqrt (cd[1,1]**2 + cd[2,1]**2) * 3600.0d0 yscale = sqrt (cd[1,2]**2 + cd[2,2]**2) * 3600.0d0 xrot = 0.0 yrot = 0.0 if (!fp_equald (cd[1,1], 0.0d0)) xrot = DRADTODEG(atan ( cd[2,1] / cd[1,1])) if (!fp_equald (cd[2,2], 0.0d0)) yrot = DRADTODEG(atan (-cd[1,2] / cd[2,2])) } else { ndim = 2 xscale = 1.0 yscale = 1.0 xrot = 0.0 yrot = 0.0 } if (IMG_DEBUG) { call printf("WCS Info:\n=========\n") call printf("R term: %g %g\n"); call pargd(r[1]); call pargd(r[2]) call printf("W term: %g %g\n"); call pargd(w[1]); call pargd(w[2]) call printf(" cd: %g %g\n %g %g\n") call pargd(cd[1,1]); call pargd(cd[1,2]) call pargd(cd[2,1]); call pargd(cd[2,2]) call printf(" scale: %g %g\n");call pargd(xscale);call pargd(yscale) call printf(" rot: %g %g\n");call pargd(xrot);call pargd(yrot) } IMG_SCALE(img) = (xscale + yscale) / 2.0d0 #IMG_ROT(img) = (xrot + yrot) / 2.0d0 IMG_ROT(img) = xrot # Now format a WCS text panel such as # # Projection: TAN System: Equatorial FK5 # Ra/Dec axes: 1/2 Dimensions: 512 x 512 # # Center Pos: RA: 13:29:52.856 Dec: +47:11:40.39 # Reference Pos: RA: 13:29:52.856 Dec: +47:11:40.39 # Ref pixel coord: X: 250.256 Y: 266.309 # Plate Scale: 0.765194 Rot Angle: 1.02939 # Equinox: J2000.000 Epoch: J1987.25775240 # MJD: 46890.39406 # Get some preliminary values. if (idxstr (radecsys, Memc[radecstr], SZ_FNAME, EQTYPE_LIST) <= 0) call strcpy ("FK5", Memc[radecstr], SZ_FNAME) call strupr (Memc[radecstr]) if (idxstr (wtype, Memc[proj], SZ_FNAME, WTYPE_LIST) <= 0) call strcpy ("logical", Memc[proj], SZ_FNAME) call strupr (Memc[proj]) call fprintf (fd, "wcsinfo {\n") call fprintf (fd, " Projection: %-6s\t System: %s %s\n") call pargstr (Memc[proj]) switch (ctype) { case CTYPE_EQUATORIAL: call pargstr ("Equatorial") call pargstr (Memc[radecstr]) case CTYPE_ECLIPTIC: call pargstr ("Ecliptic") call pargstr ("") case CTYPE_GALACTIC: call pargstr ("Galactic") call pargstr ("") case CTYPE_SUPERGALACTIC: call pargstr ("SuperGalactic") call pargstr ("") default: call pargstr ("Linear") call pargstr ("") } call fprintf (fd, " Ra/Dec axes: %d/%d") call pargi (sk_stati (co, S_PLNGAX)) call pargi (sk_stati (co, S_PLATAX)) call fprintf (fd, " Dimensions: %d x %d\n\n") call pargi (IM_LEN(im,1)) call pargi (IM_LEN(im,2)) call fprintf (fd, " Center Pos: %3s: %-12H %3s: %-12h\n") if (ctype == CTYPE_EQUATORIAL) call pargstr (" RA") else call pargstr ("Lon") call pargd (cval1) if (ctype == CTYPE_EQUATORIAL) call pargstr ("Dec") else call pargstr ("Lat") call pargd (cval2) call fprintf (fd, " Reference Pos: %3s: %-12H %3s: %-12h\n") if (ctype == CTYPE_EQUATORIAL) call pargstr (" RA") else call pargstr ("Lon") call pargd (crval1) if (ctype == CTYPE_EQUATORIAL) call pargstr ("Dec") else call pargstr ("Lat") call pargd (crval2) call fprintf (fd, " Reference Pixel: X: %-9.4f Y: %-9.4f\n") call pargd (crpix1) call pargd (crpix2) call fprintf (fd, " Plate Scale: %-8f Rot Angle: %-8f\n") call pargr (IMG_SCALE(img)) call pargr (IMG_ROT(img)) call fprintf (fd, " Equinox: %s%8f Epoch: %s%.6f\n") switch (radecsys) { case EQTYPE_FK5, EQTYPE_ICRS: call pargstr ("J") ; call pargd (sk_statd(co,S_EQUINOX)) call pargstr ("J") ; call pargd (sl_epj(sk_statd(co,S_EPOCH))) default: if (IMG_LINEAR(img) == YES) { call pargstr (" ") ; call pargd (INDEFD) call pargstr (" ") ; call pargd (INDEFD) } else { call pargstr ("B") call pargd (sk_statd(co,S_EQUINOX)) call pargstr ("B") call pargd (sl_epb(sk_statd(co,S_EPOCH))) } } call fprintf (fd, " MJD: %.6f\n") call pargd (sk_statd(co,S_EPOCH)) call fprintf (fd, "}\n \n \n") # Close the formatted string and send the message. call close (fd) call wcspix_message (Memc[buf]) call sfree (sp) end # IMG_SEND_PIXTAB -- Send a 'pixtab' message. Format of the message is # # pixtab { # { {pix} {pix} ... } # pixel table values # { {x1} {x2} ... } # column label values # { {y1} {y2} ... } # row label values # { } # pixtab statistics # } # procedure img_send_pixtab (pixtab, size, x1, x2, y1, y2) real pixtab[ARB] #i pixtab array int size #i pixtab size int x1, x2, y1, y2 #i raster boundaries pointer sp, buf, el int i, j, npix real pix, sum, sum2, mean, var, stdev, x, y define SZ_PIXTAB (6*SZ_LINE) begin call smark (sp) call salloc (buf, SZ_PIXTAB, TY_CHAR) call salloc (el, SZ_FNAME, TY_CHAR) # Begin the pixtab message. call strcpy ("pixtab {\n{\ntable {\n", Memc[buf], SZ_PIXTAB) # Format the pixels into a table for presentation. Do the y-flip # here so the pixels are in order for the List widget in the GUI. # Accumulate the pixel statistics so we don't have to do it in the # GUI where it's slower. sum = 0.0 sum2 = 0.0 npix = size * size for (i=size - 1; i >= 0; i=i-1) { for (j=1; j <= size; j=j+1) { pix = pixtab[(i * size) + j] sum = sum + pix sum2 = sum2 + (pix * pix) call sprintf (Memc[el], SZ_FNAME, " {%10.1f}") call pargr (pix) call strcat (Memc[el], Memc[buf], SZ_PIXTAB) } call strcat ("\n", Memc[buf], SZ_PIXTAB) } call strcat ("}\n}\n", Memc[buf], SZ_PIXTAB) # Do the row and column label parts of the message. call strcat ("{", Memc[buf], SZ_PIXTAB) for (x = x1; x <= x2; x = x + 1.) { call sprintf (Memc[el], SZ_FNAME, " {%10.1f}") call pargr (x) call strcat (Memc[el], Memc[buf], SZ_PIXTAB) } call strcat ("}\n", Memc[buf], SZ_PIXTAB) call strcat ("{", Memc[buf], SZ_PIXTAB) for (y = y2; y >= y1; y = y - 1.) { call sprintf (Memc[el], SZ_FNAME, " {%10.1f}") call pargr (y) call strcat (Memc[el], Memc[buf], SZ_PIXTAB) } call strcat ("}\n", Memc[buf], SZ_PIXTAB) # Compute the statistics for the raster. mean = sum / real(npix) var = (sum2 - sum * mean) / real(npix - 1) if (var <= 0) stdev = 0.0 else stdev = sqrt (var) call sprintf (Memc[el], SZ_FNAME, " { %10.2f %10.4f }\n") call pargr (mean) call pargr (stdev) call strcat (Memc[el], Memc[buf], SZ_PIXTAB) # Close the message. call strcat ("}", Memc[buf], SZ_PIXTAB) # Send the formatted message. call wcspix_message (Memc[buf]) call sfree (sp) end # IMG_AMP_WCS -- Create a WCS transformation for the amplifier coordinates. pointer procedure img_amp_wcs (im, mw) pointer im #i image pointer pointer mw #i MWCS descriptor pointer ct double r[IM_MAXDIM], w[IM_MAXDIM], cd[IM_MAXDIM,IM_MAXDIM] double imgetd() pointer mw_sctran() begin r[1] = 0.0d0 r[2] = 0.0d0 w[1] = imgetd (im, "ATV1") w[2] = imgetd (im, "ATV2") cd[1,1] = imgetd (im, "ATM1_1") cd[1,2] = 0.0d0 cd[2,1] = 0.0d0 cd[2,2] = imgetd (im, "ATM2_2") # Create a new named system. call mw_newsystem (mw, "amplifier", 2) # Set the new Wterm for the system. call mw_swtermd (mw, r, w, cd, 2) # Set up the transform. ct = mw_sctran (mw, "logical", "amplifier", 03B) # Reset the default world system. call mw_sdefwcs (mw) return (ct) end # IMG_DET_WCS -- Create a WCS transformation for the detector coordinates. pointer procedure img_det_wcs (im, mw) pointer im #i image pointer pointer mw #i MWCS descriptor pointer ct double r[IM_MAXDIM], w[IM_MAXDIM], cd[IM_MAXDIM,IM_MAXDIM] double imgetd() pointer mw_sctran() begin r[1] = 0.0d0 r[2] = 0.0d0 w[1] = imgetd (im, "DTV1") w[2] = imgetd (im, "DTV2") cd[1,1] = imgetd (im, "DTM1_1") cd[1,2] = 0.0d0 cd[2,1] = 0.0d0 cd[2,2] = imgetd (im, "DTM2_2") # Create a new named system. call mw_newsystem (mw, "detector", 2) # Set the new Wterm for the system. call mw_swtermd (mw, r, w, cd, 2) # Set up the transform. ct = mw_sctran (mw, "logical", "detector", 03B) # Reset the default world system. call mw_sdefwcs (mw) return (ct) end # IMG_COORD_LABELS -- Get the WCS name, coord labels and format strings for # the specified object. procedure img_coord_labels (cp, line, wcsname, xunits, yunits) pointer cp #i cache pointer pointer line #i WCS output line char wcsname[ARB] #o WCS name string char xunits[ARB], yunits[ARB] #o WCS coord labels pointer img, co, wp pointer sp, proj, radecstr int strcmp(), sk_stati(), idxstr() begin img = C_DATA(cp) # initialize ptrs co = IMG_CO(img) wp = IMG_WP(img) if (SYSTEMS(wp,line) == SYS_WORLD) { switch (sk_stati(co,S_CTYPE)) { case CTYPE_EQUATORIAL: call strcpy (" RA", xunits, LEN_WCSNAME) call strcpy (" Dec", yunits, LEN_WCSNAME) case CTYPE_ECLIPTIC: call strcpy ("ELon", xunits, LEN_WCSNAME) call strcpy ("ELat", yunits, LEN_WCSNAME) case CTYPE_GALACTIC: call strcpy ("GLon", xunits, LEN_WCSNAME) call strcpy ("GLat", yunits, LEN_WCSNAME) case CTYPE_SUPERGALACTIC: call strcpy ("SLon", xunits, LEN_WCSNAME) call strcpy ("SLat", yunits, LEN_WCSNAME) } } else if (SYSTEMS(wp,line) == SYS_SKY) { call strcpy (WCSNAME(wp,line), wcsname, LEN_WCSNAME) call strlwr (wcsname) if (strcmp (wcsname,"ecliptic") == 0) { call strcpy ("ELon", xunits, LEN_WCSNAME) call strcpy ("ELat", yunits, LEN_WCSNAME) } else if (strcmp (wcsname,"galactic") == 0) { call strcpy ("GLon", xunits, LEN_WCSNAME) call strcpy ("GLat", yunits, LEN_WCSNAME) } else if (strcmp (wcsname,"supergalactic") == 0) { call strcpy ("SLon", xunits, LEN_WCSNAME) call strcpy ("SLat", yunits, LEN_WCSNAME) } else { call strcpy (" RA", xunits, LEN_WCSNAME) call strcpy (" Dec", yunits, LEN_WCSNAME) } } else { call strcpy ("X", xunits, LEN_WCSNAME) call strcpy ("Y", yunits, LEN_WCSNAME) } # Now get the format strings. For systems other than the image # default just use the WCS string as the name, otherwise format a # string giving more information about the system. if (SYSTEMS(wp,line) != SYS_WORLD) call strcpy (WCSNAME(wp,line), wcsname, LEN_WCSNAME) else { call smark (sp) call salloc (radecstr, SZ_FNAME, TY_CHAR) call salloc (proj, SZ_FNAME, TY_CHAR) call sprintf (wcsname, LEN_WCSNAME, "%s-%s-%s") switch (sk_stati(co,S_CTYPE)) { case CTYPE_EQUATORIAL: call pargstr ("EQ") case CTYPE_ECLIPTIC: call pargstr ("ECL") case CTYPE_GALACTIC: call pargstr ("GAL") case CTYPE_SUPERGALACTIC: call pargstr ("SGAL") default: call pargstr ("UNKN") } if (sk_stati(co,S_CTYPE) == CTYPE_EQUATORIAL) { if (idxstr(sk_stati(co,S_RADECSYS), Memc[radecstr], SZ_FNAME, EQTYPE_LIST) <= 0) call strcpy ("FK5", Memc[radecstr], SZ_FNAME) call strupr (Memc[radecstr]) call pargstr (Memc[radecstr]) } else { if (sk_stati(co,S_CTYPE) == CTYPE_SUPERGALACTIC) call pargstr ("-") else call pargstr ("--") } if (idxstr(sk_stati(co,S_WTYPE), Memc[proj], SZ_FNAME, WTYPE_LIST) <= 0) call strcpy ("linear", Memc[proj], SZ_FNAME) call strupr (Memc[proj]) call pargstr (Memc[proj]) call sfree (sp) } # Now fix up the WCS system name. if (strcmp (wcsname, "fk4") == 0 || strcmp (wcsname, "fk5") == 0 || strcmp (wcsname, "icrs") == 0 || strcmp (wcsname, "gappt") == 0 || strcmp (wcsname, "fk4-no-e") == 0) { call strupr (wcsname) } else if (IS_LOWER(wcsname[1])) wcsname[1] = TO_UPPER(wcsname[1]) end # IMG_COORD_FMT -- Format the coordinate strings. procedure img_coord_fmt (cp, line, xval, yval, xc, yc) pointer cp #i object cache pointer int line #i output line number double xval, yval #i input coords char xc[ARB], yc[ARB] #o formatted coord strings pointer img, co, wp char xfmt[LEN_WCSNAME], yfmt[LEN_WCSNAME] int sk_stati() bool streq() begin img = C_DATA(cp) # initialize ptrs co = IMG_CO(img) wp = IMG_WP(img) # Convert coords to the requested format. if (FORMATS(wp,line) == FMT_DEFAULT) { if (IMG_MW(img) == NULL) { call strcpy ("%10.2f", xfmt, LEN_WCSNAME) call strcpy ("%10.2f", yfmt, LEN_WCSNAME) } else { if (SYSTEMS(wp,line) == SYS_WORLD || SYSTEMS(wp,line) == SYS_SKY) { if (streq(WCSNAME(wp,line),"ecliptic") || streq(WCSNAME(wp,line),"galactic") || streq(WCSNAME(wp,line),"supergalactic")) call strcpy ("%h", xfmt, LEN_WCSNAME) else call strcpy ("%.2H", xfmt, LEN_WCSNAME) call strcpy ("%.1h", yfmt, LEN_WCSNAME) } else { call strcpy ("%10.2f", xfmt, LEN_WCSNAME) call strcpy ("%10.2f", yfmt, LEN_WCSNAME) } } } else if (FORMATS(wp,line) == FMT_HMS) { if (sk_stati(co, S_CTYPE) == CTYPE_EQUATORIAL) call strcpy ("%.2H", xfmt, LEN_WCSNAME) else call strcpy ("%.1h", xfmt, LEN_WCSNAME) call strcpy ("%h", yfmt, LEN_WCSNAME) } else { call strcpy ("%10.2f", xfmt, LEN_WCSNAME) call strcpy ("%10.2f", yfmt, LEN_WCSNAME) } # Convert the value to the requested format call sprintf (xc, LEN_WCSNAME, xfmt) if (FORMATS(wp,line) != FMT_RAD) call pargd (xval) else call pargd (DEGTORAD(xval)) call sprintf (yc, LEN_WCSNAME, yfmt) if (FORMATS(wp,line) != FMT_RAD) call pargd (yval) else call pargd (DEGTORAD(yval)) end # IMG_GET_COORD -- Given an x,y position in the image return the coordinate in # the given system. procedure img_get_coord (img, x, y, system, wcsname, wx, wy) pointer img #i IMG struct pointer double x, y #i input image position int system #i coordinate system requested char wcsname[ARB] #i desired WCS name double wx, wy #o output coordinates double ox, oy real epoch pointer im, co, nco char buf[SZ_LINE] int stat real imgetr() int imaccf(), sk_stati(), sk_decwstr() bool streq() errchk imgetr begin im = IMG_IM(img) co = IMG_CO(img) wx = x # fallback values wy = y switch (system) { case SYS_NONE: wx = x wy = y case SYS_PHYSICAL: if (IMG_CTP(img) != NULL) call mw_c2trand (IMG_CTP(img), x, y, wx, wy) case SYS_WORLD: if (IMG_CTW(img) != NULL) call mw_c2trand (IMG_CTW(img), x, y, wx, wy) case SYS_AMP: if (IMG_CTA(img) != NULL) call mw_c2trand (IMG_CTA(img), x, y, wx, wy) case SYS_CCD: ; # TBD case SYS_DETECTOR: if (IMG_CTD(img) != NULL) call mw_c2trand (IMG_CTD(img), x, y, wx, wy) case SYS_SKY: # Note Ecliptic/GAPPT coords need an epoch value. if (streq (wcsname, "ecliptic") || streq (wcsname, "gappt")) { if (imaccf (im, "EPOCH") == YES) { epoch = imgetr (im, "EPOCH") if (epoch == 0.0 || IS_INDEFR(epoch)) epoch = 1950.0 } else epoch = 1950.0 call sprintf (buf, SZ_LINE, "%s %.1f") call pargstr (wcsname) call pargr (epoch) } else call strcpy (wcsname, buf, SZ_LINE) stat = sk_decwstr (buf, nco, co) if (stat != ERR) { if (IMG_CTW(img) != NULL) call mw_c2trand (IMG_CTW(img), x, y, ox, oy) call sk_lltran (co, nco, DEGTORAD(ox), DEGTORAD(oy), INDEFD, INDEFD, 0.0d0, 0.0d0, wx, wy) if (sk_stati(co,S_PLATAX) < sk_stati(co,S_PLNGAX)) { wx = RADTODEG(wy) # transposed image wy = RADTODEG(wx) } else { wx = RADTODEG(wx) # regular image wy = RADTODEG(wy) } } case SYS_OTHER: ; # TBD default: # default coords wx = x wy = y } end