Galaxy.py 19.6 KB
Newer Older
Fang Yuedong's avatar
Fang Yuedong committed
1
2
3
4
5
6
7
import numpy as np
import galsim
import os, sys
import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate

8
9
10
11
from ObservationSim.MockObject._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders
from ObservationSim.MockObject.SpecDisperser import SpecDisperser
from ObservationSim.MockObject.MockObject import MockObject

Fang Yuedong's avatar
Fang Yuedong committed
12
class Galaxy(MockObject):
13
14
    def __init__(self, param, rotation=None, logger=None):
        super().__init__(param, logger=logger)
Fang Yuedong's avatar
Fang Yuedong committed
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
        self.thetaR = self.param["theta"]
        self.bfrac = self.param["bfrac"]
        self.hlr_disk = self.param["hlr_disk"]
        self.hlr_bulge = self.param["hlr_bulge"]

        # Extract ellipticity components
        self.e_disk = galsim.Shear(g=self.param["ell_disk"], beta=self.thetaR*galsim.degrees)
        self.e_bulge = galsim.Shear(g=self.param["ell_bulge"], beta=self.thetaR*galsim.degrees)
        self.e_total = galsim.Shear(g=self.param["ell_tot"], beta=self.thetaR*galsim.degrees)
        self.e1_disk, self.e2_disk = self.e_disk.g1, self.e_disk.g2
        self.e1_bulge, self.e2_bulge = self.e_bulge.g1, self.e_bulge.g2
        self.e1_total, self.e2_total = self.e_total.g1, self.e_total.g2

        if rotation is not None:
            self.rotateEllipticity(rotation)

    def unload_SED(self):
        """(Test) free up SED memory
        """
        del self.sed

    def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.):
        if len(psf_list) != len(bandpass_list):
            raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.")
        objs = []
        if nphotons_tot == None:
            nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
        # print("nphotons_tot = ", nphotons_tot)

        try:
            full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
        except Exception as e:
            print(e)
48
            self.logger.error(e)
Fang Yuedong's avatar
Fang Yuedong committed
49
50
51
52
53
54
55
            return -1
        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]
            try:
                sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
            except Exception as e:
                print(e)
56
                self.logger.error(e)
Fang Yuedong's avatar
Fang Yuedong committed
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
                return -1
            
            ratio = sub/full
            if not (ratio == -1 or (ratio != ratio)):
                nphotons = ratio * nphotons_tot
            else:
                return -1

            psf = psf_list[i]
            disk = galsim.Sersic(n=1.0, half_light_radius=self.hlr_disk, flux=1.0)
            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
            disk = disk.shear(disk_shape)
            bulge = galsim.Sersic(n=4.0, half_light_radius=self.hlr_bulge, flux=1.0)
            bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
            bulge = bulge.shear(bulge_shape)

            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
            gal = gal.withFlux(nphotons)
            gal_shear = galsim.Shear(g1=g1, g2=g2)
            gal = gal.shear(gal_shear)
            gal = galsim.Convolve(psf, gal)
            objs.append(gal)
        final = galsim.Sum(objs)
        return final

    def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150.):
        if nphotons_tot == None:
            nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
        # print("nphotons_tot = ", nphotons_tot)

        try:
            full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
        except Exception as e:
            print(e)
91
            self.logger.error(e)
Fang Yuedong's avatar
Fang Yuedong committed
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
            return False

        nphotons_sum = 0
        photons_list = []
        xmax, ymax = 0, 0

        # print('hlr_disk = %.4f, hlr_bulge = %.4f'%(self.hlr_disk, self.hlr_bulge))
        big_galaxy = False
        if self.hlr_disk > 3.0: # Very big galaxy
            big_galaxy = True

        # (TEST) Galsim Parameters
        if self.getMagFilter(filt) <= 15 and (not big_galaxy):
            folding_threshold = 5.e-4
        else:
            folding_threshold = 5.e-3
        gsp = galsim.GSParams(folding_threshold=folding_threshold)

xin's avatar
xin committed
110
111
112
113
114
115
116
117
118
119
120
121
        self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
                                        img_real_wcs=self.real_wcs)

        x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
        x_nominal = int(np.floor(x + 0.5))
        y_nominal = int(np.floor(y + 0.5))
        dx = x - x_nominal
        dy = y - y_nominal
        offset = galsim.PositionD(dx, dy)

        real_wcs_local = self.real_wcs.local(self.real_pos)

Fang Yuedong's avatar
Fang Yuedong committed
122
123
124
125
126
127
128
        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]

            try:
                sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
            except Exception as e:
                print(e)
