Python instrument.CameraGeometry类代码示例

def test_equals():
    cam1 = CameraGeometry.from_name("LSTCam")
    cam2 = CameraGeometry.from_name("LSTCam")
    cam3 = CameraGeometry.from_name("ASTRICam")

    assert cam1 is not cam2
    assert cam1 == cam2
    assert cam1 != cam3

def test_neighbor_pixels():
    hexgeom = CameraGeometry.from_name("LSTCam")
    recgeom = CameraGeometry.make_rectangular()

    # most pixels should have 4 neighbors for rectangular geometry and 6 for
    # hexagonal
    assert int(median(recgeom.neighbor_matrix.sum(axis=1))) == 4
    assert int(median(hexgeom.neighbor_matrix.sum(axis=1))) == 6

def test_known_camera_names():
    cams = CameraGeometry.get_known_camera_names()
    assert len(cams) > 4
    assert 'FlashCam' in cams
    assert 'NectarCam' in cams

    for cam in cams:
        geom = CameraGeometry.from_name(cam)
        geom.info()

def test_calc_pixel_neighbors_square_diagonal():
    x, y = np.meshgrid(np.arange(20), np.arange(20))

    cam = CameraGeometry(
        cam_id='test',
        pix_id=np.arange(400),
        pix_type='rectangular',
        pix_x=u.Quantity(x.ravel(), u.cm),
        pix_y=u.Quantity(y.ravel(), u.cm),
        pix_area=u.Quantity(np.ones(400), u.cm**2),
    )

    cam._neighbors = cam.calc_pixel_neighbors(diagonal=True)
    assert set(cam.neighbors[21]) == {0, 1, 2, 20, 22, 40, 41, 42}

def test_border_pixels():
    from ctapipe.instrument.camera import CameraGeometry

    cam = CameraGeometry.from_name("LSTCam")

    assert np.sum(cam.get_border_pixel_mask(1)) == 168
    assert np.sum(cam.get_border_pixel_mask(2)) == 330

    cam = CameraGeometry.from_name("ASTRICam")
    assert np.sum(cam.get_border_pixel_mask(1)) == 212
    assert np.sum(cam.get_border_pixel_mask(2)) == 408

    assert cam.get_border_pixel_mask(1)[0]
    assert cam.get_border_pixel_mask(1)[2351]
    assert not cam.get_border_pixel_mask(1)[521]

def camera_waveforms():
    camera = CameraGeometry.from_name("CHEC")

    n_pixels = camera.n_pixels
    n_samples = 96
    mid = n_samples // 2
    pulse_sigma = 6
    r_hi = np.random.RandomState(1)
    r_lo = np.random.RandomState(2)

    x = np.arange(n_samples)

    # Randomize times
    t_pulse_hi = r_hi.uniform(mid - 10, mid + 10, n_pixels)[:, np.newaxis]
    t_pulse_lo = r_lo.uniform(mid + 10, mid + 20, n_pixels)[:, np.newaxis]

    # Create pulses
    y_hi = norm.pdf(x, t_pulse_hi, pulse_sigma)
    y_lo = norm.pdf(x, t_pulse_lo, pulse_sigma)

    # Randomize amplitudes
    y_hi *= r_hi.uniform(100, 1000, n_pixels)[:, np.newaxis]
    y_lo *= r_lo.uniform(100, 1000, n_pixels)[:, np.newaxis]

    y = np.stack([y_hi, y_lo])

    return y, camera

