aboutsummaryrefslogtreecommitdiff
path: root/pkg/plot/t_pradprof.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/plot/t_pradprof.x')
-rw-r--r--pkg/plot/t_pradprof.x548
1 files changed, 548 insertions, 0 deletions
diff --git a/pkg/plot/t_pradprof.x b/pkg/plot/t_pradprof.x
new file mode 100644
index 00000000..dbf7aae1
--- /dev/null
+++ b/pkg/plot/t_pradprof.x
@@ -0,0 +1,548 @@
+include <imhdr.h>
+include <gset.h>
+include <math.h>
+
+# T_ PRADPROF -- Compute a radial profile using user specified coordinates
+# and plot or list the result.
+
+procedure t_pradprof()
+
+int images # the list of images
+real xinit, yinit # the initial guess for the profile center
+real pradius # the plotting radius
+real paz1, paz2 # azimuth limits
+bool center # center the object before computing profile
+int cboxsize # the centering box width
+bool list # list output instead of plot output
+
+int rboxsize, npts
+pointer sp, imname, im, radius, azimuth, intensity
+real xcntr, ycntr
+
+bool clgetb()
+int imtopenp(), imtgetim(), clgeti(), rp_radius()
+pointer immap()
+real clgetr()
+
+begin
+ # Allocate stack space.
+ call smark (sp)
+ call salloc (imname, SZ_FNAME, TY_CHAR)
+
+ # Get the radial profiling parameters. The width of the extraction
+ # box mut be odd.
+ images = imtopenp ("input")
+ xinit = clgetr ("xinit")
+ yinit = clgetr ("yinit")
+ pradius = clgetr ("radius")
+ paz1 = clgetr ("az1")
+ paz2 = clgetr ("az2")
+ rboxsize = 2 * (nint (pradius + 1.0)) + 1
+
+ # Get the centering parameters. The centering box must be odd.
+ center = clgetb ("center")
+ if (center) {
+ cboxsize = clgeti ("cboxsize")
+ if (mod (cboxsize, 2) == 0)
+ cboxsize = cboxsize + 1
+ }
+
+ # List the radial profile instead of plotting it?
+ list = clgetb ("list")
+
+ # Allocate memory for vectors.
+ call malloc (radius, rboxsize * rboxsize, TY_REAL)
+ call malloc (azimuth, rboxsize * rboxsize, TY_REAL)
+ call malloc (intensity, rboxsize * rboxsize, TY_REAL)
+
+ # Loop over all images
+ while (imtgetim (images, Memc[imname], SZ_FNAME) != EOF) {
+
+ # Open the image.
+ iferr (im = immap (Memc[imname], READ_ONLY, 0)) {
+ call eprintf ("Image %s not found\n")
+ call pargstr (Memc[imname])
+ next
+ }
+
+ # Find the star center, if center=yes.
+ if (center)
+ call rp_cntr (im, xinit, yinit, cboxsize, xcntr, ycntr)
+ else {
+ xcntr = xinit
+ ycntr = yinit
+ }
+
+ # Get the radius and intensity vectors.
+ npts = rp_radius (im, xcntr, ycntr, rboxsize, pradius, paz1, paz2,
+ Memr[radius], Memr[azimuth], Memr[intensity])
+
+ # Make list of the radial profile if list=yes or plot if list=no.
+ if (list)
+ call rp_rlist (Memc[imname], xcntr, ycntr, paz1, paz2,
+ Memr[radius], Memr[azimuth], Memr[intensity], npts)
+ else
+ call rp_rplot (Memc[imname], xcntr, ycntr, paz1, paz2,
+ Memr[radius], Memr[azimuth], Memr[intensity], npts)
+
+ call imunmap (im)
+ }
+
+ call mfree (radius, TY_REAL)
+ call mfree (intensity, TY_REAL)
+ call imtclose (images)
+ call sfree (sp)
+end
+
+
+# RP_CNTR -- Compute the star center using a simple 1D centroiding algorithm
+# on the x and y marginals, after thresholding at the mean.
+
+procedure rp_cntr (im, xstart, ystart, boxsize, xcntr, ycntr)
+
+pointer im # pointer to the input image
+real xstart, ystart # starting coordinates
+int boxsize # width of the centering box
+real xcntr, ycntr # centered coordinates
+
+int half_box, x1, x2, y1, y2
+int ncols, nrows, nx, ny, try
+pointer bufptr, sp, x_vect, y_vect
+real xinit, yinit
+pointer imgs2r()
+
+begin
+ # Initialize.
+ half_box = (boxsize - 1) / 2
+ xinit = xstart
+ yinit = ystart
+ ncols = IM_LEN (im, 1)
+ nrows = IM_LEN (im, 2)
+
+ try = 0
+ repeat {
+
+ # Compute the extraction region.
+ x1 = max (xinit - half_box, 1.0) + 0.5
+ x2 = min (xinit + half_box, real (ncols)) + 0.5
+ y1 = max (yinit - half_box, 1.0) + 0.5
+ y2 = min (yinit + half_box, real (nrows)) + 0.5
+ nx = x2 - x1 + 1
+ ny = y2 - y1 + 1
+
+ # Get the data.
+ bufptr = imgs2r (im, x1, x2, y1, y2)
+
+ call smark (sp)
+ call salloc (x_vect, nx, TY_REAL)
+ call salloc (y_vect, ny, TY_REAL)
+
+ # Compute the marginals.
+ call aclrr (Memr[x_vect], nx)
+ call aclrr (Memr[y_vect], ny)
+ call rp_rowsum (Memr[bufptr], Memr[x_vect], nx, ny)
+ call rp_colsum (Memr[bufptr], Memr[y_vect], nx, ny)
+
+ # Compute the centers.
+ call rp_getcenter (Memr[x_vect], nx, xcntr)
+ call rp_getcenter (Memr[y_vect], ny, ycntr)
+
+ # Add in offsets to image coordinate system.
+ xcntr = xcntr + x1
+ ycntr = ycntr + y1
+
+ call sfree (sp)
+
+ # If the shifts are greater than a pixel to 1 more iteration.
+ try = try + 1
+ if (try == 1) {
+ if ((abs (xcntr-xinit) > 1.0) || (abs (ycntr-yinit) > 1.0)) {
+ xinit = xcntr
+ yinit = ycntr
+ }
+ } else
+ break
+ }
+end
+
+
+# RP_RADIUS -- Get the data and compute the radius and intensity vectors.
+
+int procedure rp_radius (im, xcntr, ycntr, rboxsize, pradius, paz1, paz2,
+ radius, azimuth, intensity)
+
+pointer im # pointer to the input image
+real xcntr, ycntr # the center of the extraction box
+int rboxsize # the width of the extraction box
+real pradius # the plotting radius
+real paz1, paz2 # the azimuth limits
+real radius[ARB] # the output radius vector
+real azimuth[ARB] # the output azimuth vector
+real intensity[ARB] # the output intensity vector
+
+int half_box, ncols, nrows, x1, x2, y1, y2, nx, ny, npts
+pointer bufptr
+real xinit, yinit
+int rp_vectors()
+pointer imgs2r()
+
+begin
+ # Initialize.
+ half_box = (rboxsize - 1) / 2
+ xinit = xcntr
+ yinit = ycntr
+ ncols = IM_LEN(im,1)
+ nrows = IM_LEN(im,2)
+
+ # Get the data.
+ x1 = max (xinit - half_box, 1.0) + 0.5
+ x2 = min (xinit + half_box, real (ncols)) + 0.5
+ y1 = max (yinit - half_box, 1.0) + 0.5
+ y2 = min (yinit + half_box, real (nrows)) + 0.5
+ nx = x2 - x1 + 1
+ ny = y2 - y1 + 1
+ bufptr = imgs2r (im, x1, x2, y1, y2)
+
+ # Compute the radius and intensity vectors.
+ npts = rp_vectors (Memr[bufptr], nx, ny, x1, y1, xcntr, ycntr,
+ pradius, paz1, paz2, radius, azimuth, intensity)
+
+ return (npts)
+end
+
+
+# RP_RLIST -- Print the intensity as a function of radial distance on the
+# standard output.
+
+procedure rp_rlist (imname, xcntr, ycntr, paz1, paz2, radius, azimuth,
+ intensity, npts)
+
+char imname[ARB] # the name of the input image
+real xcntr, ycntr # the center of the radial profile
+real paz1, paz2 # the azimuth limits
+real radius[npts] # the radius vector
+real azimuth[npts] # the azimuth vector
+real intensity[npts] # the intensity vector
+int npts # the number of points
+
+int i
+
+begin
+ call printf ("# [%s] xcntr:%7.2f ycntr:%7.2f\n")
+ call pargstr (imname)
+ call pargr (xcntr)
+ call pargr (ycntr)
+ call printf ("# az1:%7.2f az2:%7.2f\n")
+ call pargr (min(paz1,paz2))
+ call pargr (max(paz1,paz2))
+
+ do i = 1, npts {
+ call printf ("%7.2f %g\n")
+ call pargr (radius[i])
+ call pargr (intensity[i])
+ }
+end
+
+
+# RP_RPLOT -- Plot intensity as a function of radial distance.
+
+procedure rp_rplot (imname, xcntr, ycntr, az1, az2,
+ radius, azimuth, intensity, npts)
+
+char imname[ARB]
+int npts
+real xcntr, ycntr, az1, az2
+real radius[npts], azimuth[npts], intensity[npts]
+
+char device[SZ_LINE]
+int mode
+pointer gp
+
+bool clgetb()
+pointer gopen()
+
+begin
+ call clgstr ("graphics", device, SZ_LINE)
+
+ if (clgetb("append"))
+ mode = APPEND
+ else
+ mode = NEW_FILE
+
+ gp = gopen (device, mode, STDGRAPH)
+ call rp_graph (gp, imname, xcntr, ycntr, az1, az2,
+ mode, radius, intensity, npts)
+ call gclose (gp)
+end
+
+
+# RP_VECTORS -- Compute the radius and intensity vectors.
+
+int procedure rp_vectors (a, nx, ny, x1, y1, xcntr, ycntr, pradius, paz1, paz2,
+ radius, azimuth, intensity)
+
+real a[nx,ny] # the input data array
+int nx, ny # dimensions of the input array
+int x1, y1 # lower left corner of input array
+real xcntr, ycntr # coordinates of center pixel
+real pradius # the plotting radius
+real paz1, paz2 # the azimuth limits
+real radius[ARB] # the output radius vector
+real azimuth[ARB] # the output azimuth vector
+real intensity[ARB] # the output intensity vector
+
+int i, j, npts
+real dx, dy, az, pr2, az1, az2, r2, dy2
+
+begin
+ az1 = DEGTORAD (min (paz1, paz2))
+ az2 = DEGTORAD (max (paz1, paz2))
+ while (az1 < 0.) {
+ az1 = az1 + TWOPI
+ az2 = az2 + TWOPI
+ }
+ while (az1 > TWOPI) {
+ az1 = az1 - TWOPI
+ az2 = az2 - TWOPI
+ }
+ pr2 = pradius * pradius
+
+ npts = 0
+ do i = 1, ny {
+ dy = (ycntr - y1 + 1 - i)
+ dy2 = dy ** 2
+ do j = 1, nx {
+ dx = (xcntr - x1 + 1 - j)
+ r2 = dx ** 2 + dy2
+ if (r2 > pr2)
+ next
+ az = atan2 (dy, dx)
+ if (az < 0.)
+ az = az + TWOPI
+ if (az < az1 || az > az2)
+ next
+ npts = npts + 1
+ radius[npts] = sqrt (r2)
+ azimuth[npts] = RADTODEG (az)
+ intensity[npts] = a[j,i]
+ }
+ }
+
+ return (npts)
+end
+
+
+# RP_ROWSUM -- Sum all the rows in a raster.
+
+procedure rp_rowsum (v, row, nx, ny)
+
+real v[nx,ny] # the input subraster
+real row[ARB] # the output summed row
+int nx, ny # the dimensions of the input subraster
+
+int i, j
+
+begin
+ do i = 1, ny
+ do j = 1, nx
+ row[j] = row[j] + v[j,i]
+end
+
+
+# RP_COLSUM -- Sum all the columns in a raster.
+
+procedure rp_colsum (v, col, nx, ny)
+
+real v[nx,ny] # the input subraster
+real col[ARB] # the output summed column
+int nx, ny # the dimensions of the input subraster
+
+int i, j
+
+begin
+ do i = 1, ny
+ do j = 1, nx
+ col[j] = col[j] + v[i,j]
+end
+
+
+# RP_GETCENTER -- Compute the centroid of an array.
+
+procedure rp_getcenter (v, nv, vc)
+
+real v[ARB] # the input array
+int nv # length of the input array
+real vc # the output centroid
+
+int i
+real sum1, sum2, sigma, cont
+
+begin
+ sum1 = 0.0
+ sum2 = 0.0
+
+ call aavgr (v, nv, cont, sigma)
+ do i = 1, nv
+ if (v[i] > cont) {
+ sum1 = sum1 + (i-1) * (v[i] - cont)
+ sum2 = sum2 + (v[i] - cont)
+ }
+
+ vc = sum1 / sum2
+end
+
+
+define MTYPES "|point|box|plus|cross|circle|hebar|vebar|hline|vline|diamond|"
+define RP_GBUF 0.10
+define RP_SZTITLE 512
+define DEF_IMTITLE "imtitle"
+
+# RP_GRAPH -- Graph the radial profile.
+
+procedure rp_graph (gp, imname, xcntr, ycntr, az1, az2, mode, x, y, npts)
+
+pointer gp # GIO pointer
+char imname[ARB] # image name
+real xcntr # starting x coordinate
+real ycntr # starting y coordinate
+real az1, az2 # azimuth limits
+int mode # Mode
+real x[npts] # X data
+real y[npts] # Y data
+int npts # Number of points
+
+int i, marks[10], linepattern, patterns[4], clgeti(), btoi(), strdic()
+pointer sp, marker, title, xlabel, ylabel
+real x1, x2, y1, y2, wx1, wx2, wy1, wy2, vx1, vx2, vy1,vy2, temp,
+ szmarker, clgetr()
+bool clgetb(), streq()
+
+data patterns/GL_SOLID, GL_DASHED, GL_DOTTED, GL_DOTDASH/
+data marks/GM_POINT, GM_BOX, GM_PLUS, GM_CROSS, GM_CIRCLE, GM_HEBAR,
+ GM_VEBAR, GM_HLINE, GM_VLINE, GM_DIAMOND/
+
+begin
+ call smark (sp)
+ call salloc (marker, SZ_LINE, TY_CHAR)
+
+ # If a new graph setup all the axes and labeling options and then
+ # make the graph.
+
+ if (mode == NEW_FILE) {
+
+ call gclear (gp)
+
+ linepattern = 0
+
+ x1 = clgetr ("wx1")
+ x2 = clgetr ("wx2")
+ y1 = clgetr ("wy1")
+ y2 = clgetr ("wy2")
+
+ if (IS_INDEF (x1) || IS_INDEF (x2))
+ call gascale (gp, x, npts, 1)
+ if (IS_INDEF (y1) || IS_INDEF (y2))
+ call gascale (gp, y, npts, 2)
+
+ call gswind (gp, x1, x2, y1, y2)
+ call ggwind (gp, wx1, wx2, wy1, wy2)
+
+ temp = wx2 - wx1
+ if (IS_INDEF (x1))
+ wx1 = wx1 - RP_GBUF * temp
+ if (IS_INDEF (x2))
+ wx2 = wx2 + RP_GBUF * temp
+
+ temp = wy2 - wy1
+ if (IS_INDEF (y1))
+ wy1 = wy1 - RP_GBUF * temp
+ if (IS_INDEF (y2))
+ wy2 = wy2 + RP_GBUF * temp
+
+ call gswind (gp, wx1, wx2, wy1, wy2)
+ call gsetr (gp, G_ASPECT, 0.)
+ call gseti (gp, G_ROUND, btoi (clgetb ("round")))
+
+ if (clgetb("fill"))
+ call gsetr (gp, G_ASPECT, 0.0)
+ else
+ call gsetr (gp, G_ASPECT, 1.0)
+
+ i = GW_LINEAR
+ if (clgetb ("logx"))
+ i = GW_LOG
+ call gseti (gp, G_XTRAN, i)
+ i = GW_LINEAR
+ if (clgetb ("logy"))
+ i = GW_LOG
+ call gseti (gp, G_YTRAN, i)
+
+ # Set the view port
+ vx1 = clgetr ("vx1")
+ vx2 = clgetr ("vx2")
+ vy1 = clgetr ("vy1")
+ vy2 = clgetr ("vy2")
+ call gsview (gp, vx1, vx2, vy1, vy2)
+
+ if (clgetb ("box")) {
+
+ # Get number of major and minor tick marks.
+ call gseti (gp, G_XNMAJOR, clgeti ("majrx"))
+ call gseti (gp, G_XNMINOR, clgeti ("minrx"))
+ call gseti (gp, G_YNMAJOR, clgeti ("majry"))
+ call gseti (gp, G_YNMINOR, clgeti ("minry"))
+
+ # Label tick marks on axes?
+ call gseti (gp, G_LABELTICKS,
+ btoi (clgetb ("ticklabels")))
+
+ # Fetch labels and plot title string.
+ call salloc (title, RP_SZTITLE, TY_CHAR)
+ call salloc (xlabel, SZ_LINE, TY_CHAR)
+ call salloc (ylabel, SZ_LINE, TY_CHAR)
+
+ # Build system info string
+ call sysid (Memc[title], SZ_LINE)
+ call strcat ("\n", Memc[title], RP_SZTITLE)
+
+ # Build the title string
+ call clgstr ("title", Memc[marker], SZ_LINE)
+ if (streq (Memc[marker], DEF_IMTITLE)) {
+ call sprintf (Memc[marker], SZ_LINE,
+ "Radial Plot of %s at [%0.2f,%0.2f] az=[%.1f,%.1f]")
+ call pargstr (imname)
+ call pargr (xcntr)
+ call pargr (ycntr)
+ call pargr (az1)
+ call pargr (az2)
+ }
+ call strcat (Memc[marker], Memc[title], RP_SZTITLE)
+
+ call clgstr ("xlabel", Memc[xlabel], SZ_LINE)
+ call clgstr ("ylabel", Memc[ylabel], SZ_LINE)
+
+ call glabax (gp, Memc[title], Memc[xlabel], Memc[ylabel])
+ }
+ }
+
+ # Draw the data.
+ if (clgetb ("pointmode")) {
+ call clgstr ("marker", Memc[marker], SZ_LINE)
+ i = strdic (Memc[marker], Memc[marker], SZ_LINE, MTYPES)
+ if (i == 0)
+ i = 2
+ if (marks[i] == GM_POINT)
+ szmarker = 0.0
+ else
+ szmarker = clgetr ("szmarker")
+ call gpmark (gp, x, y, npts, marks[i], szmarker, szmarker)
+ }
+ else {
+ linepattern = min (4, linepattern + 1)
+ call gseti (gp, G_PLTYPE, patterns[linepattern])
+ call gpline (gp, x, y, npts)
+ }
+ call gflush (gp)
+
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 05:11:38 -0400