Initial commit

README.md

@@ -1,2 +1,9 @@
 # idl_mommaps
 IDL package for generating moment maps from radio data cubes
+
+Requires: IDL Astronomy Library <http://idlastro.gsfc.nasa.gov>
+
+Contributions from Tony Wong, Rui Xue, Annie Hughes, and Erik Rosolowsky.
+
+* v1 - 27may2015 - initial release, still cleaning up documentation.
+* v2 - 04jun2015 - include plotting routines

cmreplicate.pro

@@ -0,0 +1,151 @@
+;+
+; NAME:
+;   CMREPLICATE
+;
+; AUTHOR:
+;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
+;   craigm@lheamail.gsfc.nasa.gov
+;
+; PURPOSE:
+;   Replicates an array or scalar into a larger array, as REPLICATE does.
+;
+; CALLING SEQUENCE:
+;   ARRAY = CMREPLICATE(VALUE, DIMS)
+;
+; DESCRIPTION: 
+;
+;   The CMREPLICATE function constructs an array, which is filled with
+;   the specified VALUE template.  CMREPLICATE is very similar to the
+;   built-in IDL function REPLICATE.  However there are two
+;   differences:
+;
+;      * the VALUE can be either scalar or an ARRAY.
+;
+;      * the dimensions are specified as a single vector rather than
+;        individual function arguments.
+;
+;   For example, if VALUE is a 2x2 array, and DIMS is [3,4], then the
+;   resulting array will be 2x2x3x4.
+;
+; INPUTS:
+;
+;   VALUE - a scalar or array template of any type, to be replicated.
+;           NOTE: These two calls do not produce the same result:
+;                  ARRAY = CMREPLICATE( 1,  DIMS)
+;                  ARRAY = CMREPLICATE([1], DIMS)
+;           In the first case the output dimensions will be DIMS and
+;           in the second case the output dimensions will be 1xDIMS
+;           (except for structures).  That is, a vector of length 1 is
+;           considered to be different from a scalar.
+;
+;   DIMS - Dimensions of output array (which are combined with the
+;          dimensions of the input VALUE template).  If DIMS is not
+;          specified then VALUE is returned unchanged.
+;
+; RETURNS:
+;   The resulting replicated array.  
+;
+; EXAMPLE:
+;   x = [0,1,2]
+;   help, cmreplicate(x, [2,2])
+;     <Expression>    INT       = Array[3, 2, 2]
+;   Explanation: The 3-vector x is replicated 2x2 times.
+;
+;   x = 5L
+;   help, cmreplicate(x, [2,2])
+;     <Expression>    LONG      = Array[2, 2]
+;   Explanation: The scalar x is replicated 2x2 times.
+;
+; SEE ALSO:
+;
+;   REPLICATE
+;
+; MODIFICATION HISTORY:
+;   Written, CM, 11 Feb 2000
+;   Fixed case when ARRAY is array of structs, CM, 23 Feb 2000 
+;   Apparently IDL 5.3 can't return from execute().  Fixed, CM, 24 Feb
+;     2000
+;   Corrected small typos in documentation, CM, 22 Jun 2000
+;   Removed EXECUTE() call by using feature of REBIN() new in IDL 5.6,
+;     (thanks to Dick Jackson) CM, 24 Apr 2009 
+;   Remove some legacy code no longer needed after above change
+;     (RETVAL variable no longer defined; thanks to A. van Engelen),
+;     CM, 08 Jul 2009
+;   Change to square bracket array index notation; there were reports
+;     of funny business with parenthesis indexing (thanks Jenny Lovell),
+;     CM, 2012-08-16
+;
+;-
+; Copyright (C) 2000, 2009, 2012, Craig Markwardt
+; This software is provided as is without any warranty whatsoever.
+; Permission to use, copy, modify, and distribute modified or
+; unmodified copies is granted, provided this copyright and disclaimer
+; are included unchanged.
+;-
+function cmreplicate, array, dims
+
+  COMPILE_OPT strictarr
+  if n_params() EQ 0 then begin
+      message, 'RARRAY = CMREPLICATE(ARRAY, DIMS)', /info
+      return, 0L
+  endif
+  on_error, 2
+
+  if n_elements(dims) EQ 0 then return, array
+  if n_elements(array) EQ 0 then $
+    message, 'ERROR: ARRAY must have at least one element'
+
+  ;; Construct new dimensions, being careful about scalars
+  sz = size(array)
+  type = sz[sz[0]+1]
+  if sz[0] EQ 0 then return, make_array(value=array, dimension=dims)
+  onedims = [sz[1:sz[0]], dims*0+1] ;; For REFORM, to extend # of dims.
+  newdims = [sz[1:sz[0]], dims    ] ;; For REBIN, to enlarge # of dims.
+  nnewdims = n_elements(newdims)
+
+  if nnewdims GT 8 then $
+    message, 'ERROR: resulting array would have too many dimensions.'
+
+  if type NE 7 AND type NE 8 AND type NE 10 AND type NE 11 then begin
+      ;; Handle numeric types
+
+      ;; Passing REBIN(...,ARRAY) is a feature introduced in
+      ;; IDL 5.6, and will be incompatible with previous versions.
+      array1 = rebin(reform([array], onedims), newdims)
+
+      return, array1
+
+  endif else begin
+      ;; Handle strings, structures, pointers, and objects separately
+
+      ;; Handle structures, which are never scalars
+      if type EQ 8 AND sz[0] EQ 1 AND n_elements(array) EQ 1 then $
+        return, reform(make_array(value=array, dimension=dims), dims, /over)
+
+      nold = n_elements(array)
+      nadd = 1L
+      for i = 0L, n_elements(dims)-1 do nadd = nadd * round(dims[i])
+      if type EQ 8 then $
+        array1 = make_array(value=array[0], dimension=[nold,nadd]) $
+      else $
+        array1 = make_array(type=type, dimension=[nold,nadd], /nozero)
+      array1 = reform(array1, [nold,nadd], /overwrite)
+      array2 = reform([array], n_elements(array))
+
+      ;; Efficient copying, done by powers of two
+      array1[0,0] = temporary(array2)
+      stride = 1L   ;; stride increase by a factor of two in each iteration
+      i = 1L & nleft = nadd - 1
+      while nleft GT stride do begin
+          array1[0,i] = array1[*,0:stride-1]  ;; Note sneaky IDL optimization
+          i = i + stride & nleft = nleft - stride
+          stride = stride * 2
+      endwhile
+      if nleft GT 0 then array1[0,i] = array1[*,0:nleft-1]
+
+      return, reform(array1, newdims, /overwrite)
+  endelse
+
+end
+
+