def display_event(event):
    """an extremely inefficient display. It creates new instances of
    CameraDisplay for every event and every camera, and also new axes
    for each event. It's hacked, but it works
    """
    print("Displaying... please wait (this is an inefficient implementation)")
    global fig
    ntels = len(event.r0.tels_with_data)
    fig.clear()

    plt.suptitle("EVENT {}".format(event.r0.event_id))

    disps = []

    for ii, tel_id in enumerate(event.r0.tels_with_data):
        print("\t draw cam {}...".format(tel_id))
        nn = int(ceil(sqrt(ntels)))
        ax = plt.subplot(nn, nn, ii + 1)

        x, y = event.inst.pixel_pos[tel_id]
        geom = CameraGeometry.guess(x, y, event.inst.optical_foclen[tel_id])
        disp = CameraDisplay(geom, ax=ax, title="CT{0}".format(tel_id))
        disp.pixels.set_antialiaseds(False)
        disp.autoupdate = False
        disp.cmap = random.choice(cmaps)
        chan = 0
        signals = event.r0.tel[tel_id].adc_sums[chan].astype(float)
        signals -= signals.mean()
        disp.image = signals
        disp.set_limits_percent(95)
        disp.add_colorbar()
        disps.append(disp)

    return disps

def test_intensity():
    from ctapipe.image.toymodel import Gaussian

    np.random.seed(0)
    geom = CameraGeometry.from_name('LSTCam')

    x, y = u.Quantity([0.2, 0.3], u.m)
    width = 0.05 * u.m
    length = 0.15 * u.m
    intensity = 50
    psi = '30d'

    # make a toymodel shower model
    model = Gaussian(x=x, y=y, width=width, length=length, psi=psi)

    image, signal, noise = model.generate_image(
        geom, intensity=intensity, nsb_level_pe=5,
    )

    # test if signal reproduces given cog values
    assert np.average(geom.pix_x.to_value(u.m), weights=signal) == approx(0.2, rel=0.15)
    assert np.average(geom.pix_y.to_value(u.m), weights=signal) == approx(0.3, rel=0.15)

    # test if signal reproduces given width/length values
    cov = np.cov(geom.pix_x.value, geom.pix_y.value, aweights=signal)
    eigvals, eigvecs = np.linalg.eigh(cov)

    assert np.sqrt(eigvals[0]) == approx(width.to_value(u.m), rel=0.15)
    assert np.sqrt(eigvals[1]) == approx(length.to_value(u.m), rel=0.15)

    # test if total intensity is inside in 99 percent confidence interval
    assert poisson(intensity).ppf(0.05) <= signal.sum() <= poisson(intensity).ppf(0.95)

def test_skewed():
    from ctapipe.image.toymodel import SkewedGaussian
    # test if the parameters we calculated for the skew normal
    # distribution produce the correct moments
    np.random.seed(0)
    geom = CameraGeometry.from_name('LSTCam')

    x, y = u.Quantity([0.2, 0.3], u.m)
    width = 0.05 * u.m
    length = 0.15 * u.m
    intensity = 50
    psi = '30d'
    skewness = 0.3

    model = SkewedGaussian(
        x=x, y=y, width=width,
        length=length, psi=psi, skewness=skewness
    )
    image, signal, _ = model.generate_image(
        geom, intensity=intensity, nsb_level_pe=5,
    )

    a, loc, scale = model._moments_to_parameters()
    mean, var, skew = skewnorm(a=a, loc=loc, scale=scale).stats(moments='mvs')

    assert np.isclose(mean, 0)
    assert np.isclose(var, length.to_value(u.m)**2)
    assert np.isclose(skew, skewness)

def test_convert_geometry(cam_id, rot):

    geom = CameraGeometry.from_name(cam_id)

    if geom.pix_type=='rectangular':
        return  # skip non-hexagonal cameras, since they don't need conversion

    model = generate_2d_shower_model(centroid=(0.4, 0), width=0.01, length=0.03,
                                     psi="25d")

    _,image,_ = make_toymodel_shower_image(geom, model.pdf,
                                           intensity=50,
                                           nsb_level_pe=100)

    hillas_0 = hillas_parameters(geom.pix_x, geom.pix_y, image)

    geom2d, image2d = convert_geometry_1d_to_2d(geom, image,
                                                 geom.cam_id+str(rot),
                                                 add_rot=-2)
    geom1d, image1d = convert_geometry_back(geom2d, image2d,
                                            geom.cam_id+str(rot),
                                            add_rot=rot)
    hillas_1 = hillas_parameters(geom1d.pix_x, geom1d.pix_y, image1d)

    if __name__ == "__main__":
        plot_cam(geom, geom2d, geom1d, image, image2d, image1d)

    assert np.abs(hillas_1.phi - hillas_0.phi).deg < 1.0

