updated the measurement and evaluation parts

# Pipeline for source injection on CSST images

## Installation

Create or update the `csst-injection` Conda environment, install SExtractor,

PSFEx, and this package with:

```

    source install.sh

    ./install.sh

 ```

The native tools are installed from conda-forge as

`astromatic-source-extractor` and `astromatic-psfex`. To update the environment

without reinstalling the editable project manually, rerun the same command.

The installer also migrates environments created by older versions from

`cfitsio=3.*` to the compatible `cfitsio=4.6.4` release.

## Usage

### Run L1 Detrending pipeline (optional)

* If the L1 images are not available, we need first to run the CSST detrending pipeline on the raw images to get the L1 calibrated products. Otherwise this step can be ignored.

    * Detect and do photometry measurements on the L1 calibrated images

    * Re-detect and redo photometry measurements on the injected images

* modify the corresponding ```config_photometry.yaml``` file to configure all running options

    * Injected multi-extension L1 products are consumed directly from their

      `IMAGE`, `IVAR`, and `FLAG` extensions.

    * Use `input_image_list` for flat injection outputs, or the pointing and

      chip settings for products stored below pointing directories.

```

    run_csst_photometry /path/to/config_photometry.yaml

```

    * Injected objects

    * Detection completeness

    * Photometric accuracy

* ```evaluation/measurement_catalog_evaluation.ipynb``` demonstrates one-to-one

  cross-matching and completeness/photometry evaluation for newly generated

  ```measurement_pipeline``` catalogs.

# input_dir: "/public/home/fangyuedong/project/demo_csst_injection/50sqDeg_Photo_W2/"

input_dir: "/public/home/fangyuedong/project/demo_csst_injection/test/"

flag_weight_dir: "/public/home/fangyuedong/project/demo_csst_injection/L1_products/50sqDeg_Photo_W2/"

output_dir: "/public/home/fangyuedong/project/demo_csst_injection/test/"

pointing_label_list: ["MSC_0000000"]

chip_label_list: ["07"]

calib_data_path: "/public/home/fangyuedong/project/calib_data/"

 No newline at end of file

input_dir: "/home/yuedong/Work/demo_csst_injection/new_L1/"

# For flat injection outputs, provide their complete paths here and omit the

# pointing traversal below. Multi-extension products use IMAGE/IVAR/FLAG.

input_image_list:

  - "/home/yuedong/Work/demo_csst_injection/new_L1/CSST_MSC_MS_WIDE_20260820131350_20260820131620_10100328199_09_L1_V01_injected.fits"

# Optional legacy directory containing separate wht/flg files. It is ignored

# when the input product contains IVAR and FLAG extensions.

flag_weight_dir: null

output_dir: "/home/yuedong/Work/demo_csst_injection/photometry/"

pointing_label_list: null

chip_label_list: null

calib_data_path: null

name: csst-injection

channels:

  - conda-forge

  - nodefaults

dependencies:

  - python=3.11

  - pip

  - setuptools<82

  - wheel

  - cython=3.0.6

  - numpy=1.26.4

  - cfitsio=4.6.4

  - astromatic-source-extractor=2.28.2

  - astromatic-psfex=3.24.2

import numpy as np

"""Catalog cross-matching helpers used by the evaluation notebooks.

The legacy :func:`match_catalogs_img` API is retained for older notebooks.

New code should use :func:`match_catalogs_nearest`, which records separations,

invalid coordinates, ambiguous neighborhoods, and one-to-one assignments.

"""

import astropy.units as u

import numpy as np

from astropy.coordinates import SkyCoord

from astropy.io import fits

from astropy.io import ascii

from sklearn.neighbors import BallTree

from astropy.io import ascii, fits

from astropy.table import Table

from scipy.spatial import cKDTree

TU_catalog = "injected_bkgsub_img.cat"

source_catalog = "extracted_injected_bkgsub_img.fits"

def convert_catalog(catname):

