Commit 35600e7a authored by yuedong0607's avatar yuedong0607
Browse files

remove hard-coded focal plane definition. add channel and layout parameters to apply_crosstalk

parent f2e7f6de
Loading
Loading
Loading
Loading
+39 −20
Original line number Diff line number Diff line
@@ -239,27 +239,40 @@ class Observation(object):
                overall_config=self.config, copy_obs_config=True
            )

            self.focal_plane = FocalPlane(chip_list=pointing.obs_param["run_chips"])
            survey_type = pointing.obs_param.get(
                "survey_type",
                self.config.get("obs_setting", {}).get(
                    "survey_type", "Photometric"
                ),
            )
            self.focal_plane = FocalPlane(
                chip_list=pointing.obs_param["run_chips"],
                survey_type=survey_type,
                bad_chips=pointing.obs_param.get("bad_chips"),
                instrument_repository=self.instrument_repository,
            )
            # Make Chip & Filter lists
            self.chip_list = []
            self.filter_list = []
            self.all_filters = []
            for i in range(self.focal_plane.nchips):
                chipID = i + 1
                chip = Chip(
                    chipID=chipID,
                    config=self.config,
                    instrument_repository=self.instrument_repository,
                )
                filter_id, filter_type = chip.getChipFilter()
            selected_chip_ids = set(self.focal_plane.selected_chip_ids)
            for detector in self.focal_plane.iter_detectors():
                chip_id = detector.id
                filter_type = detector.instance.optical_role.filter
                filter_id = self.instrument_repository.filter_id(filter_type)
                filt = Filter(
                    filter_id=filter_id,
                    filter_type=filter_type,
                    filter_param=self.filter_param,
                    detector_id=chip.chipID,
                    detector_id=chip_id,
                    instrument_repository=self.instrument_repository,
                )
                if chip_id in selected_chip_ids:
                    chip = Chip(
                        chipID=chip_id,
                        config=self.config,
                        instrument_repository=self.instrument_repository,
                    )
                if not self.focal_plane.isIgnored(chipID=chipID):
                    self.chip_list.append(chip)
                    self.filter_list.append(filt)
                self.all_filters.append(filt)
@@ -271,15 +284,21 @@ class Observation(object):
                nchips_per_fp = len(self.chip_list)
            else:
                # Only run a particular set of chips
                run_chips = []
                run_filts = []
                for ichip in range(len(self.chip_list)):
                    chip = self.chip_list[ichip]
                    filt = self.filter_list[ichip]
                    if chip.chipID in chips:
                        run_chips.append(chip)
                        run_filts.append(filt)
                nchips_per_fp = len(chips)
                requested_ids = set(self.focal_plane.canonicalize_chip_ids(chips))
                unavailable_ids = requested_ids - selected_chip_ids
                if unavailable_ids:
                    raise ValueError(
                        "requested chips are not enabled for this pointing: "
                        f"{sorted(unavailable_ids)}"
                    )
                selected_pairs = [
                    (chip, filt)
                    for chip, filt in zip(self.chip_list, self.filter_list)
                    if str(chip.chipID).zfill(2) in requested_ids
                ]
                run_chips = [chip for chip, _ in selected_pairs]
                run_filts = [filt for _, filt in selected_pairs]
                nchips_per_fp = len(run_chips)

            for ichip in range(nchips_per_fp):
                i_process = process_counter + ichip
+92 −36
Original line number Diff line number Diff line
import galsim
import numpy as np

from csst_msc_instrument import get_builtin_repository

class FocalPlane(object):
    def __init__(self, chip_list=None, survey_type="Photometric", bad_chips=None):
        """Get the focal plane layout"""
        self.nchips = 42
        self.ignore_chips = []

        if bad_chips is None:
            self.bad_chips = []
        else:
            self.bad_chips = bad_chips
            for chip_id in bad_chips:
                self.ignore_chips.append(chip_id)

