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



contains

logical function select_pts(pts_field,vel_s,var_s,var,env_s,env,coh_s,coh)
type(sep3d)   :: coh_s,vel_s,var_s,env_s
real, pointer :: pts_field(:,:),temp(:,:),temp3(:,:,:)
real          :: coh(:,:,:),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,br3,er3,nbuf,ibvel(3)
logical,allocatable:: free(:,:,:)


select_pts=.false.
call set_n_o_d_3d(vel_s,nvel,ovel,dvel)
call set_n_o_d_3d(coh_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
nbuf=1000
allocate(temp(3,nbuf))

!variance contains true image space limits
!only allow pts within this regio 
ibvel(1)=max(1,ceiling(((ovar(1)-ovel(1))/dvel(1))+1.))
ibvel(2)=max(1,ceiling(((ovar(2)-ovel(2))/dvel(2))+1.))
nvel(1)=min(nvel(1),floor(((ovar(1)+dvar(1)*(nvar(1)-1)-ovel(1))/dvel(1))+1.))
nvel(2)=min(nvel(2),floor(((ovar(2)+dvar(2)*(nvar(2)-1)-ovel(2))/dvel(2))+1.))
if(nvel(3) >1) then
  ibvel(3)=max(1,ceiling(((ovar(3)-ovel(3))/dvel(3))+1.))
  nvel(3)=min(nvel(3),floor(((ovar(3)+dvar(3)*(nvar(3)-1)-ovel(2))/dvel(2))+1.))
else 
  ibvel(3)=1
end if



!LOOP IN SECTIONS
do i=1,nloop_pts_select
  !I AM X PCT WAY DONE

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

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

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


  do i3=ibvel(3),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=ibvel(2),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=ibvel(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(coh(ic1,ic2,ic3) > my_coh .and. var(iv1,iv2,iv3) < my_var &
          .and. env(ie1,ie2,ie3) > my_env .and. free(i1,i2,i3)) then

          br3=max(1,i3-min_dist(3))
          er3=min(nvel(3),i3+min_dist(3))
          br2=max(1,i2-min_dist(2))
          er2=min(nvel(2),i2+min_dist(2))
          br1=max(1,i1-min_dist(1))
          er1=min(nvel(1),i1+min_dist(1))
          free(br1:er1,br2:er2,br3:er3)=.false.
          npts=npts+1
          temp(1:3,npts)=(/pos(1),pos(2),pos(3)/)
          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
      end do
    end do
    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
  allocate(pts_field(3,npts))
  pts_field=temp(:,1:npts)
  deallocate(temp)
  deallocate(free)
  select_pts=.true.
end function

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

 parse_pts=.false.
 allocate(temp(3,size(pts,2)))
 npts=0
 do  i=1,size(pts,2)
   good=.true.
   do ireg=1,nkill_reg
     if(pts(1,i) >= beg1(ireg)  .and. pts(1,i) <= end1(ireg)) then
      if(pts(2,i) >= beg2(ireg)  .and. pts(2,i) <= end2(ireg)) then
         if(pts(3,i) >= beg3(ireg)  .and. pts(3,i) <= end3(ireg))  good=.false.
       end if
     end if
   end do
   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)=min(n(1),max(1,nint((locs(1,i)-o(1))/d(1)+1.5)))
  iloc(2)=min(n(2),max(1,nint((locs(2,i)-o(2))/d(2)+1.5)))
  iloc(3)=min(n(3),max(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

logical function grab_pts_select_pars()
  real :: first,last,pct(20)
  integer :: i
  grab_pts_select_pars=.false.
  call from_param("verb",verb,.false.)
  call from_param("nloop_pts_select",nloop_pts_select,5)
  do i=1,nloop_pts_select; pct(i)=real(i-1)/real(nloop_pts_select);end do
  call from_param("min_dist",min_dist,(/4,4,4/))
  call from_param("first_image_pct",first,85.)
  call from_param("last_image_pct",last,40.)
  do i=1,nloop_pts_select 
    image_pct(i)=(1.-pct(i))*first+pct(i)*last
  end do
  call from_param("first_dip_coh_pct",first,80.)
  call from_param("last_dip_coh_pct",last,30.)
  do i=1,nloop_pts_select 
    coh_pct(i)=(1.-pct(i))*first+pct(i)*last
  end do
  call from_param("first_var_pct",first,20.)
  call from_param("last_var_pct",last,90.)
  do i=1,nloop_pts_select 
    var_pct(i)=(1.-pct(i))*first+pct(i)*last
  end do
  grab_pts_select_pars=.true.
end function
                                                                                
                                                                                
logical function grab_parse_pts_pars(nvel,ovel,dvel)
integer :: nvel(3)
real    :: ovel(3)
real    :: dvel(3)
integer :: i
logical :: found
real    :: default(20)
                                                                                
  grab_parse_pts_pars=.false.
  default=-999999
  beg1=default
  end1=default
  beg2=default
  end2=default
  beg3=default
  end3=default
  call from_param("verb",verb,.false.)
  call from_param("beg1_kill",beg1,default)
  call from_param("beg2_kill",beg2,default)
  call from_param("beg3_kill",beg3,default)
  call from_param("end1_kill",end1,default)
  call from_param("end2_kill",end2,default)
  call from_param("end3_kill",end3,default)
  found=.false.
  i=1
  nkill_reg=0
  do while (.not. found .and. i < 20)
    if (beg1(i) /=default(i) .or.  beg2(i) /=default(i) .or.  &
      beg3(i) /=default(i) .or.  end1(i) /=default(i) .or. &
      end2(i) /=default(i) .or.  end3(i) /=default(i) ) then
      nkill_reg=nkill_reg+1
      if(beg1(i)==default(i)) beg1(i)=ovel(1)-dvel(1)
      if(end1(i)==default(i)) end1(i)=ovel(1)+dvel(1)*nvel(1)
      if(beg2(i)==default(i)) beg2(i)=ovel(2)-dvel(2)
      if(end2(i)==default(i)) end2(i)=ovel(2)+dvel(2)*nvel(2)
      if(beg3(i)==default(i)) beg3(i)=ovel(3)-dvel(3)
      if(end3(i)==default(i)) end3(i)=ovel(3)+dvel(3)*nvel(3)
      if(verb) write(0,*) "kill region axis1 ",i,beg1(i),end1(i)
      if(verb) write(0,*) "kill region axis2",i,beg2(i),end2(i)
      if(verb) write(0,*) "kill region axis3",i,beg3(i),end3(i)
      i=i+1
    else
      found=.true.
    end if
  end do
  grab_parse_pts_pars=.true.
end function


end module