def create_sample_image(
        psi='-30d',
        x=0.2 * u.m,
        y=0.3 * u.m,
        width=0.05 * u.m,
        length=0.15 * u.m,
        intensity=1500
):
    seed(10)

    geom = CameraGeometry.from_name('LSTCam')

    # make a toymodel shower model
    model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi)

    # generate toymodel image in camera for this shower model.
    image, _, _ = model.generate_image(
        geom,
        intensity=1500,
        nsb_level_pe=3,
    )

    # calculate pixels likely containing signal
    clean_mask = tailcuts_clean(geom, image, 10, 5)

    return geom, image, clean_mask

def test_convert_geometry_mock(cam_id, rot):
    """here we use a different key for the back conversion to trigger the mock conversion
    """

    geom = CameraGeometry.from_name(cam_id)
    image = create_mock_image(geom)
    hillas_0 = hillas_parameters(geom, image)

    if geom.pix_type == 'hexagonal':
        convert_geometry_1d_to_2d = convert_geometry_hex1d_to_rect2d
        convert_geometry_back = convert_geometry_rect2d_back_to_hexe1d

        geom2d, image2d = convert_geometry_1d_to_2d(geom, image, key=None,
                                                    add_rot=rot)
        geom1d, image1d = convert_geometry_back(geom2d, image2d,
                                                "_".join([geom.cam_id,
                                                          str(rot), "mock"]),
                                                add_rot=rot)
    else:
        # originally rectangular geometries don't need a buffer and therefore no mock
        # conversion
        return

    hillas_1 = hillas_parameters(geom, image1d)
    assert np.abs(hillas_1.phi - hillas_0.phi).deg < 1.0

def create_sample_image(psi='-30d'):

    seed(10)

    # set up the sample image using a HESS camera geometry (since it's easy
    # to load)
    geom = CameraGeometry.from_name("LSTCam")

    # make a toymodel shower model
    model = toymodel.generate_2d_shower_model(centroid=(0.2, 0.3),
                                              width=0.001, length=0.01,
                                              psi=psi)

    # generate toymodel image in camera for this shower model.
    image, signal, noise = toymodel.make_toymodel_shower_image(geom, model.pdf,
                                                               intensity=50,
                                                               nsb_level_pe=100)

    # denoise the image, so we can calculate hillas params
    clean_mask = tailcuts_clean(geom, image, 10,
                                5)  # pedvars = 1 and core and boundary
    # threshold in pe
    image[~clean_mask] = 0

    # Pixel values in the camera
    pix_x = geom.pix_x.value
    pix_y = geom.pix_y.value

    return pix_x, pix_y, image

def test_neighbor_pixels(cam_id):
    """
    test if each camera has a reasonable number of neighbor pixels (4 for
    rectangular, and 6 for hexagonal.  Other than edge pixels, the majority
    should have the same value
    """

    geom = CameraGeometry.from_name(cam_id)
    n_pix = len(geom.pix_id)
    n_neighbors = [len(x) for x in geom.neighbors]


    if geom.pix_type.startswith('hex'):
        assert n_neighbors.count(6) > 0.5 * n_pix
        assert n_neighbors.count(6) > n_neighbors.count(4)

    if geom.pix_type.startswith('rect'):
        assert n_neighbors.count(4) > 0.5 * n_pix
        assert n_neighbors.count(5) == 0
        assert n_neighbors.count(6) == 0

    # whipple has inhomogenious pixels that mess with pixel neighborhood
    # calculation
    if cam_id != 'Whipple490':
        assert np.all(geom.neighbor_matrix == geom.neighbor_matrix.T)
        assert n_neighbors.count(1) == 0  # no pixel should have a single neighbor

def test_slicing_rotation(cam_id):
    cam = CameraGeometry.from_name(cam_id)
    cam.rotate('25d')

    sliced1 = cam[5:10]

    assert sliced1.pix_x[0] == cam.pix_x[5]