#<;

module pts_select_mod{ 
use quantile2_mod
use sep3d_struct_mod
use basic_util_mod
implicit none
logical, private, save :: verb
real, private, save    :: beg1(99),beg2(99),beg3(99)
real, private, save    :: end1(99),end2(99),end3(99)
integer,private,save   :: nkill_reg,nloop_pts_select,min_dist(3)
real, private, save    :: image_pct(99),coh_pct(99),var_pct(99)



contains

logical function select_pts(pts_field,vel_s,var_s,var,env_s,env,coh_s,coh)
type(sep3d)   :: coher_s,vel_s
real, pointer :: pts_field(:,:),temp(:,:)
real          :: coher(:,:,:),var(:,:,:),env(:,:,:)
integer       :: ncoh(3),nvar(3),nenv(3),nvel(3)
real          :: ocoh(3),dcoh(3),ovar(3),dvar(3),oenv(3),denv(3),ovel(3),dvel(3)
integer       :: i,i1,i2,i3
real          :: my_coh,my_var,my_env,pos(3)
integer       :: iv1,iv2,iv3,ie1,ie2,ie3,ic1,ic2,ic3,npts
integer       :: br1,br2,er1,er2,b3,er3
logical,allocatable:: free(:,:,:)


select_pts=.false.
call set_n_o_d_3d(vel_s,nvel,ovel,dvel)
call set_n_o_d_3d(coher_s,ncoh,ocoh,dcoh)
call set_n_o_d_3d(var_s,nvar,ovar,dvar)
call set_n_o_d_3d(env_s,nenv,oenv,denv)
allocate(free(nvel(1),nvel(2),nvel(3)))
free=.true.
npts=0
allocate(temp(3,1000))


#LOOP IN SECTIONS
do i=1,nloop_pts_select
  #I AM X PCT WAY DONE
  ib=0;ia=0; id=0;ie=0

  #ACCEPT COHERENCE AT THIS LEVEL
  i1=min(size(coher),max(1,coh_pct(i)*size(coher)/100.))
  my_coh=quantile(i1,reshape(coher,(/size(coher)/)))

  #ACCEPT IMAGE AMPLITUDE AT THIS LEVEL
  i1=min(size(env),max(1,env_pct(i)*size(env)/100.))
  my_env=quantile(i1,reshape(env,(/size(env)/)))

  #ACCEPT VARIANCE AMPLITUDE AT THIS LEVEL
  i1=min(size(coher),max(1,var_pct(i)*size(var)/100.))
  my_var=quantile(i1,reshape(var,(/size(var)/)))


  do i3=1,nvel(3)
    pos(3)=(i3-1)*dvel(3)+ovel(3)
    iv3=min(nvar(3),max(1,nint((pos(3)-ovar(3))/dvar(3)+1.5)))
    ic3=min(ncoh(3),max(1,nint((pos(3)-ocoh(3))/dcoh(3)+1.5)))
    ie3=min(nenv(3),max(1,nint((pos(3)-oenv(3))/denv(3)+1.5)))
  
    do i2=1,nvel(2)
      pos(2)=(i2-1)*dvel(2)+ovel(2)
      iv2=min(nvar(2),max(1,nint((pos(2)-ovar(2))/dvar(2)+1.5)))
      ic2=min(ncoh(2),max(1,nint((pos(2)-ocoh(2))/dcoh(2)+1.5)))
      ie2=min(nenv(2),max(1,nint((pos(2)-oenv(2))/denv(2)+1.5)))
      do i1=1,nvel(1)
        pos(1)=(i1-1)*dvel(1)+ovel(1)
        iv1=min(nvar(1),max(1,nint((pos(1)-ovar(1))/dvar(1)+1.5)))
        ic1=min(ncoh(1),max(1,nint((pos(1)-ocoh(1))/dcoh(1)+1.5)))
        ie1=min(nenv(1),max(1,nint((pos(1)-oenv(1))/denv(1)+1.5)))

        if(coher(ic1,ic2,ic3) > my_coher .and. var(iv1,iv2,iv3) < my_var &
          .and. env(ie1,ie,i3) > my_env .and. free(i1,i2,i3)) then

          br3=max(1,i3-max_pt_dist(3))
          er3=min(nvel(3),i3+max_pt_dist(3))
          br2=max(1,i2-max_pt_dist(2))
          er2=min(nvel(2),i2+max_pt_dist(2))
          br1=max(1,i1-max_pt_dist(1))
          er1=min(nvel(1),i1+max_pt_dist(1))
          free(br1:er1,br2:er2,br3:er3)=.false.
          temp(1:3,npts)=(/pos(1),pos(2),pos(3)/)
          npts=npts+1
          if(npts+1==size(temp,2)) then
             allocate(pts_field(size(temp,1),nbuf)); pts_field=temp
             deallocate(temp); nbuf=nbuf*1.5
             allocate(temp(size(pts_field,1),nbuf)); 
             temp(:,:size(pts_field,2))=pts_field
             deallocate(pts_field)
          end if
        end if
        if (verb) write(0,*) "loop coh, var, env, npts"
        if (verb) write(0,*) i,my_coh,my_var,my_env,npts
      end do
    end do
  end do
  allocate(pts_field(3,npts))
  pts_field=temp(:,1:npts)
  deallocate(temp)
  select_pts=.true.
end function

logical function parse_pts(pts)
 real, pointer :: pts(:,:),temp(:,:)
 integer :: i,npts
 logical :: good

 parse_pts=.false.
 allocate(temp(3,size(pts,2)))
 npts=0
 do  i=1,size(pts,2)
   good=.true.
   if(pts(1,i) >= beg(1)  .and. pts(1,i) <= end(1)) then
     if(pts(2,i) >= beg(2)  .and. pts(2,i) <= end(2)) then
       if(pts(3,i) >= beg(3)  .and. pts(3,i) <= end(3))  good=.false.
     end if
   end if
   if(good) then 
     npts=npts+1
     temp(:,npts)=pts(:,i)
   end if
 end do
 deallocate(pts)
 allocate(pts(3,npts))
 pts=temp(:,:npts)
 parse_pts=.true.
end function
   
subroutine extract_vals(tag,locs,vals){
integer :: n(3)
real    :: o(3),d(3),locs(:,:),vals(:)
character(len=*) :: tag
type(sep3d) :: struct
integer     :: i,iloc(3)
integer,external :: sreed
real, allocatable :: buf(:,:,:)

call init_sep3d(tag,struct,"INPUT")
call set_n_o_d_3d(struct,n,o,d)
allocate(buf(n(1),n(2),n(3)))
if (4*size(buf)/=sreed(tag,buf,size(buf)*4)) &
  call srite("trouble reading in buffer")

do i=1,size(vals)
  iloc(1)=max(n(1),min(1,nint((locs(1,i)-o(1))/d(1)+1.5)))
  iloc(2)=max(n(2),min(1,nint((locs(2,i)-o(2))/d(2)+1.5)))
  iloc(3)=max(n(3),min(1,nint((locs(3,i)-o(3))/d(3)+1.5)))
  vals(i)=buf(iloc(1),iloc(2),iloc(3))
end do
  
deallocate(buf)

end subroutine

end module
