diff options
Diffstat (limited to 'pkg/xtools/xtstat.x')
| -rw-r--r-- | pkg/xtools/xtstat.x | 337 |
1 files changed, 337 insertions, 0 deletions
diff --git a/pkg/xtools/xtstat.x b/pkg/xtools/xtstat.x new file mode 100644 index 00000000..1979fddf --- /dev/null +++ b/pkg/xtools/xtstat.x @@ -0,0 +1,337 @@ +# XT_STAT -- Compute statistics from a sample. +# +# The sample array will be sorted. + +define NMIN 10 # Minimum number of pixels for mode calculation +define ZSTEP 0.01 # Step size for search for mode +define ZBIN 0.1 # Bin size for mode. + + +procedure xt_stats (sample, nsample, frac, mean, sigma, median, mode) + +short sample[nsample] #I Sample +int nsample #I Number of sample pixels +real frac #I Fraction of data to use +real mean, sigma, median, mode #O Statistics + +int i, j, k, nmax +real z1, z2, zstep, zbin +bool fp_equalr() + +begin + # Sort the sample. + call asrts (sample, sample, nsample) + + # Set fraction to use. + i = max (1, 1 + nsample * (1. - frac) / 2.) + j = min (nsample, 1 + nsample * (1. + frac) / 2.) + z1 = sample[i] + z2 = sample[j] + + # Compute the mean and sigma. + call aavgs (sample[i], j-i+1, mean, sigma) + + # Compute the median. + median = sample[nsample/2] + + z1 = median - 2 * sigma + if (z1 < sample[1]) + i = 1 + else { + k = i + do i = k, 2, -1 { + if (sample[i] <= z1) + break + } + } + z1 = sample[i] + + z2 = median + 2 * sigma + if (z2 > sample[nsample]) + i = nsample + else { + k = j + do j = k, nsample-1 { + if (sample[j] >= z1) + break + } + } + z2 = sample[j] + + # Compute the mode. + + if (nsample < NMIN) + mode = median + + else if (fp_equalr (z1, z2)) + mode = z1 + + else { + zstep = ZSTEP * sigma + zbin = ZBIN * sigma + zstep = max (1., zstep) + zbin = max (1., zbin) + + z1 = z1 - zstep + k = i + nmax = 0 + repeat { + z1 = z1 + zstep + z2 = z1 + zbin + for (; i < j && sample[i] < z1; i=i+1) + ; + for (; k < j && sample[k] < z2; k=k+1) + ; + if (k - i > nmax) { + nmax = k - i + mode = sample[(i+k)/2] + } + } until (k >= j) + } +end + +procedure xt_stati (sample, nsample, frac, mean, sigma, median, mode) + +int sample[nsample] #I Sample +int nsample #I Number of sample pixels +real frac #I Fraction of data to use +real mean, sigma, median, mode #O Statistics + +int i, j, k, nmax +real z1, z2, zstep, zbin +bool fp_equalr() + +begin + # Sort the sample. + call asrti (sample, sample, nsample) + + # Set fraction to use. + i = max (1, 1 + nsample * (1. - frac) / 2.) + j = min (nsample, 1 + nsample * (1. + frac) / 2.) + z1 = sample[i] + z2 = sample[j] + + # Compute the mean and sigma. + call aavgi (sample[i], j-i+1, mean, sigma) + + # Compute the median. + median = sample[nsample/2] + + z1 = median - 2 * sigma + if (z1 < sample[1]) + i = 1 + else { + k = i + do i = k, 2, -1 { + if (sample[i] <= z1) + break + } + } + z1 = sample[i] + + z2 = median + 2 * sigma + if (z2 > sample[nsample]) + i = nsample + else { + k = j + do j = k, nsample-1 { + if (sample[j] >= z1) + break + } + } + z2 = sample[j] + + # Compute the mode. + + if (nsample < NMIN) + mode = median + + else if (fp_equalr (z1, z2)) + mode = z1 + + else { + zstep = ZSTEP * sigma + zbin = ZBIN * sigma + zstep = max (1., zstep) + zbin = max (1., zbin) + + z1 = z1 - zstep + k = i + nmax = 0 + repeat { + z1 = z1 + zstep + z2 = z1 + zbin + for (; i < j && sample[i] < z1; i=i+1) + ; + for (; k < j && sample[k] < z2; k=k+1) + ; + if (k - i > nmax) { + nmax = k - i + mode = sample[(i+k)/2] + } + } until (k >= j) + } +end + +procedure xt_statr (sample, nsample, frac, mean, sigma, median, mode) + +real sample[nsample] #I Sample +int nsample #I Number of sample pixels +real frac #I Fraction of data to use +real mean, sigma, median, mode #O Statistics + +int i, j, k, nmax +real z1, z2, zstep, zbin +bool fp_equalr() + +begin + # Sort the sample. + call asrtr (sample, sample, nsample) + + # Set fraction to use. + i = max (1, 1 + nsample * (1. - frac) / 2.) + j = min (nsample, 1 + nsample * (1. + frac) / 2.) + z1 = sample[i] + z2 = sample[j] + + # Compute the mean and sigma. + call aavgr (sample[i], j-i+1, mean, sigma) + + # Compute the median. + median = sample[nsample/2] + + z1 = median - 2 * sigma + if (z1 < sample[1]) + i = 1 + else { + k = i + do i = k, 2, -1 { + if (sample[i] <= z1) + break + } + } + z1 = sample[i] + + z2 = median + 2 * sigma + if (z2 > sample[nsample]) + i = nsample + else { + k = j + do j = k, nsample-1 { + if (sample[j] >= z1) + break + } + } + z2 = sample[j] + + # Compute the mode. + + if (nsample < NMIN) + mode = median + + else if (fp_equalr (z1, z2)) + mode = z1 + + else { + zstep = ZSTEP * sigma + zbin = ZBIN * sigma + + z1 = z1 - zstep + k = i + nmax = 0 + repeat { + z1 = z1 + zstep + z2 = z1 + zbin + for (; i < j && sample[i] < z1; i=i+1) + ; + for (; k < j && sample[k] < z2; k=k+1) + ; + if (k - i > nmax) { + nmax = k - i + mode = sample[(i+k)/2] + } + } until (k >= j) + } +end + +procedure xt_statd (sample, nsample, frac, mean, sigma, median, mode) + +double sample[nsample] #I Sample +int nsample #I Number of sample pixels +real frac #I Fraction of data to use +double mean, sigma, median, mode #O Statistics + +int i, j, k, nmax +double z1, z2, zstep, zbin +bool fp_equald() + +begin + # Sort the sample. + call asrtd (sample, sample, nsample) + + # Set fraction to use. + i = max (1, 1 + nsample * (1. - frac) / 2.) + j = min (nsample, 1 + nsample * (1. + frac) / 2.) + z1 = sample[i] + z2 = sample[j] + + # Compute the mean and sigma. + call aavgd (sample[i], j-i+1, mean, sigma) + + # Compute the median. + median = sample[nsample/2] + + z1 = median - 2 * sigma + if (z1 < sample[1]) + i = 1 + else { + k = i + do i = k, 2, -1 { + if (sample[i] <= z1) + break + } + } + z1 = sample[i] + + z2 = median + 2 * sigma + if (z2 > sample[nsample]) + i = nsample + else { + k = j + do j = k, nsample-1 { + if (sample[j] >= z1) + break + } + } + z2 = sample[j] + + # Compute the mode. + + if (nsample < NMIN) + mode = median + + else if (fp_equald (z1, z2)) + mode = z1 + + else { + zstep = ZSTEP * sigma + zbin = ZBIN * sigma + + z1 = z1 - zstep + k = i + nmax = 0 + repeat { + z1 = z1 + zstep + z2 = z1 + zbin + for (; i < j && sample[i] < z1; i=i+1) + ; + for (; k < j && sample[k] < z2; k=k+1) + ; + if (k - i > nmax) { + nmax = k - i + mode = sample[(i+k)/2] + } + } until (k >= j) + } +end + |