FieldDistortion.py 4.92 KB
Newer Older
Fang Yuedong's avatar
Fang Yuedong committed
1
2
import galsim
import numpy as np
3
import cmath
Fang Yuedong's avatar
Fang Yuedong committed
4
5
6

class FieldDistortion(object):

7
    def __init__(self, chip, fdModel=None, fdModel_path=None, img_rot=0.):
Fang Yuedong's avatar
Fang Yuedong committed
8
        if fdModel is None:
Fang Yuedong's avatar
Fang Yuedong committed
9
10
11
12
            if hasattr(chip, 'fdModel'):
                self.fdModel = chip.fdModel
            elif fdModel_path is not None:
                import pickle
Fang Yuedong's avatar
Fang Yuedong committed
13
14
15
                with open(fdModel_path, "rb") as f:
                    self.fdModel = pickle.load(f)
            else:
Fang Yuedong's avatar
Fang Yuedong committed
16
                raise ValueError("Error: no field distortion model has been specified!")
Fang Yuedong's avatar
Fang Yuedong committed
17
18
        else:
            self.fdModel = fdModel
19
        self.img_rot = img_rot
Fang Yuedong's avatar
Fang Yuedong committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
        self.ifdModel = self.fdModel["wave1"]
        self.ixfdModel = self.ifdModel["xImagePos"]
        self.iyfdModel = self.ifdModel["yImagePos"]
        # first-order derivatives of the global field distortion model
        self.ifx_dx    = self.ixfdModel.partial_derivative(1,0)
        self.ifx_dy    = self.ixfdModel.partial_derivative(0,1)
        self.ify_dx    = self.iyfdModel.partial_derivative(1,0)
        self.ify_dy    = self.iyfdModel.partial_derivative(0,1)
        if "residual" in self.fdModel["wave1"]:
            self.irsModel = self.fdModel["wave1"]["residual"]["ccd" + chip.getChipLabel(chipID=chip.chipID)]
            self.ixrsModel = self.irsModel["xResidual"]
            self.iyrsModel = self.irsModel["yResidual"]
            # first-order derivatives of the residual field distortion model
            self.irx_dx    = self.ixrsModel.partial_derivative(1,0)
            self.irx_dy    = self.ixrsModel.partial_derivative(0,1)
            self.iry_dx    = self.iyrsModel.partial_derivative(1,0)
            self.iry_dy    = self.iyrsModel.partial_derivative(0,1)
        else:
            self.irsModel = None
Fang Yuedong's avatar
Fang Yuedong committed
39
40

    def isContainObj_FD(self, chip, pos_img):
Fang Yuedong's avatar
Fang Yuedong committed
41
        xLowI, xUpI, yLowI, yUpI = self.ifdModel["interpLimit"]
Fang Yuedong's avatar
Fang Yuedong committed
42
        x, y = pos_img.x, pos_img.y
Fang Yuedong's avatar
Fang Yuedong committed
43
        if (xLowI - x) * (xUpI - x) > 0 or (yLowI - y) * (yUpI - y) > 0:
Fang Yuedong's avatar
Fang Yuedong committed
44
45
46
            return False
        return True

47
    def get_distorted(self, chip, pos_img, bandpass=None, img_rot=None):
Fang Yuedong's avatar
Fang Yuedong committed
48
49
50
51
52
53
54
55
56
57
58
59
60
61
        """ Get the distored position for an undistorted image position

        Parameters:
            chip:       A 'Chip' object representing the
                        chip we want to extract PSF from.
            pos_img:    A 'galsim.Position' object representing
                        the image position.
            bandpass:   A 'galsim.Bandpass' object representing
                        the wavelength range.
        Returns:
            pos_distorted:      A 'galsim.Position' object representing
                                the distored position.
        """
        if not self.isContainObj_FD(chip=chip, pos_img=pos_img):
Fang Yuedong's avatar
Fang Yuedong committed
62
            return galsim.PositionD(-1, -1), None
63
64
65
66
        if not img_rot:
            img_rot = np.radians(self.img_rot)
        else:
            img_rot = np.radians(img_rot)
Fang Yuedong's avatar
Fang Yuedong committed
67
        x, y = pos_img.x, pos_img.y
Fang Yuedong's avatar
Fang Yuedong committed
68
69
        x = self.ixfdModel(x, y)[0][0]
        y = self.iyfdModel(x, y)[0][0]
70
        if self.irsModel:
Fang Yuedong's avatar
Fang Yuedong committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
            # x1LowI, x1UpI, y1LowI, y1UpI = self.irsModel["interpLimit"]
            # if (x1LowI-x)*(x1UpI-x) <=0 and (y1LowI-y)*(y1UpI-y)<=0:
            #     dx = self.ixrsModel(x, y)[0][0]
            #     dy = self.iyrsModel(x, y)[0][0]
            #     x += dx
            #     y += dy
            #     # field distortion induced ellipticity
            #     ix_dx = self.ifx_dx(x, y) + self.irx_dx(x, y)
            #     ix_dy = self.ifx_dy(x, y) + self.irx_dy(x, y)
            #     iy_dx = self.ify_dx(x, y) + self.iry_dx(x, y)
            #     iy_dy = self.ify_dy(x, y) + self.iry_dy(x, y)
            # else:
            #     return galsim.PositionD(-1, -1), None
            dx = self.ixrsModel(x, y)[0][0]
            dy = self.iyrsModel(x, y)[0][0]
86
87
            x += dx
            y += dy
Fang Yuedong's avatar
Fang Yuedong committed
88
89
90
91
92
            # field distortion induced ellipticity
            ix_dx = self.ifx_dx(x, y) + self.irx_dx(x, y)
            ix_dy = self.ifx_dy(x, y) + self.irx_dy(x, y)
            iy_dx = self.ify_dx(x, y) + self.iry_dx(x, y)
            iy_dy = self.ify_dy(x, y) + self.iry_dy(x, y)
93
94
95
96
97
        else:
            ix_dx = self.ifx_dx(x, y)
            ix_dy = self.ifx_dy(x, y)
            iy_dx = self.ify_dx(x, y)
            iy_dy = self.ify_dy(x, y)
Fang Yuedong's avatar
Fang Yuedong committed
98
99
        g1k_fd   = 0.0 + (iy_dy - ix_dx) / (iy_dy + ix_dx)
        g2k_fd   = 0.0 - (iy_dx + ix_dy) / (iy_dy + ix_dx)
100
101
102
103
104
105
106
107
108
109

        # [TODO] [TESTING] Rotate the shear:
        g_abs = np.sqrt(g1k_fd**2 + g2k_fd**2)
        phi = cmath.phase(complex(g1k_fd, g2k_fd))
        # g_abs = 0.7
        g1k_fd = g_abs * np.cos(phi - 2*img_rot)
        g2k_fd = g_abs * np.sin(phi - 2*img_rot)
        # g1k_fd = g_abs * np.cos(phi - 2*img_rot)
        # g2k_fd = -g_abs * np.sin(phi - 2*img_rot)
        # g1k_fd = g_abs * np.cos(0. - 2*img_rot)
110
        # g2k_fd = g_abs * np.sin(0. - 2*img_rot)
111

Fang Yuedong's avatar
Fang Yuedong committed
112
113
114
        fd_shear = galsim.Shear(g1=g1k_fd, g2=g2k_fd)
        return galsim.PositionD(x, y), fd_shear