c********************************************************************* c Filename: wv.f c Version: 2.0 c Description: transform an image into wavelet space, filter, come back c Author: Pia Mukherjee c Created at: Fri Oct 31 09:24:59 2003 c Modified at: Fri Oct 31 09:35:09 2003 c Modified by: Eddie Baron c********************************************************************* c c f77 -o wv wv.f -L/home1/andromeda/pia/PGPLOT c -R/home1/andromeda/pia/PGPLOT -lpgplot c c testing the reconstruction of P*(k) with P_j program testprog implicit none real max,min,mean,rms real lena(512,512),lenacopy(512,512) real b(0:8,0:8,0:255,0:255),x real test(512),xarr(512) integer i,j,j1,j2,l1,l2,nn(2),num,jbig character*80 infile external pwt ! read in and plot the image file infile='/home2/andromeda/pia/lena.dat' open(1,file=infile,form='formatted',status='old') do i=1,512 do j=1,512 read(1,*) lena(i,j) lenacopy(512-i,512-j)=lena(i,j) c lenacopy(i,j)=lena(i,j) end do end do close(1) open(1,file='original.bin',form='unformatted') write(1)lenacopy close(1) call plotimage(lenacopy,512,512) ! choose wavelet, example the Daubachies 4 ! can also choose d6,d8,d20 (from the Daubachies series) c2,c3 ! (from the Coiflet series) s6,s8 (from the symmlet series) or h (haar) call pwtset('d', 4) ! say that it is a 2d transform nn(1)=512 nn(2)=512 ! transform call wtn(lenacopy,nn,2,1,pwt) ! image of wavelet coeffs, mostly 0 call plotimage(lenacopy,512,512) max=-1.e10 min=1.e10 mean=0. rms=0. do i=1,512 do j=1,512 max=amax1(max,lenacopy(i,j)) min=amin1(min,lenacopy(i,j)) mean=mean+lenacopy(i,j) end do end do mean=mean/(512.**2.) do i=1,512 do j=1,512 rms=rms+(lenacopy(i,j)-mean)**2. end do end do rms=sqrt(rms/(512.**2.)) write(*,*) 'min,max,mean,rms' write(*,*) min,max,mean,rms ! just to get an idea of the min,max,mean,rms of coeffs; ! mean is small, even though the range is large write(*,*) write(*,*) & 'the scales and positions of the largest coefficients ..' ! the arrangement of the coeffs is as follows jbig=9 do j1=0,jbig-1 do j2=0,jbig-1 do l1=0,(2**j1)-1 do l2=0,(2**j2)-1 i=(2**j1)+1+l1 j=(2**j2)+1+l2 b(j1,j2,l1,l2)=lenacopy(i,j) ! b coefficients ! to get a feel for which are the largest coefficients if(lenacopy(i,j).gt.(0.5*max) .or. & lenacopy(i,j).lt.(0.5*min)) then write(*,*) j1,j2,l1,l2,b(j1,j2,l1,l2) end if end do end do end do end do ! where for 512 by 512 pixels, 512=2^jbig implies that jbig=9 ! b subscript j1,j2,l1,l2 are the wavelet coefficients for ! different scales and positions along the 2 axis (for a 2d image) write(*,*) ! some form of thresholding x=0.3 ! you could simply do this num=0 do i=1,512 do j=1,512 if(abs(lenacopy(i,j)).lt.(x*rms)) then num=num+1 lenacopy(i,j)=0. end if end do end do write(*,*) 'discarded ',num,' coeffs' write(*,*) 'or ', num/(512.**2), ' of them' ! you could have thresholded wavelet coefficients of ! different scales separately, finding the rms at each scale; ! done something more sofisticated; ! maybe sorted the coefficients and retained the highest 10%; ! but this is fine for now; can try different x ! call plotimage(lenacopy,512,512) ! inverse wavelet transform the remaining wavelet coeffs to ! reconstruct and plot image open(1,file='compressed.bin',form='unformatted') write(1)lenacopy close(1) call wtn(lenacopy,nn,2,-1,pwt) call plotimage(lenacopy,512,512) ! to see how each basis function looks for any j and l, set to unity the ! corresponding element in wavelet space and inverse transform ! the corresponding element can be found from knowing that in 1d the ! coefficients are arranged in order of increasing j (0 to 8 here) and ! within each j they are arranged in order of increasing l (0 to 2^j-1); so ! the arrangement is a00, b00, b10, b11, b20, b21, b22, b23, b30 .... do i=1,512 test(i)=0. xarr(i)=real(i) end do test(250)=1. call pwtset('d', 4) call wtn(test,512,1,-1,pwt) min=1.e10 max=-1.e10 do i=1,512 min=amin1(min,test(i)) max=amax1(max,test(i)) end do min=min-0.1 max=max+0.1 call pgbegin(0,'?',1,1) call pgenv(0.,512.,min,max,0,1) call pgline(512,xarr,test) call pgend stop end include '/home1/andromeda/pia/wavelets/wvtrans/nrwave.f' subroutine plotimage(dummy,nx,ny) IMPLICIT NONE integer nx,ny real dummy(nx,ny) INTEGER nxd,nyd,imin,imax,i real cellx,celly INTEGER nc,count,iposmin,jposmin,iposmax,jposmax PARAMETER (nc=5) REAL c(nc),r(nc),g(nc),b(nc) REAL min, max, immean, imsigma, trans(6),smin,smax REAL xmin,xmax,ymin,ymax nxd=nx nyd=ny cellx=1. celly=1. call stats1(dummy,nx,ny,nxd,nyd,min,max,immean,imsigma,count, & iposmin,jposmin,iposmax,jposmax) write(*,*) write(*,*) ' Good pixels = ', count write(*,*) ' Posn of maximum = ', iposmax,jposmax write(*,*) ' Posn of minimum = ', iposmin,jposmin write(*,*) ' Maximum of image = ', max write(*,*) ' Minimum of image = ', min write(*,*) ' Mean of image = ', immean write(*,*) ' RMS of image = ', imsigma write(*,*) call viewpt(nxd,nyd,cellx,celly,trans,xmin,xmax,ymin,ymax) c write(*,*) xmin,xmax,ymin,ymax xmin=0. xmax=cellx*nx c xmin=0. c xmax=360. ymin=0. ymax=celly*ny+0. call pgbegin(0.,'?',1,1) call pgenv(xmin,xmax,ymin,ymax,0,-1) call ctab (c,r,g,b,nc) call pgctab(c,r,g,b,nc,1.,.5) call pggray(dummy,nx,ny,1,nx,1,ny,trans) call pgimag(dummy,nx,ny,1,nx,1,ny,max,min,trans) !reversed call pgbox('ABCTSN',0.0,0,'ABCTSN',0.0,0) call pgwedg('TI',1.,4.,max,min,' ') !reversed call pgend return end c................................................................. c subroutine to define pgplot viewport subroutine viewpt(nx,ny,cellx,celly,trans,xmin,xmax,ymin,ymax) implicit none integer nx,ny,centrex,centrey real cellx,celly real trans(6) real xmin,xmax,ymin,ymax centrex=nx/2+1 centrey=ny/2+1 c trans(1)=-float(centrex)*cellx c trans(1)=30.0 trans(2)=cellx trans(3)=0. c trans(4)=133.39 c trans(4)=-float(centrey)*celly trans(5)=0. trans(6)=celly c xmin=-(float(centrex-1)+0.5)*cellx c xmax=(float(centrex-1)-0.5)*cellx c ymin=-(float(centrey-1)+0.5)*celly c ymax=(float(centrey-1)-0.5)*celly end c-----*----------------------------------------------------------------- c subroutine to find useful statistics of a real image subroutine stats1(image,nx,ny,nxd,nyd,min,max,mean,sigma,count, & imin,jmin,imax,jmax) implicit none integer*4 i,j,nx,ny,nxd,nyd,count,imin,jmin,imax,jmax real*4 image(nx,ny) real*4 min,max real*4 mean,sigma real*4 sum,sum2,bad count=0 sum=0. sum2=0. max=-1e32 min=1e32 imin=1 jmin=1 imax=1 jmax=1 bad=-999999. do i=1,nxd do j=1,nyd if (image(i,j).ne.bad) then if(image(i,j).ne.0.) then count=count+1 sum=sum+image(i,j) sum2=sum2+(image(i,j)**2) if (image(i,j).gt.max) then max=image(i,j) imax=i jmax=j end if if (image(i,j).lt.min) then min=image(i,j) imin=i jmin=j end if end if end if end do end do mean=sum/float(count) if ((sum2/float(count)-mean**2).lt.0.) then sigma=0. else sigma=sqrt(sum2/float(count)-mean**2) end if c sigma=sqrt(sum2/float(count)) c if (count.eq.1) count=nx*ny end c................................................................... subroutine ctab(x,r,g,b,nc) real x(nc),r(nc),g(nc),b(nc) c colours set as black body orange spectrum x(1)=0.00 r(1)=0.00 g(1)=0.00 b(1)=0.00 x(2)=0.25 r(2)=0.5 g(2)=0.0 b(2)=0.0 x(3)=0.5 r(3)=1.0 g(3)=0.5 b(3)=0.0 x(4)=0.75 r(4)=1.0 g(4)=1.0 b(4)=0.5 x(5)=1.0 r(5)=1.0 g(5)=1.0 b(5)=1.0 return end