129
                self.logger.error(e)
Fang Yuedong's avatar
Fang Yuedong committed
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
                # return False
                continue
            
            ratio = sub/full
            if not (ratio == -1 or (ratio != ratio)):
                nphotons = ratio * nphotons_tot
            else:
                # return False
                continue
            nphotons_sum += nphotons
            # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))

            psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
            disk = galsim.Sersic(n=1.0, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
            disk = disk.shear(disk_shape)
            bulge = galsim.Sersic(n=4.0, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
            bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
            bulge = bulge.shear(bulge_shape)

150
            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
Fang Yuedong's avatar
Fang Yuedong committed
151

152
153
154
155
            # # (TEST) Random knots
            # knots = galsim.RandomKnots(npoints=100, profile=disk)
            # kfrac = np.random.random()*(1.0 - self.bfrac)
            # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
Fang Yuedong's avatar
Fang Yuedong committed
156

Fang Yuedong's avatar
Fang Yuedong committed
157
158
159
160
161
162
163
164
            gal = gal.withFlux(nphotons)
            gal_shear = galsim.Shear(g1=g1, g2=g2)
            gal = gal.shear(gal_shear)

            if self.hlr_disk < 10.0: # Not apply PSF for very big galaxy
                gal = galsim.Convolve(psf, gal)

            # Use (explicit) stamps to draw
xin's avatar
xin committed
165
166
167
168
169
170
171
172
173
            # stamp = gal.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True)
            # xmax = max(xmax, stamp.xmax)
            # ymax = max(ymax, stamp.ymax)
            # photons = stamp.photons
            # photons.x += self.x_nominal
            # photons.y += self.y_nominal
            # photons_list.append(photons)

            stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
Fang Yuedong's avatar
Fang Yuedong committed
174
175
176
            xmax = max(xmax, stamp.xmax)
            ymax = max(ymax, stamp.ymax)
            photons = stamp.photons
xin's avatar
xin committed
177
178
            photons.x += x_nominal
            photons.y += y_nominal
Fang Yuedong's avatar
Fang Yuedong committed
179
180
181
182
            photons_list.append(photons)

        # print('xmax = %d, ymax = %d '%(xmax, ymax))

xin's avatar
xin committed
183
184
185
186
187
        stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1))
        stamp.wcs = real_wcs_local
        stamp.setCenter(x_nominal, y_nominal)
        bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
        chip.img.setOrigin(0, 0)
Fang Yuedong's avatar
Fang Yuedong committed
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
        stamp[bounds] = chip.img[bounds]

        if not big_galaxy:
            for i in range(len(photons_list)):
                if i == 0:
                    chip.sensor.accumulate(photons_list[i], stamp)
                else:
                    chip.sensor.accumulate(photons_list[i], stamp, resume=True)
        else:
            sensor = galsim.Sensor()
            for i in range(len(photons_list)):
                if i == 0:
                    sensor.accumulate(photons_list[i], stamp)
                else:
                    sensor.accumulate(photons_list[i], stamp, resume=True)

        chip.img[bounds] = stamp[bounds]
xin's avatar
xin committed
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233

        chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)



        # stamp = galsim.ImageF(int(xmax*1.1), int(ymax*1.1))
        # stamp.wcs = self.localWCS
        # stamp.setCenter(self.x_nominal, self.y_nominal)
        # bounds = stamp.bounds & chip.img.bounds
        # stamp[bounds] = chip.img[bounds]
        #
        # if not big_galaxy:
        #     for i in range(len(photons_list)):
        #         if i == 0:
        #             chip.sensor.accumulate(photons_list[i], stamp)
        #         else:
        #             chip.sensor.accumulate(photons_list[i], stamp, resume=True)
        # else:
        #     sensor = galsim.Sensor()
        #     for i in range(len(photons_list)):
        #         if i == 0:
        #             sensor.accumulate(photons_list[i], stamp)
        #         else:
        #             sensor.accumulate(photons_list[i], stamp, resume=True)
        #
        # # print(stamp.array.sum())
        # # chip.img[bounds] += stamp[bounds]
        # chip.img[bounds] = stamp[bounds]
        # # print("nphotons_sum = ", nphotons_sum)
Fang Yuedong's avatar
Fang Yuedong committed
234
235
236
237
238
        del photons_list
        del stamp
        return True, pos_shear

    def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
Xin Zhang's avatar
Xin Zhang committed
239
                         exptime=150., normFilter=None, grating_split_pos=3685):