def convert_catalog(catname, output_path=None):

    """Convert an ASCII catalog to FITS and return the output path."""

    output_path = output_path or "test_ascii_to_fits.fits"

    text_file = ascii.read(catname)

    text_file.write('test_ascii_to_fits.fits', overwrite=True)

    text_file.write(output_path, overwrite=True)

    return output_path

def read_catalog(catname, ext_num=1, ra_name='ra', dec_name='dec', col_list=[]):

    hdu = fits.open(catname)

    hdr = hdu[ext_num].header

    # print(hdr)

def read_catalog(catname, ext_num=1, ra_name="ra", dec_name="dec", col_list=None):

    """Read selected columns from a FITS catalog (legacy API)."""

    col_list = [] if col_list is None else col_list

    with fits.open(catname, memmap=False) as hdu:

        data = hdu[ext_num].data

    ra = data[ra_name]

    dec = data[dec_name]

    col_other = []

    if len(col_list) > 0:

        for col in col_list:

            col_other.append(data[col])

        ra = np.asarray(data[ra_name]).copy()

        dec = np.asarray(data[dec_name]).copy()

        col_other = [np.asarray(data[col]).copy() for col in col_list]

    return ra, dec, col_other

def match_catalogs_sky(ra1, dec1, ra2, dec2, max_dist=0.6, others1=[], others2=[], thresh=[]):

    cat1 = SkyCoord(ra=ra1 * u.degree, dec=dec1 * u.degree)

    cat2 = SkyCoord(ra=ra2 * u.degree, dec=dec2 * u.degree)

    idx1, idx2, d2d, d3d = cat2.search_around_sky(cat1, 1 * u.deg)

    # print(idx1)

    # print(idx2)

    # print(np.shape(idx1))

    # print(np.shape(idx2))

    # TODO

def match_catalogs_img(x1, y1, x2, y2, max_dist=2, others1=[], others2=[], thresh=[]):

    cat1 = np.array([(x, y) for x, y in zip(x1, y1)])

    cat2 = np.array([(x, y) for x, y in zip(x2, y2)])

    tree = BallTree(cat2)

    idx1 = tree.query_radius(cat1, r=max_dist)

    tree = BallTree(cat1)

    idx2 = tree.query_radius(cat2, r=max_dist)

    tot = 0

    for idx in idx1:

        if len(idx) == 0:

def _empty_match_table(size):

    return Table(

        {

            "truth_index": np.arange(size, dtype=np.int64),

            "measured_index": np.full(size, -1, dtype=np.int64),

            "separation": np.full(size, np.nan, dtype=float),

            "matched": np.zeros(size, dtype=bool),

            "valid_truth_position": np.zeros(size, dtype=bool),

            "candidate_count": np.zeros(size, dtype=np.int32),

            "ambiguous": np.zeros(size, dtype=bool),

            "duplicate_rejected": np.zeros(size, dtype=bool),

        }

    )

def _assign_pairs(size, truth_valid, measured_valid_indices, pair_truth,

                  pair_measured_local, distances, candidate_count,

                  one_to_one):

    result = _empty_match_table(size)

    result["valid_truth_position"] = truth_valid

    result["candidate_count"] = candidate_count

    result["ambiguous"] = candidate_count > 1

    if len(distances) == 0:

        return result

    used_truth = set()

    used_measured = set()

    rejected_truth = set()

    for pair_index in np.argsort(distances, kind="stable"):

        truth_index = int(pair_truth[pair_index])

        measured_index = int(measured_valid_indices[pair_measured_local[pair_index]])

        if truth_index in used_truth:

            continue

        if one_to_one and measured_index in used_measured:

            rejected_truth.add(truth_index)

            continue

        if len(idx) > 1: print(len(idx))

        tot += 1

    print("number of matched sources = ", tot)

        result["measured_index"][truth_index] = measured_index

        result["separation"][truth_index] = float(distances[pair_index])

        result["matched"][truth_index] = True

        used_truth.add(truth_index)

        used_measured.add(measured_index)

    if rejected_truth:

        result["duplicate_rejected"][list(rejected_truth)] = True

    return result

