aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/ccdred/src/t_skyreplace.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/imred/ccdred/src/t_skyreplace.x')
-rw-r--r--noao/imred/ccdred/src/t_skyreplace.x301
1 files changed, 301 insertions, 0 deletions
diff --git a/noao/imred/ccdred/src/t_skyreplace.x b/noao/imred/ccdred/src/t_skyreplace.x
new file mode 100644
index 00000000..9bd2e9d0
--- /dev/null
+++ b/noao/imred/ccdred/src/t_skyreplace.x
@@ -0,0 +1,301 @@
+include <imhdr.h>
+
+
+# T_SKYREPLACE -- Replace objects by sky. This development code as is not
+# used in the package. It is here to be worked on further when an image
+# display interface is added.
+
+procedure t_skyreplace ()
+
+char image[SZ_FNAME] # Image to be modified
+
+char graph[SZ_LINE], display[SZ_LINE], cmd[SZ_LINE]
+pointer im, immap()
+int clgeti(), wcs, key, clgcur(), nrep, skyreplace()
+real wx, wy, xc, yc, r, s
+
+begin
+ call clgstr ("image", image, SZ_FNAME)
+ call sprintf (graph, SZ_LINE, "contour %s")
+ call pargstr (image)
+ call sprintf (display, SZ_LINE, "display %s %d")
+ call pargstr (image)
+ call pargi (clgeti ("frame"))
+
+ im = immap (image, READ_WRITE, 0)
+ while (clgcur ("cursor",wx, wy, wcs, key, cmd, SZ_LINE) != EOF) {
+ switch (key) {
+ case 'a':
+ r = sqrt ((wx - xc) ** 2 + (wy - yc) ** 2)
+ s = 2 * r
+ case 'b':
+ nrep = skyreplace (im, xc, yc, r, s)
+ case 'c':
+ xc = wx
+ yc = wy
+ case 'd':
+ call imunmap (im)
+ call clcmdw (display)
+ im = immap (image, READ_WRITE, 0)
+ case 'g':
+ call imunmap (im)
+ call clcmdw (graph)
+ im = immap (image, READ_WRITE, 0)
+ case 'q':
+ break
+ default:
+ call printf ("\007")
+ }
+ }
+
+ call imunmap (im)
+end
+
+
+define NSKY 100 # Minimum number of sky points
+
+int procedure skyreplace (im, xc, yc, r, s)
+
+pointer im # IMIO pointer
+real xc, yc # Object center
+real r # Object aperture radius
+real s # Sky aperture radius
+
+real avg, sigma, urand(), mode, find_mode()
+long seed
+int xlen, ylen, nx, nx1, nx2, ny, ny1, ny2, ntotal, nobj, nallsky, nsky[4]
+int i, j, x1, x2, x3, x4, y1, y2, y3, y4, y
+pointer sp, allsky, sky[4], ptr1, ptr2
+pointer datain, dataout, imgs2r(), imps2r()
+
+begin
+ xlen = IM_LEN(im,1)
+ ylen = IM_LEN(im,2)
+ x1 = max (1, int (xc - s))
+ x4 = min (xlen, int (xc + s + 0.5))
+ y1 = max (1, int (yc - s))
+ y4 = min (ylen, int (yc + s + 0.5))
+ nx = x4 - x1 + 1
+ ny = y4 - y1 + 1
+ ntotal = nx * ny
+
+ x2 = max (1, int (xc - r))
+ x3 = min (xlen, int (xc + r + 0.5))
+ y2 = max (1, int (yc - r))
+ y3 = min (xlen, int (yc + r + 0.5))
+ nx1 = (x3 - x2 + 1)
+ ny1 = (y3 - y2 + 1)
+ nobj = nx1 * ny1
+ nallsky = ntotal - nobj
+
+ if ((nallsky < NSKY) || (nobj < 1))
+ return (0)
+
+ call smark (sp)
+ call salloc (allsky, nallsky, TY_REAL)
+ datain = imgs2r (im, x1, x4, y1, y4)
+ dataout = imps2r (im, x2, x3, y2, y3)
+ ptr2 = allsky
+
+ # First quadrant
+ x2 = max (1, int (xc - r))
+ x3 = min (xlen, int (xc + 0.5))
+ y2 = max (1, int (yc - r))
+ y3 = min (xlen, int (yc + 0.5))
+ nx1 = x3 - x1 + 1
+ nx2 = x3 - x2
+ ny1 = y2 - y1
+ ny2 = y3 - y2 + 1
+ nsky[1] = nx1 * ny1 + nx2 * ny2
+ sky[1] = ptr2
+
+ if (nsky[1] > 0) {
+ ptr1 = datain
+ for (y=y1; y<y2; y=y+1) {
+ call amovr (Memr[ptr1], Memr[ptr2], nx1)
+ ptr1 = ptr1 + nx
+ ptr2 = ptr2 + nx1
+ }
+ for (; y<=y3; y=y+1) {
+ call amovr (Memr[ptr1], Memr[ptr2], nx2)
+ ptr1 = ptr1 + nx
+ ptr2 = ptr2 + nx2
+ }
+ }
+
+ # Second quadrant
+ x2 = max (1, int (xc + 1.5))
+ x3 = min (xlen, int (xc + r + 0.5))
+ y2 = max (1, int (yc - r))
+ y3 = min (xlen, int (yc + 0.5))
+ nx1 = x4 - x2 + 1
+ nx2 = x4 - x3
+ ny1 = y2 - y1
+ ny2 = y3 - y2 + 1
+ nsky[2] = nx1 * ny1 + nx2 * ny2
+ sky[2] = ptr2
+
+ if (nsky[2] > 0) {
+ ptr1 = datain + x2 - x1
+ for (y=y1; y<y2; y=y+1) {
+ call amovr (Memr[ptr1], Memr[ptr2], nx1)
+ ptr1 = ptr1 + nx
+ ptr2 = ptr2 + nx1
+ }
+ ptr1 = ptr1 + x3 - x2 + 1
+ for (; y<=y3; y=y+1) {
+ call amovr (Memr[ptr1], Memr[ptr2], nx2)
+ ptr1 = ptr1 + nx
+ ptr2 = ptr2 + nx2
+ }
+ }
+
+ # Third quadrant
+ x2 = max (1, int (xc - r))
+ x3 = min (xlen, int (xc + 0.5))
+ y2 = max (1, int (yc + 1.5))
+ y3 = min (xlen, int (yc + r + 0.5))
+ nx1 = x3 - x2
+ nx2 = x3 - x1 + 1
+ ny1 = y3 - y2 + 1
+ ny2 = y4 - y3
+ nsky[3] = nx1 * ny1 + nx2 * ny2
+ sky[3] = ptr2
+
+ if (nsky[3] > 0) {
+ ptr1 = datain + (y2 - y1) * nx
+ for (y=y2; y<=y3; y=y+1) {
+ call amovr (Memr[ptr1], Memr[ptr2], nx1)
+ ptr1 = ptr1 + nx
+ ptr2 = ptr2 + nx1
+ }
+ for (; y<=y4; y=y+1) {
+ call amovr (Memr[ptr1], Memr[ptr2], nx2)
+ ptr1 = ptr1 + nx
+ ptr2 = ptr2 + nx2
+ }
+ }
+
+ # Fourth quadrant
+ x2 = max (1, int (xc + 1.5))
+ x3 = min (xlen, int (xc + r + 0.5))
+ y2 = max (1, int (yc + 1.5))
+ y3 = min (xlen, int (yc + r + 0.5))
+ nx1 = x4 - x3
+ nx2 = x4 - x2 + 1
+ ny1 = y3 - y2 + 1
+ ny2 = y4 - y3
+ nsky[4] = ny1 * nx1 + ny2 * nx2
+ sky[4] = ptr2
+
+ if (nsky[4] > 0) {
+ ptr1 = datain + (y2 - y1) * nx + x3 - x1 + 1
+ for (y=y2; y<=y3; y=y+1) {
+ call amovr (Memr[ptr1], Memr[ptr2], nx1)
+ ptr1 = ptr1 + nx
+ ptr2 = ptr2 + nx1
+ }
+ ptr1 = ptr1 - (x3 - x2 + 1)
+ for (; y<=y4; y=y+1) {
+ call amovr (Memr[ptr1], Memr[ptr2], nx2)
+ ptr1 = ptr1 + nx
+ ptr2 = ptr2 + nx2
+ }
+ }
+
+ # This part is for doing a gradient correction. It is not implemented.
+# if ((nsky[1]>NSKY)&&(nsky[2]>NSKY)&&(nsky[3]>NSKY)&&(nsky[4]>NSKY)) {
+# call asrtr (Memr[sky[1]], Memr[sky[1]], nsky[1])
+# call asrtr (Memr[sky[2]], Memr[sky[2]], nsky[2])
+# call asrtr (Memr[sky[3]], Memr[sky[3]], nsky[3])
+# call asrtr (Memr[sky[4]], Memr[sky[4]], nsky[4])
+
+ # Add a gradient correction here.
+
+# seed = dataout
+# do i = dataout, dataout+nobj-1 {
+# j = 4 * urand (seed) + 1
+# k = 0.95 * nsky[j] * urand (seed)
+# Memr[i] = Memr[sky[j]+k]
+# }
+# } else {
+ call asrtr (Memr[allsky], Memr[allsky], nallsky)
+
+ # Find the mean and sigma excluding the outer 20%
+ x1 = 0.1 * nallsky
+ x2 = 0.9 * nallsky
+ call aavgr (Memr[allsky+x1-1], x2-x1+1, avg, sigma)
+ mode = find_mode (Memr[allsky], nallsky, nallsky / 20)
+ call printf ("Mean = %g, Median = %g, Mode = %g\n")
+ call pargr (avg)
+ call pargr (Memr[allsky+nallsky/2-1])
+ call pargr (mode)
+ for (x1=0; (x1<nallsky)&&(Memr[allsky+x1]<avg-3*sigma); x1=x1+1)
+ ;
+ for (x2=nallsky-1; (x2>0)&&(Memr[allsky+x2]>avg+3*sigma); x2=x2-1)
+ ;
+ nx = x2 - x1 - 1
+
+ seed = dataout
+ do i = dataout, dataout+nobj-1 {
+ j = nx * urand (seed) + x1
+ Memr[i] = Memr[allsky+j]
+ }
+# }
+
+ call sfree (sp)
+ return (nobj)
+end
+
+real procedure find_mode (data, npts, n)
+
+real data[npts] # Data
+int npts # Number of data points
+int n # Bin size
+
+int x, xlast, xmin
+real sumx, sumy, sumxx, sumxy, a, amin
+pointer sp, slope
+
+begin
+ call smark (sp)
+ call salloc (slope, npts - n, TY_REAL)
+
+ sumx = 0.
+ sumy = 0.
+ sumxx = 0.
+ sumxy = 0.
+
+ x = 0
+ xlast = 0
+ while (x < n) {
+ x = x + 1
+ sumx = sumx + x
+ sumy = sumy + data[x]
+ sumxx = sumxx + x ** 2
+ sumxy = sumxy + x * data[x]
+ }
+ amin = (n * sumxy - sumx * sumy) / (n * sumxx - sumx ** 2)
+ xmin = (x + xlast) / 2
+ Memr[slope] = amin
+
+ while (x < npts - n) {
+ x = x + 1
+ xlast = xlast + 1
+ sumx = sumx + x - xlast
+ sumy = sumy + data[x] - data[xlast]
+ sumxx = sumxx + x * x - xlast * xlast
+ sumxy = sumxy + x * data[x] - xlast * data[xlast]
+
+ a = (n * sumxy - sumx * sumy) / (n * sumxx - sumx ** 2)
+ if (a < amin) {
+ amin = a
+ xmin = (x + xlast) / 2
+ }
+ Memr[slope+xlast] = a
+ }
+
+ call gplotv (Memr[slope+11], npts-2*n-22, 1., real (npts-2*n-22), "")
+ call sfree (sp)
+ return (data[xmin])
+end
generated by cgit (git 2.34.1) at 2025-09-20 08:50:30 -0400