! makeXcorr compute the correlation functions
! Copyright (C) 1999,2010 Gauthier Delerce
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2, or (at your option)
! any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software Foundation,
! Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
!
! Gauthier Delerce <gauthier@delerce.fr>
subroutine make_Xcorrelation
use mod_parameters
!use modvsv
!use modimage
!use modgrille
!use modcorr
implicit none
integer,parameter::border=0
integer::i,j,ii,jj,nbnodes,k,j1
real::cm,cr
integer,dimension(ib(1,1),ib(1,2))::at
integer,dimension(is(1,1),is(1,2))::bt
real,dimension(is(1,1)-ib(1,1)+1,is(1,2)-ib(1,2)+1)::ct
real,dimension(:),allocatable::avg,arms
real,dimension(:,:,:),allocatable::bvg,brms
integer,dimension(:),allocatable::it,jt,tot,iit,jjt
logical,dimension(:,:),allocatable::mask
!*** FFT stuff
real(8),dimension(:,:),allocatable::ma,a2,ma2,b2,iffmab,iffmb,iffmb2,p
complex(8),dimension(:,:),allocatable::fma,fb_ron,fm_ron,fb2,ffmab,ffmb,ffmb2
real(8),dimension(:,:),allocatable::a_ron,b_ron,m_ron
!INTEGER FORWARD,BACKWARD,FFTW_PATIENT,FFTW_ESTIMATE
integer,parameter::FORWARD=-1,BACKWARD=1,FFTW_PATIENT=32,FFTW_ESTIMATE=64
integer,dimension(2)::maxfft,maxcor
integer(8):: plan,plan_ma,plan_m_ron,plan_b_ron,plan_b2,plan_ffmab,plan_ffmb,plan_ffmb2
integer::nxf,nyf,xm,ym
real::maxc
real(8)::sum_ma,sum_ma2,sum2_ma
real(4),dimension(:,:),allocatable::r_tmp
!***
print*,'Begin XCORR',mg*ng,' nodes'
allocate(ncor(mg*ng,2))
ncor(:,1)=(is(:,1)-ib(:,1)+1)
ncor(:,2)=(is(:,2)-ib(:,2)+1)
!print*,' ncor : ',ncor
nbnodes=mg*ng
allocate(corr_tab(nbnodes,maxval(ncor(:,1)),maxval(ncor(:,2))))
if(test_seuil) allocate(corr_tab2(nbnodes,maxval(ncor(:,1)),maxval(ncor(:,2))))
allocate(it(nbnodes),jt(nbnodes),tot(nbnodes),avg(nbnodes),arms(nbnodes))
allocate(bvg(nbnodes,maxval(ncor(:,1)),maxval(ncor(:,2))),brms(nbnodes,maxval(ncor(:,1)),maxval(ncor(:,2))))
allocate(mask(-border+1:nx+border,-border+1:ny+border))
allocate(iit(maxval(ncor(:,1))),jjt(maxval(ncor(:,2))))
!*** FFT stuff
nxf=is(1,1)+1
nyf=is(1,2)+1
allocate(a_ron(nxf,nyf),b_ron(nxf,nyf),m_ron(nxf,nyf),r_tmp(nxf,nyf))
allocate(ma(nxf,nyf),a2(nxf,nyf),ma2(nxf,nyf),b2(nxf,nyf))
allocate(iffmab(nxf,nyf),iffmb(nxf,nyf),iffmb2(nxf,nyf),p(nxf,nyf),R(nbnodes,nxf,nyf))
allocate(fma(nxf/2+1,nyf),fb_ron(nxf/2+1,nyf),fm_ron(nxf/2+1,nyf))
allocate(fb2(nxf/2+1,nyf),ffmab(nxf/2+1,nyf),ffmb(nxf/2+1,nyf),ffmb2(nxf/2+1,nyf))
!****
mask=.true.
!forall (i=1:nx,j=1:ny)
do j=1,ny
do i=1,nx
mask(i,j)=.false.
enddo
enddo
!endforall
print*,' All arrays allocated, using RMSmin = '
print*,'test seuil',test_seuil
corr_tab=1.
do i=1,nbnodes
! print*,i,nbnodes
if( modulo(i,nbnodes/5)==0 ) print*,' correlation completed at' ,nint(float(i)/float(nbnodes)*100.),'%'
if(grille(i,3).eq.0.) then
corr_tab(i,:,:)=1.
! print*,i,'masked'
goto 4
else
if(corrMod/=1) then
it(i)=nint(grille(i,1)-(ib(i,1)-1)/2)
jt(i)=nint(grille(i,2)-(ib(i,2)-1)/2)
tot(i)=product(ib(i,:))
! print*,'it,jt',it(i),jt(i),ib(i,1),ib(i,2)
avg(i)=sum(imagea( it(i):(it(i)+ib(i,1)-1) , jt(i):(jt(i)+ib(i,2)-1) ) )/float(tot(i))
arms(i)=sqrt(sum((imagea (it(i):it(i)+ib(i,1)-1 , jt(i):jt(i)+ib(i,2)-1 ) -avg(i))**2) /float(tot(i)))
if(test_seuil) at=imagea (it(i):it(i)+ib(i,1)-1 , jt(i):jt(i)+ib(i,2)-1 )
! print*,arms(i)
if(avg(i)<=intensity_min) then
print*,i,' AVG in the box is less than ',intensity_min
grille(i,3)=3.
goto 4
endif
! if(sum(imagea (it(i):it(i)+ib(i,1)-1 , jt(i):jt(i)+ib(i,2)-1 ))/=0.) then
! print*,'***'
! print*,grille(i,1),grille(i,2)
! do j=jt(i),jt(i)+ib(i,2)-1
! write(*,200) (imagea(k,j),k=it(i),it(i)+ib(i,1)-1)
! write(*,200) (k,k=it(i),it(i)+ib(i,1)-1)
! enddo
! print*,''
200 format(20(2x,i7))
do jj=1,(is(i,2)-ib(i,2)+1)
jjt(jj)=jj-1+nint(grille(i,2)-(is(i,2)-1)/2) -nint(grille(i,5))
do ii=1,(is(i,1)-ib(i,1)+1)
iit(ii)=ii-1+nint(grille(i,1)-(is(i,1)-1)/2) +nint(grille(i,4))
bvg(i,ii,jj)=sum( imageb(iit(ii):(iit(ii)+ib(i,1)-1),jjt(jj):(jjt(jj)+ib(i,2)-1)))/float(tot(i))
brms(i,ii,jj)=sqrt(sum((imageb (iit(ii):iit(ii)+ib(i,1)-1 , jjt(jj):jjt(jj)+ib(i,2)-1 ) -bvg(i,ii,jj))**2)/float(tot(i)))
! print*,'bvg,brms,1',bvg(i,ii,jj),brms(i,ii,jj)
! if(brms(i,ii,jj)==0.) then
! print*,'Zone d intensite nulle '
! grille(i,3)=3.1
! goto 4
! endif
Corr_tab(i,ii,jj)=sum( ( imagea ( it(i):(it(i)+ib(i,1)-1) , jt(i):(jt(i)+ib(i,2)-1) ) -avg(i) )* &
( imageb ( iit(ii):(iit(ii)+ib(i,1)-1) , jjt(jj):(jjt(jj)+ib(i,2)-1) ) - bvg(i,ii,jj) ))/ &
arms(i)/brms(i,ii,jj) /float(tot(i))
enddo
enddo
if(test_seuil) then
bt=imageb (nint(grille(i,1)-(is(i,1)-1)/2) +nint(grille(i,4)):nint(grille(i,1)-(is(i,1)-1)/2) &
+nint(grille(i,4))+is(i,1)-ib(i,1),nint(grille(i,2)-(is(i,2)-1)/2) -nint(grille(i,5)):nint(grille(i,2)-(is(i,2)-1)/2)&
-nint(grille(i,5))+is(i,2)-ib(i,2))
call XcorrMASK(level1,level2,at,ib(1,1),ib(1,2),bt,is(1,1),is(1,2),ct)
corr_tab2(i,:,:)=ct
endif
! print*,maxval(corr_tab(i,:,:))
! print*,maxval(corr_tab2(i,:,:))
! print*,'******'
if( (avg(i)<=intensity_min).and.(maxval(corr_tab(i,:,:))<= (correlation_min/100.))) then
print*,i,' intensity and correlation are less than ',intensity_min,correlation_min
grille(i,3)=3.2
goto 4
endif
endif
if(corrMod/=0)then
it(i)=nint(grille(i,1)-(is(i,1)+1-1)/2)
jt(i)=nint(grille(i,2)-(is(i,2)+1-1)/2)
tot(i)=(is(i,1)+1)*(is(i,2)+1)
xm=(is(1,1)-ib(1,1))/2+1
ym=(is(1,2)-ib(1,2))/2+1
a_ron=0.
b_ron=0.
a_ron(1:is(i,1)+1,1:is(i,2)+1)=real(imagea( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2))))! - avg(i)
b_ron(1:is(i,1)+1,1:is(i,2)+1)=real(imageb( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2))))! - bvg(i,1,1)
m_ron=0.
m_ron(xm:xm+ib(1,1),ym:ym+ib(1,2))=1./float((ib(i,1)+1)*(ib(i,2)+1))
! ma(:,:)=m_ron(:,:)*a_ron(:,:)
b2(:,:)=b_ron(:,:)*b_ron(:,:)
if(i==1) then
call dfftw_plan_dft_r2c_2d (plan_ma,nxf,nyf,ma,fma,FFTW_PATIENT)
a_ron(1:is(i,1)+1,1:is(i,2)+1)=real(imagea( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2))))! - avg(i)
ma(:,:)=m_ron(:,:)*a_ron(:,:)
call dfftw_execute(plan_ma)
call dfftw_plan_dft_r2c_2d (plan_b_ron,nxf,nyf,b_ron,fb_ron,FFTW_PATIENT)
b_ron(1:is(i,1)+1,1:is(i,2)+1)=real(imageb( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2))))
call dfftw_execute(plan_b_ron)
call dfftw_plan_dft_r2c_2d (plan_m_ron,nxf,nyf,m_ron,fm_ron,FFTW_PATIENT)
m_ron=0.
m_ron(xm:xm+ib(1,1),ym:ym+ib(1,2))=1./real((ib(i,1)+1)*(ib(i,2)+1))
call dfftw_execute(plan_m_ron)
call dfftw_plan_dft_r2c_2d (plan_b2,nxf,nyf,b2,fb2,FFTW_PATIENT)
b_ron(1:is(i,1)+1,1:is(i,2)+1)=real(imageb( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2))))
b2(:,:)=b_ron(:,:)*b_ron(:,:)
call dfftw_execute(plan_b2)
else
! a_ron(1:is(i,1)+1,1:is(i,2)+1)=real(imagea( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2))))! - avg(i)
ma(:,:)=m_ron(:,:)*real(imagea( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2))))
call dfftw_execute(plan_ma)
b_ron(1:is(i,1)+1,1:is(i,2)+1)=real(imageb( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2))))
call dfftw_execute(plan_b_ron)
m_ron=0.
m_ron(xm:xm+ib(1,1),ym:ym+ib(1,2))=1./real((ib(i,1)+1)*(ib(i,2)+1))
call dfftw_execute(plan_m_ron)
b2=(real(imageb( it(i):(it(i)+is(i,1)) , jt(i):(jt(i)+is(i,2)))))**2
! b2(:,:)=b_ron(:,:)*b_ron(:,:)
call dfftw_execute(plan_b2)
endif
!****************** FFT inverse****************
if(i==1) then
call dfftw_plan_dft_c2r_2d(plan_ffmab,nxf,nyf,ffmab,iffmab, FFTW_PATIENT)!+backward)
ffmab(:,:)=(DCONJG(fma(:,:)))*(fb_ron(:,:))
call dfftw_execute(plan_ffmab)
call dfftw_plan_dft_c2r_2d(plan_ffmb,nxf,nyf,ffmb,iffmb, FFTW_PATIENT)!+backward)
ffmb(:,:)= (DCONJG(fm_ron(:,:)))*(fb_ron(:,:))
call dfftw_execute(plan_ffmb)
call dfftw_plan_dft_c2r_2d(plan_ffmb2,nxf,nyf,ffmb2, iffmb2, FFTW_PATIENT)!+backward)
ffmb2(:,:)= (DCONJG(fm_ron(:,:)))*(fb2(:,:))
call dfftw_execute(plan_ffmb2)
else
ffmab(:,:)=(DCONJG(fma(:,:)))*(fb_ron(:,:))
call dfftw_execute(plan_ffmab)
ffmb(:,:)= (DCONJG(fm_ron(:,:)))*(fb_ron(:,:))
call dfftw_execute(plan_ffmb)
ffmb2(:,:)= (DCONJG(fm_ron(:,:)))*(fb2(:,:))
call dfftw_execute(plan_ffmb2)
endif
a2(:,:)=a_ron(:,:)*a_ron(:,:)
ma2(:,:)=m_ron(:,:)*a2(:,:)
p(:,:)=iffmb(:,:)*iffmb(:,:)
m_ron=0.
sum_ma=sum(ma(:,:))
sum_ma2=sum(ma2(:,:))
sum2_ma=sum_ma**2
R(i,:,:)=(iffmab(:,:)-(sum_ma*(iffmb(:,:))))/sqrt(abs((sum_ma2-sum2_ma)*(iffmb2(:,:)-p(:,:))))
!maxfft(:) =maxloc(R(:,:))
!m_ron=0.
! R(i,:,:)=r_tmp(:,:)
!******** BE CAREFUL ; R need to circular translation ******
goto 115
m_ron(:,:)=R(i,:,:)
r_tmp(1:nxf/2,1:nyf/2)=m_ron(nxf/2+1:nxf,nyf/2+1:nyf)
r_tmp(nxf/2+1:nxf,nyf/2+1:nyf)=m_ron(1:nxf/2,1:nyf/2)
r_tmp(nxf/2+1:nxf,1:nyf/2)=m_ron(1:nxf/2,nyf/2+1:nyf)
r_tmp(1:nxf/2,nyf/2+1:nyf)=m_ron(nxf/2+1:nxf,1:nyf/2)
! R(i,1:nxf/2,1:nyf/2)=r_tmp(nxf/2+1:nxf,nyf/2+1:nyf)
! R(i,nxf/2+1:nxf,nyf/2+1:nyf)=r_tmp(1:nxf/2,1:nyf/2)
! R(i,nxf/2+1:nxf,1:nyf/2)=r_tmp(1:nxf/2,nyf/2+1:nyf)
! R(i,1:nxf/2,nyf/2+1:nyf)=r_tmp(nxf/2+1:nxf,1:nyf/2)
R(i,:,:)=r_tmp(:,:)
maxc=maxval(R)
115 continue
endif
! do j=1,(is(i,2)-ib(i,2)+1)
! write(*,201) (corr_tab(i,k,j),k=1,(is(i,1)-ib(i,1)+1))
! enddo
! print*,'is was corr'
! endif
endif
if(corrMod==2) then !Test mode
print*,"Xcorr ",maxloc(Corr_tab(i,:,:))
print*,"FFT ",maxloc(R(i,:,:))
endif
4 continue
!do jj=1,ncor(i,2)
! write(*,111) (corr_tab(i,ii,jj),ii=1,ncor(i,1))
!enddo
!print*,grille(i,4),grille(i,5)
!print*,''
enddo
!111 format (20(f6.2,2x))
deallocate(it,jt,tot,avg,arms,iit,jjt)
deallocate(bvg,brms)
deallocate(mask)
goto 1009
print*,'Valeurs globales: ',nbnodes
cm=sum(corr_tab(:,:,:)) /float(nbnodes)/float((is(i,2)-ib(i,2)+1))/float((is(i,1)-ib(i,1)+1))
cr=sqrt ( (sum( (corr_tab(:,:,:)-cm)**2 ) )/ float(nbnodes)/float((is(i,2)-ib(i,2)+1))/float((is(i,1)-ib(i,1)+1)))
print*, cm,cr,minval(corr_tab(:,:,:)),maxval(corr_tab(:,:,:))
print*,''
201 format(5(2x,f6.3))
1009 continue
!**** FFT stuff
deallocate(a_ron,b_ron,m_ron)
deallocate(ma,a2,ma2,b2)
deallocate(iffmab,iffmb,iffmb2,p)
deallocate(fma,fb_ron,fm_ron)
deallocate(fb2,ffmab,ffmb,ffmb2)
!*****
print*,'END XCORR'
return
end