1248 correlationStrength=0.1, maxLag=10, nSigmaClip=5, border=10):
1249 """Calculate the bias induced when sigma-clipping non-Gassian distributions. 1251 Fill image-pairs of the specified size with Poisson-distributed values, 1252 adding correlations as necessary. Then calculate the cross correlation, 1253 and calculate the bias induced using the cross-correlation image 1254 and the image means. 1258 fluxLevels : `list` of `int` 1259 The mean flux levels at which to simiulate. 1260 Nominal values might be something like [70000, 90000, 110000] 1261 imageShape : `tuple` of `int` 1262 The shape of the image array to simulate, nx by ny pixels. 1263 repeats : `int`, optional 1264 Number of repeats to perform so that results 1265 can be averaged to improve SNR. 1266 seed : `int`, optional 1267 The random seed to use for the Poisson points. 1268 addCorrelations : `bool`, optional 1269 Whether to add brighter-fatter-like correlations to the simulated images 1270 If true, a correlation between x_{i,j} and x_{i+1,j+1} is introduced 1271 by adding a*x_{i,j} to x_{i+1,j+1} 1272 correlationStrength : `float`, optional 1273 The strength of the correlations. 1274 This is the value of the coefficient `a` in the above definition. 1275 maxLag : `int`, optional 1276 The maximum lag to work to in pixels 1277 nSigmaClip : `float`, optional 1278 Number of sigma to clip to when calculating the sigma-clipped mean. 1279 border : `int`, optional 1280 Number of border pixels to mask 1284 biases : `dict` of `list` of `float` 1285 A dictionary, keyed by flux level, containing a list of the biases 1286 for each repeat at that flux level 1287 means : `dict` of `list` of `float` 1288 A dictionary, keyed by flux level, containing a list of the average mean 1289 fluxes (average of the mean of the two images) 1290 for the image pairs at that flux level 1291 xcorrs : `dict` of `list` of `np.ndarray` 1292 A dictionary, keyed by flux level, containing a list of the xcorr 1293 images for the image pairs at that flux level 1295 means = {f: []
for f
in fluxLevels}
1296 xcorrs = {f: []
for f
in fluxLevels}
1297 biases = {f: []
for f
in fluxLevels}
1299 config = MakeBrighterFatterKernelTaskConfig()
1300 config.isrMandatorySteps = []
1301 config.isrForbiddenSteps = []
1302 config.nSigmaClipXCorr = nSigmaClip
1303 config.nPixBorderXCorr = border
1304 config.maxLag = maxLag
1305 task = MakeBrighterFatterKernelTask(config=config)
1307 im0 = afwImage.maskedImage.MaskedImageF(imageShape[1], imageShape[0])
1308 im1 = afwImage.maskedImage.MaskedImageF(imageShape[1], imageShape[0])
1310 random = np.random.RandomState(seed)
1312 for rep
in range(repeats):
1313 for flux
in fluxLevels:
1314 data0 = random.poisson(flux, (imageShape)).astype(float)
1315 data1 = random.poisson(flux, (imageShape)).astype(float)
1317 data0[1:, 1:] += correlationStrength*data0[: -1, : -1]
1318 data1[1:, 1:] += correlationStrength*data1[: -1, : -1]
1319 im0.image.array[:, :] = data0
1320 im1.image.array[:, :] = data1
1322 _xcorr, _means = task._crossCorrelate(im0, im1, runningBiasCorrSim=
True)
1324 means[flux].
append(_means)
1325 xcorrs[flux].
append(_xcorr)
1327 bias = xcorrs[flux][-1][1, 1]/means[flux][-1]*(1 + correlationStrength)/correlationStrength
1328 print(
"Simulated/expected avg. flux: %.1f, %.1f" % (flux, means[flux][-1]/2))
1329 print(
"Bias: %.6f" % bias)
1331 bias = xcorrs[flux][-1][0, 0]/means[flux][-1]
1332 print(
"Simulated/expected avg. flux: %.1f, %.1f" % (flux, means[flux][-1]/2))
1333 print(
"Bias: %.6f" % bias)
1334 biases[flux].
append(bias)
1336 return biases, means, xcorrs
1337 std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.