# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<;
!dip_field
!
!Usage:
!dip_field.x <in.H >out.H
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author: Robert Clapp, ESMB 463, 7230253
!
!Date Created:Mon May 21 15:10:48 PDT 2001
!
!Purpose:
!
! Module for parsing a set of pts
!
module parse_region_mod
  use basic_util_mod
  use quantile2_mod
  use param_check_mod
  use sort_subs_mod
  implicit none
  real, private, save :: last_image_pct,last_dip_pct,last_var_pct
  real, private, save :: first_image_pct,first_dip_pct,first_var_pct
  real, private, save :: end1(20),end2(20),end3(20)
  real, private, save :: beg1(20),beg2(20),beg3(20)
  integer, private, save :: min_pt_dist(3)
  integer, private, save :: iref_use(30),nloop_pts_select
  logical, private, save :: verb
  interface select_ref_locations
    module procedure select_ref_locations_2d
    module procedure select_ref_locations_3d
  end interface
  contains
  logical function select_pt_locations(output_space,image,image_space&
    &,coher,coher_space,var,var_space,pts_field)
    type(sep3d) :: output_space,image_space,coher_space,var_space
    real :: image(:,:,:),coher(:,:,:),var(:,:,:)
    real,pointer:: pts_field(:,:),temp(:,:),locs(:,:)
    logical,pointer:: have_pt(:,:,:)
    integer :: io(3),ii(3),ic(3),iv(3),io1,io2,io3
    integer :: no(3),ni(3),nc(3),nv(3)
    real :: oo(3),oi(3),oc(3),ov(3)
    real :: do(3),di(3),dc(3),dv(3)
    integer :: nbuf,npts,mind(3)
    integer :: i1,i,b(3),e(3)
    integer :: iloop_beg(3),iloop_end(3)
    integer :: ia,ib,id,ie,nb(3),nreal(3),cyclera
    real :: pct,my_coh,my_amp,loc(3),my_var
    select_pt_locations=.false.
    call srite("aa1.H",image,size(coher)*4)
!call srite("aa1.H",coher,size(coher)*4)
    call set_n_o_d_3d(output_space,no,oo,do)
    call set_n_o_d_3d(coher_space,nc,oc,dc)
    call set_n_o_d_3d(image_space,ni,oi,di)
    call set_n_o_d_3d(var_space,nv,ov,dv)
    do i=1,3
      nreal(i)=dv(i)*nv(i)/di(i)
    end do
    if (nv(3).eq.1) then
      nreal(3)=1
    end if
!write(0,*) "CHECK ov",ov
!write(0,*) "CHECK dv",dv
!write(0,*) "CHECK nv",nv
!write(0,*) "CHECK oi",oi
!write(0,*) "CHECK di",di
!write(0,*) "CHECK ni",ni
!write(0,*) "CHECK 1",nreal
!write(0,*) "CHECK 2",nv
!write(0,*) "CHECK 3",ni
!write(0,*) "CHECK 4",nc
    cyclera=product(ni)-product(nreal)
    mind=min_pt_dist
    nbuf=1000
    npts=0
    allocate(pts_field(4,1000))
    allocate(have_pt(no(1),no(2),no(3)))
    have_pt=.false.
!FIGURE OUT THE MINIMUM AREA OF THE 3 RELEVANT CUBES
    ii=(oi-oo)/do+1.501
    ic=(oc-oo)/do+1.501
    iv=(ov-oo)/do+1.501
    iloop_beg(1)=maxval((/1,ii(1),ic(1),iv(1)/))
    iloop_beg(2)=maxval((/1,ii(2),ic(2),iv(2)/))
    iloop_beg(3)=maxval((/1,ii(3),ic(3),iv(3)/))
    ii=(oi+di*(ni-1)-oo)/do+1.499
    ic=(oc+dc*(nc-1)-oo)/do+1.499
    iv=(ov+dv*(nv-1)-oo)/do+1.499
    iloop_end(1)=minval((/no(1),ii(1),ic(1),iv(1)/))
    iloop_end(2)=minval((/no(2),ii(2),ic(2),iv(2)/))
    iloop_end(3)=minval((/no(3),ii(3),ic(3),iv(3)/))
    iloop_end(3)=max(iloop_end(3),1)
!LOOP IN SECTIONS
    do i=1,nloop_pts_select
!I AM X PCT WAY DONE
      pct=real(i-1.)/real(nloop_pts_select-1.)
      ib=0
      ia=0
      id=0
      ie=0
!ACCEPT COHERENCE AT THIS LEVEL
      my_coh=(1.-pct)*first_dip_pct + pct *last_dip_pct
      i1=(size(coher)-cyclera)*my_coh/100+cyclera
      my_coh=quantile(i1,reshape(coher,(/size(coher)/)))
!ACCEPT IMAGE AMPLITUDE AT THIS LEVEL
      my_amp=(1.-pct)*first_image_pct + pct *last_image_pct
      i1=(size(image)-cyclera)*my_amp/100+cyclera
      my_amp=quantile(i1,reshape(image,(/size(image)/)))
!ACCEPT VARIANCE AMPLITUDE AT THIS LEVEL
      my_var=(1.-pct)*first_var_pct + pct *last_var_pct
      i1=size(var)*my_var/100
      my_var=quantile(i1,reshape(var,(/size(var)/)))
! write(0,*) "pct C A V",pct,(1.-pct)*first_dip_pct + pct *last_dip_pct,
!(1.-pct)*first_image_pct + pct *last_image_pct,(1.-pct)*first_var_pct + pct *last_var_pct
! write(0,*) "pct C A V",my_coh,my_amp,my_var,nextra
!LOOP OVER OUTPUT SPACE
      do io3=iloop_beg(3),iloop_end(3)
        io(3)=io3
        loc(3)=(io(3)-1)*do(3)+oo(3)
        ii(3)=nint((loc(3)-oi(3))/di(3)+1)
        ic(3)=nint((loc(3)-oc(3))/dc(3)+1)
        iv(3)=nint((loc(3)-ov(3))/dv(3)+1)
        do io2=iloop_beg(2),iloop_end(2)
          io(2)=io2
          loc(2)=(io(2)-1)*do(2)+oo(2)
          ii(2)=nint((loc(2)-oi(2))/di(2)+1)
          ic(2)=nint((loc(2)-oc(2))/dc(2)+1)
          iv(2)=nint((loc(2)-ov(2))/dv(2)+1)
          do io1=iloop_beg(1),iloop_end(1)
            io(1)=io1
            loc(1)=(io(1)-1)*do(1)+oo(1)
            ii(1)=nint((loc(1)-oi(1))/di(1)+1)
            ic(1)=nint((loc(1)-oc(1))/dc(1)+1)
            iv(1)=nint((loc(1)-ov(1))/dv(1)+1)
! if(coher(ic(1),ic(2),ic(3))>my_coh ) ia+=1
! if(image(ii(1),ii(2),ii(3)) > my_amp) ib+=1
! if(var(iv(1),iv(2),iv(3)) < my_var ) id+=1
! if(.not.have_pt(io1,io2,io3)) ie+=1
            if (coher(ic(1),ic(2),ic(3))>my_coh .and. image(ii(1),ii&
              &(2),ii(3)) > my_amp .and. var(iv(1),iv(2),iv(3)) <&
              & my_var .and. .not. have_pt(io1,io2,io3)) then
              b(3)=max(1,io(3)-mind(3))
              e(3)=min(no(3),io(3)+mind(3))
              b(2)=max(1,io(2)-mind(2))
              e(2)=min(no(2),io(2)+mind(2))
              b(1)=max(1,io(1)-mind(1))
              e(1)=min(no(1),io(1)+mind(1))
              have_pt(b(1):e(1),b(2):e(2),b(3):e(3))=.true.
              npts= npts+1
              pts_field(1:3,npts)=(/loc(2),loc(3),loc(1)/)
              if (npts+1.eq.size(pts_field,2)) then
                allocate(temp(size(pts_field,1),nbuf))
                temp=pts_field
                deallocate(pts_field)
                nbuf=nbuf*1.5
                allocate(pts_field(size(temp,1),nbuf))
                pts_field(:,:size(temp,2))=temp
                deallocate(temp)
!write(0,*) "CHECK LOC",npts,pts_field(:,1:2)
              end if
            end if
          end do
        end do
      end do
      if (verb) then
        write(0,*) "      Have ",npts, " on pass ",i
      end if
    end do
    if (.not. parse_sort_points(pts_field,npts,no,oo,do)) then
      write(0,*) "trouble parsing and sorting points"
      return
    end if
    call srite("pts.H",pts_field,16*npts)
    if (npts.eq.0) then
      write(0,*) "Didn't find  any acceptable back projection pts"
      return
    end if
    select_pt_locations=.true.
  end function
  logical function parse_sort_points(pts,npts,ns,os,ds)
    real, pointer :: pts(:,:),temp(:,:)
    integer, pointer ::ilocs(:,:)
    real :: os(3),ds(3)
    integer :: ns(3),npts,i2
    if (verb) then
      write(0,*) "      Before removing points from selected regions"
    end if
    if (.not. parse_regioned_pts(pts,npts,ns,os,ds)) then
      write(0,*) "trouble parsing points"
      return
    end if
!SORT THE POINTS
    allocate(ilocs(4,npts))
    do i2=1,npts
      ilocs(:,i2) =(/nint(((/pts(3,i2),pts(1:2,i2)/)-(/os(1),os(2),os(3&
        &)/))/(/ds(1),ds(2),ds(3)/) +1.),i2/)
    end do
    call sort_it(ilocs,3,npts)
    allocate(temp(4,npts))
    do i2=1,npts
      temp(:,i2) =(/ ((/ilocs(2:3,i2),ilocs(1,i2)/)-1)*(/ds(2:3),ds(1)&
        &/)+(/os(2:3),os(1)/),pts(4,i2)/)
    end do
    deallocate(ilocs)
    allocate(pts(size(temp,1),size(temp,2)))
    pts=temp
    deallocate(temp)
    parse_sort_points=.true.
  end function
  logical function parse_regioned_pts(pts_field,npts,ns,os,ds)
    integer :: npts,ipts,nout,ireg
    integer :: nreg,ns(3),i(6),i2,ia
    real :: os(3),ds(3),beg(3),end(3)
    logical :: endit
    real, pointer :: pts_field(:,:),temp(:,:)
    logical,allocatable :: kill_pts(:)
    integer, external :: getch
    parse_regioned_pts=.false.
    beg=os-ds
    end=os+ds*(ns)
    i2=1
    endit=.false.
    do while (i2 < size(beg1) .and. .not. endit)
      endit=.true.
      if (beg1(i2).eq.default_real) then
        beg1(i2)=beg(1)
      else
        endit=.false.
      end if
      if (beg2(i2).eq.default_real) then
        beg2(i2)=beg(2)
      else
        endit=.false.
      end if
      if (beg3(i2).eq.default_real) then
        beg3(i2)=beg(3)
      else
        endit=.false.
      end if
      if (end1(i2).eq.default_real) then
        end1(i2)=end(1)
      else
        endit=.false.
      end if
      if (end2(i2).eq.default_real) then
        end1(i2)=end(2)
      else
        endit=.false.
      end if
      if (end3(i2).eq.default_real) then
        end1(i2)=end(3)
      else
        endit=.false.
      end if
      if (.not. endit) then
        write(0,*) "REGION",beg1(i2),end1(i2),beg2(i2),end2(i2),beg3(i2&
          &),end3(i2)
      end if
      i2=i2+1
    end do
    nreg=i2-1
    write(0,*) "NREG",nreg
    nout=npts
    allocate(kill_pts(npts))
    kill_pts=.false.
    do ipts=1,npts
!write(0,*) "DIE 2",size(pts_field,2),npts
      do ireg=1,nreg