Fang Yuedong's avatar
Fang Yuedong committed
240
241
242
243
244
245
246
247
248
249
250

        norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
        sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
                                                      norm_thr=normFilter,
                                                      sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
                                                      eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
        if sedNormFactor == 0:
            return False
        normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
                          names=('WAVELENGTH', 'FLUX'))

xin's avatar
xin committed
251
252
253
254
255
256
257
258
259
260
261
262
263
        self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
                                        img_real_wcs=self.real_wcs)

        x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
        x_nominal = int(np.floor(x + 0.5))
        y_nominal = int(np.floor(y + 0.5))
        dx = x - x_nominal
        dy = y - y_nominal
        offset = galsim.PositionD(dx, dy)

        real_wcs_local = self.real_wcs.local(self.real_pos)


Fang Yuedong's avatar
Fang Yuedong committed
264
265
266
267
268
269
270
271
272
273
        big_galaxy = False
        if self.hlr_disk > 3.0: # Very big galaxy
            big_galaxy = True

        if self.getMagFilter(filt) <= 15 and (not big_galaxy):
            folding_threshold = 5.e-4
        else:
            folding_threshold = 5.e-3
        gsp = galsim.GSParams(folding_threshold=folding_threshold)
        # nphotons_sum = 0
274
275
276

        flat_cube = chip.flat_cube

Zhang Xin's avatar
Zhang Xin committed
277
        xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
xin's avatar
xin committed
278
        grating_split_pos_chip = 0 + grating_split_pos
Fang Yuedong's avatar
Fang Yuedong committed
279
280
281
282
283
284
285
286
287
288
289
        for i in range(len(bandpass_list)):
            bandpass = bandpass_list[i]

            psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
            disk = galsim.Sersic(n=1.0, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
            disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
            disk = disk.shear(disk_shape)
            bulge = galsim.Sersic(n=4.0, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
            bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
            bulge = bulge.shear(bulge_shape)

290
            gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
Fang Yuedong's avatar
Fang Yuedong committed
291
292

            # (TEST) Random knots
293
294
295
            # knots = galsim.RandomKnots(npoints=100, profile=disk)
            # kfrac = np.random.random()*(1.0 - self.bfrac)
            # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
Fang Yuedong's avatar
Fang Yuedong committed
296

Fang Yuedong's avatar
Fang Yuedong committed
297
298
299
300
301
            gal = gal.withFlux(tel.pupil_area * exptime)
            gal_shear = galsim.Shear(g1=g1, g2=g2)
            gal = gal.shear(gal_shear)
            gal = galsim.Convolve(psf, gal)

xin's avatar
xin committed
302
            starImg = gal.drawImage(wcs=real_wcs_local)
Fang Yuedong's avatar
Fang Yuedong committed
303

xin's avatar
xin committed
304
305
            origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
                           x_nominal - (starImg.center.x - starImg.xmin)]
Fang Yuedong's avatar
Fang Yuedong committed
306

Xin Zhang's avatar
Xin Zhang committed
307
308
309
310
311
312
313
314
315
316
            gal_origin = [origin_star[0], origin_star[1]]
            gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]

            if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
                subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
                ## part img disperse

                subImg_p1 = starImg.array[:, 0:subSlitPos]
                star_p1 = galsim.Image(subImg_p1)
                origin_p1 = origin_star
xin's avatar
xin committed
317
318
                xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0
                ycenter_p1 = y_nominal-0
Xin Zhang's avatar
Xin Zhang committed
319
320
321
322
323
324

                sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
                                    ycenter=ycenter_p1, origin=origin_p1,
                                    tar_spec=normalSED,
                                    band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
                                    conf=chip.sls_conf[0],
325
326
                                    isAlongY=0,
                                    flat_cube=flat_cube)
Xin Zhang's avatar
Xin Zhang committed
327

xin's avatar
xin committed
328
                self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
Xin Zhang's avatar
Xin Zhang committed
329
330
331
332

                subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
                star_p2 = galsim.Image(subImg_p2)
                origin_p2 = [origin_star[0], grating_split_pos_chip]
xin's avatar
xin committed
333
334
                xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
                ycenter_p2 = y_nominal - 0
Xin Zhang's avatar
Xin Zhang committed
335
336
337
338
339
340

                sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
                                       ycenter=ycenter_p2, origin=origin_p2,
                                       tar_spec=normalSED,
                                       band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
                                       conf=chip.sls_conf[1],
341
342
                                       isAlongY=0,
                                       flat_cube=flat_cube)
Xin Zhang's avatar
Xin Zhang committed
343

xin's avatar
xin committed
344
                self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
