[GRASS-dev] anyone interested in helping to port an algorithm
from matlab?
Dylan Beaudette
dylan.beaudette at gmail.com
Wed Nov 5 12:53:21 EST 2008
On Wednesday 05 November 2008, Dylan Beaudette wrote:
> Hi,
>
> I was taking another look on a paper published in Computers and GeoSciences
> on something called a 'Variance Quadtree' algorithm[1]. The main idea is to
> recursively partition an image into smaller and smaller rectangles, based
> on the amount of variance within each rectangle. Thus, regions of more
> variation get split into smaller nested rectangles.
>
> The source code for the algorithm is given in Matlab source[2]. I have
> tried unsuccessfully to port the code to R, mostly because I do not
> completely understand several of the matlab matrix idioms used in the code.
>
> I think that a C version of the algorithm, along with some minor
> modifications to convert the quadtree structure into GRASS vectors (similar
> to what v.surf.rst can output) would be a nice addition to imagery analysis
> in GRASS.
>
> Anyone want to give some pointers / help on porting this to C? The
> algorithm is relatively simple, however some of the matlab-specific matrix
> operations may be a challenge.
>
> 1. http://dx.doi.org/10.1016/j.cageo.2006.08.009
> 2. http://www.usyd.edu.au/su/agric/acpa/software/matlab/vqt_matlab.zip
>
>
> Cheers,
>
> Dylan
Just in case this help, attached is a working copy of a FORTRAN95 version of
the same algorithm.
Dylan
--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341
-------------- next part --------------
! compile like this:
! f95 vqt2c.for -o vqt2c
program variance_quadtree
c varaince quadtree
c (c) 2005 Australian Centre for precision Agriculture
c The University of Sydney
c version May 2006
! use dflib
include 'par.f'
real z(MaxDat,MaxDat)
real*8 qvar(MaxVar), qs2(MaxVar),qh, qmax, sum,sum2, s2
integer Indx(4),Icord(4,4),Icd(MaxVar,4),ids(MaxDat,MaxDat)
character*100 infile,datafile,outfile,quadfile
character*8 dum
logical same,iplot
! TYPE (QWINFO) winfo
namelist / vqtin / datafile,outfile,quadfile,niter,iplot,iform
! idum = ABOUTBOXQQ('VQT -Variance Quadtree')
! open (6, file='USER',title='VQT')
! winfo.H = 40
! winfo.W = 120
! winfo.TYPE = QWIN$SET
! status = SETWSIZEQQ(6, winfo)
infile='vqt_input.txt'
open(1,file=infile,status='old')
read(1, NML = vqtin)
close(1)
iplot=.false.
open (8, file= outfile, status='unknown')
open (5, file= quadfile, status='unknown')
open (4, file='all_vq.txt', status='unknown')
open (2, file='vqt_iter.txt', status='unknown')
open (1, file='vqt_o.txt', status='unknown')
call random_seed()
open(3,file=datafile,status='OLD')
read(3,*) dum,Ncol
read(3,*) dum,Nrow
read(3,*) dum,Xmin
read(3,*) dum,Ymin
read(3,*) dum,dis
read(3,*) dum,Valmis
Nr=Nrow
Nc=Ncol
if(mod(Ncol,2).ne.0) Ncol=Ncol+1
if(mod(Nrow,2).ne.0) Nrow=Nrow+1
Ncb=0
Nce=0
Nrb=0
Nre=0
Nradd=0
Ncadd=0
if(Ncol.gt.Nrow) then
Nradd=(Ncol-Nrow)
Nrb=Nradd/2
Nre=Nradd-Nrb
Ymin=Ymin-dis*Nrb
Nrow=Nrow+Nradd
elseif(Nrow.gt.Ncol) then
Ncadd=(Nrow-Ncol)
Ncb=Ncadd/2
Nce=Ncadd-Ncb
Xmin=Xmin-dis*Ncb
Ncol=Ncol+Ncadd
end if
do i=1,MaxDat
do j=1,MaxDat
z(i,j)=ValMis
end do
end do
do i=Nr+Nrb,1+Nrb, -1
read(3,*) (z(i,j), j=1+Ncb, Nc+Ncb)
end do
close(3)
Indx(1) = 1
Indx(2) = Ncol
Indx(3) = Nrow
Indx(4) = 1
! if (iplot) call DrawAxis(Indx)
nvar=0
do i=1,4
icord(i,1)=Indx(i)
end do
iter=0
if (iform.eq.1) write(5,130) 'Iter', 'x', 'y'
write(2,120) 'Iter','Max_Qh','Mean_Qh','Sd_Qh'
do iter=1, niter
write (6,50) iter
c divide into 4 rectangle
call DivideIndex(Indx,Icord)
! if (iplot) call plot4(Icord)
c calculate variance at each rect
do i=1, 4
call semivar(Z,Icord,i,qh,s2,ncol,same,ValMis)
if (same) then
c if rect reached min, find another rect
if(maxV.lt.nvar) then
do j=1, 4
Icd(maxV,j)=Icd(nvar,j)
qvar(maxV)=qvar(nvar)
qs2(maxV)=qs2(nvar)
end do
end if
nvar=nvar-1
goto 10
else
c accumulate variance list
nvar=nvar+1
qvar(nvar)=qh
qs2(nvar)=s2
do j=1, 4
Icd(nvar,j)=Icord(i,j)
end do
end if
cLeft = Xmin+dis*(Icord(i,1)-1)
cRight = Xmin+dis*(Icord(i,2)-1)
cTop = Ymin+dis*(Icord(i,3)-1)
cBottom= Ymin+dis*(Icord(i,4)-1)
write(4,80) iter,cLeft,cRight,cTop,cBottom,qh,s2
end do
c find which rectangle has the largest variance
10 qmax=-1.d0
sum=0.d0
sum2=0.d0
do iv = 1, nvar
sum=sum+qvar(iv)
sum2=sum2+qvar(iv)*qvar(iv)
if(qvar(iv).gt.qmax) then
qmax=qvar(iv)
maxV=iv
end if
end do
qn=float(nvar)
sum2=(qn*sum2-sum*sum)/(qn*(qn-1.0))
sum2=dsqrt(sum2)
sum=sum/qn
c rect with max qh is chosen
do i=1,4
Indx(i) = Icd(maxV,i)
end do
! if (iplot) call plot1(Indx)
write(2,160) iter,qMax,sum,sum2
if(iter<niter) then
c remove the chosen rect from variance list
nvar=nvar-1
do iv = maxV, nvar
qvar(iv)=qvar(iv+1)
qs2(iv)=qs2(iv+1)
do j=1, 4
Icd(iv,j)=Icd(iv+1,j)
end do
end do
end if
end do
! select a sample from all quadrants
write(8,120) 'Sample','x_Center','y_Center','x_random','y_random',
* ' Qh','s2'
ids=-9
do is=1,nvar
call alloc(Icd,is,ids)
cLeft = Xmin+dis*(Icd(is,1)-1)
cRight = Xmin+dis*(Icd(is,2)-1)
cTop = Ymin+dis*(Icd(is,3)-1)
cBottom= Ymin+dis*(Icd(is,4)-1)
c write the coordinates of quadrant
if (iform.eq.1) then
write(5,150) is,cLeft, cBottom ! txt
write(5,150) is,cRight,cBottom
write(5,150) is,cRight,cTop
write(5,150) is,cLeft, cTop
write(5,150) is,cLeft, cBottom
write(5,*)
else
write(5,*)is ! Arc format
write(5,140) cLeft,',', cBottom
write(5,140) cRight,',',cBottom
write(5,140) cRight,',',cTop
write(5,140) cLeft,',', cTop
write(5,140) cLeft,',', cBottom
write(5,*)'END'
end if
c calculate center of rect
xCenter= cLeft +(cRight-cLeft)/2
yCenter= cBottom+(cTop-cBottom)/2
c pick a random coordinate from rect
ipick=0
20 call pick(Icd(is,1),Icd(is,2),ipix)
call pick(Icd(is,4),Icd(is,3),ipiy)
if(Z(ipiy,ipix).eq.-9999) then
ipick=ipick+1
if(ipick.lt.50) goto 20
end if
xpick=Xmin+dis*(ipix-1)
ypick=Ymin+dis*(ipiy-1)
write(8,160) is,xCenter,yCenter,xpick,ypick,qvar(is),qs2(is)
end do
do i=1,Nrow
yc=Ymin+dis*(i-1)
do j=1,Ncol
xc=Xmin+dis*(j-1)
if(z(i,j).gt.valmis) write(1,*) xc,yc,z(i,j),ids(i,j)
end do
end do
if(iform.ne.1) write(5,*)'END'
50 format ('+',' Iter ', i5)
80 format ((i5,x),10(f14.2,2x))
100 format (4(i5,x),10(f15.4,2x))
120 format (a6,4x,10(3x,a10,4x))
130 format (a4,x,2(2x,a10,2x))
140 format (f12.1,a,f12.1)
150 format (i5,x,5(f15.4,2x))
160 format (i5,x,10(f15.2,2x))
close(2)
close(5)
close(8)
end program variance_quadtree
subroutine DivideIndex(Indx,Icord)
c divide rect Indx into 4 rect Icord
dimension Icord(4,4), Indx(4)
IdxLeft = Indx(1)
IdxRight = Indx(2)
IdxTop = Indx(3)
IdxBottom= Indx(4)
iwidth = IdxRight - IdxLeft
iheight= IdxTop - IdxBottom
icx = IdxLeft + iwidth/2
icy = IdxBottom+ iheight/2
Icord(1,1) = IdxLeft
Icord(1,2) = icx
Icord(1,3) = icy
Icord(1,4) = IdxBottom
Icord(2,1) = icx
Icord(2,2) = IdxRight
Icord(2,3) = icy
Icord(2,4) = IdxBottom
Icord(3,1) = IdxLeft
Icord(3,2) = icx
Icord(3,3) = Idxtop
Icord(3,4) = icy
Icord(4,1) = icx
Icord(4,2) = IdxRight
Icord(4,3) = IdxTop
Icord(4,4) = icy
return
end
subroutine semivar(Z,Icord,ir,Qh,s2,ncol,same,ValMis)
c calculate semivariance of a quadrant
include 'par.f'
dimension Z(MaxDat,MaxDat)
dimension ident(MaxGrd,2),Icord(4,4)
real*8 Qh,semv,sum,sum2,s2
logical same
integer*8 npoint,ndata
same=.false.
IdxLeft = Icord(ir,1)
IdxRight = Icord(ir,2)
IdxTop = Icord(ir,3)
IdxBottom= Icord(ir,4)
iwidth = IdxRight - IdxLeft
iheight= IdxTop - IdxBottom
if(iwidth.le.0.or.iheight.le.0) then
same=.true.
return
end if
ndata=0
c from bottom to top
do i= IdxBottom, IdxTop
c from left to right find data that belongs to this rectangle
do j= IdxLeft, IdxRight
ndata=ndata+1
ident(ndata,1)=i
ident(ndata,2)=j
end do
end do
npoint=0
nd=0
Qh=0.d0
sum=0.d0
sum2=0.d0
do i = 1, ndata
it=ident(i,1)
jt=ident(i,2)
Zi=Z(it,jt)
if(Zi.gt.ValMis) then
nd = nd+1
sum=sum+zi
sum2=sum2+zi*zi
end if
do j = i+1, ndata
it=ident(j,1)
jt=ident(j,2)
Zj=Z(it,jt)
if(Zi.gt.ValMis.and.Zj.gt.ValMis) then
npoint = npoint+1
difz = Zi - Zj
semv = difz*difz
Qh = Qh+semv
end if
end do
end do
if(npoint.le.0) then
Qh=0.d0
s2=0.d0
else
Qh=dsqrt(Qh)
s2=(float(nd)*sum2-sum*sum)/float(nd*(nd-1))
end if
return
end
subroutine pick(Idx1, Idx2, ipick)
include 'par.f'
integer idx(MaxDat)
nt=Idx2-Idx1
do i=1,nt
ij=Idx1+i
idx(i)=ij
end do
do j=nt, 1, -1
! generate uniform rand no.
call random_number(rn)
rn=float(j)*rn
k=floor(rn)+1
! exchange idj(k) <-> idx(j)
itemp=idx(k)
idx(k)=idx(j)
idx(j)=itemp
end do
ipick=idx(1)
return
end
subroutine alloc(Icord,is,ids)
c determine which strata
include 'par.f'
dimension ids(MaxDat,MaxDat),Icord(MaxVar,4)
integer*8 ndata
IdxLeft = Icord(is,1)
IdxRight = Icord(is,2)
IdxTop = Icord(is,3)
IdxBottom= Icord(is,4)
c from bottom to top
do i= IdxBottom, IdxTop
c from left to right find data that belongs to this rectangle
do j= IdxLeft, IdxRight
ids(i,j)=is
end do
end do
return
end
More information about the grass-dev
mailing list