Add files via upload

Author: yusuke1999

README.md

@@ -0,0 +1,125 @@
+# PDRL: Position-Dependent Richardson-Lucy deconvolution
+Code for this paper
+
+**Please note that the provided code is not a finished product, and we do not guarantee its functionality or suitability for any specific purpose. Use it at your own risk.**
+
+## Overview
+Position-Dependent Richardson-Lucy (PDRL) deconvolution is an extension of the Richardson-Lucy deconvolution method, specifically designed for images with multiple positional point spread functions (PSFs).
+
+![001](figures/casA_PDRL.jpg)<br>
+Figure 2. (a) X-ray image in the 0.5\UTF{2013}7.0 keV band of Cas A obtained with Chandra. (a-1, -2, -3) Enlarged images specified by the colored frames in (a). (b) Same as (a), but for the RLsv-deconvolved results. The unit of flux in the images is $\rm{photons~cm^{-2}~s^{-1}}$.
+
+The 35 x 35 pixel spacing of PSFs used in the above PDRL image.
+![002](figures/PSFs_35bin.jpg)
+Cas A image (Obs. ID = 4636) and the two-dimensional probabil- ities of the PSFs. The integral of each PSF is normalized to be 1. The PSF color scale is a fixed range. The location of the optical axis is indicated with a green cross.
+
+## Our Environments
+* DISTRIB_ID=Ubuntu
+* DISTRIB_RELEASE=20.04
+* DISTRIB_CODENAME=focal
+* DISTRIB_DESCRIPTION="Ubuntu 20.04.3 LTS"
+* Chandra CIAO-4.13
+* marx-5.5.1
+* python=3.8.2
+* astropy==5.1
+* numpy==1.22.4
+* tqdm==4.64.0
+
+## Description of Each File
+Our code is optimized for Chandra X-ray observed images and utilizes [CIAO](https://cxc.cfa.harvard.edu/ciao/ahelp/merge_obs.html) and [MARX](https://space.mit.edu/ASC/MARX/) for PSF simulation.
+
+### 0. Install MARX and build environment
+Please refer to the CIAO [simulate_psf](https://cxc.cfa.harvard.edu/ciao/ahelp/simulate_psf.html) documentation for instructions on installing MARX and setting up the environment. (If you have not installed CIAO, please refer to the official CIAO documentation.)
+
+### 1. Prepare count map and exposure map
+Prepare the count map and exposure map files required for the PDRL method.
+
+### 2. Obtain ra, dec coordinates for PSF simulation
+Use the "pixel2skycoord.py" script to obtain the sky coordinates for each position and save them in a .txt file for use with CIAO simulate_psf.
+```
+python pixel2skycoord.py [-h] infile outfile psf_bins
+```
+infile: Input counts map file<br>
+outfile: Output .txt file containing image_y, image_x, ra, dec<br>
+psf_bins: Pixel spacing for creating PSF (Smaller values result in better accuracy but require more computation time)
+
+### 3. Generate PSF from coordinates file created in Step 2
+Use CIAO simulate_psf to generate the PSF for each position using the ra, dec .txt file. 
+
+**It is recommended to run the code in the background (e.g., using the "screen" command) as simulating PSFs can take a significant amount of time.**
+```
+bash simulate_psf.sh [-h] infile_txt infile_fits monoenergy outdir
+```
+infile_txt: Input .txt file containing image_y, image_x, ra, dec<br>
+infile_fits: Input event FITS file (Note: not a counts map file)<br>
+monoenergy: Set the PSF monoenergy<br>
+outdir: Output directory for each position's PSF
+
+### 4. PSF adjustment
+Move the center of the PSFs to the image center and adjust all PSFs to the same size.
+```
+repro_psfs.py [-h] infile psfs_dir outfile psfs_bins
+```
+infile: input counts map file for made psfs (Note: not evt file)<br>
+psfs_dir: input raw psf dir for each location<br>
+outfile: output all psf .npz file<br>
+psfs_bins: spacing of pixels you made each position psf<br>
+
+### 5. Run PDRL
+This is PDRL method. Results for all iterations are saved.
+```
+python LDRL.py [-h] thresh_img expmap psfs outdir [--num_iter NUM_ITER] [--lambda_tv LAMBDA_TV]
+               [--data_type DATA_TYPE] [--im_deconv_0_flat] [--poisson_err]
+               [--boundary_px BOUNDARY_PX]
+```
+
+thresh_img: Input counts map file for PDRL<br>
+expmap: Input exposure map file for PDRL<br>
+psfs: Input PSF file for each infile position (.npz)<br>
+outdir: Output directory for each iteration of PDRL results<br>
+--num_iter NUM_ITER: Number of iterations for PDRL (default: 200)<br>
+--lambda_tv LAMBDA_TV: TV (Total Variation) regularization lambda for PDRL (default: 0.002)<br>
+--data_type DATA_TYPE: Numpy data type (default: float64). If "killed" is returned, using float32 can help reduce memory usage.<br>
+--im_deconv_0_flat: Initial value for the 0th iteration of the LDRL method (default: False, meaning the input file)<br>
+--poisson_err: Add a Poisson distribution random number according to the input count file for each iteration (default: False)<br>
+--boundary_px BOUNDARY_PX: Randomly select a PSF from boundary_px pixels near the boundary of the PSFs (default: 1)
+
+### 6. Run PDRL error analysis
+Calculate the uncertainty of the PDRL method using error propagation.
+```
+python PDRL_err.py [-h] thresh_img expmap psfs im_deconv outfile [--data_type DATA_TYPE]
+```
+
+thresh_img: Input counts map file for PDRL
+expmap: Input exposure map file for PDRL
+psfs: Input PSF file for each infile position (.npz)
+im_deconv: Input PDRL result to calculate the next iteration's PDRL uncertainty
+outfile: Output file for PDRL error result
+--data_type DATA_TYPE: Default: float64. Numpy data type. (Note: If `killed` is returned, using float32 can be useful to reduce memory usage)
+
+## Demonstration
+The demonstration code provides the results mentioned above. After completing step 0, run the following code:
+```
+bash demo.sh
+```
+
+## Quick Demonstration
+Place the following files directly in the demo directory:
+* [merged_4636_4637_4639_5319](https://drive.google.com/file/d/1jNSmAmFOXCZtuikxnf2gb5QXl2mmt0Gg/view?usp=share_link)(merged Obs. ID=4636, 4637, 4639, and 5319 counts map and exposure map files)
+* [repro_psfs_35bin.npz](https://drive.google.com/file/d/1qa6gsLghqQbbsElNKirOgW-UdVKlTzIv/view?usp=share_link)(35 × 35 pixels Obs. ID 4636 PSFs)
+
+Then, in the demo directory, run the PDRL method with the following command:
+```
+python ../pyscript/PDRL.py merged_4636_4637_4639_5319/broad_thresh.img merged_4636_4637_4639_5319/broad_thresh.expmap \
+       repro_psfs_35bin.npz PDRL_results --lambda_tv 0 --boundary_px 0
+```
+
+After completion, to obtain the uncertainty of the 200th iteration based on the 199th iteration of the PDRL method, use the following command:
+
+```
+python ../pyscript/PDRL_err.py merged_4636_4637_4639_5319/broad_thresh.img merged_4636_4637_4639_5319/broad_thresh.expmap \
+       repro_psfs_35bin.npz PDRL_results/iter_0199.fits iter_0200_err.fits
+```
+
+## Citing
+If you find our code helpful, we would appreciate it if you could cite the following reference as a token of acknowledgment:

demo/demo.sh

@@ -0,0 +1,25 @@
+#!/bin/bash
+# Merge the data.
+merge_obs 4636/repro/acisf04636_repro_evt2.fits,\
+4637/repro/acisf04637_repro_evt2.fits,\
+4639/repro/acisf04639_repro_evt2.fits,\
+5319/repro/acisf05319_repro_evt2.fits \
+merged_4636_4637_4639_5319/ bands=broad binsize=1
+
+# Create a txt file converted from pixel coordinates to radec coordinates.
+python ../pyscript/pixel2skycoord.py merged_4636_4637_4639_5319/broad_thresh.img pixel2sky_35bin.txt 35
+
+# Simulation of psf by marx using Obs. ID 4636 as a representative.
+bash ../shscript/simulate_psfs.sh pixel2sky_35bin.txt 4636/repro/acisf04636_repro_evt2.fits 2.3 psfs_35bin
+
+# Adjust the psf so that the center of the psf is at the center of the image.
+# Save an array with psf for each location.
+python ../pyscript/repro_psfs.py merged_4636_4637_4639_5319/broad_thresh.img psfs_35bin repro_psfs_35bin.npz 35
+
+# Run PDRL method.
+python ../pyscript/PDRL.py merged_4636_4637_4639_5319/broad_thresh.img merged_4636_4637_4639_5319/broad_thresh.expmap \
+       repro_psfs_35bin.npz PDRL_results --lambda_tv 0 --boundary_px 0
+
+# Run the error of PDRL by the law of error propagation. 
+python ../pyscript/PDRL_err.py merged_4636_4637_4639_5319/broad_thresh.img merged_4636_4637_4639_5319/broad_thresh.expmap \
+       repro_psfs_35bin.npz PDRL_results/iter_0199.fits iter_0200_err.fits

figures/PSFs_35bin.jpg

@@ -0,0 +1,133 @@
+import os
+import random
+from tqdm import tqdm
+import numpy as np
+from scipy.signal import convolve
+import astropy.io.fits as fits
+import argparse
+import lib_fits
+
+
+def main():
+    # Parameters
+    parser = argparse.ArgumentParser(description='This is PDRL method. Results for all iterations are saved.')
+    parser.add_argument('thresh_img', type=str, help='input counts map file for PDRL')
+    parser.add_argument('expmap', type=str, help='input exposure map file for PDRL')
+    parser.add_argument('psfs', type=str, help='input psf file for each infile position (.npz)')
+    parser.add_argument('outdir', type=str, help='output dir for each iteration of PDRL results')
+    parser.add_argument('--num_iter', type=int, default=200, help='default:200. the number of iterations for PDRL')
+    parser.add_argument('--lambda_tv', type=float, default=0.002, help='default:0.002. the TV regularization lambda')
+    parser.add_argument('--data_type', type=str, default='float64', help='default:float64. numpy data type (Note:If `killed` is returned, `float32` is useful to lower the memory.)')
+    parser.add_argument("--im_deconv_0_flat", action='store_true', help='default:False. initial value for the 0th iteration of the ldrl method (`False` mean the input file)')
+    parser.add_argument("--poisson_err", action='store_true', help='default:False. add a Poisson distribution random number according to the input count file for each iteration')
+    parser.add_argument("--boundary_px", type=int, default=1, help='default:1. Randomly select a PSF from `boundary_px` pixels near the boundary of the PSFs')
+    args = parser.parse_args()
+
+    # Loads the thresh.img and expmap.
+    thresh_img, wcs = lib_fits.load_fits(infile=args.thresh_img, data_type=args.data_type)
+    expmap, _ = lib_fits.load_fits(infile=args.expmap, data_type=args.data_type)
+
+    # Loads all PSFs as numpy array.
+    print('PSF loading in progress.')
+    psfs, psfs_bins = load_all_psf(args.psfs, args.data_type)
+
+    # PDRL method.
+    np.random.seed(seed=2023)
+    print('PDRL method is start.')
+    os.makedirs(args.outdir, exist_ok=True)
+    pdrl = Position_Dependent_Richardson_Lucy(args.outdir, thresh_img, expmap, psfs, psfs_bins, args.num_iter, args.lambda_tv, wcs, args.data_type, args.im_deconv_0_flat, args.poisson_err, args.boundary_px)
+    deconvs = pdrl.position_dependent_richardson_lucy()
+
+
+def load_all_psf(psfs_infile, data_type):
+    repro_psfs = np.load(psfs_infile)
+    psfs = np.array(repro_psfs['repro_psfs'], dtype=data_type)
+    psfs_bins = repro_psfs['psfs_bins'].tolist()
+    return psfs, psfs_bins
+
+
+class Position_Dependent_Richardson_Lucy:
+    def __init__(self, outdir, thresh_img, expmap, psfs, psfs_bins, num_iter=200, lambda_tv=0.002, wcs=None, data_type='float64', im_deconv_0_flat=False, poisson_err=False, boundary_px=1):
+        self.outdir = outdir
+        self.thresh_img = thresh_img
+        self.expmap = expmap + 1e-15  # Used to avoid 0 divisions
+        self.psfs = psfs
+        self.psfs_bins = psfs_bins
+        self.num_iter = num_iter
+        self.lambda_tv = lambda_tv
+        self.wcs = wcs
+        self.data_type = data_type
+        self.im_deconv_0_flat = im_deconv_0_flat
+        self.poisson_err = poisson_err
+        self.boundary_px = boundary_px
+
+
+    def _convolve_each_kernel(self, image, kernels, kernels_xgrid, kernels_ygrid):
+        y_shape, x_shape = image.shape
+        *_, k_y_shape, k_x_shape = kernels.shape
+        pad_y, pad_x = (k_y_shape // 2, k_x_shape // 2)
+        pad_conv = np.zeros([y_shape + pad_y*2, x_shape + pad_x*2], dtype=self.data_type)  # For convolution edge works
+        y_range = range(self.boundary_px, y_shape-self.boundary_px)
+        x_range = range(self.boundary_px, x_shape-self.boundary_px)
+        for y in tqdm(y_range):
+            for x in x_range:
+                pad_conv[y:y+k_y_shape, x:x+k_x_shape] += image[y, x] * kernels[(y+kernels_ygrid[y, x])//self.psfs_bins, (x+kernels_xgrid[y, x])//self.psfs_bins]
+        conv = pad_conv[pad_y:-pad_y, pad_x:-pad_x]  # Make array size equal to input size
+        return conv
+
+
+    def _convolve_each_kernel_dependent(self, image, kernels, kernels_xgrid, kernels_ygrid):
+        y_shape, x_shape = image.shape
+        *_, k_y_shape, k_x_shape = kernels.shape
+        k_y_half, k_x_half = (k_y_shape // 2, k_x_shape // 2)
+        conv = np.zeros_like(image, dtype=self.data_type)
+        y_range = range(k_y_half+1+self.boundary_px, y_shape-k_y_half-1-self.boundary_px)
+        x_range = range(k_x_half+1+self.boundary_px, x_shape-k_x_half-1-self.boundary_px)
+        for y in tqdm(y_range):
+            for x in x_range:
+                conv[y, x] = np.sum(image[y-k_y_half:y+k_y_half+1, x-k_x_half:x+k_x_half+1] * kernels[(y+kernels_ygrid[y, x])//self.psfs_bins, (x+kernels_xgrid[y, x])//self.psfs_bins])
+        return conv
+
+    def _generate_poisson_img(self, thresh_img):
+        poisson_func = lambda x: np.random.poisson(x)
+        poisson_thresh_img = poisson_func(thresh_img)
+        return poisson_thresh_img
+
+    def position_dependent_richardson_lucy(self):
+        if self.im_deconv_0_flat:
+            im_deconv = np.full(self.thresh_img.shape, 0.5, dtype=self.data_type)
+        else:
+            im_deconv = np.copy(self.thresh_img) / self.expmap
+
+        eps = 1e-15  # Small regularization parameter used to avoid 0 divisions
+
+        for _iter in range(1, self.num_iter+1):
+            print(f'iter {_iter} start.')
+            print('Bottom calculation in progress.')
+            psfs_xgrid = np.random.randint(-self.boundary_px, self.boundary_px+1, size=im_deconv.shape)
+            psfs_ygrid = np.random.randint(-self.boundary_px, self.boundary_px+1, size=im_deconv.shape)
+            conv = self._convolve_each_kernel(im_deconv, self.psfs, psfs_xgrid, psfs_ygrid) + eps
+
+            if self.poisson_err:
+                relative_blur = self._generate_poisson_img(self.thresh_img) / self.expmap / conv
+            else:
+                relative_blur = self.thresh_img / self.expmap / conv
+
+            # TVregularization
+            grad_x = np.gradient(im_deconv, axis=1)
+            grad_y = np.gradient(im_deconv, axis=0)
+            grad_norm = np.sqrt(grad_x**2+grad_y**2) + eps
+            tv_reg = np.gradient(grad_x / grad_norm, axis=1) + np.gradient(grad_y / grad_norm, axis=0)
+
+            print('Top calculation in progress.')
+            im_deconv = im_deconv/(1-self.lambda_tv*tv_reg)*self._convolve_each_kernel_dependent(relative_blur, self.psfs, psfs_xgrid, psfs_ygrid)
+
+            # Save the results of deconvolution for each iteration
+            outfile = os.path.join(self.outdir, f'iter_{_iter:0>4}.fits')
+            lib_fits.np2fits(np_array=im_deconv, outfile=outfile, wcs=self.wcs)
+        
+        return
+
+
+if __name__ == '__main__':
+    main()

figures/casA_PDRL.jpg

@@ -0,0 +1,104 @@
+import os
+from tqdm import tqdm
+import numpy as np
+import astropy.io.fits as fits
+import argparse
+import lib_fits
+
+
+def main():
+    # Parameters
+    parser = argparse.ArgumentParser(description='Calculate the uncertainty of PDRL method by the law of error propagation.')
+    parser.add_argument('thresh_img', type=str, help='input counts map file for PDRL')
+    parser.add_argument('expmap', type=str, help='input exposure map file for PDRL')
+    parser.add_argument('psfs', type=str, help='input psf file for each infile position (.npz)')
+    parser.add_argument('im_deconv', type=str, help='input a PDRL result to get the next iteration PDRL uncertainty')
+    parser.add_argument('outfile', type=str, help='output file for PDRL error result')
+    parser.add_argument('--data_type', type=str, default='float64', help='default:float64. numpy data type (Note:If `killed` is returned, `float32` is useful to lower the memory.)')
+    args = parser.parse_args()
+
+    # Loads the thresh.img, expmap and im_deconv.
+    thresh_img, wcs = lib_fits.load_fits(infile=args.thresh_img, data_type=args.data_type)
+    expmap, _ = lib_fits.load_fits(infile=args.expmap, data_type=args.data_type)
+    im_deconv, _ = lib_fits.load_fits(infile=args.im_deconv, data_type=args.data_type)
+
+    # Loads all PSFs as numpy array.
+    print('PSF loading in progress.')
+    psfs, psfs_bins = load_all_psf(args.psfs, args.data_type)
+
+    # PDRL error method.
+    print('Estimate the PDRL error method is start.')
+    outdir = os.path.dirname(args.outfile)
+    os.makedirs(outdir, exist_ok=True)
+    pdrl_err = Position_Dependent_Richardson_Lucy_Err(args.outfile, thresh_img, expmap, psfs, psfs_bins, im_deconv, wcs, args.data_type)
+    deconvs_err = pdrl_err.position_dependent_richardson_lucy_err()
+    print('Done.')
+
+
+def load_all_psf(psfs_infile, data_type):
+    repro_psfs = np.load(psfs_infile)
+    psfs = np.array(repro_psfs['repro_psfs'], dtype=data_type)
+    psfs_bins = repro_psfs['psfs_bins'].tolist()
+    return psfs, psfs_bins
+
+
+class Position_Dependent_Richardson_Lucy_Err:
+    def __init__(self, outfile, thresh_img, expmap, psfs, psfs_bins, im_deconv, wcs=None, data_type='float64'):
+        self.outfile = outfile
+        self.thresh_img = thresh_img
+        self.expmap = expmap + 1e-15  # Used to avoid 0 divisions
+        self.psfs = psfs
+        self.psfs_bins = psfs_bins
+        self.im_deconv = im_deconv
+        self.wcs = wcs
+        self.data_type = data_type
+
+
+    def _convolve_each_kernel(self, image, kernels):
+        y_shape, x_shape = image.shape
+        *_, k_y_shape, k_x_shape = kernels.shape
+        pad_y, pad_x = (k_y_shape // 2, k_x_shape // 2)
+        pad_conv = np.zeros([y_shape + pad_y*2, x_shape + pad_x*2], dtype=self.data_type)  # For convolution edge works
+        y_range = range(y_shape)
+        x_range = range(x_shape)
+        for y in tqdm(y_range):
+            for x in x_range:
+                pad_conv[y:y+k_y_shape, x:x+k_x_shape] += image[y, x] * kernels[y//self.psfs_bins, x//self.psfs_bins]
+        conv = pad_conv[pad_y:-pad_y, pad_x:-pad_x]  # Make array size equal to input size
+        return conv
+
+
+    def _convolve_each_kernel_dependent(self, image, kernels):
+        y_shape, x_shape = image.shape
+        *_, k_y_shape, k_x_shape = kernels.shape
+        k_y_half, k_x_half = (k_y_shape // 2, k_x_shape // 2)
+        conv = np.zeros_like(image, dtype=self.data_type)
+        y_range = range(k_y_half+1, y_shape-k_y_half-1, 1)
+        x_range = range(k_x_half+1, x_shape-k_x_half-1, 1)
+        for y in tqdm(y_range):
+            for x in x_range:
+                conv[y, x] = np.sum(image[y-k_y_half:y+k_y_half+1, x-k_x_half:x+k_x_half+1] * kernels[y//self.psfs_bins, x//self.psfs_bins])
+        return conv
+
+    
+    def position_dependent_richardson_lucy_err(self):
+        psfs_sq = self.psfs ** 2
+        expmap_sq = self.expmap ** 2
+
+        eps = 1e-15  # Small regularization parameter used to avoid 0 divisions
+        
+        print('Bottom calculation in progress.')
+        conv = self._convolve_each_kernel(self.im_deconv, self.psfs) + eps
+        relative_blur = self.thresh_img / expmap_sq / conv**2
+
+        print('Top calculation in progress.')
+        im_deconv_err = self.im_deconv*np.sqrt(self._convolve_each_kernel_dependent(relative_blur, psfs_sq))
+
+        # Save the results of deconvolution err
+        lib_fits.np2fits(np_array=im_deconv_err, outfile=self.outfile, wcs=self.wcs)
+        
+        return
+
+
+if __name__ == '__main__':
+    main()