1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
|
tol.asy.fix <- function(f,g,pooled,verbose=FALSE){
if(verbose) cat("Initialising window...")
if(!compatible(f$z,g$z,pooled$z)) stop("incompatible images 'f', 'g', 'pooled'... kernel estimates must be evaluated on identical domains")
M <- Window(pooled$z)
inside <- M$m
pres <- nrow(inside)
res2 <- 2*pres
resseq <- 1:pres
xcol.ker <- M$xstep*c(0:(pres-1),-rev(resseq))
yrow.ker <- M$ystep*c(0:(pres-1),-rev(resseq))
kerpixarea <- M$xstep*M$ystep
len.pad <- res2^2
Mpad <- matrix(0, ncol=2*pres, nrow=2*pres)
Mpad[1:pres, 1:pres] <- inside
fM <- fft(Mpad)
fg.h0 <- mean(c(f$h0,g$h0))
qb <- function(oo,hfac){
hfp <- hfac*fg.h0
densX.ker <- dnorm(xcol.ker,sd=hfp)
densY.ker <- dnorm(yrow.ker,sd=hfp)
Kern <- outer(densY.ker,densX.ker,"*")*kerpixarea
con <- fft(fM*fft(Kern), inverse=TRUE)/len.pad
qhz <- im(Mod(con[1:pres,1:pres]),xcol=oo$z$xcol,yrow=oo$z$yrow)
qhz[qhz>1] <- 1
qhz[!inside] <- NA
return(as.matrix(qhz))
}
if(verbose) cat("Done.\nPerforming kernel*window convolution(s)...")
qq <- pooled$q
if(is.null(qq)||is.vector(qq)){
qq <- qb(pooled,1)
}
qq <- as.matrix(qq)
qq2 <- (1/(4*pi))*qb(pooled,sqrt(0.5))
if(verbose) cat("Done.\n")
RrzK <- qq2/qq^2
denominator <- sqrt(RrzK*(1/length(f$h)+1/length(g$h)))/(fg.h0*sqrt(as.matrix(pooled$z)))
suppressWarnings(numerator <- as.matrix(log(f$z)-log(g$z)))
zstandard <- numerator/denominator
P <- pnorm(zstandard,lower.tail=FALSE)
return(list(num=numerator,den=denominator,zst=zstandard,p=im(P,xcol=pooled$z$xcol,yrow=pooled$z$yrow),rz=RrzK))
}
|
