UNDER DEVELOPMENT

The software package is currently being dramatically improved and updated.
This is a (very) old version, potentially with bugs.
  Shapelets web page     Shapelets IDL code     Installation     Help pages     Simulated images  
List of routines
Instrutions for general code
Instructions for image simulation code

Source code for simage_add_noise.pro:

You can also view the help page for this routine.

pro simage_add_noise, filename, time=time, IMAGE=image, HEADER=header, $ SEED=seed, HDF=hdf, SMOOTH=smooth, NOSMOOTH=nosmooth, $ COSMIC_RAYS=cosmic_rays ;$Id: simage_add_noise.pro, v1.0$ ; ; Copyright © 2004 Richard Massey and Alexandre Refregier. ; ; This file is a part of the Shapelets analysis code. ; www.ast.cam.ac.uk/~rjm/shapelets/ ; ; The Shapelets code is free software; you can redistribute it and/or ; modify it under the terms of the GNU General Public License as published ; by the Free Software Foundation; either version 2 of the License, or ; (at your option) any later version. ; ; The Shapelets code is distributed in the hope that it will be useful, ; but WITHOUT ANY WARRANTY; without even the implied warranty of ; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ; GNU General Public License for more details. ; ; You should have received a copy of the GNU General Public License ; along with the Shapelets code; if not, write to the Free Software ; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ; ;+ ; NAME: ; SIM_ADDNOISE ; ; PURPOSE: ; Adds a model of observational noise to a simulated image. ; ; CATEGORY: ; Astronomical image simulation. ; ; INPUTS: ; filename - Name of input image, with path but w/o extension. ; time - Equivalent exposure time on WFPC2, in ks. ; Ought to be < 100, since that is the depth of the HDF! ; ; OPTIONAL INPUTS: ; [image] - Image array, to save time loading it from disc. ; [header] - Image FITS header ; [seed] - Random number seed. ; ; KEYWORD PARAMETERS: ; [/HDF] - Cuts out WFPC L shape, for fun. ; [/NOSMOOTH] - Explicitly set to -1 to force correlated background noise. ; DEFAULT is now unsmoothed, uncorrelated pixels. ; [/SMOOTH] - Set in order to enable correlated background noise. ; [/COSMIC_RAYS] - Add cosmic rays, using model by Jodi Lamoureux. ; ; OUTPUTS: ; Image with noise added is written to file image_t100.fits. ; The updated image, updated header and the inverse variance map ; are also returned in variables image, header and image_noise ; respectively. ; ; NOTES: ; The current implementation consists of a noise model including ; photon counting noise plus Gaussian background noise, with their ; levels set to roughly match those in the HDF. ; ; MODIFICATION HISTORY: ; Oct 03 - Jodi Lamoureux's cosmic ray model incorporated. ; Sep 03 - Inverse variance map output created by RM. ; May 03 - Default behaviour for correlated bground noise ; changed to unsmoothed by RM. ; Mar 03 - Program name changed to protect the innocent. ; Jul 02 - Exposure time glitch corrected by RM. ; Feb 02 - add_noise.pro written by Richard Massey. ;- ; ; Parse input requirements ; if n_elements(filename) eq 0 then filename=get_path(1)+'hst/simulation/synth2' if not keyword_set(seed) then seed=101L if not keyword_set(time) then time=100. if not keyword_set(nosmooth) then nosmooth=1 ; ; Set default noise levels ; background = 0. ; HDF-N is roughly 1.33150e-05 shot_noise = 4e-3/2. ; rms of photon shot noise (estimated from HDF) gauss_noise = 7e-5/2. ; rms of background (calculated from mean of inv. variance map) if keyword_set(smooth) then gauss_noise = 12e-5 ; compensate for smoothed noise gain = 14.0 ; [DN/e-] ? sat_level = 4096.0 ; [DN] saturation level ? read_noise = 5.0 ; [e-] ? DQE = 0.40 ; 40% maximum at 7000A ? e- to photon ratio? bias = 300.0 ; [DN] ? Fuck knows! hdftime = 100. ; [ks] (HDF-N is 100.3ks, HDF-S is 123.6ks) CCDwidth = 10.5 ; CCD width in microns CCDthickness= 200.0 ; CCD thickness in microns CRdiffusion = 0 ; Cosmic rays diffusion 0=off, 1=on CRflux = 4.5 ; Cosmic rays per second per cm^2 CRn_elec = 70.0 ; Electrons emitted per micron of silicon CRsigma = 3.5 ; Microns diffusion for 200 micron drift kernel = 2.0 ; [pixels] to simulate DRIZZLE and seeing (sky background is colvolved with this) ; The pixel scale will affect the noise smoothing kernel ; - or if images are to be drizzled, there should be no kernel at all! ; ; Read in image (NB: rd_shapecathdf has normalised things to HDF-S) ; if not keyword_set(image) or not keyword_set(header) then begin message,'Reading in noise-free image '+filename+'.fits',/info fits_read,filename+'.fits',image,header endif ; ; Determine key image properties ; im_size=size(image,/dimensions) pixel_scale_line=where(strmid(header[*],0,9) eq "PIXSIZE =") if pixel_scale_line[0] eq -1 then pixel_scale=1. $ else pixel_scale=float(strmid(header[pixel_scale_line],9,22)) ; ; Update header with parameters of the noise model ; message,'Adapting .fits header to incorporate noise parameters',/info for i=0,n_elements(header)-1 do begin ; Could use sxpar.pro in astlib case strmid(header[i],0,9) of "EXPTIME =": eline=i "PHOTOZP =": header[i] = "PHOTOZP = "+string(22.08-2.5*alog10(hdftime/time))+" /" "DATAMAX =": header[i] = "DATAMAX = "+string(float(max(image)))+" /" "DATAMIN =": header[i] = "DATAMIN = "+string(float(min(image)))+" /" else: endcase endfor if keyword_set(eline) then begin header[eline] = "EXPTIME = "+string(float(time)*1000.)+" / s Normalised to HST throughput" cortext="NOISECOR= 0 / Uncorrelated between adjacent pixels" if nosmooth eq -1 then cortext="NOISECOR= "+string(float(kernel))+" / pixels Smoothing length of background noise" header=[header[0:eline],$ ;"COMMENT = / Parameters of applied noise model",$ "NOISEPOS= 1 / Positive-definite constraint imposed",$ "NOISELEV= "+string(float(background))+" / Sky background level",$ "NOISESKY= "+string(float(gauss_noise))+" / Gaussian sky noise rms",cortext,$ "NOISEPHT= "+string(float(shot_noise))+" / Photon counting noise rms x sqrt(flux)",$ "NOISECR = 0 / No cosmic rays",$ "NOISESAT= 0 / No CCD saturation limit",$ header[eline+1:*]] endif ; Positive-definite constraint message,'Imposing positive-definite constraint',/info ;image = abs(temporary(image)) ; Impose constraint but conserve clumpiness image = temporary(image) > 0.; (truncation was used for Massey etal 2004) ; but both of these ideas might affect shearing if time ne hdftime then $ ; renormalise level for exposure time image=temporary(image)*time/hdftime ; Create inverse variance noise map image_noise=image*(shot_noise)^2+(gauss_noise)^2 ; Construct header for inverse variance noise map ;image_noise = temporary(image)*shot_noise^2+gauss_noise^2 header_noise = header for i=0,n_elements(header_noise)-1 do begin case strmid(header_noise[i],0,9) of "OBJECT =": header_noise[i] = "OBJECT = 'Inverse variance map'" "BUNIT =": header_noise[i] = "BUNIT = 'photons^(-2)' / '" "DATAMAX =": header_noise[i] = "DATAMAX = "+string(float(max(image_noise)))+" /" "DATAMIN =": header_noise[i] = "DATAMIN = "+string(float(min(image_noise)))+" /" else: endcase endfor ; Write inverse variance noise map to disc and free space in memory message,'Writing inverse variance noise map to file '+$ filename+'_t'+strmid(strtrim(time,1),0,5)+'_noise.fits',/info writefits,filename+'_t'+strmid(strtrim(time,1),0,5)+'_noise.fits',$ 1./image_noise, header_noise delvarx,image_noise ; Free up memory for later calculations delvarx,header_noise ; Add Photon counting noise to original image (stdev = sqrt(n)) message,'Adding Gaussian photon counting noise (should be Poisson)',/info ;b=sqrt(image)*randomn(seed,im_size[0],im_size[1])*shot_noise image=image+sqrt(image)*randomn(seed,im_size[0],im_size[1])*shot_noise ; ; *** Slow but memory-efficient method *** ;for y=0,im_size(2)-1 do for x=0,im_size(1)-1 do $ ; image(x,y) = image(x,y)+sqrt(image(x,y))*randomn(seed)/400. ; Getting y,x the right way round speeds it up a bit (see RSI website) ; However, 500 is a fudge factor because I don`t know the absolute ; background level in the HDF! LATER NOTE: SHOULDN`T DEPEND ON THIS! ; Add sky background THE SMOOTHING TAKES LOADS OF MEMORY! message,'Adding Gaussian sky background',/info ; NB: This really should be Gaussian. The background level is very high, but ; pre-subtracted from my image. s_noise = randomn(seed,im_size[0],im_size[1])*gauss_noise ; Scale to level observe in HDF if nosmooth eq -1 then s_noise = smooth(s_noise,kernel,/edge) ; Only top-hat, but quick (actual kernel produced by dithering is on web) s_noise = temporary(s_noise)*sqrt(time/hdftime) ; Exposure time calibration image = temporary(image)+s_noise delvarx,snoise ; Add a constant background level if background ne 0 then begin message,'Adding constant background level',/info image=temporary(image)+background endif ; Add cosmic rays (NB: should also add these to the inverse variance map) if keyword_set(cosmic_rays) then begin message,'Superimposing cosmic rays',/info cnoise = fltarr(im_size[0],im_size[1]) dummy = call_external(shapelets_paths(9,/silent)+'simages/cosmic_rays/CosmicRay.so',$ 'cosmicRayIDL', cnoise, im_size[0], im_size[1],$ CCDwidth, CCDthickness, time, $ CRdiffusion, CRflux, CRn_elec, CRsigma) image = temporary(image)+cnoise/time delvarx,cnoise endif ; Chop out a HST-type mask for fun! if keyword_set(hdf) then begin image[0:80,*]=0. ; left image[4015:4095,*]=0. ; right image[*,0:80]=0. ; bottom image[*,4015:4095]=0. ; top image[0:1000,2049:4095]=0. ; L1 image[0:2048,3050:4095]=0. ; L2 endif ; Write image to file message,'Writing noisy image to file '+$ filename+'_t'+strmid(strtrim(time,1),0,5)+'.fits',/info writefits, filename+'_t'+strmid(strtrim(time,1),0,5)+'.fits', image, header end

View the help page for this routine, return to the shapelets web page or return to the code help menu.


Last modified on 02nd May 2008 by Richard Massey.

Valid HTML 4.01!