aboutsummaryrefslogtreecommitdiff
path: root/pkg/xtools/rmmed.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/xtools/rmmed.x')
-rw-r--r--pkg/xtools/rmmed.x446
1 files changed, 446 insertions, 0 deletions
diff --git a/pkg/xtools/rmmed.x b/pkg/xtools/rmmed.x
new file mode 100644
index 00000000..940fda9f
--- /dev/null
+++ b/pkg/xtools/rmmed.x
@@ -0,0 +1,446 @@
+include <mach.h>
+include <pkg/rmsorted.h>
+
+# RM_MED -- Running median/maximum/minimum library.
+#
+# This is a layer over the sorted running routines.
+# This layer provides:
+#
+# 1. Support for multiple datasets (e.g. pixels in running image)
+# 2. Support for an interior average
+# 3. Support for masks
+# 4. Support for excluded index (e.g. image)
+
+
+# Method object structure.
+define RM_LEN 26 # Structure size
+define RM_RMS Memi[$1] # Pointer to RMEDSRT method
+define RM_BOX Memi[$1+1] # Box size
+define RM_TYPE Memi[$1+2] # Type of output
+define RM_NDATA Memi[$1+3] # Number of datasets
+define RM_PIXTYPE Memi[$1+4] # Internal storage type
+define RM_GOOD Memi[$1+5] # Ptr to good array (box)
+define RM_MASK Memi[$1+6] # Ptr to mask array
+define RM_PWIN Memi[$1+7] # Ptr to packed window data (n*box)
+define RM_POUT Memi[$1+8] # Ptr to packed outlist data (n*box)
+define RM_PMASK Memi[$1+9] # Ptr to mask data (n*(box+15)/16)
+define RM_SETMASK P2S($1+10) # Ptr to set mask array (16)
+define RM_UNSETMASK P2S($1+18) # Ptr to unset mask array (16)
+
+define GOOD Memr[RM_GOOD($1)+$2]
+define MASK Mems[RM_MASK($1)+$2/16]
+define SETMASK Mems[RM_SETMASK($1)+mod($2,16)]
+define UNSETMASK Mems[RM_UNSETMASK($1)+mod($2,16)]
+
+define PWINR Memr[RM_PWIN($1)+RM_BOX($1)*($2-1)]
+define PWINS Mems[RM_PWIN($1)+RM_BOX($1)*($2-1)]
+define POUT Mems[RM_POUT($1)+(RM_BOX($1)+1)/2*($2-1)]
+
+define RM_TYPES "|median|maximum|minimum|"
+define RM_TYMED 1 # Medain
+define RM_TYMAX 2 # Maximum
+define RM_TYMIN 3 # Maximum
+
+
+# RM_MED -- Compute next running value.
+
+real procedure rm_med (rm, nclip, navg, blank, exclude, index, in, mask, nused)
+
+pointer rm #I RM pointer
+real nclip #I Clipping factor
+int navg #I Number of central values to average
+real blank #I Blank value
+int exclude #I Index of excluded data (one indexed)
+int index #I Index of new data (one indexed)
+real in #I Input data value
+short mask #I Input mask value
+short nused #O Number of values in calculated result
+real val #R Return value
+
+int i, j, iexclude
+short s1, s2, ors(), ands()
+pointer rms
+real clip, rmsorted()
+
+begin
+ # Call sorted running routine.
+ rms = RM_RMS(rm)
+ val = rmsorted (rms, nclip, index, in)
+
+ # Set mask if needed.
+ s2 = mod (index-1, RM_BOX(rm))
+ s1 = MASK(rm,s2)
+ if (mask != 0 || s1 != 0) {
+ if (mask != 0)
+ MASK(rm,s2) = ors (s1, SETMASK(rm,s2))
+ else
+ MASK(rm,s2) = ands (s1, UNSETMASK(rm,s2))
+ s1 = MASK(rm,s2)
+ }
+
+ # Recompute value if there are masks or an excluded value.
+ iexclude = mod (exclude-1, RM_BOX(rm))
+ if (s1 == 0 && iexclude < 0) {
+ do s2 = 0, RM_BOX(rm)-1, 16 {
+ s1 = MASK(rm,s2)
+ if (s1 != 0)
+ break
+ }
+ }
+ if (s1 != 0 || iexclude >= 0) {
+ nused = 0
+ do i = 0, RM_BOX(rm)-1 {
+ s2 = IN(rms,i)
+ if (s2 != iexclude) {
+ s1 = MASK(rm,s2)
+ if (s1 == 0) {
+ GOOD(rm,nused) = DATA(rms,i)
+ nused = nused + 1
+ } else if (ands (s1, SETMASK(rm,s2)) == 0) {
+ GOOD(rm,nused) = DATA(rms,i)
+ nused = nused + 1
+ }
+ }
+ }
+
+ if (nused > 2 && nclip > 0.) {
+ i = nused / 2
+ if (mod (nused, 2) == 0)
+ val = (GOOD(rm,i) + GOOD(rm,i-1)) / 2
+ else
+ val = GOOD(rm,i)
+ clip = val + nclip * (val - GOOD(rm,0))
+ do i = nused, 1, -1 {
+ if (GOOD(rm,i-1) < clip)
+ break
+ }
+ nused = i
+ }
+
+ switch (RM_TYPE(rm)) {
+ case RM_TYMED:
+ switch (nused) {
+ case 0:
+ val = blank
+ case 1:
+ val = GOOD(rm,0)
+ case 2:
+ val = (GOOD(rm,0) + GOOD(rm,1)) / 2.
+ default:
+ for (i = 0; nused-2*i>max(2,navg); i=i+1)
+ ;
+ val = GOOD(rm,i)
+ do j = i+1, nused-i-1 {
+ val = val + GOOD(rm,j)
+ }
+ nused = nused - 2 * i
+ val = val / nused
+ }
+ case RM_TYMAX:
+ switch (nused) {
+ case 0:
+ val = blank
+ default:
+ val = GOOD(rm,nused-1)
+ }
+ case RM_TYMIN:
+ switch (nused) {
+ case 0:
+ val = blank
+ default:
+ val = GOOD(rm,0)
+ }
+ }
+ } else
+ nused = min (navg, RM_BOX(rm))
+
+ return (val)
+end
+
+
+# RM_GMED -- Running sorted value.
+
+real procedure rm_gmed (rm, nclip, navg, blank, exclude, nused)
+
+pointer rm #I RM pointer
+real nclip #I Clipping factor
+int navg #I Number of central values to average
+real blank #I Blank value
+int exclude #I Index of excluded data (one indexed)
+short nused #O Number of values in calculated result
+real val #R Return value
+
+int i, j, iexclude
+short mask, ands()
+real clip
+pointer rms
+
+begin
+ rms = RM_RMS(rm)
+ iexclude = mod (exclude-1, RM_BOX(rm))
+
+ # Extract good values to use.
+ nused = 0
+ do i = 0, RM_BOX(rm)-1 {
+ j = IN(rms,i)
+ mask = MASK(rm,j)
+ if (j != iexclude) {
+ if (mask == 0) {
+ GOOD(rm,nused) = DATA(rms,i)
+ nused = nused + 1
+ } else if (ands (mask, SETMASK(rm,j)) == 0) {
+ GOOD(rm,nused) = DATA(rms,i)
+ nused = nused + 1
+ }
+ }
+ }
+
+ if (nused > 2 && nclip > 0.) {
+ i = nused / 2
+ if (mod (nused, 2) == 0)
+ val = (GOOD(rm,i) + GOOD(rm,i-1)) / 2
+ else
+ val = GOOD(rm,i)
+ clip = val + nclip * (val - GOOD(rm,0))
+ do i = nused, 1, -1 {
+ if (GOOD(rm,i-1) < clip)
+ break
+ }
+ nused = i
+ }
+
+ switch (RM_TYPE(rm)) {
+ case RM_TYMED:
+ switch (nused) {
+ case 0:
+ val = blank
+ case 1:
+ val = GOOD(rm,0)
+ case 2:
+ val = (GOOD(rm,0) + GOOD(rm,1)) / 2.
+ default:
+ for (i = 0; nused-2*i>max(2,navg); i=i+1)
+ ;
+ val = GOOD(rm,i)
+ do j = i+1, nused-i-1 {
+ val = val + GOOD(rm,j)
+ }
+ nused = nused - 2 * i
+ val = val / nused
+ }
+ case RM_TYMAX:
+ switch (nused) {
+ case 0:
+ val = blank
+ default:
+ val = GOOD(rm,nused-1)
+ }
+ case RM_TYMIN:
+ switch (nused) {
+ case 0:
+ val = blank
+ default:
+ val = GOOD(rm,0)
+ }
+ }
+
+ return (val)
+end
+
+
+# RM_GDATA -- Get data value for specified index
+
+real procedure rm_gdata (rm, index)
+
+pointer rm #I RM pointer
+int index #I Index of new data (one indexed)
+
+int i, j
+pointer rms
+
+begin
+ rms = RM_RMS(rm)
+ i = mod (index-1, RM_BOX(rm))
+ do j = 0, RM_BOX(rm)-1 {
+ if (IN(rms,j) == i)
+ return (DATA(rms,j))
+ }
+end
+
+
+
+# RM_OPEN -- Open running sorted package.
+#
+# This is called once to allocate memory and initialize the algorithms.
+
+pointer procedure rm_open (box, type, ndatasets, pixtype)
+
+int box #I Median box size (<= 128)
+char type[ARB] #I Output type
+int ndatasets #I Number of datasets
+int pixtype #I Internal storage type
+pointer rm #O RM pointer
+
+char str[8]
+int i, j, strdic()
+short s, nots(), shifts()
+real val
+pointer rms, rms_open()
+
+begin
+ # Set internal storage type.
+ if (pixtype == TY_SHORT)
+ i = TY_SHORT
+ else
+ i = TY_REAL
+
+ # Set the output type.
+ j = strdic (type, str, 8, RM_TYPES)
+ switch (j) {
+ case RM_TYMED:
+ val = 0.
+ rms = rms_open (box, RMS_TYMED, val)
+ case RM_TYMAX:
+ switch (i) {
+ case TY_SHORT:
+ val = -MAX_SHORT
+ rms = rms_open (box, RMS_TYMAX, val)
+ case TY_REAL:
+ val = -MAX_REAL
+ rms = rms_open (box, RMS_TYMAX, val)
+ }
+ case RM_TYMIN:
+ switch (i) {
+ case TY_SHORT:
+ val = MAX_SHORT
+ rms = rms_open (box, RMS_TYMIN, val)
+ case TY_REAL:
+ val = MAX_REAL
+ rms = rms_open (box, RMS_TYMIN, val)
+ }
+ default:
+ call error (1, "Unknown running type")
+ }
+
+ call calloc (rm, RM_LEN, TY_STRUCT)
+ call calloc (RM_GOOD(rm), box, TY_REAL)
+ call calloc (RM_PWIN(rm), box*ndatasets, i)
+ call calloc (RM_POUT(rm), ndatasets*(box+1)/2, TY_SHORT)
+ call calloc (RM_PMASK(rm), ndatasets*(box+15)/16, TY_SHORT)
+
+ RM_RMS(rm) = rms
+ RM_BOX(rm) = box
+ RM_TYPE(rm) = j
+ RM_NDATA(rm) = ndatasets
+ RM_PIXTYPE(rm) = i
+ RM_MASK(rm) = RM_PMASK(rm)
+
+ # Set mask flags.
+ s = 1
+ do i = 0, 15 {
+ SETMASK(rm,i) = s
+ UNSETMASK(rm,i) = nots (s)
+ s = shifts (s, short(1))
+ }
+
+ do i = 1, ndatasets
+ call rm_pack (rm, i)
+
+ return (rm)
+end
+
+
+# RM_CLOSE -- Close running sorted package.
+
+procedure rm_close (rm)
+
+pointer rm #I RM pointer
+
+begin
+ call rms_close (RM_RMS(rm))
+
+ call mfree (RM_GOOD(rm), TY_REAL)
+ call mfree (RM_PWIN(rm), RM_PIXTYPE(rm))
+ call mfree (RM_POUT(rm), TY_SHORT)
+ call mfree (RM_PMASK(rm), TY_SHORT)
+ call mfree (rm, TY_STRUCT)
+end
+
+
+# RM_PACK -- Pack data.
+
+procedure rm_pack (rm, dataset)
+
+pointer rm #I RM pointer
+int dataset #I Data set
+
+pointer rms
+
+begin
+ rms = RM_RMS(rm)
+ if (RM_PIXTYPE(rm) == TY_SHORT)
+ call anirs (DATA(rms,0), PWINS(rm,dataset), RM_BOX(rm))
+# else
+# call amovr (DATA(rms,0), PWINR(rm,dataset), RM_BOX(rm))
+ call achtsb (OUT(rms,0), POUT(rm,dataset), RM_BOX(rm))
+end
+
+
+# RM_UNPACK -- Unpack data.
+
+procedure rm_unpack (rm, dataset)
+
+pointer rm #I RM pointer
+int dataset #I Data set
+
+int i, j, box
+pointer rms
+
+begin
+ rms = RM_RMS(rm)
+ box = RM_BOX(rm)
+
+ if (RM_PIXTYPE(rm) == TY_SHORT)
+ call achtsr (PWINS(rm,dataset), DATA(rms,0), box)
+ else
+ RMS_DATA(rms) = RM_PWIN(rm) + box * (dataset - 1)
+# call amovr (PWINR(rm,dataset), DATA(rms,0), box)
+ call achtbs (POUT(rm,dataset), OUT(rms,0), box)
+ RM_MASK(rm) = RM_PMASK(rm) + (box + 15) / 16 * (dataset - 1)
+
+ do i = 0, box-1 {
+ j = OUT(rms,i)
+ IN(rms,j) = i
+ }
+end
+
+
+# ANIRS -- Convert real to short using nearest integer.
+
+procedure anirs (a, b, n)
+
+real a[n] #I Input real array
+short b[n] #O Output short array
+int n #I Number of array values
+
+int i
+
+begin
+ do i = 1, n
+ b[i] = a[i] + 0.5
+end
+
+
+# RM_DUMP -- Dump data.
+
+procedure rm_dump (rm, unsorted, sorted, in, out)
+
+pointer rm #I Method pointer
+bool unsorted #I Dump data in unsorted order?
+bool sorted #I Dump data in sorted order?
+bool in #I Dump in list?
+bool out #I Dump out list?
+
+begin
+ call rms_dump (RM_RMS(rm), unsorted, sorted, in, out)
+end
generated by cgit (git 2.34.1) at 2025-09-20 05:14:02 -0400