library(mps,lib.loc="./mps.Rcheck") timageA <- mps.cutoff(mps.pgmload("iamgT1.pgm"),128) timageB <- mps.cutoff(mps.pgmload("iamgT2.pgm"),128) observed <- mps.ExampleOps(d=c(200,200),100) iiA <- mps.whole(timageA,observed,scale=1,rot=0,replica=3) iiB <- mps.whole(timageB,observed,scale=1,rot=0,replica=3) plot(iiA$image) plot(iiB$image) postscript(file="stripIn.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(0,0,1,0),xaxt="n",yaxt="n",bty="n") plot(iiA2r$replica[[2]]) lines(c(0,200),c(0,200),lwd=3,col="red");lines(c(0,200),c(200,0),lwd=3,col="red"); title(main="True Structure",line=0) plot(observed,xlim,ylim);title(main="Observation",line=0); plot(timageA,xlim=c(0,200),ylim=c(200,0));title(main="Training Image A",line=0); plot(iiA$image,xlim=c(0,200),ylim=c(200,0));title(main="A MPS-simulation",line=0); dev.off() G.timageA <- rbind(cbind(timageA,timageA*0+1),cbind(timageA*0+1,timageA*0+1)) class(G.timageA)<-class(timageA) S.observed <- observed[1:100,1:100] class(S.observed)<-class(observed) postscript(file="stripCond.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(0,0,1,0),xaxt="n",yaxt="n",bty="n") plot(S.observed);title(main="Observations",line=0) plot(timageA);title(main="Training Image",line=0) plot(timageA);title(main="Training Image",line=0) plot(timageA);title(main="Training Image",line=0) dev.off() postscript(file="strip1.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(0,0,1,0),xaxt="n",yaxt="n",bty="n") plot(timageA,xlim,ylim);title(main="A Training 100x100 ",line=0) plot(iiA$image);title(main="A Interp. 200x200 ",line=0) #plot(iiA$replica[[1]],main="A"); plot(timageB,xlim,ylim);title(main="B Training 100x100 ",line=0) plot(iiB$image);title(main="B Interp. 200x200 ",line=0) dev.off() iiA1 <- mps.whole(timageA,observed,scale=1,rot=0,replica=3) iiA2 <- mps.whole(timageA,observed,scale=2,rot=0,replica=3) iiA3 <- mps.whole(timageA,observed,scale=4,rot=0,replica=3) iiA4 <- mps.whole(timageA,observed,scale=c(0.25,0.5,1,2,3,4),rot=0,replica=3) postscript(file="strip2.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(0,0,1,0),xaxt="n",yaxt="n",bty="n") plot(iiA$image);title(main="A) Scale 1 ",line=0) plot(iiA2$image);title(main="A) Scale 2 ",line=0) plot(iiA3$image);title(main="A) Scale 4 ",line=0) plot(iiA4$image);title(main="A) self similar ",line=0) dev.off() postscript(file="timageA.eps",height=1.5,width=1.5,horizontal=FALSE) par(mfrow=c(1,1),mar=c(0,0,0,0),xaxt="n",yaxt="n",bty="n") plot(timageA) dev.off() iiA1r <- mps.whole(timageA,observed,scale=2,rot=0,replica=3) iiA2r <- mps.whole(timageA,observed,scale=2,rot=45,replica=3) iiA3r <- mps.whole(timageA,observed,scale=2,rot=c(0,90),replica=3) iiA4r <- mps.whole(timageA,observed,scale=2,rot=mps.isorot(13),replica=3) postscript(file="strip3.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(0,0,1,0),xaxt="n",yaxt="n",bty="n") plot(iiA1r$image);title(main="A) rot=0 ",line=0) plot(iiA2r$image);title(main="A) rot=45 ",line=0) plot(iiA3r$image);title(main="A) rot=0 and 90 ",line=0) plot(iiA4r$image);title(main="A) isotropy ",line=0) dev.off() iiA1b <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,weights=c(0.7,1)) iiA2b <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,weights=c(1,1)) iiA3b <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,weights=c(1.4,1)) iiA4b <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,weights=c(1.6,1)) postscript(file="strip4.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(0,0,1,0),xaxt="n",yaxt="n",bty="n") plot(iiA1b$image);title(main="A) 0.7:1 ",line=0) plot(iiA2b$image);title(main="A) 1:1 ",line=0) plot(iiA3b$image);title(main="A) 1.4:1 ",line=0) plot(iiA4b$image);title(main="A) 1.6:1 ",line=0) dev.off() iiA1n <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,indep=0) iiA2n <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,indep=0.01) iiA3n <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,indep=0.1) iiA4n <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,indep=1) postscript(file="strip5.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(0,0,1,0),xaxt="n",yaxt="n",bty="n") plot(iiA1n$image);title(main="A) no nugget ",line=0) plot(iiA2n$image);title(main="A) 1% nugget ",line=0) plot(iiA3n$image);title(main="A) 10% nugget ",line=0) plot(iiA4n$image);title(main="A) pure nugget ",line=0) dev.off() iiA1s <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,specificy=-3) iiA2s <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,specificy=-1) iiA3s <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,specificy=0) iiA4s <- mps.whole(timageA,observed,scale=1,rot=0,replica=3,specificy=0.6) postscript(file="strip6.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(0,0,1,0),xaxt="n",yaxt="n",bty="n") plot(iiA1s$image);title(main="A) specificy -3 ",line=0) plot(iiA2s$image);title(main="A) specificy -1 ",line=0) plot(iiA3s$image);title(main="A) normal ",line=0) plot(iiA4s$image);title(main="A) specificy +0.6 ",line=0) dev.off() plot(iiA$replica[[1]]) plot(iiA$replica[[2]]) plot(iiA$replica[[3]]) plot(iiA4$replica[[1]]) ############### Up to here it should work XXobs <- ceiling(runif(400)*200^2) newObserved <- iiA$replica[[1]] newObserved[-XXobs]<- -1 newOps <- newObserved iiA1P <- mps.whole(myTimage,newOps,scale=1,replica=15,alttimage=altTimages,debug=0) iiA2P <- mps.whole(myTimage,newOps,scale=1,replica=15,alttimage=altTimages,minTotal=10) iiA3P <- mps.whole(myTimage,newOps,scale=1,replica=15,alttimage=altTimages,debug=0,minTotal=100) iiA4P <- mps.whole(myTimage,newOps,scale=1,replica=15,alttimage=altTimages,debug=0,minTotal=1000) postscript(file="strip7.eps",height=1.1*1.5,width=4*1.5,horizontal=FALSE) par(mfrow=c(1,4),mar=c(2,2,1,1),bty="n",pch=".",pty="s") plot(port(iiA1P$report),port(iiA1P$altreport[[1]]));title(main="Nmin =1 ",line=0) plot(port(iiA2P$report),port(iiA2P$altreport[[1]]));title(main="Nmin =10 ",line=0) plot(port(iiA3P$report),port(iiA3P$altreport[[1]]));title(main="Nmin =100 ",line=0) plot(port(iiA4P$report),port(iiA4P$altreport[[1]]));title(main="Nmin =1000",line=0) dev.off() XX <- matrix(c(12,2,3,17),ncol=2) W <- diag(c(1,1/2)) XX W%*% XX %*% W exp(log(XX)*2)