class FocalPlane(object):
    _SURVEY_ALIASES = {
        "all": None,
        "photometric": "photometric",
        "spectroscopic": "spectroscopic",
        "fgs": "guiding",
        "guiding": "guiding",
    }

    def __init__(
        self,
        chip_list=None,
        survey_type="Photometric",
        bad_chips=None,
        instrument_repository=None,
    ):
        """Get the focal plane layout."""
        self.instrument_repository = (
            instrument_repository or get_builtin_repository()
        )
        self._detectors = tuple(self.instrument_repository.iter_detectors())
        self.chip_ids = tuple(detector.id for detector in self._detectors)
        self._all_chip_ids = frozenset(self.chip_ids)
        self.nchips = len(self.chip_ids)

        bad_chip_ids = () if bad_chips is None else bad_chips
        self.bad_chips = frozenset(
            self._normalize_chip_ids(bad_chip_ids, "bad_chips")
        )
        if chip_list is not None:
            for i in range(42):
                if i + 1 not in chip_list:
                    self.ignore_chips.append(i + 1)
        elif survey_type == "Photometric":
            for i in range(5):
                self.ignore_chips.append(i + 1)
                self.ignore_chips.append(i + 26)
            self.ignore_chips.append(10)
            self.ignore_chips.append(21)
            for i in range(31, 43):
                self.ignore_chips.append(i)
        elif survey_type == "Spectroscopic":
            for i in range(6, 26):
                if i == 10 or i == 21:
                    continue
            selected = self._normalize_chip_ids(chip_list, "chip_list")
        else:
                    self.ignore_chips.append(i)
            for i in range(31, 43):
                self.ignore_chips.append(i)
        elif survey_type == "FGS":
            for i in range(1, 31):
                self.ignore_chips.append(i)
            selected = self._select_survey_chips(survey_type)

        selected -= self.bad_chips
        self.selected_chip_ids = tuple(
            chip_id for chip_id in self.chip_ids if chip_id in selected
        )
        self.nselected = len(self.selected_chip_ids)
        self.ignore_chips = frozenset(self._all_chip_ids - selected)

        self.nchip_x = 6
        self.nchip_y = 5
@@ -55,12 +61,62 @@ class FocalPlane(object):
    def getChipLabel(self, chipID):
        return str("0%d" % chipID)[-2:]

    def _normalize_chip_ids(self, chip_ids, parameter_name):
        if isinstance(chip_ids, (str, bytes, int)):
            raise TypeError(f"{parameter_name} must be an iterable of chip IDs")

        normalized = set()
        for chip_id in chip_ids:
            try:
                normalized.add(self.instrument_repository.detector(chip_id).id)
            except KeyError as error:
                raise ValueError(
                    f"unknown chip ID {chip_id!r} in {parameter_name}"
                ) from error
        return normalized

    def _select_survey_chips(self, survey_type):
        if survey_type is None:
            survey = None
        else:
            key = str(survey_type).strip().lower()
            if key not in self._SURVEY_ALIASES:
                valid = ", ".join(sorted(self._SURVEY_ALIASES))
                raise ValueError(
                    f"unknown survey type {survey_type!r}; expected one of {valid}"
                )
            survey = self._SURVEY_ALIASES[key]

        if survey is None:
            return set(self._all_chip_ids)
        return {
            detector.id
            for detector in self._detectors
            if detector.instance.optical_role.survey == survey
        }

    def canonicalize_chip_ids(self, chip_ids):
        """Validate IDs and return them in focal-plane order."""
        normalized = self._normalize_chip_ids(chip_ids, "chip_ids")
        return tuple(chip_id for chip_id in self.chip_ids if chip_id in normalized)

    def iter_detectors(self):
        return iter(self._detectors)

    def iter_selected_detectors(self):
        selected = set(self.selected_chip_ids)
        return (
            detector for detector in self._detectors if detector.id in selected
        )

    def isBadChip(self, chipID):
        """Check if chip #(chipID) on the focal plane is bad or not"""
        return chipID in self.bad_chips
        """Check if chip #(chipID) on the focal plane is bad or not."""
        chip_id = self.canonicalize_chip_ids([chipID])[0]
        return chip_id in self.bad_chips

    def isIgnored(self, chipID):
        return chipID in self.ignore_chips
        chip_id = self.canonicalize_chip_ids([chipID])[0]
        return chip_id in self.ignore_chips

    def getTanWCS(self, ra, dec, img_rot, pix_scale, xcen=None, ycen=None, logger=None):
        """Get the WCS of the image mosaic using Gnomonic/TAN projection
+10 −5
Original line number Diff line number Diff line
@@ -23,8 +23,15 @@ class Chip(FocalPlane):
        logger=None,
        instrument_repository=None,
    ):
        # Get focal plane (instance of paraent class) info
        super().__init__()
        # Get focal plane (instance of parent class) info
        settings = (config or {}).get("instrument_data", {})
        repository = instrument_repository or InstrumentRepository.from_settings(
            settings
        )
        super().__init__(
            survey_type=None,
            instrument_repository=repository,
        )
        self._set_attributes_from_config(config)

        self.logger = logger
@@ -36,9 +43,7 @@ class Chip(FocalPlane):
        # Resolve hardware, readout, calibration, and placement data through
        # the validated repository. Scalar attributes remain as a compatibility
        # facade for simulation code that has not yet migrated.
        self.detector_model = self._resolve_detector_model(
            config, instrument_repository
        )
        self.detector_model = self._resolve_detector_model(config, repository)
        for key, value in self.detector_model.legacy_attributes.items():
            setattr(self, key, value)
        self.image_layout = "active"
+54 −10
Original line number Diff line number Diff line
@@ -46,6 +46,7 @@ def add_channel_read_noise(image, channels, channel_slices, seed: int):
        view += generator.normal(0.0, channel.electronics.read_noise_e, view.shape)
    return image


def apply_channel_gain(
    image,
    channels,
@@ -74,24 +75,67 @@ def apply_channel_gain(
    return image, realized_gains


def apply_crosstalk(image, data_slices, matrix: np.ndarray):
    """Apply a channel-to-channel coupling matrix to active channel data."""
def _raw_to_detector(array: np.ndarray, transform: str) -> np.ndarray:
    """Undo a channel's detector-to-raw orientation transform."""

    inverse_transforms = {
        "identity": lambda value: value,
        "flip_x": lambda value: value[:, ::-1],
        "flip_y": lambda value: value[::-1, :],
        "rotate_180": lambda value: value[::-1, ::-1],
        "transpose": lambda value: value.T,
        "transpose_flip_x": lambda value: value.T[::-1, :],
        "transpose_flip_y": lambda value: value.T[:, ::-1],
        "anti_transpose": lambda value: value.T[::-1, ::-1],
    }
    return inverse_transforms[transform](array)


    data_slices = tuple(data_slices)