def match_catalogs_nearest(x_truth, y_truth, x_measured, y_measured,

                           max_dist=2.0, one_to_one=True):

    """Match image coordinates and return one row for every truth object.

    Candidate pairs are considered in order of increasing distance. With

    ``one_to_one=True``, each measured source can be assigned to at most one

    truth object. ``separation`` is in pixels.

    """

    x_truth = np.asarray(x_truth, dtype=float)

    y_truth = np.asarray(y_truth, dtype=float)

    x_measured = np.asarray(x_measured, dtype=float)

    y_measured = np.asarray(y_measured, dtype=float)

    if x_truth.shape != y_truth.shape or x_measured.shape != y_measured.shape:

        raise ValueError("x and y arrays must have matching shapes")

    if max_dist <= 0:

        raise ValueError("max_dist must be positive")

    truth_valid = np.isfinite(x_truth) & np.isfinite(y_truth)

    measured_valid = np.isfinite(x_measured) & np.isfinite(y_measured)

    measured_valid_indices = np.flatnonzero(measured_valid)

    result = _empty_match_table(len(x_truth))

    result["valid_truth_position"] = truth_valid

    if not np.any(truth_valid) or not np.any(measured_valid):

        return result

    truth_valid_indices = np.flatnonzero(truth_valid)

    truth_xy = np.column_stack((x_truth[truth_valid], y_truth[truth_valid]))

    measured_xy = np.column_stack(

        (x_measured[measured_valid], y_measured[measured_valid])

    )

    tree = cKDTree(measured_xy)

    neighbor_indices = tree.query_ball_point(truth_xy, r=max_dist)

    neighbor_distances = []

    for item_index, (point, indices) in enumerate(zip(truth_xy, neighbor_indices)):

        indices = np.asarray(indices, dtype=int)

        separations = np.linalg.norm(measured_xy[indices] - point, axis=1)

        order = np.argsort(separations, kind="stable")

        neighbor_indices[item_index] = indices[order]

        neighbor_distances.append(separations[order])

    candidate_count = np.zeros(len(x_truth), dtype=np.int32)

    candidate_count[truth_valid_indices] = [len(item) for item in neighbor_indices]

    pair_truth = []

    pair_measured_local = []

    distances = []

    for local_truth, (indices, separations) in enumerate(

        zip(neighbor_indices, neighbor_distances)

    ):

        pair_truth.extend([truth_valid_indices[local_truth]] * len(indices))

        pair_measured_local.extend(indices)

        distances.extend(separations)

    return _assign_pairs(

        len(x_truth), truth_valid, measured_valid_indices,

        np.asarray(pair_truth, dtype=np.int64),

        np.asarray(pair_measured_local, dtype=np.int64),

        np.asarray(distances, dtype=float), candidate_count, one_to_one,

    )

