include <mach.h>

procedure qdimdn( pxv, px_n, tv, t_n, median, adev)

# determine the robust median and average deviation of pixel vector PXV

real    pxv[ARB]
int	px_n
real    tv[ARB]
double  median
double  adev

int     t_n
int     t_n_old
double  t_mdn
double  t_mdn_old
double  t_adv
double  t_adv_old
double  t_limit
double	t_px
double	t_mean
double	t_min

int     j
int     jj
int     m
int     k

begin

t_n = px_n
t_n_old = t_n

t_mdn_old = pxv[1]
t_adv_old = MAX_DOUBLE

do j = 1, t_n {
   tv[j] = pxv[j]
}

t_mdn = 0
t_limit = MAX_DOUBLE
do k = 0, 5 {

   jj = t_n   
   t_n = 0
   t_mean = 0 
   do j = 1, jj {
      if ( abs(tv[j] - t_mdn) < t_limit ) {
         t_px = tv[j]
         t_mean = t_mean + t_px
         t_n = t_n + 1
         tv[t_n] = t_px
      }
   }
   t_mean = t_mean / t_n

   # sort tv
   call asrtr (tv, tv, t_n)

   # minimum value is the 1st element
   t_min = tv[1]

   # determine the median value of tv
   m = t_n / 2
   if ((m*2)==t_n) {  # t_n is even
      t_mdn = (tv[m] + tv[m+1])/2d0
   } else {            # t_n is odd
      t_mdn = tv[m]
   }

   # determine the average (mean) deviation of the bb vector
   t_adv = 0
   do j = 1, t_n {
      t_adv = t_adv + abs(tv[j] - t_mean)
   }
   t_adv = t_adv / t_n

   # set t_limit for the next iteration
   t_limit = 3*t_adv

   if ( 
      (t_n == t_n_old) 
      && (t_mdn == t_mdn_old) 
      && (t_adv == t_adv_old) 
      ) { 
      break
   } else {
      t_n_old = t_n
      t_mdn_old = t_mdn
      t_adv_old = t_adv
   }

}

median = t_mdn
adev = t_adv

end