def apply_crosstalk(image, channels, layout: str, matrix: np.ndarray):
    """Apply channel coupling to active data aligned in raw-readout coordinates.

    Channels with different active shapes are mixed only over their common
    readout-coordinate region. Pixels without a corresponding source pixel
    retain only contributions from channels that cover that location.
    """

    channels = tuple(channels)
    matrix = np.asarray(matrix)
    channel_count = len(data_slices)
    channel_count = len(channels)
    if matrix.shape != (channel_count, channel_count):
        raise ValueError(
            f"crosstalk matrix shape {matrix.shape} does not match "
            f"{channel_count} channels"
        )

    if layout == "active":
        data_slices = tuple(channel.detector_slices for channel in channels)
        source = [
            channel.detector_to_raw(image.array[channel_slice].copy())
            for channel, channel_slice in zip(channels, data_slices)
        ]
    elif layout == "raw_mosaic":
        data_slices = tuple(channel.raw_mosaic_data_slices for channel in channels)
        source = [image.array[channel_slice].copy() for channel_slice in data_slices]
    shapes = {value.shape for value in source}
    if len(shapes) != 1:
        raise ValueError("crosstalk requires channels with equal active shapes")
    for target_index, target_slice in enumerate(data_slices):
    else:
        raise ValueError(f"unsupported image layout {layout!r}")

    for target_index, (channel, target_slice) in enumerate(
        zip(channels, data_slices)
    ):
        target = np.zeros_like(source[target_index])
        for source_index, source_array in enumerate(source):
            target += matrix[target_index, source_index] * source_array
            coefficient = matrix[target_index, source_index]
            if coefficient == 0:
                continue
            rows = min(target.shape[0], source_array.shape[0])
            columns = min(target.shape[1], source_array.shape[1])
            target[:rows, :columns] += (
                coefficient * source_array[:rows, :columns]
            )
        if layout == "active":
            target = _raw_to_detector(
                target, channel.geometry.detector_to_raw
            )
        image.array[target_slice] = target
    return image
+3 −2
Original line number Diff line number Diff line
@@ -318,10 +318,11 @@ def add_crosstalk(self, chip, filt, tel, pointing, catalog, obs_param):
        1.0,
    ])

    data_slices = chip.detector_model.channel_data_slices(chip.image_layout)
    channels = tuple(chip.detector_model.iter_channels())
    chip.img = channel_effects.apply_crosstalk(
        chip.img,
        data_slices=data_slices,
        channels=channels,
        layout=chip.image_layout,
        matrix=crosstalk,
    )

Loading