!      A.Bijaoui 18-07-01
!      Haar transform
! Version MASMAI 04/11/09

      Implicit None

      Character(Len=80) Filein,FileOut,Mouch,String
      Character(Len=60) Comment
      Integer Status,BP,Nc,Nl,Np,Ip,err,i,nkeys,Ns,Nodef
      Real, allocatable :: Im1(:),Im2(:)
      Real Xd,Yd,Dx,Dy
      
      Character(Len=100) paramoption(50)

      FileIn="dummy"
      FileOut="dummy"
	  Ns=3
      mouch='mouch.lis'
      call getline(nkeys,paramoption)

    if (nkeys.eq.0) then
      Write(*,*) 'Syntax : -i Image  -o T_Haar -s Number_Scale -l logfile'
      Write(*,*) 'Haar Transform'
      Write(*,*) '-i Image: Input Image (required)'
      Write(*,*) '-o T_Haar: Haar transform (required)'
	  Write(*,*) '-s Number_Scale for the transform (def. 3)'
      Write(*,*) '-l logfile : name of the log file (def. mouch.lis)'
      Stop
     endif
       
     do i=1,Nkeys,2
       if(paramoption(i).eq.'-i') Filein    =paramoption(i+1) 
       if(paramoption(i).eq.'-s') Read(paramoption(i+1),*) Ns
       if(paramoption(i).eq.'-o') FileOut   =paramoption(i+1)
       if(paramoption(i).eq.'-l') Mouch      =paramoption(i+1)
    enddo
    
    If(FileIn .eq."dummy")    Stop "No name for input image"
    If(FileOut.eq."dummy")   Stop "No name for resulting image"
      
    Call Debut('OVN_thaar',Mouch)

    Call WC('Haar transform')

    Call WC('Input Image '//Trim(Filein))
   	Call OFits(21,FileIn,Bp,String,Nc,Nl,Dx,Dy,Xd,Yd)
    If(BP.ne.-32) Stop 'This program processes only real data'
	Np=Nc*Nl
    Allocate(Im1(Np),stat=err) ; If(err/=0) Stop "Allocation problem"
	Status=0
    Call FTGPVE(21,0,1,Np,0,Im1,Nodef,Status)
    Call FtClos(21,Status)

    Call Wi('Number of scales',Ns,1)
    Allocate(Im2(Np),stat=err) ; If(err/=0) Stop "Allocation problem"

	Call TH2D(Nc,Nl,Ns,Im1,Im2)

	Call WC('Resulting Image '//Trim(FileOut))
	Comment='Haar transform'
    Call CFits(FileOut,21,Comment,Nc,Nl,Dx,Dy,Xd,Yd,BP)
    Call FtpprE(21,0,1,Np,Im1,Status)
    Call FTClos(21,Status)

    Call Lafin('OVN_thaar')
    End

	Subroutine TH2D(Nc,Nl,Ns,W,C)
	Implicit None
	Integer*4 Nc,Nl,Is,Ncs,Nls,Ns,Ics,Ils
	Real*4 C(Nc,Nl),W(Nc,Nl),X1,X2,X3,X4,Rac05

	Ncs=Nc
	Nls=Nl
	Do Is=1,Ns
	  Do Ils=1,Nls
		Do Ics=1,Ncs
		  C(Ics,Ils)=W(Ics,Ils)
		Enddo
	  Enddo
	  Ncs=Ncs/2
	  Nls=Nls/2
	  Do Ils=1,Nls
	    Do Ics=1,Ncs
	      X1=C(2*Ics-1,2*Ils-1)
	      X2=C(2*Ics,2*Ils-1)
	      X3=C(2*Ics-1,2*Ils)
	      X4=C(2*Ics,2*Ils)
	      W(Ics,Ils)=0.5*(X1+X2+X3+X4)
	      W(Ics+Ncs,Ils)=0.5*(X1-X2+X3-X4)
	      W(Ics,Ils+Nls)=0.5*(X1+X2-X3-X4)
	      W(Ics+Ncs,Ils+Nls)=0.5*(X1-X2-X3+X4)
	    Enddo
	  Enddo
	Enddo
	Return
	End