From b83c30b594f45615d5274a754171981fd7830406 Mon Sep 17 00:00:00 2001
From: Mike Birchall <mike.birchall@mq.edu.au>
Date: Thu, 17 Jun 2010 11:24:10 +1000
Subject: [PATCH] Added RunMedan directory 	Added RunMedian directory to
 contain all development and testing code 	for efficient running median
 algorithm

---
 RunMedian/Makefile                 |  18 +
 RunMedian/MedianFilter2D.f90       | 153 ++++++++
 RunMedian/RunMedian.f90            | 574 +++++++++++++++++++++++++++++
 RunMedian/RunMedianBad.f90         | 112 ++++++
 RunMedian/bench_test.f90           | 257 +++++++++++++
 RunMedian/quicksort.c              |  47 +++
 RunMedian/runmed_common.mod        |  65 ++++
 RunMedian/runmedianbad_globals.mod |  27 ++
 RunMedian/shell.f90                |  45 +++
 RunMedian/testFilter.f90           |  44 +++
 RunMedian/testmain.f90             | 141 +++++++
 11 files changed, 1483 insertions(+)
 create mode 100644 RunMedian/Makefile
 create mode 100644 RunMedian/MedianFilter2D.f90
 create mode 100644 RunMedian/RunMedian.f90
 create mode 100644 RunMedian/RunMedianBad.f90
 create mode 100644 RunMedian/bench_test.f90
 create mode 100644 RunMedian/quicksort.c
 create mode 100644 RunMedian/runmed_common.mod
 create mode 100644 RunMedian/runmedianbad_globals.mod
 create mode 100644 RunMedian/shell.f90
 create mode 100644 RunMedian/testFilter.f90
 create mode 100644 RunMedian/testmain.f90

