program OW_tstdso

use sep 

implicit none

real                                          :: thresh,o1,o2,o3,d1,d2,d3
real, allocatable                             :: R_image(:,:,:),B_image(:,:,:),h(:),a(:)
integer				      :: n1,n2,n3,ish,ix,iz,i1,i2,i3,wsmoo

 call sep_init()
 call sep_begin_prog()

          call from_history("n1",n1); call from_history("d1",d1); call from_history("o1",o1)  
          call from_history("n2",n2) ; call from_history("d2",d2) ; call from_history("o2",o2)
          call from_history("n3",n3) ; call from_history("d3",d3) ; call from_history("o3",o3)
allocate(R_image(n1,n2,n3));R_image=0.
allocate(B_image(n1,n2,n3));B_image=0.
          call from_param("thresh",thresh,0.5)
          call from_param("wsmoo",wsmoo,21)
!!!!!!!!!!!!!!!!!!!!!!!!!
! thresholding
!!
!!	< const.bkim.hh.H Transp plane=13 | Gpow gpow=2 | Gpow gpow=0.5 | Scale scale_to=1 axis=3 > x1.H
!!	<x1.H Smooth rect1=35 rect2=1 rect3=1|Scale scale_to=1|Clip clip=.5 chop=greater to=1|Clip clip=0.5 chop=less to=0 >x2.H
!!	< const.bkim.hh.H Transp plane=13 > x3.H
!!	Math file1=x2.H file2=x3.H exp="file1*file2" | Transp plane=13 > $@
!
!
!          call scale(R_image)
!          call owsmooth(.false., R_image)
!          call scale(R_image)
!          do i3=1,img%z%n
!          do i2=1,img%th%n
!          do i1=1,img%xm%n
!             if (R_image(i1,i2,i3) > thresh) then
!                 R_image(i1,i2,i3) = 1.
!             else
!                 R_image(i1,i2,i3) = 0.
!             end if
!          end do
!          end do
!          end do
!!          where(R_image > thresh)
!!              R_image = 1.
!!          elsewhere
!!              R_image = 0.
!!          end where
!          B_image=B_image*R_image
!          call owdso(.false. , B_image)!, R_image1)
!          call owdso(.true.  , B_image)!, R_image)

          call to_history("n1",n1,"Dimage") ; call to_history("d1",d1,"Dimage") ; call to_history("o1",o1,"Dimage")  
          call to_history("n2",n2,"Dimage") ; call to_history("d2",d2,"Dimage") ; call to_history("o2",o2,"Dimage")
          call to_history("n3",n3,"Dimage") ; call to_history("d3",d3,"Dimage") ; call to_history("o3",o3,"Dimage")
          call sreed("in",B_image,4*n1*n2*n3)
          R_image=abs(B_image)
! Scaling
          allocate(a(n1))
          a=-9999.
          do i1=1,n1;do i2=1,n3;do i3=1,n2
             if (abs(R_image(i1,i3,i2)) > a(i1)) a(i1)=abs(R_image(i1,i3,i2))
          end do;end do;end do
          do i1=1,n1
             if(a(i1).ne.0)  R_image(i1,:,:)=R_image(i1,:,:)/a(i1)
          end do
          deallocate(a)
! Smoothing
          allocate(a(n3))
          a=0.
          do ish=1,n2;do ix=1,n1
             a=abs(R_image(ix,ish,:))
             do iz=int(wsmoo/2.)+1,n3-int(wsmoo/2.)
                a(iz)=sum(a(iz-int(wsmoo/2.):iz+int(wsmoo/2.)))/wsmoo
             end do
             R_image(ix,ish,:)=a
          end do;end do
          deallocate(a)
! Scaling
          allocate(a(n1))
          a=-9999.
          do i1=1,n1;do i2=1,n3;do i3=1,n2
             if (abs(R_image(i1,i3,i2)) > a(i1)) a(i1)=abs(R_image(i1,i3,i2))
          end do;end do;end do
          do i1=1,n1
             if(a(i1).ne.0)  R_image(i1,:,:)=R_image(i1,:,:)/a(i1)
          end do
          deallocate(a)
! Thresholding
          where(R_image > thresh)
              R_image = 1.
          elsewhere
              R_image = 0.
          end where
! Masking
          B_image=B_image*R_image
! DSO
          allocate(h(n2))
          do ish=1,n2
             h(ish)=(ish-1)*d2+o2
          end do
          h = (abs(h / maxval(h)))
 
          do ish=1,n2;do ix=1,n1
             B_image(ix,ish,:)=B_image(ix,ish,:)*h(ish)**2
          end do;end do

          call srite("Dimage",B_image,4*n1*n2*n3)

 if (allocated(R_image)) deallocate(R_image)
 if (allocated(B_image)) deallocate(B_image)
 call sep_end_prog()
end program
