make.data.array.maps <-
function(outenv=parent.env(environment()), env=NULL){
#Procedure creates data mask arrays, using input images and
#location data from the catalogue
#Procedure is parallelised to allow scaleability
#
# LIMITS NOMENCLATURE:
# > ap.lims.data.map is the limits of the aperture stamp in the image.env$im space
# > ap.lims.data.stamp is the limits of the aperture stamp in the data.stamp space
# > data.stamp.lims is the limits of the data.stamp in the image.env$im space
#
# > ap.lims.mask.map is the limits of the aperture stamp in the image.env$imm space
# > ap.lims.mask.stamp is the limits of the aperture stamp in the mask.stamp space
# > mask.stamp.lims is the limits of the mask.stamp in the image.env$imm space
#
# > ap.lims.error.map is the limits of the aperture stamp in the image.env$ime space
# > ap.lims.error.stamp is the limits of the aperture stamp in the error.stamp space
# > error.stamp.lims is the limits of the error.stamp in the image.env$ime space
#
#
if (!quiet) { cat('Make_Data_Masks ') }
message('--------------------------Make_Data_Mask---------------------------------')
# Load Parameter Space {{{
if(!is.null(env)) {
attach(env, warn.conflicts=FALSE)
}
if(is.null(outenv)&!is.null(env)) { outenv<-env }
else if (is.null(outenv)) {
warning("Output Environment cannot be NULL; using parent env")
outenv<-parent.env(environment())
}
#}}}
#Setup sizes {{{
npos<-length(cat.ra)
if (length(image.env$imm)!=1) {
immxpixmax<-length(image.env$imm[,1])
immypixmax<-length(image.env$imm[1,])
} else {
immxpixmax<-1
immypixmax<-1
}
if (length(image.env$ime)!=1) {
imexpixmax<-length(image.env$ime[,1])
imeypixmax<-length(image.env$ime[1,])
} else {
imexpixmax<-1
imeypixmax<-1
}
#}}}
#Check to make sure pixel values are finite {{{
if (length(which(!(is.finite(x.pix) | !(is.finite(y.pix)))))!=0) {
sink(type="message") ; stop(paste("x.pix[",which(!(is.finite(x.pix))),"] or y.pix[",which(!(is.finite(x.pix))),"] are NA/NaN/Inf in Make_Data_Mask()", sep=""))
}#}}}
#Set Stamp pixel limits {{{
message("Setting Stamp Limits")
#Details {{{
#Stampwidths: aperture width multiplied by a buffer, plus the
#PSF WIDHT determined by the desired confidence. Total axis length MUST be odd
#If stamplen==1, then make the stamplen == smallest nonzero aperture by default
#This happens when we have a point source and are not convolving with PSF
#If ALL apertures are point sources, then this will have length 5}}}
if (psf.filt) {
if (aperture.type==2) {
if (confidence==1) {
confidence=1-1E-16
}
nsig<-sqrt(qchisq(confidence,df=2))
stamplen<-(floor(ceiling(def.buff*nsig*cat.a/arcsec.per.pix)+(ceiling(psf.clip)/2))*2+5)
} else {
stamplen<-(floor(ceiling(def.buff*cat.a/arcsec.per.pix)+(ceiling(psf.clip)/2))*2+5)
}
} else {
if (aperture.type==2) {
if (confidence==1) {
confidence=1-1E-16
}
nsig<-sqrt(qchisq(confidence,df=2))
stamplen<-(floor(ceiling(def.buff*nsig*cat.a/arcsec.per.pix))*2+5)
} else {
stamplen<-(floor(ceiling(def.buff*cat.a/arcsec.per.pix))*2+5)
}
}
#Deal with Point Sources {{{
if (all(stamplen==5)) { stamplen[which(stamplen==5)]<-15 }
else if (any(stamplen==5)) { stamplen[which(stamplen==5)]<-floor(min(cat.a[which(cat.a>0)],na.rm=TRUE))*2+15 }
#}}}
#Calculate Stamp limits in image-pixel space {{{
ap.lims.data.stamp<-cbind(x.pix-floor(stamplen/2), x.pix+floor(stamplen/2), y.pix-floor(stamplen/2), y.pix+floor(stamplen/2))
#}}}
#}}}
#Check that stamp limits are within image {{{
message("Removing Stamps that lie outside the image limits")
cat.len<-length(cat.x)
inside.mask<-!((ap.lims.data.stamp[,1]<1)|(ap.lims.data.stamp[,2]>length(image.env$im[,1]))|(ap.lims.data.stamp[,3]<1)|(ap.lims.data.stamp[,4]>length(image.env$im[1,])))
#}}}
#If not being careful {{{
if (!exists("careful")) { careful<-TRUE }
if (any(!inside.mask)& !careful) {
message(paste("There are originally ",length(which(!inside.mask)),"sources to remove because we would have been careful"))
#Trim down the offending aperture stamps
ap.lims.data.stamp[,1]<-rowMaxs(cbind(ap.lims.data.stamp[,1],rep(1,cat.len)))
ap.lims.data.stamp[,3]<-rowMaxs(cbind(ap.lims.data.stamp[,3],rep(1,cat.len)))
ap.lims.data.stamp[,2]<-rowMins(cbind(ap.lims.data.stamp[,2],rep(length(image.env$im[,1]),cat.len)))
ap.lims.data.stamp[,4]<-rowMins(cbind(ap.lims.data.stamp[,4],rep(length(image.env$im[1,]),cat.len)))
#Update the stamp lengths
stamplen[which(!inside.mask)]<-(ap.lims.data.stamp[which(!inside.mask),2]-ap.lims.data.stamp[which(!inside.mask),1])+1
check.stamplen<-(ap.lims.data.stamp[,4]-ap.lims.data.stamp[,3])+1
#Check that the stamps are square
if (any(stamplen != check.stamplen)) {
for (i in which(stamplen != check.stamplen)) {
if (stamplen[i]<check.stamplen[i]) {
#push out the x-axis length
if (ap.lims.data.stamp[i,1]==1) {
ap.lims.data.stamp[i,2]<-ap.lims.data.stamp[i,1]+check.stamplen[i]-1
} else if (ap.lims.data.stamp[i,2]==length(image.env$im[,1])) {
ap.lims.data.stamp[i,1]<-ap.lims.data.stamp[i,2]-check.stamplen[i]+1
}
} else {
#push out the y-axis length
if (ap.lims.data.stamp[i,3]==1) {
ap.lims.data.stamp[i,4]<-ap.lims.data.stamp[i,3]+stamplen[i]-1
} else if (ap.lims.data.stamp[i,4]==length(image.env$im[1,])) {
ap.lims.data.stamp[i,3]<-ap.lims.data.stamp[i,4]-stamplen[i]+1
}
}
}
}
#Update the stamp lengths
stamplen[which(!inside.mask)]<-(ap.lims.data.stamp[which(!inside.mask),2]-ap.lims.data.stamp[which(!inside.mask),1])+1
check.stamplen<-(ap.lims.data.stamp[,4]-ap.lims.data.stamp[,3])+1
#Check that the stamps are odd
if (any(stamplen%%2 == 0)) {
for (i in which(stamplen%%2 == 0)) {
#push out the x-axis length
if (ap.lims.data.stamp[i,1]==1) {
ap.lims.data.stamp[i,2]<-ap.lims.data.stamp[i,2]+1
} else if (ap.lims.data.stamp[i,2]==length(image.env$im[,1])) {
ap.lims.data.stamp[i,1]<-ap.lims.data.stamp[i,1]-1
}
#push out the y-axis length
if (ap.lims.data.stamp[i,3]==1) {
ap.lims.data.stamp[i,4]<-ap.lims.data.stamp[i,4]+1
} else if (ap.lims.data.stamp[i,4]==length(image.env$im[1,])) {
ap.lims.data.stamp[i,3]<-ap.lims.data.stamp[i,3]-1
}
}
}
stamplen[which(!inside.mask)]<-floor((ap.lims.data.stamp[which(!inside.mask),2]-ap.lims.data.stamp[which(!inside.mask),1])/2)*2+1
check.stamplen<-floor((ap.lims.data.stamp[,4]-ap.lims.data.stamp[,3])/2)*2+1
#recalculate the inside.mask values
inside.mask<-!((ap.lims.data.stamp[,1]<1)|(ap.lims.data.stamp[,2]>length(image.env$im[,1]))|(ap.lims.data.stamp[,3]<1)|(ap.lims.data.stamp[,4]>length(image.env$im[1,])))
message(paste("There are now ",length(which(!inside.mask)),"sources to remove after not being careful"))
}
#}}}
#Remove any apertures whos stamps would cross the boundary {{{
if (length(which(inside.mask==TRUE))==0) {
#sink(type="message") ; stop("No Single Aperture Stamps are entirely inside the image.") }
#Return NULL
assign("cutup",TRUE,envir=outenv)
return(NULL)
}
x.pix<-x.pix[which(inside.mask)]
y.pix<-y.pix[which(inside.mask)]
cat.x<-cat.x[which(inside.mask)]
cat.y<-cat.y[which(inside.mask)]
cat.id<-cat.id[which(inside.mask)]
cat.ra<-cat.ra[which(inside.mask)]
cat.dec<-cat.dec[which(inside.mask)]
cat.theta<-cat.theta[which(inside.mask)]
cat.a<-cat.a[which(inside.mask)]
cat.b<-cat.b[which(inside.mask)]
stamplen<-stamplen[which(inside.mask)]
if (use.segmentation) { seg.x<-seg.x[which(inside.mask)] }
if (use.segmentation) { seg.y<-seg.y[which(inside.mask)] }
ap.lims.data.stamp<-rbind(ap.lims.data.stamp[which(inside.mask),])
if (length(flux.weight)!=1) { flux.weight<-flux.weight[which(inside.mask)] }
if (filt.contam) { contams<-contams[which(inside.mask)] }
if (exists("groups")) { groups<-groups[which(inside.mask)] }
if (exists("psf.id")) { psf.id<-psf.id[which(inside.mask)] }
if (exists("skyest.sources")) { skyest.sources<-skyest.sources[which(inside.mask)] }
inside.mask<-inside.mask[which(inside.mask)]
npos<-length(inside.mask)
#}}}
#Setup MPI Options {{{
chunk.size=ceiling(npos/getDoParWorkers())
mpi.opts<-list(chunkSize=chunk.size)
#}}}
#Calculate image Stamp limits in image-pixel space {{{
factor<-ifelse(floor(cat.a/arcsec.per.pix)*5>floor(psffwhm*10), floor(cat.a/arcsec.per.pix)*5, floor(psffwhm*10))
factor<-ifelse(floor(stamplen/2) >factor , floor(stamplen/2) , factor )
data.stamp.lims<-cbind(x.pix-factor, x.pix+factor, y.pix-factor, y.pix+factor)
data.stamp.lims[which(data.stamp.lims<1)]<-1
data.stamp.lims[which(data.stamp.lims[,4]>length(image.env$im[1,])),4]<-length(image.env$im[1,])
data.stamp.lims[which(data.stamp.lims[,2]>length(image.env$im[,1])),2]<-length(image.env$im[,1])
#}}}
#Construct Data Stamps {{{
data.stamp<-list(NULL)
if (mem.safe) {
#Determine if cutting up the images is worthwhile {{{
nchild<-getDoParWorkers()
#Memory of Cutups {{{
cutmem<-sum(stamplen^2)*3+sum((factor*2+1)^2)
#}}}
if ((cutmem)>(length(image.env$im[,1])*length(image.env$im[1,])*nchild)) {
cutup<-FALSE
message("Memory required by cutting up image is more than parsing whole image to children\n> USE INDIVIDUAL STAMPS: FALSE")
} else {
cutup<-TRUE
message("Memory required by parsing whole image to children is more than cutting into individual stamps\n> USE INDIVIDUAL STAMPS: TRUE")
}
#}}}
} else {
#Memory not a concern {{{
cutup<-TRUE
message("Memory considerations not required. Cutting datamap into individual stamps\n> USE INDIVIDUAL STAMPS: TRUE")
#}}}
}
if (!cutup) {
#total memory from stamps is larger than total memory of moving whole image around; Not worth cutting up {{{
data.stamp<-NULL
#}}}
#Calculate Stamp limits in image-stamp space {{{
data.stamp.lims<-cbind(rep(1,npos),rep(length(image.env$im[,1]),npos),rep(1,npos),rep(length(image.env$im[1,]),npos))
ap.lims.data.map<-ap.lims.data.stamp
ap.lims.data.stamp<-ap.lims.data.stamp
#}}}
} else {
#Total memory of stamps is less that total memory of moving whole image around; Cut up image {{{
#Create image stamps {{{
for (i in 1:npos) {
data.stamp[[i]]<-image.env$im[data.stamp.lims[i,1]:data.stamp.lims[i,2],data.stamp.lims[i,3]:data.stamp.lims[i,4]]
}
#}}}
#}}}
#Calculate Stamp limits in image-stamp space {{{
ap.lims.data.map<-ap.lims.data.stamp
ap.lims.data.stamp<-ap.lims.data.stamp-(data.stamp.lims[,c(1,1,3,3)]-1)
#}}}
}
#}}}
#Notify {{{
if (verbose) { message(paste("There are",length(cat.x),"supplied objects have stamps entirely inside the image (",
round(((cat.len-length(cat.x))/cat.len)*100, digits=2),"% of supplied lie across the edge of the image )")) }
#}}}
#Create Mask Stamps {{{
#Details {{{
#For each catalogue entry, make a stamp of the mask at that position
#If the mask is unity everywhere (or is of length 1)
#we don't need to use the imm mask at all; so skip}}}
if (((length(image.env$imm)!=1)&(length(which(image.env$imm!=1))!=0))) {
#Get Mask Limits {{{
#Get Pixel Locations in Mask {{{
if (mask.map=="NONE") {
#If mask map is NONE, it has been created by the weight.map, the error.map, or using the
#Aperture exclusion zone
if (weight.map=="NONE" & error.map=="NONE") {
#Aperture Exclusion Zone
astr.struc.mask<-read.astrometry(file.path(path.root,path.work,data.map))
} else if (weight.map=="NONE") {
#Error map
astr.struc.mask<-read.astrometry(file.path(path.root,path.work,error.map))
} else {
#Weight map
astr.struc.mask<-read.astrometry(file.path(path.root,path.work,weight.map))
}
} else {
astr.struc.mask<-read.astrometry(file.path(path.root,path.work,mask.map))
}
if (!(all(astr.struc$CRVAL==astr.struc.mask$CRVAL,na.rm=TRUE)&
all(astr.struc$CRPIX==astr.struc.mask$CRPIX,na.rm=TRUE)&
all(astr.struc$CD ==astr.struc.mask$CD ,na.rm=TRUE))){
#Get object locations in pixel space {{{
gama.pos<-ad.to.xy(cat.ra,cat.dec,astr.struc.mask)
mask.x.pix<-floor(gama.pos[,1])
mask.y.pix<-floor(gama.pos[,2])
#}}}
} else {
#Use Image pixel locations {{{
mask.x.pix<-x.pix
mask.y.pix<-y.pix
#}}}
}
#}}}
#}}}
#Calculate Mask limits in mask-pixel space {{{
mask.stamp.lims<-cbind(mask.x.pix-factor, mask.x.pix+factor, mask.y.pix-factor, mask.y.pix+factor)
mask.stamp.lims[which(mask.stamp.lims<1)]<-1
mask.stamp.lims[which(mask.stamp.lims[,4]>length(image.env$imm[1,])),4]<-length(image.env$imm[1,])
mask.stamp.lims[which(mask.stamp.lims[,2]>length(image.env$imm[,1])),2]<-length(image.env$imm[,1])
ap.lims.mask.map<-cbind(mask.x.pix-floor(stamplen/2), mask.x.pix+floor(stamplen/2), mask.y.pix-floor(stamplen/2), mask.y.pix+floor(stamplen/2))
#}}}
#Check that stamp limits are within image {{{
cat.len<-length(cat.x)
inside.mask<-!((ap.lims.mask.map[,1]<1)|(ap.lims.mask.map[,2]>length(image.env$imm[,1]))|(ap.lims.mask.map[,3]<1)|(ap.lims.mask.map[,4]>length(image.env$imm[1,])))
#}}}
#Remove any apertures whos stamps would cross the boundary {{{
if (any(!inside.mask)) {
message("Removing Stamps that lie outside the mask image limits")
if (length(which(inside.mask==TRUE))==0) {
#sink(type="message") ; stop("No Single Aperture Stamps are entirely inside the image.")
#Return NULL
assign("cutup",cutup,envir=outenv)
return(NULL)
}
factor<-factor[which(inside.mask)]
x.pix<-x.pix[which(inside.mask)]
y.pix<-y.pix[which(inside.mask)]
mask.x.pix<-mask.x.pix[which(inside.mask)]
mask.y.pix<-mask.y.pix[which(inside.mask)]
cat.x<-cat.x[which(inside.mask)]
cat.y<-cat.y[which(inside.mask)]
cat.id<-cat.id[which(inside.mask)]
cat.ra<-cat.ra[which(inside.mask)]
cat.dec<-cat.dec[which(inside.mask)]
cat.theta<-cat.theta[which(inside.mask)]
cat.a<-cat.a[which(inside.mask)]
cat.b<-cat.b[which(inside.mask)]
stamplen<-stamplen[which(inside.mask)]
if (use.segmentation) { seg.x<-seg.x[which(inside.mask)] }
if (use.segmentation) { seg.y<-seg.y[which(inside.mask)] }
if (length(data.stamp )>1) { data.stamp<-data.stamp[which(inside.mask)] }
ap.lims.data.stamp<-rbind(ap.lims.data.stamp[which(inside.mask),])
ap.lims.data.map<-rbind(ap.lims.data.map[which(inside.mask),])
data.stamp.lims<-rbind(data.stamp.lims[which(inside.mask),])
ap.lims.mask.map<-rbind(ap.lims.mask.map[which(inside.mask),])
mask.stamp.lims<-rbind(mask.stamp.lims[which(inside.mask),])
if (length(flux.weight)!=1) { flux.weight<-flux.weight[which(inside.mask)] }
if (filt.contam) { contams<-contams[which(inside.mask)] }
if (exists("groups")) { groups<-groups[which(inside.mask)] }
if (exists("psf.id")) { psf.id<-psf.id[which(inside.mask)] }
if (exists("skyest.sources")) { skyest.sources<-skyest.sources[which(inside.mask)] }
inside.mask<-inside.mask[which(inside.mask)]
npos<-length(inside.mask)
}
#}}}
#Setup MPI Options {{{
chunk.size=ceiling(npos/getDoParWorkers())
mpi.opts<-list(chunkSize=chunk.size)
#}}}
if (careful) {
#Mask the area where apertures cannot be placed {{{
minlen<-ceiling(min(stamplen)/2)
image.env$imm[1:minlen,]<-0
image.env$imm[,1:minlen]<-0
image.env$imm[dim(image.env$imm)[1]-(1:minlen-1),]<-0
image.env$imm[,dim(image.env$imm)[2]-(1:minlen-1)]<-0
#}}}
}
#Create Mask Stamps {{{
message('Creating Mask Stamps')
mask.stamp<-list(NULL)
if (!cutup) {
mask.stamp<-1
ap.lims.mask.stamp<-ap.lims.mask.map
mask.stamp.lims<-cbind(rep(1,npos),rep(length(image.env$imm[,1]),npos),rep(1,npos),rep(length(image.env$imm[1,]),npos))
} else {
for (i in 1:npos) {
mask.stamp[[i]]<-image.env$imm[mask.stamp.lims[i,1]:mask.stamp.lims[i,2],mask.stamp.lims[i,3]:mask.stamp.lims[i,4]]
}
ap.lims.mask.stamp<-ap.lims.mask.map-(mask.stamp.lims[,c(1,1,3,3)]-1)
}
} else {
#Set the mask parameters to image parameters
mask.x.pix<-x.pix
mask.y.pix<-y.pix
ap.lims.mask.stamp<-ap.lims.data.stamp
ap.lims.mask.map<-ap.lims.data.map
mask.stamp.lims<-data.stamp.lims
#Update mask for the invalid area around the image (where aps cannot be placed)
#Initialise mask as transparant everywhere
image.env$imm<-array(1,dim=dim(image.env$im))
if (careful) {
#Mask the area where apertures cannot be placed {{{
minlen<-ceiling(min(stamplen)/2)
image.env$imm[1:minlen,]<-0
image.env$imm[,1:minlen]<-0
image.env$imm[dim(image.env$imm)[1]-(1:minlen-1),]<-0
image.env$imm[,dim(image.env$imm)[2]-(1:minlen-1)]<-0
#}}}
}
#If needed, now cutup the mask
if (cutup) {
mask.stamp<-list(NULL)
for (i in 1:npos) {
mask.stamp[[i]]<-image.env$imm[mask.stamp.lims[i,1]:mask.stamp.lims[i,2],mask.stamp.lims[i,3]:mask.stamp.lims[i,4]]
}
} else {
mask.stamp<-1
}
}#}}}
#}}}
#Create Error Stamps {{{
#Details {{{
#For each catalogue entry, make a stamp of the error map at that position
#If the error map is a single value everywhere (or is of length 1)
#we don't need to use the ime mask at all; so skip}}}
if ((length(image.env$ime)!=1)&(any(image.env$ime!=image.env$ime[1]))&
(file.exists(file.path(path.root,path.work,error.map))|file.exists(file.path(path.root,path.work,weight.map)))) {
#Check Errormap Astrometry {{{
if (file.exists(file.path(path.root,path.work,error.map))) {
astr.struc.err<-read.astrometry(file.path(path.root,path.work,error.map))
} else if (file.exists(file.path(path.root,path.work,weight.map))) {
astr.struc.err<-read.astrometry(file.path(path.root,path.work,weight.map))
}
if (!(all(astr.struc$CRVAL==astr.struc.err$CRVAL,na.rm=TRUE)&
all(astr.struc$CRPIX==astr.struc.err$CRPIX,na.rm=TRUE)&
all(astr.struc$CD ==astr.struc.err$CD ,na.rm=TRUE))){
#Get object locations in pixel space {{{
gama.pos<-ad.to.xy(cat.ra,cat.dec,astr.struc.err)
error.x.pix<-floor(gama.pos[,1])
error.y.pix<-floor(gama.pos[,2])
#}}}
} else {
#Use Image pixel locations {{{
error.x.pix<-x.pix
error.y.pix<-y.pix
#}}}
}
#}}}
#Calculate Mask limits in mask-pixel space {{{
error.stamp.lims<-cbind(error.x.pix-factor, error.x.pix+factor, error.y.pix-factor, error.y.pix+factor)
error.stamp.lims[which(error.stamp.lims<1)]<-1
error.stamp.lims[which(error.stamp.lims[,4]>length(image.env$ime[1,])),4]<-length(image.env$ime[1,])
error.stamp.lims[which(error.stamp.lims[,2]>length(image.env$ime[,1])),2]<-length(image.env$ime[,1])
ap.lims.error.map<-cbind(error.x.pix-floor(stamplen/2), error.x.pix+floor(stamplen/2), error.y.pix-floor(stamplen/2), error.y.pix+floor(stamplen/2))
#}}}
#Check that stamp limits are within image {{{
cat.len<-length(cat.x)
inside.mask<-!((ap.lims.error.map[,1]<1)|(ap.lims.error.map[,2]>length(image.env$ime[,1]))|(ap.lims.error.map[,3]<1)|(ap.lims.error.map[,4]>length(image.env$ime[1,])))
#}}}
#Remove any apertures whos stamps would cross the boundary {{{
if (any(!inside.mask)) {
message("Removing Stamps that lie outside the error image limits")
if (length(which(inside.mask==TRUE))==0) {
#sink(type="message") ; stop("No Single Aperture Stamps are entirely inside the image.")
#Return NULL
assign("cutup",cutup,envir=outenv)
return(NULL)
}
factor<-factor[which(inside.mask)]
x.pix<-x.pix[which(inside.mask)]
y.pix<-y.pix[which(inside.mask)]
mask.x.pix<-mask.x.pix[which(inside.mask)]
mask.y.pix<-mask.y.pix[which(inside.mask)]
error.x.pix<-error.x.pix[which(inside.mask)]
error.y.pix<-error.y.pix[which(inside.mask)]
cat.x<-cat.x[which(inside.mask)]
cat.y<-cat.y[which(inside.mask)]
cat.id<-cat.id[which(inside.mask)]
cat.ra<-cat.ra[which(inside.mask)]
cat.dec<-cat.dec[which(inside.mask)]
cat.theta<-cat.theta[which(inside.mask)]
cat.a<-cat.a[which(inside.mask)]
cat.b<-cat.b[which(inside.mask)]
if (use.segmentation) { seg.x<-seg.x[which(inside.mask)] }
if (use.segmentation) { seg.y<-seg.y[which(inside.mask)] }
if (length(data.stamp )>1) { data.stamp<-data.stamp[which(inside.mask)] }
if (length(mask.stamp)>1) { mask.stamp<-mask.stamp[which(inside.mask)] }
stamplen<-stamplen[which(inside.mask)]
ap.lims.data.stamp<-rbind(ap.lims.data.stamp[which(inside.mask),])
ap.lims.data.map<-rbind(ap.lims.data.map[which(inside.mask),])
data.stamp.lims<-rbind(data.stamp.lims[which(inside.mask),])
ap.lims.mask.stamp<-rbind(ap.lims.mask.stamp[which(inside.mask),])
ap.lims.mask.map<-rbind(ap.lims.mask.map[which(inside.mask),])
mask.stamp.lims<-rbind(mask.stamp.lims[which(inside.mask),])
ap.lims.error.map<-rbind(ap.lims.error.map[which(inside.mask),])
error.stamp.lims<-rbind(error.stamp.lims[which(inside.mask),])
if (length(flux.weight)!=1) { flux.weight<-flux.weight[which(inside.mask)] }
if (filt.contam) { contams<-contams[which(inside.mask)] }
inside.mask<-inside.mask[which(inside.mask)]
npos<-length(inside.mask)
}
#}}}
#Setup MPI Options {{{
chunk.size=ceiling(npos/getDoParWorkers())
mpi.opts<-list(chunkSize=chunk.size)
#}}}
#Create Error Stamps {{{
message('Creating Error Stamps')
error.stamp<-list(NULL)
if (!cutup) {
error.stamp<-NULL
ap.lims.error.stamp<-ap.lims.error.map
error.stamp.lims<-cbind(rep(1,npos),rep(length(image.env$ime[,1]),npos),rep(1,npos),rep(length(image.env$ime[1,]),npos))
} else {
for (i in 1:npos) {
error.stamp[[i]]<-image.env$ime[error.stamp.lims[i,1]:error.stamp.lims[i,2],error.stamp.lims[i,3]:error.stamp.lims[i,4]]
}
ap.lims.error.stamp<-ap.lims.error.map-(error.stamp.lims[,c(1,1,3,3)]-1)
}
#}}}
} else if (error.map=="NONE"|!is.na(suppressWarnings(as.numeric(error.map)))) {
#Use Image pixel locations {{{
error.x.pix<-x.pix
error.y.pix<-y.pix
#}}}
#Calculate Mask limits in mask-pixel space {{{
error.stamp.lims<-data.stamp.lims
ap.lims.error.map<-ap.lims.data.map
ap.lims.error.stamp<-ap.lims.data.stamp
#}}}
#Create Error Stamps {{{
message('Creating Error Stamps')
error.stamp<-list(NULL)
if (!cutup) {
error.stamp<-NULL
} else {
for (i in 1:npos) {
error.stamp[[i]]<-image.env$ime[error.stamp.lims[i,1]:error.stamp.lims[i,2],error.stamp.lims[i,3]:error.stamp.lims[i,4]]
}
}
#}}}
} else if (length(image.env$ime)==1) {
error.x.pix<-x.pix
error.y.pix<-y.pix
error.stamp<-image.env$ime
error.stamp.lims<-data.stamp.lims
ap.lims.error.map<-ap.lims.data.map
ap.lims.error.stamp<-ap.lims.data.stamp
} else if (!any(image.env$ime!=image.env$ime[1])) {
error.x.pix<-x.pix
error.y.pix<-y.pix
image.env$ime<-image.env$ime[1]
error.stamp<-image.env$ime[1]
error.stamp.lims<-data.stamp.lims
ap.lims.error.map<-ap.lims.data.map
ap.lims.error.stamp<-ap.lims.data.stamp
}#}}}
#}}}
#If we've made multiple masks, Check that all arrays are conformable {{{
check<-FALSE
if (cutup) {
if ((length(mask.stamp)!=1)|(length(error.stamp)!=1)) { check<-TRUE }
} else {
l1<-length(image.env$im)!=1
l2<-length(image.env$imm)!=1
l3<-length(image.env$ime)!=1
d1<-dim(image.env$im)
d2<-dim(image.env$imm)
d3<-dim(image.env$ime)
if (l1&l2&(any(d1!=d2)|any(x.pix!=mask.x.pix)|any(y.pix!=mask.y.pix))) { check<-TRUE }
if (l1&l3&(any(d1!=d3)|any(x.pix!=error.x.pix)|any(y.pix!=error.y.pix))) { check<-TRUE }
if (l2&l3&(any(d2!=d3)|any(mask.x.pix!=error.x.pix)|any(mask.y.pix!=error.y.pix))) { check<-TRUE }
}
if (check) {
for (i in 1:npos) {
if (!cutup) {
im.mk<-matrix(0,nrow=length(data.stamp.lims[i,1]:data.stamp.lims[i,2]),ncol=length(data.stamp.lims[i,3]:data.stamp.lims[i,4]))
} else {
im.mk<-data.stamp[[i]]
}
#Get the masks' dimensions {{{
dims<-rbind(dim(im.mk))
index<-1
if (length(image.env$imm)!=1) {
if (!cutup) {
imm.mk<-matrix(0,nrow=length(mask.stamp.lims[i,1]:mask.stamp.lims[i,2]),ncol=length(mask.stamp.lims[i,3]:mask.stamp.lims[i,4]))
} else {
imm.mk<-mask.stamp[[i]]
}
#If the mask.stamp exists, get its dimension {{{
dims<-rbind(dims,dim(imm.mk))
index<-c(index,2)
#}}}
}
if (length(image.env$ime)!=1) {
if (!cutup) {
ime.mk<-matrix(0,nrow=length(error.stamp.lims[i,1]:error.stamp.lims[i,2]),ncol=length(error.stamp.lims[i,3]:error.stamp.lims[i,4]))
} else {
ime.mk<-error.stamp[[i]]
}
#If the error.stamp exists, get its dimension{{{
dims<-rbind(dims,dim(ime.mk))
index<-c(index,3)
#}}}
}#}}}
#If any x-dimensions are different {{{
if (any(dims[,1]!=dims[1,1])) {
#Determine how many dimensions are different {{{
len=length(which(dims[,1]!=min(dims[,1])))
#}}}
#For the number of different dimensions:
for (j in 1:len) {
#Get the matrix with the smallest dimension {{{
lo<-index[which.min(dims[,1])]
#}}}
#Get the matrix with the largest dimension {{{
hi<-index[which.max(dims[,1])]
#}}}
#Determine if the matrix has been cutoff at the low or high end {{{
if (lo==1) {
#Image matrix
end<-ifelse(data.stamp.lims[i,1]==1,FALSE,TRUE)
} else if (lo==2) {
#Mask matrix
end<-ifelse(mask.stamp.lims[i,1]==1,FALSE,TRUE)
} else {
#error matrix
end<-ifelse(error.stamp.lims[i,1]==1,FALSE,TRUE)
}#}}}
#Find how much they differ {{{
differ<-diff(c(dims[which.min(dims[,1]),1],dims[which.max(dims[,1]),1]))
#}}}
if (end) {
#Matrix is cutoff at the high end; modify the high matrix {{{
if (hi==1) {
#Image matrix
im.mk<-im.mk[1:(dims[which(index==hi),1]-differ),]
if (cutup) { data.stamp[[i]]<-data.stamp[[i]][1:(dims[which(index==hi),1]-differ),] }
data.stamp.lims[i,2]<-data.stamp.lims[i,2]-differ
dims[which(index==hi),]<-dim(im.mk)
} else if (hi==2) {
#Mask matrix
imm.mk<-imm.mk[1:(dims[which(index==hi),1]-differ),]
if (cutup) { mask.stamp[[i]]<-mask.stamp[[i]][1:(dims[which(index==hi),1]-differ),] }
mask.stamp.lims[i,2]<-mask.stamp.lims[i,2]-differ
dims[which(index==hi),]<-dim(imm.mk)
} else {
#error matrix
ime.mk<-ime.mk[1:(dims[which(index==hi),1]-differ),]
if (cutup) { error.stamp[[i]]<-error.stamp[[i]][1:(dims[which(index==hi),1]-differ),] }
error.stamp.lims[i,2]<-error.stamp.lims[i,2]-differ
dims[which(index==hi),]<-dim(ime.mk)
}#}}}
} else {
#Matrix is cutoff at the low end {{{
if (hi==1) {
#Image matrix
im.mk<-im.mk[(differ+1):dims[which(index==hi),1],]
if (cutup) { data.stamp[[i]]<-data.stamp[[i]][(differ+1):dims[which(index==hi),1],] }
data.stamp.lims[i,1]<-data.stamp.lims[i,1]+differ
dims[which(index==hi),]<-dim(im.mk)
} else if (hi==2) {
#Mask matrix
imm.mk<-imm.mk[(differ+1):dims[which(index==hi),1],]
if (cutup) { mask.stamp[[i]]<-mask.stamp[[i]][(differ+1):dims[which(index==hi),1],] }
mask.stamp.lims[i,1]<-mask.stamp.lims[i,1]+differ
dims[which(index==hi),]<-dim(imm.mk)
} else {
#error matrix
ime.mk<-ime.mk[(differ+1):dims[which(index==hi),1],]
if (cutup) { error.stamp[[i]]<-error.stamp[[i]][(differ+1):dims[which(index==hi),1],] }
error.stamp.lims[i,1]<-error.stamp.lims[i,1]+differ
dims[which(index==hi),]<-dim(ime.mk)
}
}#}}}
}
}#}}}
#If any y-dimensions are different {{{
if (any(dims[,2]!=dims[1,2])) {
#Determine how many dimensions are different {{{
len=length(which(dims[,2]!=min(dims[,2])))
#}}}
#For the number of different dimensions:
for (j in 1:len) {
#Get the matrix with the smallest dimension #{{{
lo<-index[which.min(dims[,2])]
#}}}
#Get the matrix with the largest dimension #{{{
hi<-index[which.max(dims[,2])]
#}}}
#Determine if the matrix has been cutoff at the low or high end {{{
if (lo==1) {
#Image matrix
end<-ifelse(data.stamp.lims[i,3]==1,FALSE,TRUE)
} else if (lo==2) {
#Mask matrix
end<-ifelse(mask.stamp.lims[i,3]==1,FALSE,TRUE)
} else {
#error matrix
end<-ifelse(error.stamp.lims[i,3]==1,FALSE,TRUE)
}#}}}
#Find how much they differ {{{
differ<-diff(c(dims[which.min(dims[,2]),2],dims[which.max(dims[,2]),2]))
#}}}
if (end) {
#Matrix is cutoff at the high end; modify the high matrix {{{
if (hi==1) {
#Image matrix
im.mk<-im.mk[,1:(dims[which(index==hi),2]-differ)]
if (cutup) { data.stamp[[i]]<-data.stamp[[i]][,1:(dims[which(index==hi),2]-differ)] }
data.stamp.lims[i,4]<-data.stamp.lims[i,4]-differ
dims[which(index==hi),]<-dim(im.mk)
} else if (hi==2) {
#Mask matrix
imm.mk<-imm.mk[,1:(dims[which(index==hi),2]-differ)]
if (cutup) { mask.stamp[[i]]<-mask.stamp[[i]][,1:(dims[which(index==hi),2]-differ)] }
mask.stamp.lims[i,4]<-mask.stamp.lims[i,4]-differ
dims[which(index==hi),]<-dim(imm.mk)
} else {
#error matrix
ime.mk<-ime.mk[,1:(dims[which(index==hi),2]-differ)]
if (cutup) { error.stamp[[i]]<-error.stamp[[i]][,1:(dims[which(index==hi),2]-differ)] }
error.stamp.lims[i,4]<-error.stamp.lims[i,4]-differ
dims[which(index==hi),]<-dim(ime.mk)
}#}}}
} else {
#Matrix is cutoff at the low end #{{{
if (hi==1) {
#Image matrix
im.mk<-im.mk[,(differ+1):dims[which(index==hi),2]]
if (cutup) { data.stamp[[i]]<-data.stamp[[i]][,(differ+1):dims[which(index==hi),2]] }
data.stamp.lims[i,3]<-data.stamp.lims[i,3]+differ
dims[which(index==hi),]<-dim(im.mk)
} else if (hi==2) {
#Mask matrix
imm.mk<-imm.mk[,(differ+1):dims[which(index==hi),2]]
if (cutup) { mask.stamp[[i]]<-mask.stamp[[i]][,(differ+1):dims[which(index==hi),2]] }
mask.stamp.lims[i,3]<-mask.stamp.lims[i,3]+differ
dims[which(index==hi),]<-dim(imm.mk)
} else {
#error matrix
ime.mk<-ime.mk[,(differ+1):dims[which(index==hi),2]]
if (cutup) { error.stamp[[i]]<-error.stamp[[i]][,(differ+1):dims[which(index==hi),2]] }
error.stamp.lims[i,3]<-error.stamp.lims[i,3]+differ
dims[which(index==hi),]<-dim(ime.mk)
}
}#}}}
}
}#}}}
if (!cutup){
#If not cutting up, will be the same for all apertures; Change all and break
data.stamp.lims<-cbind(rep(data.stamp.lims[i,1],npos),rep(data.stamp.lims[i,2],npos),rep(data.stamp.lims[i,3],npos),rep(data.stamp.lims[i,4],npos))
mask.stamp.lims<-cbind(rep(mask.stamp.lims[i,1],npos),rep(mask.stamp.lims[i,2],npos),rep(mask.stamp.lims[i,3],npos),rep(mask.stamp.lims[i,4],npos))
error.stamp.lims<-cbind(rep(error.stamp.lims[i,1],npos),rep(error.stamp.lims[i,2],npos),rep(error.stamp.lims[i,3],npos),rep(error.stamp.lims[i,4],npos))
break
}
}
#Update Stamp limits in image-stamp space {{{
ap.lims.data.stamp<-ap.lims.data.map-(data.stamp.lims[,c(1,1,3,3)]-1)
if (length(image.env$imm)!=1) {
ap.lims.mask.stamp<-ap.lims.mask.map-(mask.stamp.lims[,c(1,1,3,3)]-1)
} else {
mask.stamp.lims<-data.stamp.lims
ap.lims.mask.stamp<-ap.lims.data.stamp
}
if (length(image.env$ime)!=1) {
ap.lims.error.stamp<-ap.lims.error.map-(error.stamp.lims[,c(1,1,3,3)]-1)
} else {
error.stamp.lims<-data.stamp.lims
ap.lims.error.stamp<-ap.lims.data.stamp
}
#}}}
#Check that no stamps have become too small {{{
inside.mask<-!((ap.lims.data.stamp[,1]<1)|(ap.lims.data.stamp[,3]<1)|(ap.lims.data.stamp[,2]>as.numeric(abs(diff(t(data.stamp.lims[,c(1,2)])))+1))|(ap.lims.data.stamp[,4]>as.numeric(abs(diff(t(data.stamp.lims[,c(3,4)])))+1)))
if (length(which(!inside.mask))!=0) {
warning(paste("Mismatches in the edges of the input images (image, mask.map, error.map) means",length(which(!inside.mask)),"stamps have to be discarded"))
#Remove stamps {{{
data.stamp<-data.stamp[which(inside.mask)]
if (length(mask.stamp)>1) { mask.stamp<-mask.stamp[which(inside.mask)] }
if (length(error.stamp)>1) { error.stamp<-error.stamp[which(inside.mask)] }
x.pix<-x.pix[which(inside.mask)]
y.pix<-y.pix[which(inside.mask)]
mask.x.pix<-mask.x.pix[which(inside.mask)]
mask.y.pix<-mask.y.pix[which(inside.mask)]
error.x.pix<-error.x.pix[which(inside.mask)]
error.y.pix<-error.y.pix[which(inside.mask)]
cat.x<-cat.x[which(inside.mask)]
cat.y<-cat.y[which(inside.mask)]
cat.id<-cat.id[which(inside.mask)]
cat.ra<-cat.ra[which(inside.mask)]
cat.dec<-cat.dec[which(inside.mask)]
cat.theta<-cat.theta[which(inside.mask)]
cat.a<-cat.a[which(inside.mask)]
cat.b<-cat.b[which(inside.mask)]
stamplen<-stamplen[which(inside.mask)]
if (use.segmentation) { seg.x<-seg.x[which(inside.mask)] }
if (use.segmentation) { seg.y<-seg.y[which(inside.mask)] }
if (length(data.stamp )>1) { data.stamp<-data.stamp[which(inside.mask)] }
if (length(mask.stamp)>1) { mask.stamp<-mask.stamp[which(inside.mask)] }
if (length(error.stamp)>1) { error.stamp<-error.stamp[which(inside.mask)] }
ap.lims.data.stamp<-rbind(ap.lims.data.stamp[which(inside.mask),])
ap.lims.mask.stamp<-rbind(ap.lims.mask.stamp[which(inside.mask),])
ap.lims.error.stamp<-rbind(ap.lims.error.stamp[which(inside.mask),])
ap.lims.data.map<-rbind(ap.lims.data.map[which(inside.mask),])
ap.lims.mask.map<-rbind(ap.lims.mask.map[which(inside.mask),])
ap.lims.error.map<-rbind(ap.lims.error.map[which(inside.mask),])
data.stamp.lims<-rbind(data.stamp.lims[which(inside.mask),])
mask.stamp.lims<-rbind(mask.stamp.lims[which(inside.mask),])
error.stamp.lims<-rbind(error.stamp.lims[which(inside.mask),])
if (length(flux.weight)!=1) { flux.weight<-flux.weight[which(inside.mask)] }
if (filt.contam) { contams<-contams[which(inside.mask)] }
inside.mask<-inside.mask[which(inside.mask)]
npos<-length(inside.mask)
#}}}
}
#}}}
#Update the Masked the area for where apertures cannot be placed {{{
minlen<-ceiling(min(stamplen)/2)
min.x<-ap.lims.mask.map[which.min(ap.lims.mask.map[,1]),1]+minlen
max.x<-ap.lims.mask.map[which.max(ap.lims.mask.map[,2]),2]-minlen
min.y<-ap.lims.mask.map[which.min(ap.lims.mask.map[,3]),3]+minlen
max.y<-ap.lims.mask.map[which.max(ap.lims.mask.map[,4]),4]-minlen
image.env$imm[1:min.x,]<-0
image.env$imm[,1:min.y]<-0
image.env$imm[max.x:dim(image.env$imm)[1],]<-0
image.env$imm[,max.y:dim(image.env$imm)[2]]<-0
mask.xrange<-c(min(mask.stamp.lims[,1]),max(mask.stamp.lims[,2]))
mask.yrange<-c(min(mask.stamp.lims[,3]),max(mask.stamp.lims[,4]))
data.xrange<-c(min(data.stamp.lims[,1]),max(data.stamp.lims[,2]))
data.yrange<-c(min(data.stamp.lims[,3]),max(data.stamp.lims[,4]))
#Initilise array as opaque beyond the image limits
image.env$imm.dimim<-array(0,dim=dim(image.env$im))
image.env$imm.dimim[data.xrange[1]:data.xrange[2],data.yrange[1]:data.yrange[2]]<-image.env$imm[mask.xrange[1]:mask.xrange[2],mask.yrange[1]:mask.yrange[2]]
#}}}
} else {
image.env$imm.dimim<-image.env$imm
}
#}}}
#Parse Parameter Space {{{
if (!is.null(env)) { detach(env) }
assign("cutup",cutup,envir=outenv)
assign("x.pix",x.pix,envir=outenv)
assign("y.pix",y.pix,envir=outenv)
assign("cat.x",cat.x,envir=outenv)
assign("cat.y",cat.y,envir=outenv)
assign("cat.id",cat.id,envir=outenv)
assign("cat.ra",cat.ra,envir=outenv)
assign("cat.dec",cat.dec,envir=outenv)
assign("cat.theta",cat.theta,envir=outenv)
assign("cat.a",cat.a,envir=outenv)
assign("cat.b",cat.b,envir=outenv)
assign("stamplen",stamplen,envir=outenv)
assign("ap.lims.data.map",ap.lims.data.map,envir=outenv)
assign("ap.lims.mask.map",ap.lims.mask.map,envir=outenv)
assign("ap.lims.error.map",ap.lims.error.map,envir=outenv)
assign("ap.lims.data.stamp",ap.lims.data.stamp,envir=outenv)
assign("ap.lims.mask.stamp",ap.lims.mask.stamp,envir=outenv)
assign("ap.lims.error.stamp",ap.lims.error.stamp,envir=outenv)
assign("data.stamp.lims",data.stamp.lims,envir=outenv)
assign("mask.stamp.lims",mask.stamp.lims,envir=outenv)
assign("error.stamp.lims",error.stamp.lims,envir=outenv)
assign("mask.stamp",mask.stamp,envir=outenv)
assign("error.stamp",error.stamp,envir=outenv)
assign("flux.weight",flux.weight,envir=outenv)
if (filt.contam) { assign("contams",contams,envir=outenv) }
if (use.segmentation) { assign("seg.x",seg.x,envir=outenv) }
if (use.segmentation) { assign("seg.y",seg.y,envir=outenv) }
assign("inside.mask",inside.mask,envir=outenv)
assign("mpi.opts",mpi.opts,envir=outenv)
#}}}
message('===========END============Make_SA_MASK=============END=================\n')
#Return array of apertures {{{
return=data.stamp
#}}}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.