include include define X0 96 # input x dimension define Y0 64 # input y dimension define X1 51 # trimmed x dimension define Y1 63 # trimmed y dimension define X2 48 # x restored dimension define Y2 48 # y restored dimension define PIX1 $1[$2 + ($3 - 1)] define PIX2 $1[$2 + ($3 - 1) + ($4 - 1) * $5] # T_TAREST -- Restore 96 x 64 FOS Mode I image # # Steve Hulbert, Feb90: Original # Howard Bushouse, Mar94: Only write files for labels if requested # Howard Bushouse, Apr94: Added MWCS update for output images procedure t_tarest () char input[SZ_FNAME] # input image int trim # number of columns to trim off left # side of input image char mag_y[SZ_FNAME] # magic matrix image for y restoration char out_y[SZ_FNAME] # output image (y-restored) char mag_x[SZ_FNAME] # magic matrix image for x restoration char out_x[SZ_FNAME] # output image (x- and y-restored) char lab1[SZ_FNAME] # file containing labels for frame 1 char lab2[SZ_FNAME] # file containing labels for frame 2 char lab3[SZ_FNAME] # file containing labels for frame 3 pointer im_in, im_mag_y, im_out_y, im_mag_x, im_out_x pointer buf_trim, buf_mag_y, buf_out_y, buf_mag_x, buf_out_x int clgeti() int immap(), sp, imgl2r() int groups, i, istat, imgl1r(), imgs2r(), imps2r(), imgeti() real datamax, datamin int x, y pointer mw real shifts[IM_MAXDIM] bool envgetb() pointer mw_openim() begin # get input from the cl call clgstr ("input", input, SZ_FNAME) call clgstr ("magic_y", mag_y, SZ_FNAME) call clgstr ("output_y", out_y, SZ_FNAME) call clgstr ("magic_x", mag_x, SZ_FNAME) call clgstr ("output_x", out_x, SZ_FNAME) trim = clgeti ("trim") call clgstr ("label1", lab1, SZ_FNAME) call clgstr ("label2", lab2, SZ_FNAME) call clgstr ("label3", lab3, SZ_FNAME) # map images im_in = immap (input, READ_ONLY, 0) im_out_y = immap (out_y, NEW_COPY, im_in) im_out_x = immap (out_x, NEW_COPY, im_in) im_mag_x = immap (mag_x, READ_ONLY, 0) im_mag_y = immap (mag_y, READ_ONLY, 0) # check for correct image size allowing for raw or calibrated data groups = imgeti (im_in, "GCOUNT") x = IM_LEN (im_in, 1) if (groups == 1) { y = IM_LEN (im_in, 2) } else { y = groups } if (x != X0 || y != Y0) call error (0,"Image does not have dimensions of 96 x 64") # allocate buffer call smark (sp) call salloc (buf_trim, X1 * Y1, TY_REAL) call salloc (buf_mag_y, Y1 * Y2, TY_REAL) call salloc (buf_mag_x, X1 * X2, TY_REAL) call salloc (buf_out_y, X1 * Y2, TY_REAL) call salloc (buf_out_x, X2 * Y2, TY_REAL) # read observation do i = 1, Y1 { if (groups != 1) { call gf_opengr (im_in, i, datamin, datamax, istat) call amovr [Memr[imgl1r(im_in) + trim], Memr[buf_trim + (i - 1) * X1], X1) } else { call amovr [Memr[imgl2r(im_in, i) + trim], Memr[buf_trim + (i - 1) * X1], X1) } } # read magic y-matrix call amovr (Memr[imgs2r(im_mag_y, 1, Y1, 1, Y2)], Memr[buf_mag_y], Y1 * Y2) # restore in y call mmult (buf_mag_y, Y1, Y2, buf_trim, X1, Y1, buf_out_y) # write y-restored image IM_NDIM (im_out_y) = 2 IM_LEN (im_out_y, 1) = X1 IM_LEN (im_out_y, 2) = Y2 call amovr (Memr[buf_out_y], Memr[imps2r(im_out_y, 1, X1, 1, Y2)], X1 * Y2) # read magic x-matrix call amovr (Memr[imgs2r(im_mag_x, 1, X1, 1, X2)], Memr[buf_mag_x], X1 * X2) # restore in x call mtran (buf_out_y, X1, Y2) call mmult (buf_mag_x, X1, X2, buf_out_y, Y2, X1, buf_out_x) call mtran (buf_out_x, Y2, X2) # write x-restored image IM_NDIM (im_out_x) = 2 IM_LEN (im_out_x, 1) = X2 IM_LEN (im_out_x, 2) = Y2 call amovr (Memr[buf_out_x], Memr[imps2r(im_out_x, 1, X2, 1, Y2)], X2 * Y2) # create overlay for x-restored image if (lab1[1] != EOS) call panger (im_in, lab1, lab2, lab3) # Update the world coordinate system. if (!envgetb ("nomwcs")) { mw = mw_openim (im_out_y) shifts[1] = -trim shifts[2] = (Y2-Y1)/2.0 call mw_shift (mw, shifts, (2 ** IM_NDIM(im_out_y) - 1)) call mw_saveim (mw, im_out_y) shifts[1] = (X2-X1)/2.0 shifts[2] = 0.0 call mw_shift (mw, shifts, (2 ** IM_NDIM(im_out_x) - 1)) call mw_saveim (mw, im_out_x) call mw_close (mw) } call sfree (sp) call imunmap (im_in) call imunmap (im_mag_y) call imunmap (im_mag_x) call imunmap (im_out_x) call imunmap (im_out_y) end # MMULT -- Multiply matrices # procedure mmult (buf1, x1, y1, buf2, x2, y2, buf3) pointer buf1, buf2 # input buffers pointer buf3 # output buffer int x1, x2 # input x axis dimensions int y1, y2 # input y axis dimensions int i, j, k real r begin if (x1 != y2) call error (0,"Bad dimensions -- cannot multiply matrices.") do i = 1, x2 { do j = 1, y1 { r = 0 do k = 1, y2 { r = r + PIX2 (Memr, buf1, k, j, x1) * PIX2 (Memr, buf2, i, k, x2) } PIX2 (Memr, buf3, i, j, x2) = r } } end # MTRAN -- Transpose a small matrix # procedure mtran (buf, nx, ny) pointer buf # output buffer int nx # output x axis dimension int ny # output y axis dimension pointer sp, buf_temp int i, j begin call smark (sp) call salloc (buf_temp, nx * ny, TY_REAL) do j = 1, ny { do i = 1, nx PIX2 (Memr, buf_temp, j, i, ny) = PIX2 (Memr, buf, i, j, nx) } call amovr (Memr[buf_temp], Memr[buf], nx * ny) call sfree (sp) end define PX1 53.0 define PY1 53.0 define R 5.0 # PANGER -- Calculate aperture orientation on sky and set up # input command files for tvmark # procedure panger (im_in, lab1, lab2, lab3) pointer im_in #input image real pangaper, x, y int fd[3], open, i real imgetr, scale[3] char lab1[SZ_FNAME],lab2[SZ_FNAME],lab3[SZ_FNAME] int signr() data scale /1.5, 1.0, 1.0/ begin pangaper = imgetr(im_in, "PA_APER") if (IS_INDEFR (pangaper)) call error (0, "Cannot read keyword: PA_APER.") pangaper = pangaper * 3.1415926 / 180.0 y = R * sin (pangaper) x = R * cos (pangaper) fd[1] = open (lab1, NEW_FILE, TEXT_FILE) fd[2] = open (lab2, NEW_FILE, TEXT_FILE) fd[3] = open (lab3, NEW_FILE, TEXT_FILE) # set up to print the static labels for file #1 call fprintf (fd[1], "96.00 -3.00 101 :text RAW IMAGE\n") # set up to print the static labels for file #2 call fprintf (fd[2], "51.00 -3.00 101 :text Y-RESTORED IMAGE\n") # set up to print the static labels for file #3 call fprintf (fd[3], "48.00 -3.00 101 :text XY-RESTORED IMAGE\n") call fprintf (fd[3], "24.5 24.5 101 v\n") call fprintf (fd[3], "30.08 24.5 101 v\n") call fprintf (fd[3], "23.94 7.2 101 b\n") call fprintf (fd[3], "25.06 8.31 101 b\n") call fprintf (fd[3], "18.92 2.17 101 b\n") call fprintf (fd[3], "30.08 13.34 101 b\n") call fprintf (fd[3], "26.83 2.67 101 :text LOWER\n") call fprintf (fd[3], "23.94 40.69 101 b\n") call fprintf (fd[3], "25.06 41.80 101 b\n") call fprintf (fd[3], "18.92 46.83 101 b\n") call fprintf (fd[3], "30.08 35.66 101 b\n") call fprintf (fd[3], "26.83 45.58 101 :text UPPER\n") call fprintf (fd[3], "0.0 30.08 101 s\n") call fprintf (fd[3], "-2.0 30.08 101 s\n") call fprintf (fd[3], "-1.0 30.08 101 s\n") call fprintf (fd[3], "-1.0 18.92 101 s\n") call fprintf (fd[3], "0.0 18.92 101 s\n") call fprintf (fd[3], "-2.0 18.92 101 s\n") call fprintf (fd[3], "-2.0 24.5 101 :text 1""\n") # set up the N-E axis based on pangaper do i = 1, 3 { call fprintf (fd[i], "%6.2f %6.2f 101 s\n") call pargr (scale[i] * PX1) call pargr (scale[i] * PY1) call fprintf (fd[i], "%6.2f %6.2f 101 s\n") call pargr (scale[i] * (PX1 + x)) call pargr (scale[i] * (PY1 + y)) call fprintf (fd[i], "%6.2f %6.2f 101 :text E\n") call pargr (scale[i] * (PX1 + x) + (signr(x) * 1.5)) call pargr (scale[i] * (PY1 + y) + (signr(y) * 1.5)) call fprintf (fd[i], "%6.2f %6.2f 101 s\n") call pargr (scale[i] * PX1) call pargr (scale[i] * PY1) call fprintf (fd[i], "%6.2f %6.2f 101 s\n") call pargr (scale[i] * (PX1 - y)) call pargr (scale[i] * (PY1 + x)) call fprintf (fd[i], "%6.2f %6.2f 101 :text N\n") call pargr (scale[i] * (PX1 - y) + (signr(x) * 1.5)) call pargr (scale[i] * (PY1 + x) + (signr(y) * 1.5)) } # set up x-y axes do i = 1, 3 { call fprintf (fd[i], "%6.2f %6.2f 101 s\n") call pargr (-5.0 * scale[i]) call pargr ( 0.0 * scale[i]) call fprintf (fd[i], "%6.2f %6.2f 101 s\n") call pargr (-5.0 * scale[i]) call pargr (-5.0 * scale[i]) call fprintf (fd[i], "%6.2f %6.2f 101 s\n") call pargr (-5.0 * scale[i]) call pargr (-5.0 * scale[i]) call fprintf (fd[i], "%6.2f %6.2f 101 s\n") call pargr ( 0.0 * scale[i]) call pargr (-5.0 * scale[i]) call fprintf (fd[i], "%6.2f %6.2f 101 :text Y\n") call pargr (-4.0 * scale[i]) call pargr ( 0.5 * scale[i]) call fprintf (fd[i], "%6.2f %6.2f 101 :text X\n") call pargr ( 2.0 * scale[i]) call pargr (-5.5 * scale[i]) } do i = 1, 3 call close (fd[i]) end int procedure signr (x) real x int ix begin if (x <= 0.) ix = -1 else ix = 1 return ix end