diff --git a/RunMedian/Makefile b/RunMedian/Makefile
new file mode 100644
index 0000000..5cfa1cf
--- /dev/null
+++ b/RunMedian/Makefile
@@ -0,0 +1,18 @@
+BenchTest :
+	gcc -c quicksort.c
+#	gfortran -g -fbounds-check bench_test.f90 RunMedian.f90 quicksort.o -o benchtest
+	gfortran  bench_test.f90 RunMedian.f90 quicksort.o -o benchtest
+Test :
+	gfortran -g -fbounds-check testmain.f90 RunMedian.f90 -o test
+RunMedian :
+	gfortran -c RunMedian.f90
+
+MedianFilter :
+	gfortran -c RunMedian.f90 RunMedianBad.f90 MedianFilter2D.f90 shell.f90
+
+TestFilter :
+	gfortran -g -fbounds-check  RunMedian.f90 RunMedianBad.f90 MedianFilter2D.f90 shell.f90 testFilter.f90 -o testfilter
+	
+clean : 
+	rm -rdf *.o
+	rm -rdf test benchtest testfilter	
diff --git a/RunMedian/MedianFilter2D.f90 b/RunMedian/MedianFilter2D.f90
new file mode 100644
index 0000000..f1aadb3
--- /dev/null
+++ b/RunMedian/MedianFilter2D.f90
@@ -0,0 +1,153 @@
+Subroutine MedianFilter2D(InArray, OutArray, n1, n2, nf1, nf2, bad_val, min_val, max_val)
+
+   implicit none
+
+   integer n1,n2,nf1,nf2
+   real InArray(n1,n2)
+   real OutArray(n1,n2)
+   real bad_val, min_val, max_val
+
+   integer filter_size
+
+   filter_size=nf1*nf2
+
+   if (filter_size.ge.200) then
+      call MedianFilter2DLarge(InArray, OutArray, n1, n2, nf1, nf2, bad_val, min_val, max_val)
+   else
+      call MedianFilter2DSmall(InArray, OutArray, n1, n2, nf1, nf2, bad_val)
+   endif
+
+
+End Subroutine MedianFilter2D
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+Subroutine MedianFilter2DSmall(InArray, OutArray, n1, n2, nf1, nf2, bad_val)
+
+   implicit none
+
+   integer n1,n2,nf1,nf2
+   real InArray(n1,n2)
+   real OutArray(n1,n2)
+   real bad_val
+
+   integer i,j
+   integer isub,i1,i2,jsub,j1,j2,idx0
+   integer k1,k2
+   real med
+   real vec(100)
+   integer nvals
+
+   k1=(nf1-1)/2
+   k2=(nf2-1)/2
+
+   ! Iterate through each i,j cell of output array
+   do i=1,n1
+      ! Define the min and max rows of the filter window
+      i1=max(1,i-k1)
+      i2=min(i+k1,n1)
+      do j=1,n2
+        ! Define the min and max cols of the filter window
+         j1=max(1,j-k2)
+         j2=min(j+k2,n2)
+
+         ! Collect all good values in the filter window and determine the median
+         nvals=0
+         do isub=i1,i2
+            do jsub=j1,j2
+               if (InArray(isub,jsub).ne.bad_val) then
+                  nvals=nvals+1
+                  vec(nvals)=InArray(isub,jsub)
+               endif
+            enddo
+         enddo
+
+         ! If no good values then filter value is bad
+         if (nvals.eq.0) then
+            OutArray(i,j)=bad_val
+            cycle
+         endif
+
+         ! Determine the median via shell sort
+         call Shell(vec,nvals)
+         if(mod(nvals,2).eq.1) then
+            ! Odd no of values so median is the middle value
+            idx0=(nvals-1)/2+1
+            med=vec(idx0)
+         else
+            ! Even no of values so median is the average of the two middle values
+            idx0=(nvals)/2-1
+            med=(vec(idx0)+vec(idx0+1))*0.5
+         endif
+         OutArray(i,j)=med
+
+     enddo ! j
+   enddo ! i
+
+End Subroutine MedianFilter2DSmall
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+
+Subroutine MedianFilter2DLarge(InArray, OutArray, n1, n2, nf1, nf2, bad_val, min_val, max_val)
+
+   implicit none
+
+   integer n1,n2,nf1,nf2
+   real InArray(n1,n2)
+   real OutArray(n1,n2)
+   real bad_val, min_val, max_val
+
+   real vec(100), edge_vec(100)
+   integer i,j
+   integer row,col,row_lowbd,row_uppbd
+   integer k,k1,k2,win_size,win_step
+   real med
+   real RunMedianBadInit
+   real RunMedianBad
+
+   bad_val=1.0e08
+   win_size=nf1*nf2
+   win_step=nf1
+   k1=(nf1-1)/2
+   k2=(nf2-1)/2
+
+   ! Define an edge vector of bad values
+   do k=1, win_size
+      edge_vec(k)=bad_val
+   enddo
+
+   do row=1,n1
+
+      ! Determine the upper and lower rows we filter across
+      row_uppbd=max(1,row-k1)
+      row_lowbd=min(row+k1,n1)
+
+      ! The number of rows in the filter determine the window size and step size
+      win_step=row_lowbd-row_uppbd+1
+      win_size=win_step*nf2
+
+      ! Initialise the Bad Running Median with the edgevector (i.e. all bad)
+      med=RunMedianBadInit(edge_vec,win_size,bad_val,-1.0e08,1.0e08)
+
+      ! Now step through the columns up to n2-k2
+      do col=1, n2-k2
+         vec(1:win_step)=InArray(row_uppbd:row_lowbd,col)   ! OutArray(i,j)=RunMedian( InArray(i1:i2,j) ,win_step)
+         med=RunMedianBad(vec,win_step)
+         OutArray(row,col)=med
+      enddo ! col <=n2-k2
+
+      ! Now step through the remaining columns by adding bad vals 
+      vec(1:win_step)=edge_vec(1:win_step)
+      do col=n2-k2+1, n2
+         med=RunMedianBad(vec,win_step)
+         OutArray(row,col)=med
+      enddo ! col > n2-k2
+
+   enddo ! row
+      
+
+End Subroutine MedianFilter2DLarge
+
diff --git a/RunMedian/RunMedian.f90 b/RunMedian/RunMedian.f90
new file mode 100644
index 0000000..a7d3291
--- /dev/null
+++ b/RunMedian/RunMedian.f90
@@ -0,0 +1,574 @@
+!--------------------------------------------------------
+!--------------------------------------------------------
+!         ==================================
+!          RunMedian Common Variable Module
+!         ==================================
+
+module RUNMED_COMMON
+   integer H(81)         ! The H index array. H(i) is the index of the value at node Hi
+   integer H_parent(81)  ! The H node of the parrent of node Hi
+   integer H_child1(81)  ! The H node of the first child of node Hi (if there is one, 0 otherwise)
+   integer H_child2(81)  ! The H node of the second child of node Hi (if there is one, 0 otherwise)
+   integer WIN(81)       ! The Window indexing array
+   real X(81)            ! The vector values, circular buffer.
+   logical inMaxHeap(81) ! Specifies if node Hi is in the max heap (i.e. if node i index < node 0 index)
+   logical inMinHeap(81) ! Specifies if node Hi is in the min heap (i.e. if node i index > node 0 index)
+   logical has2Kids(81)  ! Specifies if node Hi has 2 kids
+   logical hasNoKids(81) ! Specifies if node Hi has no kids
+   integer win_k, win_size ! Key structure parameters. win_size is the size of the window=2*win_k+1
+   integer shift_h       ! The shift in the H index array to correspond to node index,eg 
+                         ! H0 =H( 0+shift_h)
+                         ! H-5=H(-5+shift_h
+   integer startx_idx    ! The next index of the circular buffer array (X) for overwritting with new data
+   integer h0_idx, min_apex, max_apex ! Key structure indicies: index of H0, the min heap apex and the max heap apex
+
+end module 
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+Real Function Median(vec,n)
+
+   implicit none
+
+   integer n
+   real vec(n)
+
+   real buffer(n)
+   real med
+   integer i
+   real RunMedian
+   real RunMedianInit
+
+   buffer=0.0
+   med=RunMedianInit(buffer,n)
+   do i=1,n
+      med=RunMedian(vec(i),1)
+   enddo
+
+   Median=med
+
+End Function Median
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+Real Function RunMedianInit(xvec,vec_size)
+
+   use RUNMED_COMMON
+   implicit none
+   integer vec_size
+   real xvec(vec_size)
+
+   integer i 
+   real med 
+
+   real RunMedian
+
+   win_size=vec_size
+   win_k=(win_size-1)/2
+   shift_h=win_k+1
+   h0_idx=shift_h
+   min_apex=h0_idx+1
+   max_apex=h0_idx-1
+
+   ! Initialise the data stack and the H index vector
+   inMaxHeap=.false.
+   inMinHeap=.false.
+   do i=1,vec_size
+      H(i)=i
+      H_parent(i)=(i-shift_h)/2+shift_h
+      H_child1(i)=(i-shift_h)*2+shift_h
+      has2Kids(i)=.true.
+      hasNoKids(i)=.false.
+      if (i.gt.h0_idx) inMinHeap(i)=.true.
+      if (i.lt.h0_idx) inMaxHeap(i)=.true.
+      if (i.lt.h0_idx) H_child2(i)=H_child1(i)-1
+      if (i.gt.h0_idx) H_child2(i)=H_child1(i)+1
+      if (i.lt.h0_idx) H_child2(i)=H_child1(i)-1
+      if (abs(H_child2(i)-shift_h).gt.win_k) then 
+         H_child2(i)=0
+         has2Kids(i)=.false.
+      endif
+      if (abs(H_child1(i)-shift_h).gt.win_k) then
+         H_child1(i)=0
+         has2Kids(i)=.false.
+      endif
+      if (H_child1(i).eq.0 .and.H_child2(i).eq.0) hasNoKids(i)=.true.
+      
+      WIN(i)=i
+   end do
+
+   ! Specify from which index we start overwritting the data stack with new data
+   startx_idx=1   
+
+   ! Now do an initial sorting of the vector data into the circular buffer.
+   X=0.0 ! Give default values to the circular buffer
+         ! (not necessary but makes me feel better)
+   ! Add the initial values to the heap system one val at a time
+   do i=1,vec_size
+      med=RunMedian(xvec(i),1)
+   enddo
+
+  RunMedianInit=X(shift_h)
+
+End Function RunMedianInit
+
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+real Function RunMedian(shift_vec,shift_size)
+   use RUNMED_COMMON
+   implicit none
+   integer shift_size
+   real shift_vec(shift_size)
+
+   integer i, x_idx,h_idx
+   real  new_x, h0_val
+
+   ! Shuffle the new set a values into the heap system one at a time.
+   do i=1,shift_size
+
+      !call CheckH ! For debugging purposes, checks that the H data 
+                   ! structure is correct. Will stop and print out 
+                   ! the faulty structure if found.
+
+      new_x=shift_vec(i)
+
+      ! Overwrite old datum to be removed with new datum.
+      x_idx=startx_idx
+      X(x_idx)=new_x
+
+      ! Move the start_xidx to point to the next datum index for overwritting
+      startx_idx=startx_idx+1
+      if (startx_idx.gt.win_size) startx_idx=1
+
+      ! Get the H node index for this datum's location in the dual heap
+      h_idx=WIN(x_idx)
+
+      ! Get the current value at H0, i.e. the current median value for
+      ! comparison purpose.
+      h0_val=X(H(h0_idx))
+
+      ! ===============
+      ! DO THE SHUFFLE:
+      ! ===============
+
+      ! So lets check if we are in the min or max heap
+      if (inMinHeap(h_idx)) then 
+
+         !-----------------------
+         ! We are in the min heap
+         !-----------------------
+
+         ! Lets check if we want to stay in the min heap
+         ! want to move to just h0 or move into the max heap
+         if(new_x.ge.X(H(h0_idx))) then
+            ! We just want to shuffle within the min heap
+            if (new_x.le.X(H(H_parent(h_idx)))) then
+               call ShuffleDownMinHeap(h_idx,new_x)
+               cycle
+            else
+               call ShuffleUpMinHeap(h_idx,new_x)
+               cycle
+            endif
+         elseif(new_x.ge.X(H(h0_idx))) then
+            call Shuffle2H0(h_idx)     ! Just shuffle to H0
+            cycle
+         else 
+            call Shuffle2H0(h_idx)               ! Shuffle to H0
+            if (new_x.ge.X(H(max_apex))) cycle   ! Cannot shuffle any further
+            call SwapNodes(h0_idx,max_apex)      ! Shuffle from H0 to max apex
+            h_idx=max_apex
+            call ShuffleDownMaxHeap(h_idx,new_x) ! Then shuffle down the max heap
+            cycle
+         endif
+      elseif (inMaxHeap(h_idx)) then
+
+         !-----------------------
+         ! We are in the max heap
+         !-----------------------
+
+         ! Lets check if we want to stay in the max heap
+         ! want to move to just h0 or move into the min heap
+         if(new_x.le.X(H(h0_idx))) then
+            ! We just want to shuffle within the max heap
+            if (new_x.ge.X(H(H_parent(h_idx)))) then
+               call ShuffleUpMaxHeap(h_idx,new_x)
+            else
+               call ShuffleDownMaxHeap(h_idx,new_x)
+            endif
+            cycle
+         elseif(new_x.le.X(H(min_apex))) then
+            call Shuffle2H0(h_idx)     ! Just shuffle to H0
+            cycle
+         else 
+            call Shuffle2H0(h_idx)              ! Shuffle to H0
+            call SwapNodes(h0_idx,min_apex)     ! Shuffle from H0 to min apex
+            h_idx=min_apex
+            call ShuffleUpMinHeap(h_idx,new_x) ! Then shuffle up the min heap
+            cycle
+         endif
+
+      else
+
+          !-------------
+          ! We are in H0
+          !-------------
+
+          ! First check if new value is going to be the new median
+          ! In which case do nothing.
+          if (new_x.ge.X(H(max_apex)).and.new_x.le.X(H(min_apex))) cycle
+
+          ! Now check if new value should go to the min heap or max heap
+          if (new_x.gt.X(H(min_apex))) then
+             ! From H0 To Min Heap
+             call SwapNodes(h0_idx,min_apex)     ! Shuffle from H0 to min apex
+             h_idx=min_apex
+             call ShuffleUpMinHeap(h_idx,new_x)  ! Shuffle up the min heap
+             cycle
+          else
+             ! From H0 To Max Heap
+             call SwapNodes(h0_idx,max_apex)      ! Shuffle from H0 to max apex
+             h_idx=max_apex
+             call ShuffleDownMaxHeap(h_idx,new_x) ! Shuffle down the min heap
+             cycle
+          endif
+
+      endif
+
+   enddo
+
+   ! Return the value at H0 as the median.
+   RunMedian=X(H(h0_idx))
+
+End Function RunMedian
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+subroutine ShuffleDownMinHeap_old(idx,xval)
+
+   use RUNMED_COMMON
+   implicit none
+   real xval
+   integer idx
+   integer x_idx0, parent_idx
+   x_idx0=H(1)
+   x_idx0=H(idx)
+   parent_idx=h_parent(idx)
+   do 
+      H(idx)=H(parent_idx)
+      WIN(H(idx))=idx
+      idx=parent_idx
+      parent_idx=h_parent(idx)
+      if (xval.ge.X(H(parent_idx))) exit
+   enddo
+
+   H(idx)=x_idx0
+   WIN(H(idx))=idx
+
+end subroutine ShuffleDownMinHeap_old
+
+subroutine ShuffleDownMinHeap(idx,xval)
+
+   use RUNMED_COMMON
+   implicit none
+   real xval
+   integer idx
+
+   do 
+      call SwapNodes(idx,h_parent(idx))
+      idx=h_parent(idx)
+      if (xval.ge.X(H(h_parent(idx)))) exit
+   enddo
+
+end subroutine ShuffleDownMinHeap
+
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+subroutine ShuffleUpMaxHeap(idx,xval)
+
+   use RUNMED_COMMON
+   implicit none
+   real xval
+   integer idx
+   integer x_idx0, parent_idx
+
+   x_idx0=H(idx)
+   parent_idx=h_parent(idx)
+   do
+      H(idx)=H(parent_idx)
+      WIN(H(idx))=idx
+      idx=parent_idx
+      parent_idx=h_parent(idx)
+      if (xval.le.X(H(parent_idx))) exit
+   enddo
+
+   H(idx)=x_idx0
+   WIN(H(idx))=idx
+
+end subroutine ShuffleUpMaxHeap
+!--------------------------------------------------------
+!--------------------------------------------------------
+subroutine ShuffleUpMinHeap(idx,xval)
+   use RUNMED_COMMON
+   implicit none
+   real xval
+   integer idx
+
+   real xval1,xval2
+
+   ! Assume that idx is in the minimum heap already
+
+   ! Now shuffle down the children of minimum value until
+   ! we cannot go any further.
+   
+   ! Keep descending whilst there are two children
+   do while (has2Kids(idx))
+   
+      xval1=X(H(h_child1(idx)))
+      xval2=X(H(h_child2(idx)))
+
+      ! Find which of the two children have the minimum value
+      ! If possible to swap then do so else we have finished.
+      if ( xval1 .le. xval2 ) then
+         ! child 1 has the minimum value
+         if (xval.gt.xval1) then
+            ! xval is greater than child 1 so swap.
+            call SwapNodes(idx,h_child1(idx))
+            idx=h_child1(idx)
+         else
+            ! We cannot shuffle any further
+            return 
+         endif
+      else 
+         ! child 2 has the minimum value 
+         if (xval.gt.xval2) then
+            ! xval is greater than child 2 so swap.
+            call SwapNodes(idx,h_child2(idx))
+            idx=h_child2(idx)
+         else
+            ! We cannot shuffle any further
+            return 
+         endif
+      endif
+   enddo
+
+   ! If there are no more children then job done
+   if (hasNoKids(idx)) return
+
+   ! Now special case of only one child, which should always be child1
+   xval1=X(H(h_child1(idx)))
+   if (xval.gt.xval1) then
+      ! xval is greater than child 1 so swap.
+      call SwapNodes(idx,h_child1(idx))
+      idx=h_child1(idx)
+   endif
+            
+   ! We cannot shuffle any further
+   return 
+
+end subroutine ShuffleUpMinHeap
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+subroutine ShuffleDownMaxHeap(idx,xval)
+
+   use RUNMED_COMMON
+   implicit none
+   real xval
+   integer idx
+   real xval1, xval2
+ 
+   ! Assume that idx is in the maximum heap already
+
+   ! Now shuffle down the children of maximum value until
+   ! we cannot go any further.
+   
+   ! Keep descending whilst there are two children
+   do while (has2Kids(idx))
+      xval1=X(H(h_child1(idx)))
+      xval2=X(H(h_child2(idx)))
+      ! Find which of the two children have the maximum value
+      ! If possible to swap then do so else we have finished.
+      if ( xval1 .ge. xval2 ) then
+         ! child 1 has the maximum value
+         if ( xval.lt. xval1 ) then
+            ! xval is smaller than child 1 so swap.
+            call SwapNodes(idx,h_child1(idx))
+            idx=h_child1(idx)
+         else
+            ! We cannot shuffle any further
+            return 
+         endif
+      else 
+         ! child 2 has the maximum value 
+         if ( xval.lt.xval2 ) then
+            ! xval is smaller than child 2 so swap.
+            call SwapNodes(idx,h_child2(idx))
+            idx=h_child2(idx)
+         else
+            ! We cannot shuffle any further
+            return 
+         endif
+      endif
+   enddo
+
+   ! If there are no more children then job done
+   if (hasNoKids(idx)) return
+
+   ! Now special case of only one child, which should always be child1
+   xval1=X(H(h_child1(idx)))
+   if (xval.lt.xval1 ) then
+      ! xval is smaller than child 1 so swap.
+      call SwapNodes(idx,h_child1(idx))
+      idx=h_child1(idx)
+   endif
+            
+   ! We cannot shuffle any further
+   return 
+
+end subroutine ShuffleDownMaxHeap
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+subroutine Shuffle2H0(idx)
+
+   use RUNMED_COMMON
+   implicit none
+   integer idx,idx0,x_idx0, parent_idx
+
+   idx0=idx
+   x_idx0=H(idx)
+   parent_idx=h_parent(idx)
+
+   do
+      H(idx)=H(parent_idx)
+      WIN(H(idx))=idx
+      idx=parent_idx
+      parent_idx=h_parent(idx)
+      if (idx.eq.h0_idx) exit
+   enddo
+
+   H(h0_idx)=x_idx0
+   WIN(x_idx0)=h0_idx
+
+end subroutine Shuffle2H0
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+Subroutine SwapNodes(i1,i2)
+
+   use RUNMED_COMMON
+   implicit none
+
+   integer i1,i2
+
+
+   integer tmp
+
+   WIN(H(i2))=i1
+   WIN(H(i1))=i2
+   tmp=H(i2)
+   H(i2)=H(i1)
+   H(i1)=tmp
+
+End Subroutine SwapNodes
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+!              ======================
+!                  DEBUG ROUTINES
+!              ======================
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+subroutine PrintH() 
+   use RUNMED_COMMON
+   character*3 HC(81)
+   integer i, shift_hc
+   HC='**'
+   shift_hc=shift_h+7
+   do i=-win_k,win_k
+      if (abs(i).le.win_k) then
+         write(HC(i+shift_hc),'(I3)'), INT(X(H(i+shift_h)))
+      endif
+    enddo
+ 
+   print *,''
+   print *,'H(4)=',HC(4+shift_hc),'     H(5)=',HC(5+shift_hc),'     H(6)=',HC(6+shift_hc),'     H(7)=',HC(7+shift_hc)
+   print *,'      \        /                 \        /   '
+   print *,'       \      /                   \      /   '
+   print *,'        \    /                     \    /   '
+   print *,'       H(2)=',HC(2+shift_hc),'                  H(3)=',HC(3+shift_hc)
+   print *,'                \              /'
+   print *,'                 \            /'
+   print *,'                  \          /'
+   print *,'                    H(1)=',HC(1+shift_hc)
+   print *,'                        |'
+   print *,'                    H(0)=',HC(shift_hc)
+   print *,'                        |'
+   print *,'                    H(-1)=',HC(shift_hc-1)
+   print *,'                  /          \'
+   print *,'                 /            \'
+   print *,'                /              \'
+   print *,'       H(-2)=',HC(shift_hc-2),'                  H(-3)=',HC(shift_hc-3)
+   print *,'        /    \                     /    \   '
+   print *,'       /      \                   /      \   '
+   print *,'      /        \                 /        \   '
+   print *,'H(-4)=',HC(shift_hc-4),'     H(-5)=',HC(shift_hc-5),'     H(-6)=',HC(shift_hc-6),'     H(-7)=',HC(shift_hc-7)
+   print *,''
+end Subroutine PrintH
+subroutine PrintH2() 
+   use RUNMED_COMMON
+   character*3 HC(256)
+   integer i, shift_hc
+   HC='**'
+   do i=1, win_size
+      print *,i, 'H(',i-shift_h,')=',H(i)
+   enddo
+end subroutine PrintH2
+subroutine CheckH 
+
+use RUNMED_COMMON
+implicit none
+integer i,j
+logical ok_status
+
+
+do i=1, win_size
+   ok_status=.false.
+   if (i==h0_idx) then
+      if ( X(H(min_apex)) .ge. X(H(h0_idx)) .and. X(H(h0_idx)) .ge. X(H(max_apex)) ) then
+         ok_status=.true.
+      endif
+   elseif(i.gt.h0_idx) then
+      if ( X(H(i)) .ge. X(H(H_parent(i))) ) ok_status=.true.
+   else
+      if ( X(H(i)) .le. X(H(H_parent(i))) ) ok_status=.true.
+   endif
+
+   if(ok_status) cycle
+
+   print *, 'ERROR !!!'
+  if (win_size.le.50) then
+     call PrintH
+  else
+     !call PrintH2
+     do j=1,win_size
+print *, 'H(',j-shift_h,')=',X(H(j)),'Parent=H(',H_parent(H(j))-shift_h,')=', X(H_parent(H(j)))
+     enddo
+  endif
+stop
+
+enddo
+print *, 'Check Ok', X(H(h0_idx))
+end subroutine CheckH
+
diff --git a/RunMedian/RunMedianBad.f90 b/RunMedian/RunMedianBad.f90
new file mode 100644
index 0000000..76df416
--- /dev/null
+++ b/RunMedian/RunMedianBad.f90
@@ -0,0 +1,112 @@
+
+!         =====================================
+!                    RunMedianBad 
+!         =====================================
+
+! Wrappers to RunMedian that allows for Bad Values
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+!         =====================================
+!          RunMedianBad Common Variable Module
+!         =====================================
+
+
+Module RunMedianBad_Globals
+   logical next_bad_use_max  ! In next bad val substition use max val
+   real bad_val              ! Value of the bad value magic number
+   real max_val              ! Value of max value substition of bad value (this must be greater than max data set value)
+   real min_val              ! Value of mix value substition of bad value (this must be smaller than min data set value)
+End Module RunMedianBad_Globals
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+Real Function RunMedianBadInit(vec,size,std_bad_val,use_bad_min_val,use_bad_max_val)
+
+   use RunMedianBad_Globals
+   implicit none
+
+   integer size
+   real vec(size)
+   real std_bad_val,use_bad_max_val,use_bad_min_val
+
+   integer i
+   real use_vec(100)
+   real med
+
+   real RunMedianInit
+
+   ! Initialise the global variables
+   bad_val=std_bad_val
+   max_val=use_bad_max_val
+   min_val=use_bad_min_val
+   next_bad_use_max=.true.
+
+   ! Create a subsitute vector alternately replacing all bad_vals with max and min vals
+   do i=1,size
+      if (vec(i).eq.bad_val) then
+         if (next_bad_use_max) then
+            use_vec(i)=max_val
+            next_bad_use_max=.false.
+          else
+            use_vec(i)=min_val
+            next_bad_use_max=.true.
+         endif
+      else
+         use_vec(i)=vec(i)
+      endif
+   enddo
+
+   ! Call the RunMedianInit with the substitute vector
+   med=RunMedianInit(use_vec,size)
+
+   ! If returned med value is either the min or max value then return bad val
+   if (med==max_val) med=bad_val
+   if (med==min_val) med=bad_val
+
+   RunMedianBadInit=med
+
+End Function RunMedianBadInit
+
+!--------------------------------------------------------
+!--------------------------------------------------------
+
+Real Function RunMedianBad(vec,size)
+
+   use RunMedianBad_Globals
+   implicit none
+   integer size
+   real vec(size)
+   integer i
+   real use_vec(100)
+   real med
+
+   real RunMedian
+
+   ! Create a subsitute vector alternately replacing all bad_vals with max and min vals
+   do i=1,size
+      if (vec(i).eq.bad_val) then
+         if (next_bad_use_max) then
+            use_vec(i)=max_val
+            next_bad_use_max=.false.
+         else
+            use_vec(i)=min_val
+            next_bad_use_max=.true.
+         endif
+      else
+         use_vec(i)=vec(i)
+      endif
+   enddo
+
+   ! Call RunMedian with the substitute vector
+   med=RunMedian(use_vec,size)
+
+   ! If returned med value is either the min or max value then return bad val
+   if (med==max_val) med=bad_val
+   if (med==min_val) med=bad_val
+
+   RunMedianBad=med
+
+End Function RunMedianBad
+
diff --git a/RunMedian/bench_test.f90 b/RunMedian/bench_test.f90
new file mode 100644
index 0000000..a3c9d72
--- /dev/null
+++ b/RunMedian/bench_test.f90
@@ -0,0 +1,257 @@
+
+program bench_test
+
+   implicit none
+
+   real rdata(8482860)
+   real RunMedian, RunMedianInit, Median
+   integer n,i,i0,i1,i2
+   integer win_size, win_step
+   real, allocatable:: window(:)
+   real, allocatable:: vec(:)
+   real xmed,rmed,xmed2
+   integer ndiff, option
+   integer count, count0,j
+   real fixed_vec(1000)
+   real fixed_vec1(1000)
+   real fixed_vec0(1000)
+   real vec2(1000)
+   logical first
+
+
+   !real quicksort
+
+   option=1
+print *, 'option=',option
+
+
+   n=8482860!2000*2000
+   win_size=9*9!15!9*9!25!9*9!11*11!9!25!81!25!81!11*11!81!25
+   win_step=9!7!9!5!9!11!3!5!9!5!9!11!9!5
+   allocate(window(win_size))
+   allocate(vec(win_step))
+
+   CALL RANDOM_NUMBER(rdata)
+   rdata=int(rdata*100)
+
+if (option.eq.1) then
+
+   print *, 'Option=1: Comparing values from RunMedian and CLA_1MED'
+   rmed= RunMedianInit(rdata(1:win_size),win_size)
+
+   ndiff=0
+   first=.true.
+   do i=win_step+1,n-win_size, win_step
+      ! A full window vector for the CLA_1MED subroutine
+      window=rdata(i:i+win_size-1)
+
+      ! Just the window step vector for the RunMedian function
+      vec=rdata(i+win_size-win_step:i+win_size-1)
+
+      ! Call RunMedian and CLA_1MED
+      rmed=RunMedian(vec,win_step)
+      CALL cla_1med(win_size,window, xmed)
+
+      ! Compare values
+      if (xmed.ne.rmed) then 
+         print *, i, xmed,rmed
+         ndiff=ndiff+1
+         print *, '=====================================================>'
+         vec2(1:win_size)=rdata(i:i+win_size-1)
+         call BubbleSort(vec2,win_size)
+         do j=1,win_size
+            print *, j, rdata(i+j-1), window(j), vec2(j)
+         enddo
+         call PrintH2
+         stop
+      endif
+
+   enddo
+
+   print *, 'Ndiff=',ndiff
+   stop
+endif
+
+if (option.eq.2) then
+
+   print *, 'Option 2: Timing RunMedian'
+   call system_clock(count0)
+   window=rdata(1:win_size)
+   rmed=RunMedianInit(window,win_size)
+   i0=win_step+1
+   i2=win_size+win_step
+   i1=i2-win_step+1
+   do while (i2.le.n)
+      do j=1,win_step      
+         vec(j)=rdata(i1+j-1)
+      enddo
+      rmed=RunMedian(vec,win_step)
+
+      i0=i0+win_step
+      i1=i1+win_step
+      i2=i2+win_step
+   enddo
+   call system_clock(count)
+   print *, 'Time 1:', count-count0, count, count0
+endif
+
+if (option.eq.3.or.option.eq.4) then
+
+   print *, 'Option 3: Timing CLA_1MED'
+
+   call system_clock(count0)
+   i0=win_step+1
+   i2=win_size+win_step
+   i1=i2-win_step+1
+
+   do while (i2.le.n)
+      do j=1,win_size      
+         fixed_vec0(j)=rdata(i0+j-1)
+         fixed_vec(j)=rdata(i0+j-1)
+         fixed_vec1(j)=rdata(i0+j-1)
+      enddo
+      if (option.eq.3) CALL cla_1med(win_size,fixed_vec, xmed)
+      if (option.eq.4) call quicksort(fixed_vec1,win_size)
+       xmed2=fixed_vec1(13)
+      i0=i0+win_step
+      i1=i1+win_step
+      i2=i2+win_step
+   enddo
+   call system_clock(count)
+   print *, 'Time 2:', count-count0, count, count0
+
+endif
+
+end Program bench_test
+
+Subroutine BubbleSort(x,n)
+integer n
+real x(n)
+integer i, maxi
+real max
+do j=n,2,-1
+   max=x(1)
+   maxi=1
+   do i=1,j
+      if (x(i).gt.max) then
+         max=x(i)
+         maxi=i
+      endif
+   enddo
+   x(maxi)=x(j)
+   x(j)=max
+enddo
+
+end subroutine BubbleSort
+
+      SUBROUTINE CLA_1MED(N,A,MEDIAN)
+!+
+!  Purpose:
+!     Find the median value of a vector
+!
+!  Description:
+!     This is done by (destructively) sorting the vector values into
+!     ascending order. The median is then either the central value or the
+!     average of the two central values (if N is even).
+!
+!     The sorting is done using what is essentially the `quicksort' algorithm
+!     which works by partitioning the vector repeatedly. For finding the
+!     median it is not necessary to fully sort the vector, but only at each
+!     step to sort the partition which contains the central value. This makes
+!     it easier to use the (essentially recursive) quicksort algorithm in
+!     Fortran since it is not necessary to keep track of the other unsorted
+!     partitions.
+!
+!     When the partition size falls to less than seven elements the remaining
+!     data is sorted by straight insertion which is generally faster for
+!     small datasets.
+! 
+!  N.B. This routine cannot cope with bad pixels, so these must be removed
+!       first if present.
+!
+!  Arguments:
+!     N = INTEGER (Given)
+!        Number of elements in vector A
+!     A = REAL VECTOR(N) (Given and Returned)
+!        Vector to find the median of.
+!        N.B. This vector will be partially reordered on return.
+!             It should not contain any bad values.
+!     MEDIAN = REAL (Returned)
+!        Median value
+!-
+      IMPLICIT NONE              ! No implicit typing
+
+! Arguments
+      INTEGER N
+      REAL A(N)
+      REAL MEDIAN
+
+! Local variables
+      INTEGER I,J,FIRST,LAST,CENT
+      REAL VMID,A1,TEMP
+      LOGICAL SORTED
+
+      FIRST = 1
+      LAST = N
+      CENT = N/2 + 1
+      SORTED = .FALSE.
+
+! Use modified quicksort algorithm for N > 7      
+      DO WHILE ((.NOT. SORTED).AND.(LAST-FIRST.GT.7))
+         I = FIRST
+         J = LAST
+         VMID = A((I+J)/2)
+         DO WHILE (I.LE.J)
+            DO WHILE (A(I).LT.VMID)
+               I = I+1
+            ENDDO
+            DO WHILE (VMID.LT.A(J))
+               J = J-1
+            ENDDO
+            IF (I.LE.J) THEN
+               TEMP = A(I)
+               A(I) = A(J)
+               A(J) = TEMP
+               I = I+1
+               J = J-1
+            ENDIF
+         ENDDO
+         
+! Set FIRST and LAST to sort whichever partition contains central value
+         IF ((FIRST.LT.J).AND.(CENT.LE.J)) THEN
+            LAST = J
+         ELSE IF ((I.LE.CENT).AND.(I.LT.LAST)) THEN
+            FIRST = I
+         ELSE
+            SORTED = .TRUE.
+         ENDIF
+      ENDDO    
+
+! When less than 7 elements sort by straight insertion
+      IF (.NOT. SORTED) THEN
+         DO J=FIRST+1,LAST
+            A1 = A(J)
+            DO I=J-1,FIRST,-1
+               IF (A(I).LE.A1) GOTO 10
+               A(I+1) = A(I)
+            ENDDO
+            I = FIRST-1
+ 10         CONTINUE
+            A(I+1) = A1
+         ENDDO
+      ENDIF  
+      
+! For odd number of elements median is the central value
+      MEDIAN = A(CENT)
+
+      IF (MOD(N,2).EQ.0) THEN
+
+! For even number have to use average of this and next highest
+         TEMP = A(CENT-1)
+          DO I=1,CENT-2
+             IF (A(I).GT.TEMP) TEMP = A(I)
+          ENDDO
+          MEDIAN = (MEDIAN+TEMP)*0.5
+       ENDIF
+       END
+
diff --git a/RunMedian/quicksort.c b/RunMedian/quicksort.c
new file mode 100644
index 0000000..ca0568e
--- /dev/null
+++ b/RunMedian/quicksort.c
@@ -0,0 +1,47 @@
+void quicksort_(float numbers[], int array_size)
+{
+int j;
+  void q_sort(float numbers[], int left, int right);
+  q_sort(numbers, 0, array_size - 1);
+  for (j=0;j<array_size;j++)
+  {
+     write("%d %f",j,numbers[j]);
+  }
+  
+}
+ 
+ 
+void q_sort(float numbers[], int left, int right)
+{
+  int pivot, l_hold, r_hold;
+ 
+  l_hold = left;
+  r_hold = right;
+  pivot = numbers[left];
+  while (left < right)
+  {
+    while ((numbers[right] >= pivot) && (left < right))
+      right--;
+    if (left != right)
+    {
+      numbers[left] = numbers[right];
+      left++;
+    }
+    while ((numbers[left] <= pivot) && (left < right))
+      left++;
+    if (left != right)
+    {
+      numbers[right] = numbers[left];
+      right--;
+    }
+  }
+  numbers[left] = pivot;
+  pivot = left;
+  left = l_hold;
+  right = r_hold;
+  if (left < pivot)
+    q_sort(numbers, left, pivot-1);
+  if (right > pivot)
+    q_sort(numbers, pivot+1, right);
+}
+
diff --git a/RunMedian/runmed_common.mod b/RunMedian/runmed_common.mod
new file mode 100644
index 0000000..fce7565
--- /dev/null
+++ b/RunMedian/runmed_common.mod
@@ -0,0 +1,65 @@
+GFORTRAN module created from RunMedian.f90 on Thu Jun 17 21:22:04 2010
+If you edit this, you'll get what you deserve.
+
+(() () () () () () () () () () () () () () () () () () () () ())
+
+()
+
+()
+
+()
+
+()
+
+(2 'h0_idx' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN) (INTEGER 4 ()) 0 0 () () 0 () ())
+3 'h' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
+DIMENSION) (INTEGER 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER 4 ()) 0
+'1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+4 'h_child1' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN DIMENSION) (INTEGER 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER
+4 ()) 0 '1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+5 'h_parent' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN DIMENSION) (INTEGER 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER
+4 ()) 0 '1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+6 'h_child2' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN DIMENSION) (INTEGER 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER
+4 ()) 0 '1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+7 'has2kids' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN DIMENSION) (LOGICAL 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER
+4 ()) 0 '1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+8 'inminheap' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN DIMENSION) (LOGICAL 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER
+4 ()) 0 '1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+9 'inmaxheap' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN DIMENSION) (LOGICAL 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER
+4 ()) 0 '1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+10 'max_apex' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN) (INTEGER 4 ()) 0 0 () () 0 () ())
+11 'shift_h' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN) (INTEGER 4 ()) 0 0 () () 0 () ())
+12 'runmed_common' 'runmed_common' 1 ((MODULE UNKNOWN-INTENT
+UNKNOWN-PROC UNKNOWN) (UNKNOWN 0 ()) 0 0 () () 0 () ())
+13 'startx_idx' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN) (INTEGER 4 ()) 0 0 () () 0 () ())
+14 'win_k' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN) (INTEGER 4 ()) 0 0 () () 0 () ())
+15 'x' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
+DIMENSION) (REAL 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER 4 ()) 0 '1')
+(CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+16 'win_size' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN) (INTEGER 4 ()) 0 0 () () 0 () ())
+17 'win' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN DIMENSION) (INTEGER 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER
+4 ()) 0 '1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+18 'min_apex' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN) (INTEGER 4 ()) 0 0 () () 0 () ())
+19 'hasnokids' 'runmed_common' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN DIMENSION) (LOGICAL 4 ()) 0 0 () (1 EXPLICIT (CONSTANT (INTEGER
+4 ()) 0 '1') (CONSTANT (INTEGER 4 ()) 0 '81')) 0 () ())
+)
+
+('hasnokids' 0 19 'has2kids' 0 7 'h_child2' 0 6 'h_child1' 0 4 'h' 0 3
+'h0_idx' 0 2 'h_parent' 0 5 'min_apex' 0 18 'max_apex' 0 10 'inmaxheap'
+0 9 'inminheap' 0 8 'win' 0 17 'startx_idx' 0 13 'runmed_common' 0 12
+'shift_h' 0 11 'win_size' 0 16 'win_k' 0 14 'x' 0 15)
diff --git a/RunMedian/runmedianbad_globals.mod b/RunMedian/runmedianbad_globals.mod
new file mode 100644
index 0000000..61b3d17
--- /dev/null
+++ b/RunMedian/runmedianbad_globals.mod
@@ -0,0 +1,27 @@
+GFORTRAN module created from RunMedianBad.f90 on Thu Jun 17 21:22:05 2010
+If you edit this, you'll get what you deserve.
+
+(() () () () () () () () () () () () () () () () () () () () ())
+
+()
+
+()
+
+()
+
+()
+
+(2 'bad_val' 'runmedianbad_globals' 1 ((VARIABLE UNKNOWN-INTENT
+UNKNOWN-PROC UNKNOWN) (REAL 4 ()) 0 0 () () 0 () ())
+3 'max_val' 'runmedianbad_globals' 1 ((VARIABLE UNKNOWN-INTENT
+UNKNOWN-PROC UNKNOWN) (REAL 4 ()) 0 0 () () 0 () ())
+4 'runmedianbad_globals' 'runmedianbad_globals' 1 ((MODULE
+UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN) (UNKNOWN 0 ()) 0 0 () () 0 () ())
+5 'next_bad_use_max' 'runmedianbad_globals' 1 ((VARIABLE UNKNOWN-INTENT
+UNKNOWN-PROC UNKNOWN) (LOGICAL 4 ()) 0 0 () () 0 () ())
+6 'min_val' 'runmedianbad_globals' 1 ((VARIABLE UNKNOWN-INTENT
+UNKNOWN-PROC UNKNOWN) (REAL 4 ()) 0 0 () () 0 () ())
+)
+
+('min_val' 0 6 'max_val' 0 3 'bad_val' 0 2 'next_bad_use_max' 0 5
+'runmedianbad_globals' 0 4)
diff --git a/RunMedian/shell.f90 b/RunMedian/shell.f90
new file mode 100644
index 0000000..1e6484f
--- /dev/null
+++ b/RunMedian/shell.f90
@@ -0,0 +1,45 @@
+
+      SUBROUTINE SHELL (A, N)
+!-----------------------------------------------------------------------
+!     THE SHELL SORTING PROCEDURE IS USED TO REORDER THE ELEMENTS OF A
+!     SO THAT A(I).LE.A(I+1) FOR I=1,...,N-1. IT IS ASSUMED THAT N.GE.1.
+!-----------------------------------------------------------------------
+      REAL A(N)
+      INTEGER K(10)
+!------------------------
+      DATA K(1)/1/, K(2)/4/, K(3)/13/, K(4)/40/, K(5)/121/, K(6)/364/, K(7)/1093/, K(8)/3280/, K(9)/9841/, K(10)/29524/
+!------------------------
+!
+!             SELECTION OF THE INCREMENTS K(I) = (3**I-1)/2
+!
+      IF (N .LT. 2) RETURN
+      IMAX = 1
+      DO 10 I = 3,10
+         IF (N .LE. K(I)) GO TO 20
+         IMAX = IMAX + 1
+   10 CONTINUE
+!
+!            STEPPING THROUGH THE INCREMENTS K(IMAX),...,K(1)
+!
+   20 I = IMAX
+      DO 40 II = 1,IMAX
+         KI = K(I)
+!
+!             SORTING ELEMENTS THAT ARE KI POSITIONS APART
+!               SO THAT A(J).LE.A(J+KI) FOR J=1,...,N-KI
+!
+         JMAX = N - KI
+         DO 31 J = 1,JMAX
+            L = J
+            LL = J + KI
+            S = A(LL)
+   30          IF (S .GE. A(L)) GO TO 31
+               A(LL) = A(L)
+               LL = L
+               L = L - KI
+               IF (L .GT. 0) GO TO 30
+   31    A(LL) = S
+!
+   40 I = I - 1
+      RETURN
+      END
diff --git a/RunMedian/testFilter.f90 b/RunMedian/testFilter.f90
new file mode 100644
index 0000000..3efd486
--- /dev/null
+++ b/RunMedian/testFilter.f90
@@ -0,0 +1,44 @@
+program testFilter
+
+integer n1,n2
+real,allocatable::ain(:,:),aout(:,:)
+real bad_val,max_val,min_val
+integer i,j,isub,jsub
+integer nval
+real vec(100),vec2(100)
+integer k, nf
+bad_val=1.0e08
+min_val=-1.0e08
+max_val=1.0e08
+
+n1=2048
+n2=2048
+
+allocate(ain(n1,n2),aout(n1,n2))
+
+
+
+CALL RANDOM_NUMBER(ain)
+ain=ain*100
+nf=7
+k=(nf-1)/2
+
+call MedianFilter2D(ain,aout,n1,n2,nf,nf,bad_val,min_val,max_val)
+
+i=203
+j=1049
+print *, aout(i,j)
+nval=0
+do isub=i-k,i+k
+   print *, ain(isub,j-3), ain(isub,j-2), ain(isub,j-1), ain(isub,j), ain(isub,j+1), ain(isub,j+2), ain(isub,j+3)
+   do jsub=j-k,j+k
+      nval=nval+1
+      vec(nval)=ain(isub,jsub)
+   enddo
+enddo
+vec2=vec
+call shell(vec2,nval)
+do i=1,nval
+   print *, i, vec(i), vec2(i)
+enddo
+end program testFilter
diff --git a/RunMedian/testmain.f90 b/RunMedian/testmain.f90
new file mode 100644
index 0000000..cfa6d5e
--- /dev/null
+++ b/RunMedian/testmain.f90
@@ -0,0 +1,141 @@
+
+program test
+   use RUNMED_GLOBALS
+   real X0(15)
+   integer H0(15)
+   integer WIN0(15)
+   real vec(15)
+   integer i
+   real median
+
+   real RunMedian
+
+   X0=(/1,5,6,8,9,11,13,15,17,20,25,30,31,33,35/)
+!   X0=(/35,24,19,1,27,8,8,13,19,21,22,11,4,30,29/)
+   WIN0=(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15/)
+   H0=(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15/)
+   verbose=.false.
+   call RunMedianInit(X0,15) !15
+
+print *, 'STARTING POINT:'
+print *, '==============='
+call PrintH
+   vec=0
+   vec(1)=42
+   vec(2)=24
+   vec(3)=99
+   vec(4)=74
+   vec(5)=50
+   vec(6)=56
+   vec(7)=48
+   vec(8)=5
+   vec(9)=3
+   vec(10)=-1
+
+median=RunMedian(vec,10)
+call PrintH
+stop
+   X(max_apex)=2
+call PrintH   
+
+i=max_apex
+!call Shuffle2H0(i)
+!call ShuffleUpMinHeap(i,X(min_apex))
+call ShuffleDownMaxHeap(i,X(max_apex))
+call PrintH 
+print *, i  
+stop
+   med=RunMedian(vec,2)
+   print *, 'med=',med
+   call PrintH
+   print *, '     ---------------'
+   vec(1)=18
+   vec(2)=17
+   med=RunMedian(vec,2)
+   print *, 'med=',med
+   call PrintH
+   print *, '     ---------------'
+   vec(1)=3
+   vec(2)=4
+   med=RunMedian(vec,2)
+   print *, 'med=',med
+   call PrintH
+   print *, '     ---------------'
+   vec(1)=1
+   vec(2)=3
+   med=RunMedian(vec,2)
+   print *, 'med=',med
+   call PrintH
+stop
+
+end program test
+subroutine test_old
+use RUNMED_GLOBALS
+   integer i0,old_i0
+   real X0(15)
+   integer H0(15)
+   integer WIN0(15)
+   win_k=6!7
+   win_size=2*win_k+1
+   shift_h=win_k+1
+   allocate(H(win_size))
+   allocate(WIN(win_size))
+   allocate(X(win_size))
+   X0=(/1,5,6,8,9,11,13,15,17,20,25,30,31,33,35/)
+!   X0=(/35,24,19,1,27,8,8,13,19,21,22,11,4,30,29/)
+   WIN0=(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15/)
+   H0=(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15/)
+   X=X0(1:win_size)
+   WIN=WIN0(1:win_size)
+   H=H0(1:win_size)
+
+print *, 'STARTING POINT:'
+print *, '==============='
+call PrintH
+stop
+do i=win_k,-win_k,-1
+   i0=i
+   call ShuffleHeap(i0)
+   print *, 'step', i
+   call PrintH
+enddo
+stop
+!call PrintH
+
+! Change H(5) to be 14
+i0=5
+X(i0+shift_h)=8
+print *, X(H(i0+shift_h))
+print *, 'STARTING POINT:'
+print *, '==============='
+print *,i0
+call PrintH
+
+call ShuffleHeap(i0)
+call PrintH
+stop
+do i=1,10
+   print *, 'Pushing ', i, 'i0=',i0
+   old_i0=i0
+   call Push2Root(i0)
+   print *,'i0 pushed to ',i0
+   call PrintH
+   if (i0.eq.old_i0) stop
+   if (abs(i0).eq.1) exit
+enddo
+print *, 'Now pushing Root i0=',i0
+call PushRoot(i0)
+print *,'i0 pushed to ',i0
+call PrintH
+do i=1,10
+   print *, 'Pushing ', i, 'i0=',i0
+   old_i0=i0
+   call Push2Leaf(i0)
+   print *,'i0 pushed to ',i0
+   call PrintH
+   if (i0.eq.old_i0) stop
+enddo
+
+
+end subroutine test_old
+
-- 
GitLab