! write(0,*) "LOCATION",pts_field(:,ipts)
        if (pts_field(1,ipts) > beg2(ireg) .and. pts_field(1,ipts) <&
          & end2(ireg) .and. pts_field(2,ipts) > beg3(ireg) .and.&
          & pts_field(2,ipts) < end3(ireg) .and. pts_field(3,ipts) >&
          & beg1(ireg) .and. pts_field(3,ipts) < end1(ireg)) then
          kill_pts(ipts)=.true.
          nout= nout-1
        end if
      end do
    end do
    write(0,*) "2A",size(kill_pts),size(temp,2),size(pts_field,2),size&
      &(temp,1),size(pts_field,1),npts
    allocate(temp(size(pts_field,1),size(pts_field,2)))
    do ipts=1,npts
      do ia=1,size(pts_field,1)
        temp(ia,ipts)=pts_field(ia,ipts)
      end do
    end do
    deallocate(pts_field)
    allocate(pts_field(size(temp,1),nout))
    write(0,*) "AA",size(kill_pts),size(temp,2),size(pts_field,2),size&
      &(temp,1),size(pts_field,1),npts,nout
    nout=0
    do ipts=1,npts
      if (.not. kill_pts(ipts) .and. nout < size(pts_field,2)) then
        nout= nout+1
        do ia=1,size(pts_field,1)
          write(0,*) "CHEC THIS",ia,ipts,nout,size(pts_field,2),size&
            &(temp,2),size(temp,1),size(pts_field,1)
          pts_field(ia,nout)=temp(ia,ipts)
        end do
      end if
    end do
    write(0,*) "A2",size(kill_pts)
    deallocate(kill_pts)
    write(0,*) "A3",size(temp)
    deallocate(temp)
    write(0,*) "A4"
    if (verb) then
      write(0,*) "        After parsing npts=",nout
    end if
    if (nout.eq.0) then
      write(0,*) "        ERROR: eliminated all pts wint&
        & kill_pt_region"
      return
    end if
!write(0,*) "A5"
    npts=nout
    parse_regioned_pts=.true.
  end function
  logical function select_ref_locations_3d(out_s,ref_s,refs,points)
    type(sep3d) :: out_s,ref_s
    real :: refs(:,:,:)
    real, pointer :: points(:,:),temp(:,:)
    integer :: i1,iref,ipts,nt,nuse_ref,npts,ns(3),i2
    real :: x,z,mymin(3),mymax(3),os(3),ds(3),y
    select_ref_locations_3d=.false.
    mymin=out_s%o
    mymax=out_s%o+out_s%d*(out_s%n-1)
    call set_n_o_d_3d(out_s,ns,os,ds)
    nt=500
    npts=0
    allocate(points(4,nt))
    if (setup_iref(size(refs,3),nuse_ref)) then
      write(0,*) "troublet setting up iref_use"
      return
    end if
    do iref=1,size(refs,3)
      do i2=1,size(refs,2)
        y=ref_s%o(2)+ref_s%d(2)*(ref_s%n(2)-1)
        do i1=1,size(refs,1)
          x=ref_s%o(1)+ref_s%d(1)*(ref_s%n(1)-1)
          z=refs(i1,i2,iref)
          if (x >= mymin(2) .and. x<=mymax(2) .and. z >= mymin(1) .and.&
            & z <= mymax(1) .and. y >= mymin(3) .and. y<=mymax(3)) then
            if (any(x-points(1,:npts) < min_pt_dist(2)) .or. any(y&
              &-points(2,:npts) < min_pt_dist(3)) .or. any(z-points&
              &(3,:npts) < min_pt_dist(1))) then
!
            else
              if (npts.eq.nt) then
                nt=nt*1.5
                allocate(temp(size(points,1),npts))
                temp=points
                deallocate(points)
                allocate(points(size(temp,1),nt))
                points(:,1:npts)=temp
                deallocate(temp)
              end if
              npts=npts+1
              points(:,npts)=(/x,y,z,1.*iref/)
            end if
          end if
        end do
      end do
    end do
    if (npts .eq.0) then
      write(0,*) "found no valid pts before parse"
      return
    end if
    if (verb) then
      write(0,*) "found ",npts, "valid points"
    end if
    if (.not. parse_sort_points(points,npts,ns,os,ds)) then
      write(0,*) "        ERROR: parsing and sorting points"
      return
    end if
    select_ref_locations_3d=.true.
  end function
  logical function select_ref_locations_2d(out_s,refs,points)
    type(sep3d) :: out_s
    complex :: refs(:,:)
    real, pointer :: points(:,:),temp(:,:)
    integer :: i1,iref,ipts,nt,nuse_ref,npts,ns(3)
    real :: x,z,mymin(2),mymax(2),os(3),ds(3)
    select_ref_locations_2d=.false.
    mymin=out_s%o
    mymax=out_s%o+out_s%d*(out_s%n-1)
    call set_n_o_d_3d(out_s,ns,os,ds)
    nt=500
    npts=0
    allocate(points(4,nt))
    if (setup_iref(size(refs,2),nuse_ref)) then
      write(0,*) "troublet setting up iref_use"
      return
    end if
    do iref=1,size(refs,2)
      do i1=1,size(refs,1)
        x=real(refs(i1,iref))
        z=aimag(refs(i1,iref))
      end do
      if (x >= mymin(2) .and. x<=mymax(2) .and. z >= mymin(1) .and. z&
        & <= mymax(1)) then
        if (any(x-points(1,:npts) < min_pt_dist(2)) .or. any(z-points&
          &(3,:npts) < min_pt_dist(1))) then