convol3d.pro

@@ -0,0 +1,100 @@
+;+
+; NAME:
+;	CONVOL3D
+; PURPOSE:
+;	Convolution or correlation of two 3D arrays (volumetric data),
+;	or autocorrelation of a single 3D array.
+;	Default is to compute using product of Fourier transforms.
+;
+; CALLING SEQUENCE:
+;
+;	imconv = convol3d( voldata, psf, FT_PSF = psf_FT )
+;  or:
+;	correl = convol3d( voldata1, voldata2, /CORREL )
+;  or:
+;	correl = convol3d( voldata, /AUTO )
+;
+; INPUTS:
+;	voldata = 3-D array to be convolved with PSF.
+;	psf = the 3-D array Point Spread Function,
+;		with size < or = to size of voldata.
+;
+; KEYWORDS:
+;
+;	FT_PSF = passes out/in the Fourier transform of the PSF,
+;		(so that it can be re-used the next time function is called).
+;
+;	FT_VOLDATA = passes out/in the Fourier transform of voldata.
+;
+;	/CORRELATE uses the conjugate of the Fourier transform of PSF,
+;		to compute the cross-correlation of voldata and PSF,
+;		(equivalent to IDL function convol() with NO rotation of PSF).
+;
+;	/AUTO_CORR computes the auto-correlation function of voldata using FFT.
+;
+;	/NO_FT overrides the use of FFT, using IDL function convol() instead.
+;		(then PSF is rotated by 180 degrees to give same result)
+; METHOD:
+;	When using FFT, PSF is centered & expanded to size of voldata.
+;	Basically a copy of function convolve modified from 2D to 3D.
+; HISTORY:
+;	written, Frank Varosi, NASA/GSFC 1997,
+;-
+
+function convol3d, voldata, psf, FT_PSF=psf_FT, FT_VOLDATA=vdFT, NO_FT=noft, $
+			CORRELATE=correlate, AUTO_CORRELATION=auto
+
+	sv = size( voldata )
+	sp = size( psf )
+	spf = size( psf_FT )
+	if (spf(0) ne 3) or (spf(spf(0)+1) ne 6) then no_psf_FT = 1
+
+	if (spf(0) eq 3) and (sv(0) eq 3) then begin
+		if max( abs( spf(1:sp(0)) - sv(1:sv(0)) ) ) then no_psf_FT = 1
+	   endif
+
+	if (sv(0) ne 3) or ((sp(0) ne 3) and $
+	    keyword_set( no_psf_FT ) and (NOT keyword_set( auto ))) then begin
+		print,"syntax:	result = convol3d( voldata, psf )
+		print,"    or:	result = convol3d( voldata, psf, FT_PSF=psf_FT )
+		print,"    or:	correl = convol3d( voldata1, voldata2, /CORREL )
+		print,"    or:	autocorr = convol3d( voldata, /AUTO )
+		return,0
+	   endif
+
+	if keyword_set( noft ) then begin
+		if keyword_set( auto ) then begin
+			message,"auto-correlation available only with FFT",/INFO
+			return, voldata
+		  endif else if keyword_set( correlate ) then $
+				return, convol( voldata, psf ) $
+			else	return, convol( voldata, rotate( psf, 2 ) )
+	   endif
+
+	sc = sv/2
+	npix = N_elements( voldata )
+	sif = size( vdFT )
+
+	if (sif(0) ne 3) or (sif(sif(0)+1) ne 6) or $
+	   (sif(1) ne sv(1)) or (sif(2) ne sv(2)) then vdFT = FFT( voldata,-1 )
+
+	if keyword_set( auto ) then $
+	 return, shift( npix*float( FFT( vdFT*conj( vdFT ), 1 ) ), $
+						sc(1), sc(2), sc(3) )
+
+	if Keyword_Set( no_psf_FT ) then begin
+		s2 = sc + (sv MOD 2)*(sp LT sv)   ;correct for odd size voldata.
+		Loc = (s2 - sp/2) > 0		;center PSF in new array,
+		s = (sp/2 - s2) > 0	;handle all cases: smaller or bigger
+		L = (s + sv-1) < (sp-1)
+		psf_FT = dcomplexarr( sv(1), sv(2), sv(3) )
+		psf_FT(Loc(1),Loc(2),Loc(3)) =psf(s(1):L(1),s(2):L(2),s(3):L(3))
+		psf_FT = FFT( psf_FT, -1, /OVERWRITE )
+	   endif
+
+	if keyword_set( correlate ) then $
+		conv = npix * float( FFT( vdFT * conj( psf_FT ), 1 ) ) $
+	  else	conv = npix * float( FFT( vdFT * psf_FT, 1 ) )
+
+return, shift( conv, sc(1), sc(2), sc(3) )
+end

