Coverage for larch/xray/background.py: 11%
88 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1"""
2Methods for fitting background in xray spectra (energy dispersive or diffraction)
4Authors/Modifications:
5----------------------
6* Mark Rivers, GSECARS
7* modified for Larch, M Newville
8* modified to work for XRD data, M. Koker
10Notes:
11------
12This function fits a background to MCA spectrum or XRD data. The background is
13fitted using an enhanced version of the algorithm published by
14Kajfosz, J. and Kwiatek, W .M. (1987) 'Non-polynomial approximation of
15background in x-ray spectra.' Nucl. Instrum. Methods B22, 78-81.
18Procedure:
201) A series of concave down polynomials are constructed. At each
21channel "i" an n'th degree polynomial which is concave up is
22fitted. The polynomial is slid up from below until it "just touches"
23some channel of the spectrum. Call this channel "i". The maximum
24height of the polynomial is thus
27 (e(i) - e(j))**n
28 height(j) = max ( b(j) + -------------- )
29 i width**n
32where width is a user_specified parameter.
343) Once the value of height(i) is known the polynomial is fitted. The
35background counts in each channel are then determined from:
38 (e(i) - e(j))**n
39 bgr(j) = max ( height(i) + -------------- )
40 i width**n
43bgr(j) is thus the maximum counts for any of the concave down
44polynomials passing though channel j.
46Before the concave-down polynomials are fitted the spectrum at each
47channel it is possible to subtract out a straight line which is
48tangent to the spectrum at that channel. Use the TANGENT qualifier to
49do this. This is equivalent to fitting a "tilted" polynomial whose
50apex is tangent to the spectrum at that channel. By fitting
51polynomials which are tangent rather than vertical the background fit
52is much improved on spectra with steep slopes.
54Input Parameter Fields:
56 width (double, variable):
57 Specifies the width of the polynomials which are concave downward.
58 The bottom_width is the full width in energy units at which the
59 magnitude of the polynomial is 100 counts. The default is 4.
61 exponent (int):
62 Specifies the power of polynomial which is used. The power must be
63 an integer. The default is 2, i.e. parabolas. Higher exponents,
64 for example EXPONENT=4, results in polynomials with flatter tops
65 and steeper sides, which can better fit spectra with steeply
66 sloping backgrounds.
68 compress (int):
69 Compression factor to apply before fitting the background.
70 Default=4, which means, for example, that a 2048 channel spectrum
71 will be rebinned to 512 channels before fitting.
72 The compression is done on a temporary copy of the input spectrum,
73 so the input spectrum itself is unchanged.
74 The algorithm works best if the spectrum is compressed before it
75 is fitted. There are two reasons for this. First, the background
76 is constrained to never be larger than the data itself. If the
77 spectrum has negative noise spikes they will cause the fit to be
78 too low. Compression will smooth out such noise spikes.
79 Second, the algorithm requires about 3*N**2 operations, so the time
80 required grows rapidly with the size of the input spectrum. On a
81 200 MHz Pentium it takes about 3 seconds to fit a 2048 channel
82 spectrum with COMPRESS=1 (no compression), but only 0.2 seconds
83 with COMPRESS=4 (the default).
85 Note - compress needs a data array that integer divisible.
87 tangent (True/False):
88 Specifies that the polynomials are to be tangent to the slope of the
89 spectrum. The default is vertical polynomials. This option works
90 best on steeply sloping spectra. It has trouble in spectra with
91 big peaks because the polynomials are very tilted up inside the
92 peaks. Default is False
94Inputs to calc()
95 data:
96 The raw data to fit the background
98 slope:
99 Slope for the conversion from channel number to energy.
100 i.e. the slope from calibration
102"""
105import numpy as np
107REFERENCE_AMPL=100.
108TINY = 1.E-20
109HUGE = 1.E20
110MAX_TANGENT=2
112def compress_array(array, compress):
113 """
114 Compresses an 1-D array by the integer factor compress.
115 near equivalent of IDL's 'rebin'....
116 """
118 if len(array) % compress != 0:
119 ## Trims array to be divisible by compress factor
120 rng_min = int( (len(array) % compress ) / 2)
121 rng_max = int( len(array) / compress ) * compress + 1
122 array = array[rng_min:rng_max]
124 nsize = int(len(array)/compress)
125 temp = np.resize(array, (nsize, compress))
126 return np.sum(temp, 1)/compress
129def expand_array(array, expand, sample=0):
130 """
131 Expands an 1-D array by the integer factor expand.
133 if 'sample' is 1 the new array is created with sampling,
134 if 0 then the new array is created via interpolation (default)
135 Temporary fix until the equivalent of IDL's 'rebin' is found.
136 """
137 if expand == 1:
138 return array
139 if sample == 1:
140 return np.repeat(array, expand)
142 kernel = np.ones(expand)/expand
143 # The following mimic the behavior of IDL's rebin when expanding
144 temp = np.convolve(np.repeat(array, expand), kernel, mode=2)
145 # Discard the first "expand-1" entries
146 temp = temp[expand-1:]
147 # Replace the last "expand" entries with the last entry of original
148 for i in range(1,expand):
149 temp[-i]=array[-1]
150 return temp
152class XrayBackground:
153 '''
154 Class defining a spectrum background
156 Attributes:
157 -----------
158 These may be set by kw argument upon initialization.
161 * width = 4.0 # Width
162 * exponent = 2 # Exponent
163 * compress = 2 # Compress
164 * tangent = False # Tangent flag
165 '''
167 def __init__(self, data=None, width=4, slope=1.0, exponent=2, compress=2,
168 tangent=False, type_int=False, data_type='xrf'):
170 if data_type == 'xrf': type_int = True
172 self.bgr = []
173 self.width = width
174 self.compress = compress
175 self.exponent = exponent
176 self.tangent = tangent
178 self.info = {'width': width, 'compress': compress,
179 'exponent': exponent, 'tangent': tangent}
181 self.data = data
182 if data is not None:
183 self.calc(data, slope=slope, type_int=type_int)
185 def calc(self, data=None, slope=1.0, type_int=False):
186 '''compute background
188 Parameters:
189 -----------
190 * data is the spectrum
191 * slope is the slope of conversion channels to energy
192 '''
194 if data is None:
195 data = self.data
197 width = self.width
198 exponent = self.exponent
199 tangent = self.tangent
200 compress = self.compress
202 nchans = len(data)
203 self.bgr = np.zeros(nchans, dtype=np.int32)
204 scratch = 1.0*data[:]
206 # Compress scratch spectrum
207 if compress > 1:
208 tmp = compress_array(scratch, compress)
209 if tmp is None:
210 compress = 1
211 else:
212 scratch = tmp
213 slope = slope * compress
214 nchans = len(scratch) #nchans / compress
216 # Copy scratch spectrum to background spectrum
217 bckgnd = 1.0*scratch[:]
219 # Find maximum counts in input spectrum. This information is used to
220 # limit the size of the function lookup table
221 max_counts = max(scratch)
223 # Fit functions which come up from below
224 bckgnd = np.arange(nchans, dtype=np.float64) - HUGE
226 denom = max(TINY, (width / (2. * slope)**exponent))
228 indices = np.arange(nchans*2+1, dtype=np.float64) - nchans
229 power_funct = indices**exponent * (REFERENCE_AMPL / denom)
230 power_funct = np.compress((power_funct <= max_counts), power_funct)
231 max_index = int(len(power_funct)/2 - 1)
232 for chan in range(nchans-1):
233 tan_slope = 0.
234 if tangent:
235 # Find slope of tangent to spectrum at this channel
236 chan0 = max((chan - MAX_TANGENT), 0)
237 chan1 = min((chan + MAX_TANGENT), (nchans-1))
238 denom = chan - np.arange(chan1 - chan0 + 1, dtype=np.float64)
239 # is this correct?
240 denom = max(max(denom), 1)
241 tan_slope = (scratch[chan] - scratch[chan0:chan1+1]) / denom
242 tan_slope = np.sum(tan_slope) / (chan1 - chan0)
244 chan0 = int(max((chan - max_index), 0))
245 chan1 = int(min((chan + max_index), (nchans-1)))
246 chan1 = max(chan1, chan0)
247 nc = chan1 - chan0 + 1
248 lin_offset = scratch[chan] + (np.arange(float(nc)) - nc/2) * tan_slope
250 # Find the maximum height of a function centered on this channel
251 # such that it is never higher than the counts in any channel
252 f = int(chan0 - chan + max_index)
253 l = int(chan1 - chan + max_index)
254 test = scratch[chan0:chan1+1] - lin_offset + power_funct[f:l+1]
255 height = min(test)
257 # We now have the function height. Set the background to the
258 # height of the maximum function amplitude at each channel
260 test = height + lin_offset - power_funct[f:l+1]
261 sub = bckgnd[chan0:chan1+1]
262 bckgnd[chan0:chan1+1] = np.maximum(sub, test)
264 # Expand spectrum
265 if compress > 1:
266 bckgnd = expand_array(bckgnd, compress)
268 ## Set background to be of type integer
269 if type_int:
270 bckgnd = bckgnd.astype(int)
272 ## No negative values in background
273 bckgnd[np.where(bckgnd <= 0)] = 0
275 self.bgr = bckgnd