!
        else
          if (npts.eq.nt) then
            nt=nt*1.5
            allocate(temp(size(points,1),npts))
            temp=points
            deallocate(points)
            allocate(points(size(temp,1),nt))
            points(:,1:npts)=temp
            deallocate(temp)
          end if
          npts=npts+1
          points(:,npts)=(/x,0.,z,iref*1./)
        end if
      end if
    end do
    if (npts .eq.0) then
      write(0,*) "found no valid pts before parse"
      return
    end if
    if (verb) then
      write(0,*) "found ",npts, "valid points"
    end if
    if (.not. parse_sort_points(points,npts,ns,os,ds)) then
      write(0,*) "trouble parsing and sorting points"
      return
    end if
    select_ref_locations_2d=.true.
  end function
  logical function setup_iref(nref,nuse_ref)
    integer :: nuse_ref
    integer nref,iref
    setup_iref=.false.
    if (all(iref_use.eq.default_int)) then
      do iref=1,nref
        iref_use(iref)=iref
      end do
      nuse_ref=nref
    else
      nuse_ref=1
      do while (iref_use(nuse_ref) .ne. default_int)
        nuse_ref= nuse_ref+1
      end do
      nuse_ref= nuse_ref-1
      if (nuse_ref.eq.0) then
        write(0,*) "found now valid iref_use"
        return
      end if
      do iref=1,nuse_ref
        if (iref_use(iref) < 1 .or. iref_use(iref) >nref) then
          write(0,*) "invalid iref_use",iref_use(iref)," must be between&
            & 1 and ",nref
          return
        end if
      end do
    end if
    setup_iref=.true.
  end function
  logical function check_pars_parse_region(tag,norerun)
    character(len=*) :: tag
    logical :: stat_ok(20) ,stat_rerun(20),norerun
    check_pars_parse_region=.false.
    stat_rerun=.true.
    stat_ok=.true.
    beg1=default_real
    beg2=default_real
    beg3=default_real
    end1=default_real
    end2=default_real
    end3=default_real
    iref_use=default_int
    stat_rerun(1)= check_par_file("kill_region_beg1",beg1,tag,beg1)
    stat_rerun(2)= check_par_file("kill_region_beg2",beg2,tag,beg2)
    stat_rerun(3)= check_par_file("kill_region_beg3",beg3,tag,beg3)
    stat_rerun(4)= check_par_file("kill_region_end1",end1,tag,end1)
    stat_rerun(5)= check_par_file("kill_region_end2",end2,tag,end2)
    stat_rerun(6)= check_par_file("kill_region_end3",end3,tag,end3)
    stat_rerun(7)=check_par_file("first_image_pct" ,first_image_pct,tag&
      &,90.)
    stat_rerun(8)=check_par_file("last_image_pct" ,last_image_pct,tag&
      &,62.)
    call from_param("first_image_pct",first_image_pct,90.)
    stat_rerun(9)=check_par_file("first_dip_coh_pct" ,first_dip_pct,tag&
      &,88.)
    stat_rerun(10)=check_par_file("last_dip_coh_pct" ,last_dip_pct,tag&
      &,55.)
    stat_rerun(11)=check_par_file("first_var_pct" ,first_var_pct,tag&
      &,90.)
    stat_rerun(12)=check_par_file("last_var_pct" ,last_var_pct,tag,40.&
      &)
    stat_rerun(13)=check_par_file("nloop_pts_select" ,nloop_pts_select&
      &,tag,4)
    stat_rerun(14)=check_par_file("min_pt_dist" ,min_pt_dist, tag,(&
      &/5,5,5/))
    stat_rerun(15)= check_par_file("iref_use",iref_use,tag,iref_use)
    call from_param("verb",verb,.false.)
    if (all(stat_ok)) then
      check_pars_parse_region=.true.
    end if
    if (all(stat_rerun)) then
      norerun=.true.
    else
      norerun=.false.
    end if
  end function
end module