dilate_mask.pro

@@ -0,0 +1,49 @@
+function dilate_mask, mask_in, constraint = constraint, $
+  all_neighbors = all_neighbors
+;+
+; NAME:
+;   DILATE_MASK
+;
+; PURPOSE:
+;   Starting with a core mask, the mask is expanded to contain all
+;   regions that have some overlap with the original mask.
+;
+; CALLING SEQUENCE:
+;   new_mask = DILATE_MASK(mask, constraint = constraint)
+;
+; INPUTS:
+;   MASK -- A binary mask for a data cube.
+;   CONSTRAINT -- A binary mask (presumably a superset of the original
+;                 mask)
+;
+; KEYWORD PARAMETERS:
+;   ALL_NEIGHBORS -- Keyword passed to LABEL_REGIONS().
+;
+; OUTPUTS:
+;   NEW_MASK -- A dilated mask.
+;
+; MODIFICATION HISTORY:
+;
+;       Fri Sep 2 14:07:52 2005, Erik Rosolowsky
+;       <erosolow@asgard.cfa.harvard.edu>
+;
+;		Written
+;
+;-
+
+; Generate a mask with a 2 everywhere there's an overlap and a 1
+; otherwise.  
+
+  mask = (mask_in*constraint)+constraint
+  mask_out = mask
+  mask_out[*] = 0b
+  regions = label_region(mask gt 0, all_neighbors = all_neighbors, /ulong)
+  if total(regions) eq 0 then return, !values.f_nan
+  h = histogram(regions, binsize = 1, min = 1, rev = ri)
+  for k = 0L, n_elements(h)-1 do begin
+    inds = ri[ri[k]:(ri[k+1])-1]
+    if total(mask[inds] eq 2) gt 1 then mask_out[inds] = 1b
+  endfor
+
+  return, mask_out
+end