Viewing contents of file '../idllib/sdss/allpro/dataimage__define.pro'
;----------------------------------------------------------------------------
; dataImage::INIT
;
; Purpose:
; Sets the parameters to be appropriate for a data image
;
FUNCTION dataImage::INIT, filename, tempname
dummyData=FLTARR(10,10)
SXHMAKE,dummyData,1,newHeader
self.dataHeader=PTR_NEW(newHeader)
self.numberOfExposures=1
IF N_PARAMS() GT 0 THEN BEGIN
self.fileName=filename
IF N_PARAMS() GT 1 THEN self.tempName=tempName
self->load
ENDIF
RETURN, 1
END
;----------------------------------------------------------------------------
; dataImage::SetProperty
;+
; NAME:
; dataImage::SetProperty
;
; PURPOSE:
; This method sets the properties of the data image.
;
; CALLING SEQUENCE:
; oDataImage->SetProperty, FILENAME=filename, TEMPNAME=tempName,$
; DATAARRAY=data, EXPOSURETIME=exposureTime,$
; NEXPOSURES=numberOfExposures, DATAHEADER=dataHeader
;
; KEYWORD PARAMETERS:
; FILENAME - a string containing the name of the file assiciate
; with the data image
; TEMPNAME - a string with the file name of a temporary file for the
; data to be saved in if onDisk is used
; DATAARRAY- a two dimensional floating point array containing
; the image data
; EXPOSURETIME - the exposure time of the image, such that the
; count rate for a given x,y position equals
; data(x,y)/exposureTime
; NEXPOSURES - the number of exposures combined in the image
; DATAHEADER - a structure containing the header information for the
; image, in the format used by the image loading and saving
; routines from the astrolib library
;
;
; EXAMPLE:
; oDataImage->setProperty,FILENAME='newImageFile.hhh',$
; TEMPNAME='/tmp/newImage.dat',$
; DATAARRAY='newArray'
; EXPOSURETIME=400
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::SetProperty, FILENAME=fileName,$
TEMPNAME=tempName, $
DATAARRAY=data,$
EXPOSURETIME=exposureTime,$
NEXPOSURES=numberOfExposures,$
DATAHEADER=dataHeader
IF KEYWORD_SET(fileName) THEN self.fileName=fileName
IF KEYWORD_SET(tempName) THEN self.tempName=tempName
IF KEYWORD_SET(data) THEN self->setData, data
IF KEYWORD_SET(exposureTime) THEN BEGIN
self.exposureTime=exposureTime
ENDIF
IF KEYWORD_SET(numberOfExposures) THEN $
self.numberOfExposures=numberOfExposures
IF KEYWORD_SET(dataHeader) THEN BEGIN
PTR_FREE, self.dataHeader
self.dataHeader=PTR_NEW(dataHeader)
ENDIF
IF(PTR_VALID(self.dataHeader)) THEN BEGIN
IF(PTR_VALID(self.data)) THEN BEGIN
sizeInfo=SIZE(*self.data)
naxis1=sizeInfo(1)
naxis2=sizeInfo(2)
SXADDPAR,*self.dataHeader,'NAXIS1',naxis1
SXADDPAR,*self.dataHeader,'NAXIS2',naxis2
ENDIF
IF KEYWORD_SET(exposureTime) THEN BEGIN
exptime=self.exposureTime
SXADDPAR,*self.dataHeader,'EXPTIME',exptime
ENDIF
ENDIF
END
;----------------------------------------------------------------------------
; dataImage::GetProperty
;+
; NAME:
; dataImage::GetProperty
;
; PURPOSE:
; This method gets some of the properties of the data image.
;
; CALLING SEQUENCE:
; oDataImage->GetProperty, FILENAME=filename, TEMPNAME=tempName,$
; DATAARRAY=data, EXPOSURETIME=exposureTime,$
; NEXPOSURES=numberOfExposures, DATAHEADER=dataHeader
;
; KEYWORD PARAMETERS:
; FILENAME - a string containing the name of the file assiciate
; with the data image
; TEMPNAME - a string with the file name of a temporary file for the
; data to be saved in if onDisk is used
; DATAARRAY- a two dimensional floating point array containing
; the image data
; EXPOSURETIME - the exposure time of the image, such that the
; count rate for a given x,y position equals
; data(x,y)/exposureTime
; NEXPOSURES - the number of exposures combined in the image
; DATAHEADER - a structure containing the header information for the
; image, in the format used by the image loading and saving
; routines from the astrolib library
;
; EXAMPLE:
; oDataImage->GetProperty,FILENAME=filename,$
; EXPOSURETIME=exposureTime
; PRINT,"Image ',filename,' has an exposure time of ',exposureTime
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::GetProperty, FILENAME=fileName,$
TEMPNAME=tempName, $
DATAARRAY =data,$
EXPOSURETIME=exposureTime,$
NEXPOSURES=numberOfExposures,$
DATAHEADER=dataHeader
fileName=self.fileName
tempName=self.tempName
IF PTR_VALID(self.data) THEN data=self->getData()
exposureTime=self.exposureTime
numberOfExposures=self.numberOfExposures
IF PTR_VALID(self.dataHeader) THEN $
dataHeader=*self.dataHeader
END
;----------------------------------------------------------------------------
; dataImage::getData
;+
; NAME:
; dataImage::getData
;
; PURPOSE:
; This methods returns an array containing the image data,
; optionally scaled to a provided exposure time, or for display.
;
; CALLING SEQUENCE:
; Result=oDataImage->GetData()
;
; KEYWORD PARAMETERS:
; TIME: If this keyword is set, the images is scaled to the
; provided exposure time.
; DISPLAY: If this keyword is set, the image is scaled to a
; 0:255 range suitable for display.
; RANGE: This keword is used to scale the display output to a
; specified range of pixel values. The return array
; will be scaled such that pixels of value range[0]
; will be 0, and those of range[1] will be 255.
; ERANGE: If this keyword is set, a value for the RANGE is estimated.
;
; OUTPUTS:
; This function returns the data array for the image, scaled as
; specified by any present keywords.
;
; EXAMPLE:
; To display an image such that pixel values between 50 and 4000
; are distinguishable,
; TV, oDataImage->GetData(/ERANGE)
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
;
FUNCTION dataImage::getData, TIME=time, DISPLAY=display, RANGE=range,$
ERANGE=erange
IF PTR_VALID(self.data) THEN BEGIN
returnData=self->astronomyImage::getData()
IF KEYWORD_SET(time) THEN BEGIN
returnData= returnData*(time/self.exposureTime)
ENDIF
IF KEYWORD_SET(display) OR KEYWORD_SET(range) THEN BEGIN
IF NOT KEYWORD_SET(range) THEN BEGIN
returnData=BYTSCL(returnData)
ENDIF ELSE BEGIN
returnDAta=BYTSCL(returnData,MIN=range(0),MAX=range(1))
ENDELSE
ENDIF
IF KEYWORD_SET(ERANGE) THEN BEGIN
range=self->estimateDisplayRange()
returnData=BYTSCL(returnData,MIN=range(0),MAX=range(1))
ENDIF
ENDIF
RETURN, returnData
END
;----------------------------------------------------------------------------
; dataImage::clone
;+
; NAME:
; dataImage::clone
;
; PURPOSE:
; This method creates a duplicate image object. If no temporary
; file name is provided, the data for the new image is stored in
; memory.
;
; CALLING SEQUENCE:
; Result = oDataImage->Clone(Tempfilename)
;
; OPTIONAL INPUTS:
; tempFileName: This string contains the file name for the file
; in which to store the image data, if any.
;
; OUTPUTS:
; A copy of the image is returned.
;
; SIDE EFFECTS:
; The clone is a new object; each must be destroyed separately.
;
; EXAMPLE:
; imagecopy = exampleimage->clone(tempfilename)
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen, Jr. 1997
;-
FUNCTION dataImage::clone, tempfilename
IF N_PARAMS() EQ 0 THEN BEGIN
imageClone=self->astronomyImage::clone()
ENDIF ELSE BEGIN
imageClone=self->astronomyImage::clone(tempfilename)
ENDELSE
imageClone->setProperty, EXPOSURETIME=self.exposureTime,$
NEXPOSURES=self.numberOfExposures
RETURN, imageClone
END
;----------------------------------------------------------------------------
; dataImage::load
;+
; NAME:
; dataImage::load
;
; PURPOSE:
; This method load an image from disk into the image object.
;
; CALLING SEQUENCE:
; oDataImage->Load, FileName
;
; OPTIONAL INPUTS:
; fileName: A string containing the name of the file from which
; to load the image. If no name is given, the current
; file name associated with the object is used.
;
; SIDE EFFECTS:
; The previous data and header information is lost.
;
; EXAMPLE:
; exampleImage->load, 'exampleImageFile.hhh'
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen, Jr. 1997
;-
;
PRO dataImage::load, filename
IF (N_PARAMS() GE 1) THEN BEGIN
self->astronomyImage::load, filename
ENDIF ELSE BEGIN
self->astronomyImage::load
ENDELSE
fileExists = N_ELEMENTS(FINDFILE(self.fileName))+$
N_ELEMENTS(FINDFILE(self.fileName+'.hhh'))
IF (fileExists GT 0) THEN BEGIN
IF(PTR_VALID(self.dataHeader)) THEN BEGIN
expTime=SXPAR(*self.dataHeader,'EXPTIME')
self.exposureTime=expTime
ENDIF
ENDIF
END
;----------------------------------------------------------------------------
; dataImage::fourierTransform
;
;+
; NAME:
; dataImage::fourierTransform
;
; PURPOSE:
; This method replaces the data associated with an image with it
; fourier transform
;
; CALLING SEQUENCE:
; oDataImage->fourierTransform
;
; KEYWORD PARAMETERS:
; INVERSE: If this keyword is set, the method takes the inverse
; fourier transform.
;
; SIDE EFFECTS:
; After the transform is taken, the image data is complex.
;
; PROCEDURE:
; This method uses the IDL FFT funtion.
;
; EXAMPLE:
; oDataImage->fourierTransform
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
;
PRO dataImage::fourierTransform, INVERSE=inverse
IF PTR_VALID(self.data) THEN BEGIN
complexFFT = FFT(self->getData(),-1,INVERSE=inverse)
self->setData, complexFFT
ENDIF
END
;----------------------------------------------------------------------------
; dataImage::powerSpectrum
;+
; NAME:
; dataImage::powerSpectrum
;
; PURPOSE:
; This method replaces the image data by its power spectrum.
;
; CALLING SEQUENCE:
; oDataImage->powerSpectrum
;
; EXAMPLE:
; To plot the power spectrum of an image
; testImage=OBJ_NEW('dataImage')
; testImage->load,'imageFile.hhh'
; testImage->powerSpectrum
; TVSCL,ALOG(testImage->getData())
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::powerSpectrum
IF PTR_VALID(self.data) THEN BEGIN
self->fourierTransform
pspec = (ABS(self->getData()))^2.0
self->setData, pspec
ENDIF
END
;----------------------------------------------------------------------------
; dataImage::frequencyFilter
;+
; NAME:
; dataImage::frequencyFilter
;
; PURPOSE:
; This method replaces the data with the data passed through a
; fourier frequency filter.
;
; CALLING SEQUENCE:
; oDataImage->frequencyFilter, Filter
;
; INPUTS:
; Filter: This parameter is an array with the frequency filter
; to be applied to the image. If the size of the filter
; is greater than that of the data array, the data is
; padded with zeros before the filtering, and cut back
; to its original size afterwards.
;
; KEYWORD PARAMETERS:
; MASK: Either a mask image or an array of the same size as
; the data image, where bad pixels are marked with ones
; and good pixels with zeros. Masked pixels are not used
; in the creation of the filtered image. See the
; PROCEDURE section of this documentation for a
; description of how masked images are handled.
;
; PROCEDURE:
; The filtered image is created by multiplying the fourier
; transform of the data image by the filter, and useing the
; inverse transform to bring it back into image space. To
; understand the handling of the mask, if provided, it is
; easiest to view the filter as a convolution. In the unmasked
; case, each pixel in the returned image is the weighted average
; of the surrounding pixels in the original image, where the
; weights are determined by the fourier transform of the
; filter. If a mask is provided, the weights of the masked
; pixels are set to zero. This procedure is not mathematically
; justified in the general case, but does provide good results
; where the filter lets through spacial frequencies that are
; large compared to the size of the masked regions.
;
; EXAMPLE:
; To apply an "ideal" filter to an image,
; idealFilter=FLTARR(256,256)
; idealFilter(WHERE(DIST(256,256) LT 100)) = 1
; oDataImage->frequencyFilter, idealFilter
;
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::frequencyFilter, filter, MASK=mask
data=self->getData()
filterSize=SIZE(filter)
dataSize=SIZE(data)
extend = NOT ((dataSize(1) EQ filterSize(1)) AND $
(dataSize(2) EQ filterSize(2)))
IF extend THEN BEGIN
newData=FLTARR(filterSize(1),filterSize(2))
newData(0:dataSize(1)-1,0:dataSize(2)-1)=TEMPORARY(data)
data=TEMPORARY(newData)
newMask=FIX(data*0)+1
newMask(0:dataSize(1)-1,0:dataSize(2)-1)=TEMPORARY(mask)
mask=TEMPORARY(newMask)
ENDIF
IF KEYWORD_SET(mask) THEN BEGIN
IF ((SIZE(mask))[0] EQ 0) THEN BEGIN
inverseMask=FLOAT(mask->getData() EQ 0)
ENDIF ELSE BEGIN
inverseMask=FLOAT(mask EQ 0)
ENDELSE
data=inverseMask*TEMPORARY(data)
inverseMask=FFT(COMPLEX(inverseMask),/OVERWRITE)
inverseMask=TEMPORARY(inverseMask)*filter
inverseMask=FLOAT(FFT(COMPLEX(inverseMask),/INVERSE,/OVERWRITE))
ENDIF
data=FFT(COMPLEX(data),/OVERWRITE)
data=TEMPORARY(data)*filter
data=FLOAT(FFT(COMPLEX(data),/INVERSE,/OVERWRITE))
IF extend THEN BEGIN
data=data(0:dataSize(1)-1,0:dataSize(2)-1)
inverseMask=inverseMask(0:dataSize(1)-1,0:dataSize(2)-1)
ENDIF
IF KEYWORD_SET(mask) THEN data=TEMPORARY(data)/inverseMask
self->setData, data
END
;----------------------------------------------------------------------------
; dataImage::pixelWavenumber
;+
; NAME:
; dataImage::pixelWavenumber
;
; PURPOSE:
; Returns an array, where the value of each pixel corresponds to
; the wavenumber represented by that pixel in a fourier
; transformed image.
;
; CALLING SEQUENCE:
; Result = aDataImage->pixelWavenumber()
;
; KEYWORD PARAMETERS:
; EXTENDEDSIZE: If the image to be considered is larger than
; the data image, this keyword can be set to the
; size desired. Smaller values are not allowed.
;
; OUTPUTS:
; The result is a two dimensional array of wavenumber values,
; suitable for use in creating frequency filters.
;
; PROCEDURE:
; This procedure is a wrapper for the DIST function in IDL.
;
; EXAMPLE:
; To apply an "ideal" filter to an image,
; idealFilter=FLTARR(256,256)
; idealFilter(oDataImage->pixelWavenumber() LT 100)) = 1
; oDataImage->frequencyFilter, idealFilter
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
; September, 1998 Modified to use DIST to calculte wavenumbers
;-
FUNCTION dataImage::pixelWavenumber, EXTENDEDSIZE=extendedSize
imageSize=self->size()
IF KEYWORD_SET(extendedSize) THEN imageSize=(imageSize > extendedSize)
wavenumber = DIST(imageSize[0],imageSize[1])
RETURN,wavenumber
END
;----------------------------------------------------------------------------
; dataImage::idealFilter
;+
; NAME:
; dataImage::idealFilter
;
; PURPOSE:
; This method applies an ideal filter to the image.
;
; CALLING SEQUENCE:
; oDataImage->idealFilter(Cutoff)
;
; INPUTS:
; Cutoff: This parameter is the value of the lowest (highest)
; spatial wavelength allowed through the filter.
;
; KEYWORD PARAMETERS:
; LOWPASS: By default, the image is highpass filtered. If this
; keyword is set, it is lowpass filtered instead. This
; correspondins to smoothing the image.
; MASK: One may provide a pixel mask to remove the influence
; unwanted pixels. See the documentation for
; dataImage::frequencyFilter for the handling of masks.
;
; PROCEDURE:
; The Cutoff us converted from a wavelength to a wavenumber, the
; ideal filter is constructed, and then applied to the image data.
;
; EXAMPLE:
; To apply an ideal smoothing filter,
; exampleImage->idealFilter(10,/LOWPASS)
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::idealFilter, cutoff, LOWPASS=lowpass, MASK=mask
dataSize=self->size()
xsize=dataSize(0)
wavenumber=xsize/cutoff
IF KEYWORD_SET(lowpass) THEN $
filter=FLOAT(self->pixelWavenumber() GT wavenumber) ELSE $
filter=FLOAT(self->pixelWavenumber() LT wavenumber)
self->frequencyFilter, filter, MASK=mask
END
;----------------------------------------------------------------------------
; dataImage::segment
;+
; NAME:
; dataImage::segment
;
; PURPOSE:
; Segment the data image. Pixels above the segmentation value
; are set to one, and those below it are set to zero. Often,
; one would want to use the maskImage::segment instead.
;
; CALLING SEQUENCE:
; oDataImage->segment, Cutoff
;
; INPUTS:
; Curoff: The cutoff value for the segmentation.
;
; EXAMPLE:
; To create a map of the saturated pixels in an image where the
; saturation value is 4000,
; exampleDataImage->segment,400
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::segment, cutoff
newData=self->getData() GT cutoff
self->setData, newData
END
;----------------------------------------------------------------------------
; dataImage::butterworthFilter
;+
; NAME:
; dataImage::butterworthFilter
;
; PURPOSE:
; This method applies a Butterworth filter to the image.
;
; CALLING SEQUENCE:
; oDataImage->butterworthFilter(Cutoff,Order)
;
; INPUTS:
; Cutoff: This parameter is the value of the lowest (highest)
; spatial wavelength allowed through the filter.
; Order: The order of the butterworth filter. Higher values
; will results in more ringing, and a sharper boundry
; between filtered and unfiltered frequencies.
;
; KEYWORD PARAMETERS:
; LOWPASS: By default, the image is highpass filtered. If this
; keyword is set, it is lowpass filtered instead. This
; correspondonds to smoothing the image.
; MASK: One may provide a pixel mask to remove the influence
; unwanted pixels. See the documentation for
; dataImage::frequencyFilter for the handling of masks.
; EXTENDEDSIZE: To reduce edge effects, the filtering may be
; performed on an image which has been padded with
; zeros, which are then masked. This parameter sets the
; size of the padded image.
;
; PROCEDURE:
; The Cutoff us converted from a wavelength to a wavenumber, the
; butterworth filter is constructed, and then applied to the image data.
;
; EXAMPLE:
; To apply an butterworth smoothing filter,
; exampleImage->butterworthFilter(10,2,/LOWPASS)
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::butterworthFilter, cutoff, order,LOWPASS=lowpass, MASK=mask,$
EXTENDEDSIZE=extendedSize
imageSize=self->size()
IF KEYWORD_SET(EXTENDEDSIZE) THEN imageSize=imageSize > extendedSize
wavenumber=imageSize(0)/cutoff
filter=self->pixelWavenumber(EXTENDEDSIZE=extendedSize)
IF NOT KEYWORD_SET(lowpass) THEN order=-1*order
filter=1.0/(1.0+(TEMPORARY(filter)/wavenumber)^(2.0*order))
self->frequencyFilter, filter, MASK=mask
END
;----------------------------------------------------------------------------
; dataImage::wienerFilter
;+
; NAME:
; dataImage::wienerFilter
;
; PURPOSE:
; This method replaces the method by a Wiener filtered version
; of the image. For information on using a Wiener filter,
; consult a image processing reference guide, _Digital Image
; Processing_ by Gonzales & Woods (1992), or _The Image
; Processing Handbook_ by Russ (1995).
;
; CALLING SEQUENCE:
; oDataImage->WienerFilter,Psf,K
;
; INPUTS:
; Psf: A dataImage object containing the point spread function
; of the image.
; k: The parameter "k" of the Wiener filter.
;
; RESTRICTIONS:
; There is no accomodation for handling noise or PSF's which
; vary over the image.
;
; PROCEDURE:
; The Wiener filter is constructed using the supplied PSF and K,
; and then applied to the image.
;
; EXAMPLE:
; To wiener filter an image,
; psfImage=OBJ_NEW('dataImage')
; psfImage->load, 'imagePSF.hhh'
; dataImage=OBJ_NEW(('dataImage')
; dataImage->load, 'starfield1.hhh'
; k=1.0e-7
; dataImage->wienerFilter,psfImage,k
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::wienerFilter, psf, k
psfPspec=OBJ_NEW("dataImage")
psfPspec->setData, psf->getData()
psfPspec->powerSpectrum
pspec=psfPspec->getData()
OBJ_DESTROY, psfPspec
filter=pspec/(pspec+k)
psfFFT=OBJ_NEW("dataImage")
psfFFT->setData, psf->getData()
psfFFT->fourierTransform
numberOfPixels=FLOAT((psfFFT->size())(0)*(psfFFT->size())(1))
filter=numberOfPixels*filter/ABS(psfFFT->getData())
OBJ_DESTROY, psfFFT
filter=FFT(filter,/OVERWRITE)
filter=filter/TOTAL(filter)
filter=FFT(filter,/INVERSE,/OVERWRITE)
self->frequencyFilter, filter
END
;----------------------------------------------------------------------------
; dataImage::psfFilter
;+
; NAME:
; dataImage::psfFilter
;
; PURPOSE:
; Frequency filters the image to emphasize point sources. It is
; intended to emphasize point source in front of extended
; sources, such as galaxies. If the galaxy has a bright, point
; source like center, it is often helpful to smooth it out
; before the appliation of this filter.
;
; CALLING SEQUENCE:
; oDataImage->psfFilter, Psf, Flux
;
; INPUTS:
; Psf: A dataImage object with the point spread function of the image.
; Flux: An estimate of the total flux of point sources in the image.
;
; KEYWORD PARAMETERS:
; MODEL: A dataImage object containing a model for the
; background on which you are detecting point sources.
; SKY: The sky value. If this is included already in the
; model, this keyword should not be set.
;
; PROCEDURE:
; This method uses frequency filtering to emphasize point
; sources. It begins by calculating the power spectrum of the
; point sources, using the psf and an estimate of their total
; flux, and the power spectrum of everything else, using the
; model (where available) or the power spectrum of the data
; itself (which assumes the background dominates the point
; sources). A frequency filter is then constructed which lets
; though only those frequencies where the power from the point
; sources is largest compared with other sources. This filter is
; essentially a Wiener filter where everything but the point
; sources is reguareded as noise.
;
; EXAMPLE:
; oDataImage->psfFilter, psf, estimatedFlux
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::psfFilter, psf, f, MODEL=model, SKY=sky
psfPspec=OBJ_NEW("dataImage")
psfPspec->setData, psf->getData()
psfPspec->powerSpectrum
pspec=psfPspec->getData()
OBJ_DESTROY, psfPspec
selfPspec=OBJ_NEW("dataImage")
IF NOT KEYWORD_SET(sky) THEN sky=0
IF KEYWORD_SET(model) THEN BEGIN
selfPspec->setData, (model->getData()+sky)
ENDIF ELSE BEGIN
selfPspec->setData, (self->getData())
ENDELSE
selfPspec->powerSpectrum
nPspec=selfPspec->getData()
OBJ_DESTROY,selfPspec
filter=pspec/(pspec+(nPspec/(f^2.0)))
psfFFT=OBJ_NEW("dataImage")
psfFFT->setData, psf->getData()
psfFFT->fourierTransform
numberOfPixels=FLOAT((psfFFT->size())(0)*(psfFFT->size())(1))
filter=numberOfPixels*filter/ABS(psfFFT->getData())
OBJ_DESTROY, psfFFT
self->frequencyFilter, filter
END
;----------------------------------------------------------------------------
; dataImage::normalize
;+
; NAME:
; dataImage::normalize
;
; PURPOSE:
; This method normalizes the total flux of the image to a set
; value. This is particularly useful for normalizing PSF images
; to one count.
;
; CALLING SEQUENCE:
; oDataImage->normalize
;
; KEYWORD PARAMETERS:
; TOTAL: Set this keyword to the desired total value of pixels
; in the image, if different from one.
;
; EXAMPLE:
; psfDataImage=OBJ_NEW('dataImage')
; psfDataImage->load, 'psfImageFile.hhh'
; psfDataImage->normalize
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::normalize, TOTAL=total
IF NOT KEYWORD_SET(total) THEN total=1.0
dataArray=self->getData()
currentTotal=TOTAL(dataArray)
dataArray=TEMPORARY(dataArray)*FLOAT(total)/FLOAT(currentTotal)
self->setData, dataArray
END
;----------------------------------------------------------------------------
; dataImage::resize
;+
; NAME:
; dataImage::resize
;
; PURPOSE:
; This method resizes the image. It does not change the pixel
; scale; it merely cuts off or pads the image to make it the
; desired size.
;
; CALLING SEQUENCE:
; oDataImage->resize, NewSize
;
; INPUTS:
; newSize: The dimensions of the new image.
;
; RESTRICTIONS:
; At present, images cannot be resized such that one dimension
; is larger than the original, while the other is smaller.
;
; EXAMPLE:
; To ensure the image is 256x256,
; oDataImage->resize, 256
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::resize, newSize
IF ((SIZE(newSize))[0] EQ 0) THEN newSize=[newSize,newSize]
newArray=FLTARR(newSize(0),newSize(1))
oldSize=self->size()
IF (oldSize(0) LT newSize(0)) AND (oldSize(1) LT newSize(1)) THEN BEGIN
starts=[0,0]
starts=(newSize-oldSize)/2
newArray[starts(0):starts(0)+oldSize(0)-1,starts(1):starts(1)$
+oldSize(1)-1]=self->getData()
ENDIF ELSE BEGIN
middle=oldSize/2
newArray=(self->getData())(middle(0)-newSize(0)/2:middle(0)+newSize(0)/2,$
middle(1)-newSize(1)/2:middle(1)+newSize(1)/2)
ENDELSE
self->setData,newArray
END
;----------------------------------------------------------------------------
; dataImage::center
;+
; NAME:
; dataImage::center
;
; PURPOSE:
; This method moves the center (or some specified coordinate) to
; the corner of the image. This is particularly useful for
; working with PSF images.
;
; CALLING SEQUENCE:
; oDataImage->center
;
; KEYWORD PARAMETERS:
; PIXELCOORDS: This keyword holds the coordiantes to be placed
; at the corner of the image. If not specified, the
; center is used.
;
; EXAMPLE:
; psfDataImage->center
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::center, PIXELCOORDS=pixelcoords
imageSize=self->size()
IF NOT KEYWORD_SET(pixelcoords) THEN BEGIN
pixelcoords=[imageSize(0)/2-1,imageSize(1)/2-1]
ENDIF
newImage=SHIFT(self->getData(),pixelcoords[0],pixelcoords[1])
self->setData,newImage
END
;----------------------------------------------------------------------------
; dataImage::copy
;+
; NAME:
; dataImage::copy
;
; PURPOSE:
; This method is identical to the dataImage::clone method
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
FUNCTION dataImage::copy
newImage=self->clone()
RETURN, newImage
END
;----------------------------------------------------------------------------
; dataImage::radius
;+
; NAME:
; dataImage::radius
;
; PURPOSE:
; This method replaces the data with an image where each pixel
; contains the distance of that pixel to a specified coordinate.
; It may be used, for exampled, to get distances from the center
; of a galaxy,
;
; CALLING SEQUENCE:
; oDataImage->radius, CenterCoord
;
; INPUTS:
; CenterCoord: The coordinate from the the distances are to be
; calculated.
;
; EXAMPLE:
; To create an image of the distances to the center of a galaxy
; centered at [410,400] in the image galaxyImage,
; radiusImage=galaxyImage->clone()
; radiusImage->radius, [410,400]
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
FUNCTION dataImage::radius, centerCoord
imageSize=self->size()
xArray=ROTATE(REPLICATE(1,imageSize(1))#FINDGEN(imageSize(0)),3)
yArray=ROTATE(xArray,1)
xArray=xArray-centerCoord(0)
yArray=yArray-centerCoord(1)
radius=SQRT(xArray*xArray+yArray*yArray)
RETURN,radius
END
;----------------------------------------------------------------------------
; dataImage::angle
;+
; NAME:
; dataImage::angle
;
; PURPOSE:
; This method replaces the data with an image where each pixel
; has the theta value of that pixels polar coordinate, where the
; origin is at the specified coordinate. Together with the
; dataImage::radius method, it can be used to get the polar
; coordinates of each pixel.
;
; CALLING SEQUENCE:
; oDataImage->angle, Pixelcoordinate
;
; INPUTS:
; Pixelcoordinate: A vector with the coordinates of the polar
; coordinate systems center.
;
; EXAMPLE:
; exampleImage->angle, [400,400]
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
FUNCTION dataImage::angle, centerCoord
imageSize=self->size()
r=self->radius(centerCoord)
xArray=ROTATE(REPLICATE(1,imageSize(1))#FINDGEN(imageSize(0)),3)
yArray=ROTATE(xArray,1)
xArray=xArray-centerCoord(0)
yArray=yArray-centerCoord(1)
fr=(r > 0.000001)
angle=ACOS(xArray/fr)
bottom=WHERE(yArray LT 0)
angle(bottom)=!PI+ACOS(-1*xArray(bottom)/fr(bottom))
RETURN,angle
END
;----------------------------------------------------------------------------
; dataImage::triangulate
;+
; NAME:
; dataImage::triangulate
;
; PURPOSE:
; This method replaces masked pixels with values calculated
; using triangulation from nearby food pixels.
;
; CALLING SEQUENCE:
; oDataImage->trinagulate, Mask
;
; INPUTS:
; Mask: The mask image (a maskImage object), where the pixels
; to be triangulated are masked.
;
; RESTRICTIONS:
; Triangulation of pixels is limited by the limits of IDL's
; trinagulation routines. The method can be quite slow.
;
; PROCEDURE:
; This pixels surrounding regions of bad pixels are used to
; create a triangulation grid using the IDL procedure
; TRIANGULATE. Replacement values for the bad pixels are then
; calculated using TRIGRID.
;
; EXAMPLE:
; To view an image with the masked pixels replaced:
; oDataImage=OBJ_NEW('dataImage')
; oDataImage->load,'galaxyData.hhh'
; oMaskImage=OBJ_NEW('maskImage')
; oMaskImage->load,'galaxyMask.hhh'
; oDataImage->tringulate, oMaskImage
; TV,oDataImage->getData(/DISPLAY)
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::triangulate, mask
imageSize=self->size()
xArray=ROTATE(REPLICATE(1,imageSize(1))#FINDGEN(imageSize(0)),3)
yArray=ROTATE(xArray,1)
dims=[MIN(xArray),MIN(yArray),MAX(xArray),MAX(yArray)]
data=self->getData()
;triangulate only using points which will be useful
kernel=FLTARR(5,5)+1
mask2=OBJ_NEW("maskImage")
mask2->setData, mask->getData()
mask2->dilate, KERNEL=kernel
goodPoints=WHERE(((mask->getData()) EQ 0) AND ((mask2->getData()) GT 0))
OBJ_DESTROY, mask2
xArray=xArray(goodPoints)
yArray=yArray(goodPoints)
goodData=data(goodPoints)
triangles=1
b=1
TRIANGULATE, xArray, yArray, triangles, b
newData=data
newData(*,*)=TRIGRID(xArray,yArray,goodData,triangles,$
[1,1],dims,EXTRAPOLATE=b)
maskedPoints=WHERE((mask->getData()) GT 0)
data(maskedPoints)=newData(maskedPoints)
self->setData,data
END
;----------------------------------------------------------------------------
; dataImage::ellipseSpace
;+
; NAME:
; dataImage::ellipseSpace
;
; PURPOSE:
; This method creates a model of the galaxy using fits in
; (intensity, radius, angle) space. The center of the galaxy
; must be provided. This method works well for elliptical
; galaxies where the center of the isophotes varies little with
; the intensity of the isophote.
;
; CALLING SEQUENCE:
; oDataImage->ellipseSpace, CenterCoord, TransformSize
;
; INPUTS:
; CenterCoord: A two element array giving the x,y coordinates
; of the center of the galaxy in the image.
; TransformSize: The image of the galaxy is often at much higher
; resolution than necessary for fitting
; isophotes; transformSize specifies has small an
; array to regrid the image to before examination
; of the isophotes.
;
; KEYWORD PARAMETERS:
; MASKIMAGE: An object of type maskImage with the pixels not
; to be used in the calculation of the isophotes
; masked.
; FITTERMS: The number of terms to use in polynomial
; fitting of log(radius) to log(intensity).
; FOURIERTERMS: The number of terms to retain the the fourier
; expansion of radius as a function of position
; angle.
; SKY: An estimate of the sky background value. This
; is important primarily when FITTERMS is also
; set.
; VERBOSE: If this keyword is set, the progress of the
; method is reported at several points in the
; process. This can be interesting for the
; impatient.
;
; RESTRICTIONS:
; The completeness of the mask is vital to the success of the
; fitting; a small number of high or low pixels may cause very
; serious corruption in the result. For this reason, it is often
; helpful to median filter (or otherwise smooth) the image
; before it is fit.
;
; PROCEDURE:
; For each unmasked pixel in the resized image, the angle,
; radius, and intensity are calculated. Triangulation is then
; used to calculate an array where the index of one dimension is
; proportional to the log of the intensity, the index of the
; other dimension is proportional to the position angle, and the
; value of each pixel is the triangulated radius. Depending of
; the keyword parameters set, the lines of constant intensity
; may be low pass filtered, the lines of constant angle may
; be replaced by best fits to a (log) polynomial, or
; both. Triangulation is then used to return to the original
; (x,y) coordinate system. For additional details, see Neilsen
; (1998).
;
; EXAMPLE:
; oDataImage->ellipseSpace, [402,410], [64,64], MASKIMAGE=mask,$
; FITTERMS=2, FOURIERTERMS=4, SKY=1.8
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
PRO dataImage::ellipseSpace, centerCoord, transformSize,$
MASKIMAGE=maskImage, FITTERMS=fitTerms,$
FOURIERTERMS=fourierTerms, SKY=sky,$
VERBOSE=verbose
imageSize=self->size()
radius=CONGRID((self->radius(centerCoord)>0.00001),transformSize(0),transformSize(1))
logRadius=ALOG(radius)
angle=CONGRID(self->angle(centerCoord),transformSize(0),transformSize(1))
intensity=CONGRID((self->getData()>0.00001),transformSize(0),transformSize(1))
logIntensity=ALOG(intensity)
IF KEYWORD_SET(maskImage) THEN BEGIN
goodPoints=$
WHERE(CONGRID(maskImage,transformSize(0),transformSize(1)) EQ 0 )
logRadius=logRadius(goodPoints)
angle=angle(goodPoints)
logIntensity=logIntensity(goodPoints)
ENDIF
profileEllipse=(1 EQ 0)
IF KEYWORD_SET(fitTerms) THEN BEGIN
profileEllipse=fitTerms
fitTerms = ABS(fitTerms)
ENDIF
smallAngle=WHERE(angle LE !PI/8.0)
largeAngle=WHERE(angle GE 15.0*!PI/8.0)
smallWrapAngle=angle(smallAngle)+2.0*!PI
largeWrapAngle=angle(largeAngle)-2.0*!PI
smallWrapLogRadius=logRadius(smallAngle)
largeWrapLogRadius=logRadius(largeAngle)
smallWrapLogIntensity=logIntensity(smallAngle)
largeWrapLogIntensity=logIntensity(largeAngle)
angle=[smallWrapAngle,angle,largeWrapAngle]
logRadius=[smallWrapLogRadius,logRadius,largeWrapLogRadius]
logIntensity=[smallWrapLogIntensity,logIntensity,largeWrapLogIntensity]
TRIANGULATE,angle,logIntensity,tr,b
gridSpacing=[2*!PI/(transformSize(0)-1),$
(MAX(logIntensity)-MIN(logIntensity))/(transformSize(1)-1)]
gridLimits=[0,MIN(logIntensity),2.0*!PI,MAX(logIntensity)]
newArray=TRIGRID(angle,logIntensity,logRadius,tr,gridSpacing,gridLimits)
originalArray=newArray
IF KEYWORD_SET(fitTerms) THEN BEGIN
IF KEYWORD_SET(VERBOSE) THEN PRINT,"Profile Fitting"
FOR count=0,transformSize(0)-1 DO BEGIN
line=newArray(count,*)
index=FINDGEN(transformSize(0))
goodPoints=WHERE((line NE 0) AND (index GE ALOG(5.0)),npoints)
IF (npoints GE 8) THEN BEGIN
fit=POLY_FIT(index(goodPoints),line(goodPoints),fitTerms)
fitArray=POLY(index,fit)
IF KEYWORD_SET(fourierTerms) THEN BEGIN
fitArray(goodPoints)=newArray(count,goodPoints)
ENDIF
newArray(count,*)=fitArray
ENDIF
ENDFOR
ENDIF
IF (KEYWORD_SET(fourierTerms) AND KEYWORD_SET(fitTerms)) THEN BEGIN
IF KEYWORD_SET(VERBOSE) THEN PRINT,"Fourier Filtering"
FOR count=0,transformSize(0)-1 DO BEGIN
line=EXP(newArray(*,count))
line=FFT(line,/OVERWRITE)
freqIndex=FINDGEN(transformSize(0)) < REVERSE(FINDGEN(transformSize(0)))
line(WHERE(freqIndex GT fourierTerms)) = 0.0
line=FFT(line,/INVERSE,/OVERWRITE)
newArray(*,count)=ALOG(line)
ENDFOR
IF KEYWORD_SET(fitTerms) THEN BEGIN
IF KEYWORD_SET(VERBOSE) THEN PRINT,"Profile Fitting II"
FOR count=0,transformSize(0)-1 DO BEGIN
line=newArray(count,*)
index=FINDGEN(transformSize(0))
goodPoints=WHERE((originalArray(count,*) NE 0),npoints)
IF (npoints GE 8) THEN BEGIN
fit=POLY_FIT(index(goodPoints),line(goodPoints),fitTerms)
fitArray=POLY(index,fit)
IF (profileEllipse GE 0) THEN BEGIN
fitArray(goodPoints)=newArray(count,goodPoints)
ENDIF
newArray(count,*)=fitArray
ENDIF
ENDFOR
ENDIF
ENDIF
newAngle=FINDGEN(transformSize(0))*2.0*!PI/transformSize(0)
newLogIntensity=MIN(logIntensity)+FINDGEN(transformSize(1))*$
(MAX(logIntensity)-MIN(logIntensity))/transformSize(1)
newAngle=newAngle#REPLICATE(1,transformSize(1))
newLogIntensity=REPLICATE(1,transformSize(0))#newLogIntensity
IF KEYWORD_SET(VERBOSE) THEN PRINT,"Calculating x, y, intensity"
x=EXP(newArray)*COS(newAngle)+centerCoord(0)
y=EXP(newArray)*SIN(newAngle)+centerCoord(1)
intensity=EXP(newLogIntensity)
x=x+(FLOAT(imageSize(0))/FLOAT(transformSize(0)))/2
y=y+(FLOAT(imageSize(1))/FLOAT(transformSize(1)))/2
IF KEYWORD_SET(VERBOSE) THEN PRINT,"finding triangulation"
TRIANGULATE,x,y,tr,b
gridSpacing=[FLOAT(imageSize(0))/FLOAT(transformSize(0)),$
FLOAT(imageSize(1))/FLOAT(transformSize(1))]
gridLimits=[0,0,imageSize(0)-1,imageSize(1)-1]
IF KEYWORD_SET(VERBOSE) THEN PRINT, "Grid Spacing ",gridSpacing
IF KEYWORD_SET(VERBOSE) THEN PRINT, "Grid Limits ",gridLimits
IF KEYWORD_SET(VERBOSE) THEN PRINT,"triangulating"
newArray=TRIGRID(x,y,intensity,tr,gridSpacing,gridLimits)
IF KEYWORD_SET(VERBOSE) THEN PRINT,"regriding"
newArray=CONGRID(newArray,imageSize(0),imageSize(1),CUBIC=-0.5)
IF KEYWORD_SET(sky) THEN newArray=newArray-sky
self->setData,newArray
END
;----------------------------------------------------------------------------
; dataImage::exponentialFit
;+
; NAME:
; dataImage::exponentialFit
;
; PURPOSE:
; Fit in exponential profile to a section of an image of a
; galaxy, and return the best fits.
;
; CALLING SEQUENCE:
; Result=oDataImage->exponentialFit(Centercoord,Anglerange,Radiusrange)
;
; INPUTS:
; Cenctercoord: A two element array containing the x,y
; coordinates of the center of the galaxy.
; Anglerange: The range of angles in which to fit the
; profile. This range should be small enough that
; the pixels at a given radius will have roughly
; the same flux. The angle is measure in radians.
; Radiusrange: This is the range of radii over which to make
; the fit. Particularly in double exponential
; profiles, better sky estimates can be made from
; large radii only.
;
; KEYWORD PARAMETERS:
; MASKIMAGE: An array or maskImage object with a mask
; covering points not to be used in the fit.
; INITIAL: A 3 element array, the initial parameters of
; the fits. The first element is an estimate of
; the sky, and is usually all that is
; necessary. At the completion of the
; procedure. it will hold the best exponential
; fit parameters.
; RANGEMASK: A mask covering all pixels not used in the
; profile fit, of the same type as MASKIMAGE.
;
; OUTPUTS:
; Fitsky: The returned value of the function is an array
; with the fits as follows:
; FitSky[0,*] = the radii of the fit points
; FitSky[1,*] = the fluxes of the fit points
; FitSky[2,*] = the fit value of the best fit
; exponential fit
; FitSky[3,*] = the fit value of the best fit
; R^(1/4)
;
; RESTRICTIONS:
;
; PROCEDURE:
; Arrays of good pixels and the corresponding radii are
; created. Least squares fits are then made at an expnential
; plus a constant (sky), and an r^(1/4) law plus a constant.
;
; EXAMPLE:
; To plot the profile and the different fits:
; skyFit=oDataImage->exponentialFit, [401,410], [0.0,0.52],$
; [30,500], MASKIMAGE=imageMask, INITIAL=[4.5,0.0,0.0]
; PLOT,skyFit[0,*],skyFit[1,*],PSYM=3
; OPLOT,skyFit[0,*],skyFit[2,*]
; OPLOT,skyFit[0,*],skyFit[3,*],LINE=2
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
FUNCTION dataImage::exponentialFit, centerCoord, angleRange, radiusRange, $
MASKIMAGE=maskImage, INITIAL=initial, $
RANGEMASK=rangeMask
imageSize=self->size()
radius=self->radius(centerCoord)>0.00001
angle=self->angle(centerCoord)
intensity=self->getData()>0.00001
IF KEYWORD_SET(maskImage) THEN BEGIN
IF (SIZE(maskImage))[0] EQ 0 THEN BEGIN
maskImageArray=maskImage->getData()
ENDIF ELSE BEGIN
maskImageArray=maskImage
ENDELSE
goodPoints = WHERE( (radius GE radiusRange(0)) $
AND (radius LE radiusRange(1)) $
AND (angle GE angleRange(0)) $
AND (angle LE angleRange(1)) $
AND (maskImageArray EQ 0) )
ENDIF ELSE BEGIN
goodPoints = WHERE( (radius GE radiusRange(0)) $
AND (radius LE radiusRange(1)) $
AND (angle GE angleRange(0)) $
AND (angle LE angleRange(1)) )
ENDELSE
rangeMaskArray=maskImageArray*0 + 1
rangeMaskArray(goodPoints)=0
IF KEYWORD_SET(rangeMask) THEN BEGIN
IF (SIZE(rangeMask))[0] EQ 0 THEN BEGIN
IF OBJ_VALID(rangeMask) THEN BEGIN
rangeMask->setData,rangeMaskArray
ENDIF ELSE BEGIN
rangeMask=rangeMaskArray
ENDELSE
ENDIF ELSE BEGIN
rangeMask=rangeMaskArray
ENDELSE
ENDIF ELSE BEGIN
IF KEYWORD_SET(maskImage) THEN BEGIN
IF (SIZE(maskImage))[0] EQ 0 THEN BEGIN
IF OBJ_VALID(rangeMask) THEN BEGIN
rangeMask->setData,rangeMaskArray
ENDIF ELSE BEGIN
rangeMask=rangeMaskArray
ENDELSE
ENDIF ELSE BEGIN
rangeMask=rangeMaskArray
ENDELSE
ENDIF ELSE BEGIN
rangeMask=OBJ_NEW('maskImage')
rangeMask->setData,rangeMaskArray
ENDELSE
ENDELSE
radius=radius(goodPoints)
angle=angle(goodPoints)
intensity=intensity(goodPoints)
TRIANGULATE,angle,radius,tr,b
angleWidth=(angleRange(1)-angleRange(0))/3.0
gridSpacing=[angleWidth,1.0]
gridLimits=[angleRange(0),radiusRange(0),$
angleRange(1),radiusRange(1)]
newArray=TRIGRID(angle,radius,intensity,tr,gridSpacing,gridLimits)
expTerms=[20.0,20000.0,1.0]
goodPoints=WHERE(newArray(1,*) GT 0,numGoodPoints)
npoints = numGoodPoints
radii=radiusRange(0)+FINDGEN(radiusRange(1)-radiusRange(0)+1)
intensities=radii
intensities(*)=newArray(1,*)
goodPoints=WHERE(intensities GT 0,numGoodPoints)
PRINT,numGoodPoints," data points"
radii=radii(goodPoints)
intensities=intensities(goodPoints)
weights=1.0/intensities
expSigma=[1.0,1.0,1.0]
expChisq=1
IF KEYWORD_SET(initial) THEN BEGIN
IF (initial(1) EQ 0) THEN BEGIN
logRadii=ALOG(radii)
logIntensities=ALOG(intensities-initial(0))
testFit=POLY_FIT(logRadii,logIntensities,1)
expTerms(0)=initial(0)
expTerms(1)=EXP(testFit(0))
expTerms(2)=-1.0*testFit(1)
ENDIF ELSE BEGIN
expTerms=initial
ENDELSE
ENDIF
expFit = CURVEFIT(radii,intensities,weights,expTerms,expSigma,$
FUNCTION_NAME='expProfile', CHISQ=expChisq,ITMAX=10000)
PRINT,"Exp. Sky:",expTerms(0),"+/-",expSigma(0), " R. Chi^2:",$
expChisq/(1.0*numGoodPoints-3)
fitSky=FLTARR(4,npoints)
fitSky(0,*)=radii
fitSky(1,*)=intensities
fitSky(2,*)=expFit
;===Estimate r^1/4 parameters from exponential parameters
r=0.5*(MAX(radii)+MIN(radii))
rqTerms=[expTerms(0),$
expTerms(1)*(r^(-1.0*expTerms(2)))*EXP(r^0.25),1.0]
rqSigma=[1.0,1.0,1.0]
rqFit = CURVEFIT(radii,intensities,weights,rqTerms,rqSigma,$
FUNCTION_NAME='rqProfile', CHISQ=rqChisq,ITMAX=10000)
PRINT,"R^(1/4) Sky:",rqTerms(0),"+/-",rqSigma(0), " R. Chi^2:",$
rqChisq/(1.0*numGoodPoints-3)
fitSky(3,*)=rqFit
RETURN,fitSky
END
;----------------------------------------------------------------------------
; dataImage::pspecFit
;+
; NAME:
; dataImage::pspecFit
;
; PURPOSE:
; Fit the power spectrum of the image to a linear function of
; the power spectrum of the PSF. This is primarily useful for
; separation of surface brightness fluctuation (SBF) variance
; from noise variance.
;
; CALLING SEQUENCE:
; Result=oDataImage->pspecFit(Psf)
;
; INPUTS:
; Psf: A dataImage object containing the PSF for the image.
;
; KEYWORD PARAMETERS:
; MASK: A maskImage object which masks all pixels not
; to be used.
; WAVELENGTHRANGE: A two element array giving the range of wavelength
; over which to fit the power spectra,
; POWERRANGE: A two element array giving the range of power
; over which to fit the power spectra,
; FITDATA: This keyword is set to an array with the
; results of the fit. Fitdata[0,*] is an array
; of wavenumbers, Fitdata[1,*] is the
; corresponding power, and Fitdata[2,*] is the
; corresponding fit to the power.
;
; OUTPUTS:
; This method returns a two element array. The first element is
; the constant term in the linear term, and the second element
; is the linear term.
;
; OPTIONAL OUTPUTS:
; FITDATA: see the keyword description
;
; PROCEDURE:
; The fit is performed with the IDL POLY_FIT routine.
;
; EXAMPLE:
;
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
FUNCTION dataImage::pspecFit, psf, MASK=mask, $
WAVELENGTHRANGE=wavelengthrange, POWERRANGE=powerrange,$
FITDATA=fitData
raiseFraction=1.5e13
;this is often required to prevent floating point underflows
dataPspec=OBJ_NEW("dataImage")
IF KEYWORD_SET(MASK) THEN BEGIN
maskedPoints=WHERE(mask->getData() GT 0, npoints)
IF (npoints GT 0) THEN BEGIN
newData = self->getData()
newData(maskedPoints) = 0
dataPspec->setData, newData
ENDIF
ENDIF
dataPspec->powerSpectrum
psfPspec=psf->clone()
psfPspec->powerSpectrum
wavenumberLimitsSet=KEYWORD_SET(WAVELENGTHRANGE)
IF wavenumberLimitsSet THEN BEGIN
wavenumberRange=[(self->size())(0)/wavelengthRange[1],$
(self->size())(0)/wavelengthRange[0]]
ENDIF
powerLimitsSet=KEYWORD_SET(POWERRANGE)
IF wavenumberLimitsSet THEN wavenumberLimitSet=$
(wavenumberRange(1) GT wavenumberRange(0))
IF powerLimitsSet THEN powerLimitSet = (powerRange(1) GT powerRange(0))
IF wavenumberLimitsSet AND powerLimitsSet THEN BEGIN
goodPoints=WHERE($
((dataPspec->pixelWavenumber()) GT wavenumberRange(0)) AND $
((dataPspec->pixelWavenumber()) LT wavenumberRange(1)) AND $
((dataPspec->getData()) GT powerRange(0)) AND $
((dataPspec->getData()) LT powerRange(1)) , npix)
ENDIF ELSE BEGIN
IF wavenumberLimitsSet THEN BEGIN
goodPoints=WHERE($
((dataPspec->pixelWavenumber()) GT wavenumberRange(0)) AND $
((dataPspec->pixelWavenumber()) LT wavenumberRange(1)),npix)
ENDIF ELSE BEGIN
IF powerLimitsSet THEN BEGIN
goodPoints=WHERE($
((dataPspec->pixelPower()) GT powerRange(0)) AND $
((dataPspec->pixelPower()) LT powerRange(1)), npix )
ENDIF ELSE goodPoints = WHERE( (dataPspec->getData())*0 EQ 0, npix )
ENDELSE
ENDELSE
IF (npix GT 0) THEN BEGIN
imageSize=psfPspec->size()
PRINT,"NPIX",npix
bestFit = POLY_FIT( raiseFraction*(psfPspec->getData())(goodPoints), $
(dataPspec->getData())(goodPoints), 1)
fitData=[[(dataPspec->pixelWavenumber())(goodPoints)],$
[(dataPspec->getData())(goodPoints)],$
[POLY(raiseFraction*(psfPspec->getData())(goodPoints),bestFit)]]
PRINT,"Best Fit:",bestFit
unmaskedFraction = 1.0
IF KEYWORD_SET(MASK) THEN unmaskedFraction=mask->fraction()
PRINT,"fit:",bestFit,"Fraction:",unmaskedFraction
poissonVariance = bestFit(0)/unmaskedFraction
PSFvariance = raiseFraction*bestFit(1)/$
(imageSize(0)*imageSize(1)*unmaskedFraction)
;Recall that convolution in Discrete FT space involved
;a normalization of the number of pixels
PRINT,"Poisson Noise = ",SQRT(poissonVariance)
PRINT,"SBF counts = ",PSFvariance
ENDIF ELSE BEGIN
PRINT,"No good pixels"
poissonVariance=1.0
PSFvariance=1.0
ENDELSE
OBJ_DESTROY, dataPspec
OBJ_DESTROY, psfPspec
result=[poissonVariance,PSFvariance]
RETURN,result
END
;----------------------------------------------------------------------------
; dataImage::estimateK
;+
; NAME:
; dataImage::estimateK
;
; PURPOSE:
; This uses a power spectrum fit to esitmate k, the parameter
; of the Wiener filter.
;
; CALLING SEQUENCE:
; k=oDataImage->EstimateK(Psf)
;
; INPUTS:
; Psf: A dataImage object with the PSF of the image
;
; KEYWORD PARAMETERS:
; MASK: A mask which covers pixels not to be use in the
; fits. Unlike the case of SBF power spectrum
; fitting, one does not necessarily want to cover
; all point sources.
; WAVELENGTHRANGE: as per dataImage::pspecFit()
; POWERRANGE: as per dataImage::pspecFit()
; FITDATA: as per dataImage::pspecFit()
; SBFCOUNTS: if this keyword is set, it the noise if fit
; to the histogram of negative pixel values,
; and the SBF noise removed to arrive at the
; final noise.
;
;
; OUTPUTS:
; This method returns an estimate of the value of "k" which is
; appropriate for use in a Wiener filter.
;
; RESTRICTIONS:
; This method assumes the image is dominated by randomly
; distributed point sources. If this is not the case, the
; returned value will not necessarily be the appropriate "k"
; value, though it is often still a good estimate in practice.
;
; PROCEDURE:
;
;
; EXAMPLE:
; To wiener filter an image with a psf in dataImage psf,
; k=exampleDataImage->estimateK(psf)
; exampleDataImage->wienerFilter,psf,k
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
; Sept. 30, 1997 EHN adds alternate estimation method triggered
; by setting the SBFCOUNTS keyword.
;-
FUNCTION dataImage::estimateK, psf, MASK=mask, $
WAVELENGTHRANGE=wavelengthrange, POWERRANGE=powerrange,$
FITDATA=fitData, SBFCOUNTS=SBFcounts
IF NOT KEYWORD_SET(SBFCOUNTS) THEN BEGIN
fit=self->pspecFit(psf,MASK=mask,$
WAVELENGTHRANGE=wavelengthrange, POWERRANGE=powerrange,$
FITDATA=fitData)
k=fit(0)/fit(1)
ENDIF ELSE BEGIN
noiseVariance=(self->estimateSigma(mask))^2.0
psfData=psf->getData()
psfSize=psf->size()
psfData=(ABS(FFT(psfData,-1,/OVERWRITE)))^2.0
psfMoment=MOMENT(psfData)
numberOfPixels=FLOAT(psfSize(0)*psfSize(1))
SBFvariance=SBFcounts*(numberOfPixels^2)*psfMoment(0)
trueNoiseVariance=noiseVariance-SBFvariance
dataArray=self->getData()
dataArray(WHERE(mask->getData() NE 0))=0
dataMoments=MOMENT(dataArray)/( mask->fraction() )
k=trueNoiseVariance/dataMoments(1)
ENDELSE
RETURN, k
END
;----------------------------------------------------------------------------
; dataImage::estimateSigma
;+
; NAME:
; dataImage::estimateSigma
;
; PURPOSE:
; This method estimates the sky sigma using the pixels with low
; pixel values. It only provides meaningful results where the
; sky and diffuse background objects have been subtracted.
;
; CALLING SEQUENCE:
; sigma=oDataImage->estimateSigma(Maskimage)
;
; OPTIONAL INPUTS:
; Maskimage: A maskImage object covering the pixels not to be
; considered.
;
; OUTPUTS:
; An estimate of the sky standard deviation is returned
;
; RESTRICTIONS:
; The image must be sky subtracted, and diffuse objects must
; either be subtracted or masked.
;
; PROCEDURE:
; The sky sigma is estimated to be the square root of the sum of
; the squares of the values of the unmasked, negative pixels.
;
; EXAMPLE:
; PRINT,exampleImage->estimateSigma()
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
FUNCTION dataImage::estimateSigma, maskImage
IF (N_PARAMS() GT 0) THEN BEGIN
lowPoints=WHERE((maskImage->getData() EQ 0) $
AND( self->getData() LT 0), npix)
ENDIF ELSE BEGIN
lowPoints=WHERE( (self->getData() LT 0), npix)
ENDELSE
sigma=SQRT(TOTAL(((self->getData())(lowPoints))^2)/npix)
RETURN, sigma
END
;----------------------------------------------------------------------------
; dataImage::estimateDisplayRange
;+
; NAME:
; dataImage::estimateDisplayRange
;
; PURPOSE:
; This method estimates a good display range for viewing the image.
;
; CALLING SEQUENCE:
; Result=oDataImage->estimateDisplayRange(Maskimage)
;
; OPTIONAL INPUTS:
; Maskimage: A maskImage mask covering bad pixels.
;
; KEYWORD PARAMETERS:
; STRETCH: This pararmeter is set to 16 by default. Larger
; values will result in a more uniform appearance in
; the image, with more extreme values unsaturated.
;
; OUTPUTS:
; An estimated range for display is returned.
;
; PROCEDURE:
;
;
; EXAMPLE:
; To display an image using a reasonable range:
; range=exampleImage->estimateDisplayRange(exampleMask)
; TV,exampleImage->getData(RANGE=range)
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
FUNCTION dataImage::estimateDisplayRange, maskImage, STRETCH=stretch
IF NOT KEYWORD_SET(stretch) THEN stretch=16.0
bufferData=self->getData()
IF (N_PARAMS() GT 0) THEN BEGIN
maskData=maskImage->getData()
bufferData=bufferData(WHERE(maskData EQ 0))
ENDIF
center=MEDIAN(bufferData)
upper=MEDIAN(bufferData(WHERE(bufferData GE center)))
lower=MEDIAN(bufferData(WHERE(bufferData LE center)))
displayRange=[MIN(bufferData),MAX(bufferData)]
displayRange(1)=(center+stretch*(upper-center) ) < MAX(bufferData)
displayRange(0)=(center+stretch*(lower-center)) > MIN(bufferData)
RETURN, DisplayRange
END
;----------------------------------------------------------------------------
; dataImage::tv
;
; Purpose: display the image
;
; Keywords:
; maskImage - the mask image to use
; stretch - the stretch of the range
PRO dataImage::tv, displayRange, INMASKIMAGE=inMaskImage, STRETCH=stretch,MAGSCALE=magscale
IF KEYWORD_SET(INMASKIMAGE) THEN BEGIN
maskImage=inMaskImage
ENDIF ELSE BEGIN
maskImage=OBJ_NEW('maskImage')
maskImage->setData,0*self->getData()
ENDELSE
IF N_PARAMS() LT 1 THEN displayRange=self->estimateDisplayrange(maskImage, STRETCH=stretch)
IF NOT KEYWORD_SET(INMASKIMAGE) THEN OBJ_DESTROY,maskImage
IF KEYWORD_SET(MAGSCALE) AND (displayRange(1) GT 0.0) THEN BEGIN
displayRange(0)= displayRange(0) > displayRange(1)*0.001
TV,BYTSCL(2.5*ALOG10(self->getData()),MIN=2.5*ALOG10(displayRange(0)),$
MAX=2.5*ALOG10(displayRange(1)),TOP=!D.TABLE_SIZE)
ENDIF ELSE BEGIN
TV,BYTSCL(self->getData(),MIN=displayRange(0),$
MAX=displayRange(1),TOP=!D.TABLE_SIZE)
ENDELSE
END
;----------------------------------------------------------------------------
; DATAIMAGE__DEFINE
;+
; NAME:
; DATAIMAGE__DEFINE
;
; PURPOSE:
; This routing defines the dataImage object. The dataImage
; object is intended to be the primary object used to store and
; manipulate astronomical image.
;
; CALLING SEQUENCE:
; oDataImage=OBJ_NEW('dataImage')
;
;
; MODIFICATION HISTORY:
; Written by: Eric H. Neilsen Jr., 1997.
;-
;
PRO DATAIMAGE__DEFINE
struct = {dataImage, $
exposureTime:1.0, $
numberOfExposures: 1, $
INHERITS astronomyImage $
}
END