Xin Zhang's avatar
Xin Zhang committed
345
346
347
348

                del sdp_p1
                del sdp_p2
            elif grating_split_pos_chip<=gal_origin[1]:
xin's avatar
xin committed
349
350
                sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
                                    ycenter=y_nominal - 0, origin=origin_star,
Xin Zhang's avatar
Xin Zhang committed
351
352
353
                                    tar_spec=normalSED,
                                    band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
                                    conf=chip.sls_conf[1],
354
355
                                    isAlongY=0,
                                    flat_cube=flat_cube)
xin's avatar
xin committed
356
                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
Xin Zhang's avatar
Xin Zhang committed
357
358
                del sdp
            elif grating_split_pos_chip>=gal_end[1]:
xin's avatar
xin committed
359
360
                sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
                                    ycenter=y_nominal - 0, origin=origin_star,
Xin Zhang's avatar
Xin Zhang committed
361
362
363
                                    tar_spec=normalSED,
                                    band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
                                    conf=chip.sls_conf[0],
364
365
                                    isAlongY=0,
                                    flat_cube=flat_cube)
xin's avatar
xin committed
366
                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
Xin Zhang's avatar
Xin Zhang committed
367
                del sdp
Fang Yuedong's avatar
Fang Yuedong committed
368

Xin Zhang's avatar
Xin Zhang committed
369
            # print(self.y_nominal, starImg.center.y, starImg.ymin)
Fang Yuedong's avatar
Fang Yuedong committed
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
            del psf
        return True, pos_shear

    def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
        if flux == None:
            flux = self.getElectronFluxFilt(filt, tel, exptime)
        disk = galsim.Sersic(n=1.0, half_light_radius=self.hlr_disk, flux=1.0)
        disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
        disk = disk.shear(disk_shape)

        bulge = galsim.Sersic(n=4.0, half_light_radius=self.hlr_bulge, flux=1.0)
        bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
        bulge = bulge.shear(bulge_shape)

        gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
        gal = gal.withFlux(flux)
        gal_shear = galsim.Shear(g1=g1, g2=g2)
        gal = gal.shear(gal_shear)
        final = galsim.Convolve(psf, gal)
        return final

    def rotateEllipticity(self, rotation):
        if rotation == 1:
            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = -self.e2_disk, self.e1_disk, -self.e2_bulge, self.e1_bulge, -self.e2_total, self.e1_total
        if rotation == 2:
            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = -self.e1_disk, -self.e2_disk, -self.e1_bulge, -self.e2_bulge, -self.e1_total, -self.e2_total
        if rotation == 3:
            self.e1_disk, self.e2_disk, self.e1_bulge, self.e2_bulge, self.e1_total, self.e2_total = self.e2_disk, -self.e1_disk, self.e2_bulge, -self.e1_bulge, self.e2_total, -self.e1_total

    def drawObject(self, img, final, noise_level=0.0, flux=None, filt=None, tel=None, exptime=150.):
        """ Override the method in parent class 
        Need to constrain the size of image stamp for extended objects
        """
        isUpdated = True
        if flux == None:
            flux = self.getElectronFluxFilt(filt, tel, exptime)
        stamp = final.drawImage(wcs=self.localWCS, offset=self.offset)
        stamp_arr = stamp.array
        mask = (stamp_arr >= 0.001*noise_level) # why 0.001?
        err = int(np.sqrt(mask.sum()))
        if np.mod(err, 2) == 1:
            err += 1
        # if err == 1:
        if err == 0:
            subSize = 16 # why 16?
        else:
            subSize = max([err, 16])
            fluxRatio = flux / stamp_arr[mask].sum()
            final = final.withScaledFlux(fluxRatio)

        imgSub = galsim.ImageF(subSize, subSize)

        # Draw with FFT
        # stamp = final.drawImage(image=imgSub, wcs=self.localWCS, offset=self.offset)

        # Draw with Photon Shoot
        stamp = final.drawImage(image=imgSub, wcs=self.localWCS, method='phot', offset=self.offset)
        
        stamp.setCenter(self.x_nominal, self.y_nominal)
        if np.sum(np.isnan(stamp.array)) >= 1:
            stamp.setZero()
        bounds = stamp.bounds & img.bounds
        if bounds.area() == 0:
            isUpdated = False
        else:
            img[bounds] += stamp[bounds]
        return img, stamp, isUpdated


    def getObservedEll(self, g1=0, g2=0):
        e1_obs, e2_obs, e_obs, theta = eObs(self.e1_total, self.e2_total, g1, g2)
Zhang Xin's avatar
Zhang Xin committed
441
        return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs