Nothing
#' iMad (optimized port)
#'
#' Perform iteratively re-weighted multivariate alteration detection.
#'
#' @param inDataSet1 A Raster* object of the first image.
#' @param inDataSet2 A Raster* object of the second image.
#' @param pos Integer vector. A vector of bands to use from each image. Default is use all bands.
#' @param mask1 (Optional) A Raster* object representing a mask to be used for inDataSet1 or a numeric value to be used as the mask value.
#' @param mask2 (Optional) A Raster* object representing a mask to be used for inDataSet2 or a numeric value to be used as the mask value.
#' @param mask1_band (Optional) The band from inDataSet1 to use for masking (only if class(mask1)=="numeric").
#' @param mask2_band (Optional) The band from inDataSet2 to use for masking (only if class(mask1)=="numeric").
#' @param output_basename Character. The basename (including path) for the output files.
#' @param format Character. The output format of the rasters (see ?writeFormats). Default is "raster".
#' @param maxiter Numeric (>= 1). The maximum number of iterations. Default is 100.
#' @param lam Numeric. The penalization function. CURRENTLY UNSUPPORTED.
#' @param corr_thresh Numeric. Used for situations where the canonical correlates are all nearly 1.0 (how close to 1.0 does it need to be to stop).
#' @param delta Numeric. The smallest change in canonical correlates to end the program.
#' @param auto_extract_overlap Logical. Extract the overlap zones between the images?
#' @param reuse_existing_raster Logical. If the algorithm detects pre-create overlaps, use them?
#' @param force_extent Logical. Attempt to force the input files (and masks) to be the same extent?
#' @param enable_snow Logical. Use clusterR to (potentially) speed up calculations on a cluster? Default=FALSE. EXPERIMENTAL.
#' @param cl Cluster. If not assigned, the program will attempt to figure it out. Use beginCluster() to create a cluster.
#' @param verbose Logical. Print out debugging information?
#' @param ... Passed to various raster functions (see writeRaster). Important ones include format= and overwrite=TRUE/FALSE. datatype should be left as 'FLT4S' for proper functioning.
#' @return Returns a RasterStack object where the first layer is the chisquare image, and the subsequent layers are the iMad layers.
#' @author Mort Canty (original code) and Jonathan A. Greenberg (R port).
#' @seealso \code{\link[raster]{writeRaster}}
# @keywords
#' @references
#' \itemize{
#' \item {Canty, M.J. and A.A. Nielsen. 2008. Automatic radiometric normalization of multitemporal satellite imagery with the iteratively re-weighted MAD transformation. Remote Sensing of Environment 112:1025-1036.}
#' \item {Nielsen, A.A. 2007. The regularized iteratively reweighted MAD method for change detection in multi- and hyperspectral data. IEEE Transactions on Image Processing 16(2):463-478.}
#' }
# @examples
#
# \dontrun{
# }
#' @export
# TODO: parallelize if/where possible
# TODO: integrate penalization function for hyperspectral imagery.
# TODO: write out full mask and check for existing masks.
# TODO: give the output stack names
# TODO: use layerStats instead of cov.wt.raster
# TODO: allow mask to be a single value
# TODO: subset bands
# TODO: save the mask step.
iMad <- function(inDataSet1,inDataSet2,pos,
mask1,mask2,mask1_band=1,mask2_band=1,
output_basename, format="raster",
maxiter=100,lam=0,delta=0.001,corr_thresh=0.001,
auto_extract_overlap=TRUE,reuse_existing_raster=TRUE,force_extent=TRUE,
enable_snow=FALSE,cl=NULL,
verbose=FALSE,
timing=FALSE,
inmemory=FALSE,
debug=FALSE,
debug_outputs=NULL,
...)
{
require("raster")
if(verbose) { print("Verbose mode enabled...")}
if(timing) {
print("Timing enabled...")
previous_time=proc.time()
}
# Do some pre-checks up here.
if(verbose) { print("Initial checks...")}
if(!inherits(inDataSet1,'Raster')) stop ("inDataSet1 must be a Raster* object...")
if(!inherits(inDataSet2,'Raster')) stop ("inDataSet2 must be a Raster* object...")
if(maxiter < 1) stop ("maxiter must be greater than 0...")
if(missing(output_basename)) stop ("output_basename must be set...")
# End pre-checks.
# if(enable_snow)
# {
# beginCluster()
# }
# Figure out if we can run this in memory.
inmemory_layers=
# inDataSet1 (subsetted)
length(pos) +
# inDataSet2 (subsetted)
length(pos) +
# mask1
1 +
# mask2
1 +
# wt
1 +
# dm
length(pos)*2 +
# U --> probably fixable
length(pos) +
# V --> probably fixable
length(pos) +
# MAD
length(pos) +
# chisqr (I don't know if this is right) --> probably fixable
1
if(canProcessInMemory(raster(inDataSet1),n=inmemory_layers) && inmemory)
{
if(verbose) { print("Calculations will be done in memory...") }
inmemory=TRUE
} else
{
inmemory=FALSE
}
# Setup output filenames.
output_inDataSet1_subset_filename=paste(output_basename,"_inDataSet1_overlap",sep="")
output_inDataSet2_subset_filename=paste(output_basename,"_inDataSet2_overlap",sep="")
output_inDataSet1_subset_mask_filename=paste(output_basename,"_inDataSet1_mask_overlap",sep="")
output_inDataSet2_subset_mask_filename=paste(output_basename,"_inDataSet2_mask_overlap",sep="")
output_inDataSet1_masked=paste(output_basename,"_inDataSet1_masked",sep="")
output_inDataSet2_masked=paste(output_basename,"_inDataSet2_masked",sep="")
output_MAD_filename=paste(output_basename,"_iMAD",sep="")
output_chisqr_filename=paste(output_basename,"_iMAD_chisqr",sep="")
##### Subset out bands if needed.
if(!missing(pos))
{
if(verbose) { print("Subsetting bands...")}
inDataSet1=spectral_subset(inDataSet1,pos)
inDataSet2=spectral_subset(inDataSet2,pos)
}
bands=nlayers(inDataSet1)
pos=0:(bands-1)
if(timing) {
new_time=proc.time()
print("Subsetting time:")
print(new_time-previous_time)
previous_time=new_time
}
##### Extract overlap regions
if(auto_extract_overlap)
{
if(verbose) { print("Extracting the overlap region...") }
if(missing(format))
{
format="raster"
}
# Check for existing files
output_inDataSet1_subset_full_filetype <-
raster:::.filetype(format=format, filename=output_inDataSet1_subset_filename)
output_inDataSet2_subset_full_filetype <-
raster:::.filetype(format=format, filename=output_inDataSet2_subset_filename)
output_inDataSet1_subset_full_filename <-
raster:::.getExtension(output_inDataSet1_subset_filename, output_inDataSet1_subset_full_filetype)
output_inDataSet2_subset_full_filename <-
raster:::.getExtension(output_inDataSet2_subset_filename, output_inDataSet2_subset_full_filetype)
if(
!reuse_existing_raster ||
!file.exists(output_inDataSet1_subset_full_filename) ||
!file.exists(output_inDataSet2_subset_full_filename)
)
{
# If the user does not want to reuse existing cropped datasets or if they don't exist...
inDataSet1=crop(inDataSet1,inDataSet2,filename=output_inDataSet1_subset_filename,
datatype='FLT4S',verbose=verbose,...)
inDataSet2=crop(inDataSet2,inDataSet1,filename=output_inDataSet2_subset_filename,
datatype='FLT4S',verbose=verbose,...)
} else
{
# Otherwise use the existing cropped datasets.
print("Reusing the existing regions...")
inDataSet1=brick(output_inDataSet1_subset_full_filename)
inDataSet2=brick(output_inDataSet2_subset_full_filename)
}
# Now crop the masks and add them together.
if((!missing(mask1) || !missing(mask2)) && auto_extract_overlap)
{
if(!missing(mask1) && !(class(mask1)=="numeric"))
{
if(verbose) { print("Cropping mask1...") }
mask1=crop(mask1,inDataSet1)
}
if(!missing(mask2) && !(class(mask2)=="numeric"))
{
if(verbose) { print("Cropping mask2...") }
mask2=crop(mask2,inDataSet2)
}
}
} else
{
if(verbose) { print("Not extracting the overlap region...") }
cols=ncol(inDataSet1)
rows=nrow(inDataSet1)
cols2=ncol(inDataSet2)
rows2=nrow(inDataSet2)
if(cols != cols2 || rows != rows2)
{ stop("Input rows and columns must be the same, try using auto_extract_overlap=TRUE...") }
}
if(inmemory)
# Read input data into memory if possible.
{
inDataSet1=readAll(inDataSet1)
inDataSet2=readAll(inDataSet2)
}
if(timing) {
new_time=proc.time()
print("Overlap extraction time:")
print(new_time-previous_time)
previous_time=new_time
}
##### Fix extents if need be...
# Fix for extents that differ. We should combine this with a col/row check.
if(force_extent)
{
if(verbose) { print("Forcing extents...") }
extent(inDataSet2)=extent(inDataSet1)
}
##### Set up the masks
mask=NA
# Check for numeric masks
if(!missing(mask1))
{
if(class(mask1)=="numeric")
{
if(verbose) { print("Creating mask1...") }
### TODO: HPC
if(enable_snow)
{
mask1=calc_hpc(x=raster(inDataSet1,layer=mask1_band),function(x,mask1) { x != mask1 },args=list(mask1=mask1))
} else
{
mask1=raster(inDataSet1,layer=mask1_band) != mask1
}
}
}
if(!missing(mask2))
{
if(class(mask2)=="numeric")
{
if(verbose) { print("Creating mask2...") }
if(enable_snow)
{
mask2=calc_hpc(x=raster(inDataSet2,layer=mask2_band),function(x,mask2) { x != mask2 },args=list(mask2=mask2))
} else
{
mask2=raster(inDataSet2,layer=mask2_band) != mask2
}
}
}
if(!missing(mask1) && !missing(mask2))
{
if(verbose) { print("Both masks present...") }
if(force_extent)
{
extent(mask1)=extent(inDataSet1)
extent(mask2)=extent(inDataSet1)
}
# This can be HPCed
mask=mask1*mask2
} else
{
if(!missing(mask1))
{
if(verbose) { print("Mask #1 present...") }
if(force_extent)
{
extent(mask1)=extent(inDataSet1)
}
mask=mask1
} else
{
if(verbose) { print("Mask #2 present...") }
if(force_extent)
{
extent(mask2)=extent(inDataSet1)
}
mask=mask2
}
# This can be HPCed
mask[mask==0] <- NA
}
if(class(mask)!="logical")
{
if(verbose) { print("Masking...") }
if(file.exists(build_raster_filename(output_inDataSet1_masked,format=format)) && reuse_existing_raster)
{
inDataSet1=brick(build_raster_filename(output_inDataSet1_masked,format=format))
} else
{
if(enable_snow)
{
if(verbose) { print("HPC masking inDataSet1...") }
inDataSet1=mask_hpc(inDataSet1,mask,verbose=verbose)
} else
{
inDataSet1=mask(x=inDataSet1,mask=mask,filename=output_inDataSet1_masked,...)
}
}
if(file.exists(build_raster_filename(output_inDataSet2_masked,format=format)) && reuse_existing_raster)
{
inDataSet2=brick(build_raster_filename(output_inDataSet2_masked,format=format))
} else
{
if(enable_snow)
{
if(verbose) { print("HPC masking inDataSet2...") }
inDataSet2=mask_hpc(inDataSet2,mask)
} else
{
inDataSet2=mask(x=inDataSet2,mask=mask,filename=output_inDataSet2_masked,...)
}
}
}
if(timing) {
new_time=proc.time()
print("Masking time:")
print(new_time-previous_time)
previous_time=new_time
}
# if(debug) { browser() }
# if(enable_snow)
# {
# endCluster()
# beginCluster()
# }
#
if(verbose) { print("Creating initial weighting raster and stacking inDataSets...")}
# IMAGE, nb=1
if(enable_snow)
{
wt = calc_hpc(raster(inDataSet1,layer=1),fun=function(x) { x*0+1 })
wt = mask_hpc(wt,mask)
} else
{
wt = raster(inDataSet1,layer=1)*0+1
}
# if(debug) { browser() }
if(timing) {
new_time=proc.time()
print("Weight creation time:")
print(new_time-previous_time)
previous_time=new_time
}
# IMAGE, nb=nlayers(inDataSet1)+nlayers(inDataSet2)
dm = stack(inDataSet1,inDataSet2)
if(inmemory) {
dm=brick(dm)
inDataSet1=0
inDataSet2=0
mask1=0
mask2=0
}
if(timing) {
new_time=proc.time()
print("dm creation time:")
print(new_time-previous_time)
previous_time=new_time
}
# if(debug) { browser() }
delta = 1.0
oldrho = array(data=0,dim=bands)
iter = 0
ab_nan=FALSE
if(verbose) { print(memory.profile()) }
# Mods to include the penalization function. Comment this out if this chokes.
# if(lam>0) { Omega_L = diag(bands) }
if(!is.null(debug_outputs))
{
setwd(debug_outputs)
writeRaster(inDataSet1,filename="inDataSet1",format="ENVI")
writeRaster(inDataSet2,filename="inDataSet2",format="ENVI")
writeRaster(dm,filename="debug_dm",format="ENVI")
writeRaster(wt,filename="debug_wt",format="ENVI")
writeRaster(mask,filename="debug_mask",format="ENVI")
}
### MAIN LOOP
while(delta > 0.001 && iter < maxiter && !(ab_nan))
{
if(verbose)
{
print(paste("Iteration:",iter))
}
# if(debug) { browser() }
if(enable_snow)
{
if(verbose) { print("HPC calculating weighted covariance and means...")}
sigma_means=layerStats_hpc(x=dm,w=wt,stat='weighted.cov',
na.rm=TRUE,enable_snow=TRUE,verbose=verbose,todisk=TRUE)
} else
{
if(verbose) { print("Calculating weighted covariance and means...")}
sigma_means=layerStats(dm,'weighted.cov',wt,na.rm=TRUE)
}
if(timing) {
new_time=proc.time()
print("sigma_means creation time:")
print(new_time-previous_time)
previous_time=new_time
}
# if(debug) { browser() }
if(verbose) { print(sigma_means) }
sigma=sigma_means[[1]]
means=sigma_means[[2]]
if(is.nan(mean(sigma)) || is.nan(mean(means)))
{
if(verbose)
{
print("Possible numerical precision problem with sigma or means, exiting loop and using the current results...")
}
ab_nan=TRUE
break
}
s11 = sigma[1:(bands),1:(bands)]
s22 = sigma[(bands+1):(2*bands),(bands+1):(2*bands)]
s12 = sigma[1:(bands),(bands+1):(2*bands)]
s21 = sigma[(bands+1):(2*bands),1:(bands)]
if(verbose) { print("Calculating generalized eigenvalues and eigenvectors...")}
lama_a=Rdggev(JOBVL=F,JOBVR=T,A=s12%*%solve(s22)%*%s21,B=s11)
if(lama_a$INFO!=0) {
print(lama_a$INFO)
print("Error encountered in lama_a, exiting while loop...")
break
}
lamb_b=Rdggev(JOBVL=F,JOBVR=T,A=s21%*%solve(s11)%*%s12,B=s22)
if(lamb_b$INFO!=0) {
print(lamb_b$INFO)
print("Error encountered in lamb_b, exiting while loop...")
break
}
if(debug) { browser() }
a=lama_a$VR
lama=lama_a$GENEIGENVALUES
b=lamb_b$VR
lamb=lamb_b$GENEIGENVALUES
# Eigenvalues, ranked largest to smallest
idx=rank(lama)
a=a[,idx]
idx=rank(lamb)
b=b[,idx]
# Penalization stuff needs to go somewhere around here
# IDL Code:
#
#; stopping criterion
#delta = max(abs(rho-old_rho))
#old_rho = rho
#
#; ensure sum of positive correlations between X and U is positive
#; their covariance matrix is S11##A
#invsig_x = diag_matrix(1/sqrt(diag_matrix(S11)))
#sum = total(invsig_x##S11##A,2)
# A = A##diag_matrix(sum/abs(sum))
#
#; ensure positive correlation between each pair of canonical variates
#cov = diag_matrix(transpose(A)##S12##B)
#B = B##diag_matrix(cov/abs(cov))
#
#if iter gt 0 and iter eq niter then goto, done
rho=sqrt(lamb[idx])
# Fix for near-perfect correlations
if(abs(sum(rho)-bands) < corr_thresh*bands)
{
print("Perfect correlation, exiting...")
ab_nan=TRUE
break
}
# normalize dispersions
tmp1=t(a)%*%s11%*%a
tmp2=1/(sqrt(diag(tmp1)))
tmp3=t(array(tmp2,c(bands,length(tmp2))))
a=a*tmp3
tmp1=t(b)%*%s22%*%b
tmp2=1/(sqrt(diag(tmp1)))
tmp3=t(array(tmp2,c(bands,length(tmp2))))
b=b*tmp3
# assure positive correlation
tmp=diag(t(a)%*%s12%*%b)
b=b%*%diag(tmp/abs(tmp))
if(is.nan(mean(a)) || is.nan(mean(b)))
{
if(verbose)
{
print("Possible numerical precision problem with a or b, exiting loop and using the current results...")
ab_nan=TRUE
}
break
}
if(timing) {
new_time=proc.time()
print("Eigenvalues time:")
print(new_time-previous_time)
previous_time=new_time
}
if(debug) { browser() }
# canonical and MAD variates
if(enable_snow)
{
if(verbose) {
print("HPC calculating canonical and MAD variates...")
print(means)
print(bands)
print(a)
print(b)
}
MAD=calc_hpc(x=dm,args=list(a=a,b=b,means=means,bands=bands),
fun=function(x,a,b,means,bands)
{
means_a=means[1:bands]
means_b=means[(bands+1):(bands*2)]
inDataSet1=getValues(spectral_subset(x,(1:bands)))
inDataSet2=getValues(spectral_subset(x,((bands+1):(bands*2))))
U=as.vector(t(a)%*%t(inDataSet1-means_a))
V=as.vector(t(b)%*%t(inDataSet2-means_b))
MAD=U-V
return(MAD)
}
)
} else
{
if(verbose) { print("Calculating canonical and MAD variates...")}
means_a=means[1:bands]
means_b=means[(bands+1):(bands*2)]
# This needs to be fixed to use dm
U=calc(inDataSet1,fun=function(x) { as.vector(t(a)%*%(x-means_a)) } )
V=calc(inDataSet2,fun=function(x) { as.vector(t(b)%*%(x-means_b)) } )
MAD = U-V
}
if(debug) { browser() }
if(timing) {
new_time=proc.time()
print("MAD time:")
print(new_time-previous_time)
previous_time=new_time
}
# new weights
if(verbose) { print("Generating new weights...")}
var_mad=2*(1-rho)
if(verbose) { print(var_mad)}
if(debug) { browser() }
if(enable_snow)
{
if(verbose) { print("HPC chisquare...")}
chisqr=calc_hpc(x=MAD,args=list(var_mad=var_mad),
fun=function(x,var_mad)
{
x=getValues(x)
out=sum(x^2/var_mad,na.rm=TRUE)
return(out)
})
if(debug) { browser() }
if(verbose) { print("HPC new wt...")}
wt=calc_hpc(x=chisqr,
fun=function(x)
{
bands=nlayers(x)
x=getValues(x)
out=1-pchisq(x,bands)
return(out)
})
if(debug) { browser() }
wt=mask_hpc(wt,mask)
} else
{
chisqr=calc(MAD,fun=function(x) { sum(x^2/var_mad) })
wt=1-calc(chisqr,fun=function(x) { pchisq(x,bands) })
}
if(timing) {
new_time=proc.time()
print("Chisqr and new weight time:")
print(new_time-previous_time)
previous_time=new_time
}
if(!is.null(debug_outputs))
{
setwd(debug_outputs)
writeRaster(MAD,filename=paste("debug_MAD_",add_leading_zeroes(iter,number_length=3),sep=""),format="ENVI")
writeRaster(chisqr,filename=paste("debug_chisqr_",add_leading_zeroes(iter,number_length=3),sep=""),format="ENVI")
writeRaster(wt,filename=paste("chisqr_wt_",add_leading_zeroes(iter,number_length=3),sep=""),format="ENVI")
}
delta = sum(abs(rho-oldrho))
oldrho = rho
if(verbose)
{
print(paste("Delta:",delta))
print("rho:")
print(rho)
print("****************")
}
iter = iter+1
if(debug) { browser() }
}
# Output results
MAD_brick=writeRaster(MAD,filename=output_MAD_filename,format=format,...)
chisqr_raster=writeRaster(chisqr,filename=output_chisqr_filename,format=format,...)
return(stack(chisqr_raster,MAD_brick))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.