R/make.gaussian.apertures.R

In AngusWright/LAMBDAR: Measure galaxy fluxes from an arbitrary FITS image using input aperture locations and dimensions

make.gaussian.apertures <-
function(outenv=parent.env(environment()), env=NULL){
#Procedure creates gaussian apertures, using input parameters
#from the catalogue, and places them (in order) onto stamps
#Procedure is parallelised to allow scaleability

  if (!quiet) { cat('Make_Gaussian_SA_Mask   ') }
  message('--------------------------Make_Gaussian_SA_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())
  }
  #}}}

  #Convert semi-major axis in arcsec to semi-major axis in pixels {{{
  cat.a.pix<-cat.a/(arcsec.per.pix)
  cat.b.pix<-cat.b/(arcsec.per.pix)
  #}}}

  #If needed, correct angluar coordinates {{{
  #Details {{{
  #Aperture function uses N0E90 angular inputs
  #If the input angles are not N0E90, correct for this offset }}}
  if (ang.offset) { theta.offset<-90-cat.theta } else {theta.offset<-cat.theta}
  #Correct for any reversal of fits image
  if (astr.struc$CD[1,1]*astr.struc$CD[2,2]<0){ theta.offset<-theta.offset*-1 }
  #}}}

  #Check if we want to resample.upres {{{
  if (!resample.aperture) { resample.iterations<-0 }
  #}}}

  gauss2D<-function(xi,yi,x0,y0,asig,bsig,theta) { 
    xy<-rotate.data.2d(xi,yi,theta,about.centre=TRUE)
    zp<-exp(-(((xy[,1]-x0)^2/(2*asig^2))+((xy[,2]-y0)^2/(2*bsig^2))))
    return=cbind(xi,yi,zp)
  }
  #Create Aperture Masks {{{
  message("Creating Aperture Masks")
  s_mask<-foreach(slen=stamplen,axa=theta.offset,maj=cat.a.pix,min=cat.b.pix,xdelt=(cat.x%%1),ydelt=(cat.y%%1),.inorder=TRUE, .options.mpi=mpi.opts) %dopar% {
      #Setup Grid for use in aperture generation {{{
      grid<-expand.grid(seq((1.5),(slen+0.5), by=1),seq((1.5),(slen+0.5), by=1))
      if (any(is.na(grid))){ stop(paste("NAs produced in Expand Grid. Stamplen=",slen)) }
      #}}}
      #For each stamp, place down the relevant aperture {{{
      expanded<-gauss2D(xi=grid[,1],yi=grid[,2],x0=ceiling(slen/2)+0.5+xdelt,y0=ceiling(slen/2)+0.5+ydelt,theta=axa,
                          asig=maj,bsig=min)
      matrix(expanded[,3],ncol=slen,byrow=FALSE)
      #}}}
  }
  message("Aperture Creation Complete")
  #}}}

  #Check that apertures do not cross image mask boundary {{{
  if ((cutup & length(mask.stamp)>1) | (!cutup & length(image.env$imm) > 1)) {
    #Check Mask stamps for Aperture Acceptance {{{
    message('Combining Aps with Mask Stamps')
    if (cutup) {
      sa_mask<-foreach(slen=stamplen, smask=s_mask,mmask=mask.stamp,mxl=ap.lims.mask.stamp[,1],mxh=ap.lims.mask.stamp[,2],myl=ap.lims.mask.stamp[,3],myh=ap.lims.mask.stamp[,4],.export=c("use.mask.lim","psf.filt"), .inorder=TRUE, .options.mpi=mpi.opts) %dopar% {
        #Check masking to determine if Aperture is acceptable {{{
        check<-sum(mmask[mxl:mxh,myl:myh]*smask)/sum(smask)
        if (check<use.mask.lim) {
          #Too much. Skip {{{
          array(0, dim=c(slen,slen))
          #}}}
        } else {
          #Not too much. Keep {{{
          if (psf.filt) {
            smask
          } else {
            smask*mmask[mxl:mxh,myl:myh]
          }
          #}}}
        }
        #}}}
      }
    } else {
      sa_mask<-foreach(slen=stamplen, smask=s_mask,mxl=ap.lims.mask.map[,1],mxh=ap.lims.mask.map[,2],myl=ap.lims.mask.map[,3],myh=ap.lims.mask.map[,4],.export=c("use.mask.lim","image.env","psf.filt"), .inorder=TRUE, .options.mpi=mpi.opts) %dopar% {
        #Check masking to determine if Aperture is acceptable {{{
        check<-sum(image.env$imm[mxl:mxh,myl:myh]*smask)/sum(smask)
        if (check<use.mask.lim) {
          #Too much. Skip {{{
          array(0, dim=c(slen,slen))
          #}}}
        } else {
          #Not too much. Keep {{{
          if (psf.filt) {
            smask
          } else {
            smask*image.env$imm[mxl:mxh,myl:myh]
          }
          #}}}
        }
        #}}}
      }
    }
    message('Mask Combine Finished.')
    #}}}
  } else {
    #No image mask, all can be kept {{{
    sa_mask<-s_mask
    #}}}
  }#}}}

  #Check for Foreach errors {{{
  if ((length(s_mask) > 0) && class(s_mask[[1]])=="try-error") {
    stop(paste("Fatal Error in Foreach's Multi-Core Application.\n",
    "This seems to at random, but is more likely\n",
    "to happen with increasing numbers of objects\n",
    "per thread. Try increasing the number of\n",
    "doParallel cores being used by the program.",sep=""))
  }
  #}}}

  #-----Diagnostic-----# {{{
  if (diagnostic) {
    message(paste("After assignment",round(length(which(is.na(sa_mask)))/length(sa_mask)*100,digits=2),"% of the sa_mask matrix are NA"))
  } #}}}

  #Parse Parameter Space {{{
  if (!is.null(env)) { detach(env) }
  assign("theta.offset",theta.offset,envir=outenv)
  assign("cat.theta",cat.theta,envir=outenv)
  assign("cat.a.pix",cat.a.pix,envir=outenv)
  #}}}

  message('===========END============Make_SA_MASK=============END=================\n')

  #Return array of apertures {{{
  return=sa_mask
  #}}}
}