1 files changed, 74 insertions, 0 deletions
diff --git a/pkg/images/tv/imedit/epgsfit.x b/pkg/images/tv/imedit/epgsfit.x
new file mode 100644
index 00000000..976af322
--- /dev/null
+++ b/pkg/images/tv/imedit/epgsfit.x
@@ -0,0 +1,74 @@
+include	<math/gsurfit.h>
+include	"epix.h"
+ 
+# EP_GSFIT -- Fit the background annulus.
+ 
+procedure ep_gsfit (ep, data, mask, x, y, w, nx, ny, gs)
+ 
+pointer	ep			# EPIX structure
+real	data[nx,ny]		# Data subraster
+int	mask[nx,ny]		# Mask subraster
+real	x[nx,ny]		# X positions
+real	y[nx,ny]		# Y positions
+real	w[nx,ny]		# Weights
+int	nx, ny			# Subraster size
+pointer	gs			# Surface pointer (returned)
+ 
+int	i, j, n, npts, xo, yo
+pointer	sp, work
+real	amedr()
+ 
+begin
+	call smark (sp)
+	call salloc (work, nx * ny, TY_REAL)
+
+	gs = NULL
+	npts = nx * ny
+
+	if (EP_XORDER(ep) == 0 || EP_YORDER(ep) == 0) {
+	    n = 0
+	    do j = 1, ny {
+		do i = 1, nx {
+		    if (mask[i,j] == 2) {
+			Memr[work+n] = data[i,j]
+			n = n + 1
+		    }
+		}
+	    }
+	    call amovkr (amedr (Memr[work], n), Memr[work], npts)
+	    xo = 1
+	    yo = 1
+	} else {
+	    call amovr (data, Memr[work], npts)
+	    xo = EP_XORDER(ep)
+	    yo = EP_YORDER(ep)
+	}
+
+	n = 0
+	do j = 1, ny {
+	    do i = 1, nx {
+		x[i,j] = i
+		y[i,j] = j
+		if (mask[i,j] == 2) {
+		    w[i,j] = 1.
+		    n = n + 1
+		} else
+		    w[i,j] = 0.
+	    }
+	}
+ 
+	if (n > 7) {
+	    repeat {
+		call gsinit (gs, GS_POLYNOMIAL, xo, yo, YES,
+		    1., real (nx), 1., real (ny))
+		call gsfit (gs, x, y, Memr[work], w, npts, WTS_USER, n)
+		if (n == OK)
+		    break
+		xo = max (1, xo - 1)
+		yo = max (1, yo - 1)
+	    }
+	} else
+	    call eprintf ("ERROR: Insufficient background points\n")
+
+	call sfree (sp)
+end