def match_catalogs_sky(ra1, dec1, ra2, dec2, max_dist=0.6,

                       one_to_one=True, **_ignored):

    """Match sky coordinates within ``max_dist`` arcseconds."""

    ra1 = np.asarray(ra1, dtype=float)

    dec1 = np.asarray(dec1, dtype=float)

    ra2 = np.asarray(ra2, dtype=float)

    dec2 = np.asarray(dec2, dtype=float)

    if ra1.shape != dec1.shape or ra2.shape != dec2.shape:

        raise ValueError("RA and DEC arrays must have matching shapes")

    truth_valid = np.isfinite(ra1) & np.isfinite(dec1)

    measured_valid = np.isfinite(ra2) & np.isfinite(dec2)

    measured_valid_indices = np.flatnonzero(measured_valid)

    result = _empty_match_table(len(ra1))

    result["valid_truth_position"] = truth_valid

    if not np.any(truth_valid) or not np.any(measured_valid):

        return result

    truth_valid_indices = np.flatnonzero(truth_valid)

    truth = SkyCoord(ra1[truth_valid] * u.deg, dec1[truth_valid] * u.deg)

    measured = SkyCoord(ra2[measured_valid] * u.deg, dec2[measured_valid] * u.deg)

    truth_local, measured_local, separation, _ = measured.search_around_sky(

        truth, max_dist * u.arcsec

    )

    pair_truth = truth_valid_indices[np.asarray(truth_local, dtype=int)]

    pair_measured_local = np.asarray(measured_local, dtype=int)

    distances = separation.to_value(u.arcsec)

    candidate_count = np.zeros(len(ra1), dtype=np.int32)

    np.add.at(candidate_count, pair_truth, 1)

    return _assign_pairs(

        len(ra1), truth_valid, measured_valid_indices, pair_truth,

        pair_measured_local, distances, candidate_count, one_to_one,

    )

def match_catalogs_img(x1, y1, x2, y2, max_dist=2, **_ignored):

    """Legacy radius-match API returning neighbors in both directions."""

    x1 = np.asarray(x1, dtype=float)

    y1 = np.asarray(y1, dtype=float)

    x2 = np.asarray(x2, dtype=float)

    y2 = np.asarray(y2, dtype=float)

    def query(x_from, y_from, x_to, y_to):

        output = np.empty(len(x_from), dtype=object)

        output[:] = [np.array([], dtype=int)]

        valid_from = np.isfinite(x_from) & np.isfinite(y_from)

        valid_to = np.isfinite(x_to) & np.isfinite(y_to)

        if not np.any(valid_from) or not np.any(valid_to):

            return output

        to_indices = np.flatnonzero(valid_to)

        tree = cKDTree(np.column_stack((x_to[valid_to], y_to[valid_to])))

        found = tree.query_ball_point(

            np.column_stack((x_from[valid_from], y_from[valid_from])), r=max_dist

        )

        for original_index, local_matches in zip(np.flatnonzero(valid_from), found):

            output[original_index] = to_indices[local_matches]

        return output

    idx1 = query(x1, y1, x2, y2)

    idx2 = query(x2, y2, x1, y1)

    print("number of matched sources = ", sum(len(item) > 0 for item in idx1))

    return idx1, idx2

if __name__ == "__main__":

    convert_catalog(TU_catalog)

    ra_TU, dec_TU, _ = read_catalog('test_ascii_to_fits.fits', ext_num=1, ra_name="ra", dec_name="dec")

    x_TU, y_TU, col_list = read_catalog('test_ascii_to_fits.fits', ext_num=1, ra_name="xImage", dec_name="yImage", col_list=["mag"])

    mag_TU = col_list[0]

    # ra_source, dec_source, _ = read_catalog(source_catalog, ext_num=2, ra_name="ALPHAPEAK_J2000", dec_name="DELTAPEAK_J2000")

    x_source, y_source, _ = read_catalog(source_catalog, ext_num=1, ra_name="X_IMAGE", dec_name="Y_IMAGE")

    # match_catalogs_sky(ra1=ra_TU, dec1=dec_TU, ra2=ra_source, dec2=dec_source)

    idx1, idx2, = match_catalogs_img(x1=x_TU, y1=y_TU, x2=x_source, y2=y_source)

    # print(ra_TU, dec_TU)

    converted = convert_catalog(TU_catalog)

    x_truth, y_truth, _ = read_catalog(

        converted, ra_name="xImage", dec_name="yImage"

    )

    x_source, y_source, _ = read_catalog(

        source_catalog, ra_name="X_IMAGE", dec_name="Y_IMAGE"

    )

    matches = match_catalogs_nearest(x_truth, y_truth, x_source, y_source)

    print(f"number of matched sources = {np.count_nonzero(matches['matched'])}")