Requested functionality - EranOfek/AstroPack GitHub Wiki

List of requested functionality

Priority options: low, medium, high, urgent

New functionality

Reference images

  • Define a sub-image size - e.g., 256x256 (300x300 after overlap)

If we believe that the ULTRASAT PSF does not change too much between each 2 points of the 25-point grid in radius, we can set a sub-image size to be 512 pixels with 45 pixels overlap: 2 * (512-45) + 9 * (512-90) = 4732, while the size of the tile is 4738, but we can just ignore the outer corner of 6 x 6 pixels. Thus we will have 11 x 11 = 121 sub-images per tile (484 sub-images per exposure).

  • Define reference images grid size - e.g., 256x256 (500x500 after overlap)
  • For each new image, check for all sections that are fully contained within the reference images sub-image grid.
  • Question: what is the efficiency of this process in terms of the fraction of images/pixels that will not be used for reference construction? Very rough estimate of 70% efficiency.
  • Have we got any tools for this?
    A straightforward solution would be taking a sufficiently fine (sub)grid of HEALPix and check which of them fall into each of the images. However, at 0.2 pix scale = 0.25 arcsec for LAST, Nside = 2^20 = 1048576 => N = 12 * (2^20)^2 ~ 10^13 points on the sphere -- way too much!
    Another approach could be first to work on a coarser grid of Nside = 10, N = 12 * (2^10)^2 ~ 10^7 points of ~ 3.5 arcmin resolution. As soon as all the refs are built, look which of the grid points fall into each of the refs and make a table where 2-4 grid points are assigned for each of them. Given a new image, first find all the refs whose grid points are the neighbors of the new's grid points (or the same). Then remap all the found refs to the pixel grid of the new and do the difference (also very slow?). Or, alternatively, if we are not willing to work with relatively small intersections, we may choose some set of points from the border of a rectangle centered of the new's center (say, 180 x 180 pix => border line = 720 pixels, take 12 of them, i.e. 3 on each side) and look if each of these points is inside one of the pre-selected refs. After that work only with the refs that contain one of such points. Their number is not likely to be more than 5-6.

Matched source analysis

Light curves Statistics

DB access

Requested by : Eran Priority : High Assigned : Sasha, Chen

DONE - see db.Db and pipeline.last

Light curves

  • Do we need a separate matlab class for light curves or they will be just a FITS product of a DB query?

BJD

Full implementation needed (UTC -> TDB), in both AstroCatalog and AstroImage and PhotonEvents.

DONE - celestial.time.barycentricJD is ready. Populated in DB level.

Multi-iteration PSF fitting

Requested by : Eran Priority : High Assigned : Sasha, Eran

A high-level function for multi-iteration PSF fitting.

Functions:

  • injectStamps
  • add bit mask for "in halo of bright star"
  • start-to-end function that populate: Back, Var, Mask, PSF, and Cat.
    • Populate background and variance (imProc.background.background)
    • Find sources with S/N>ThresholdPSF (=30) (imUtil.sources.findSources) Use Gaussian + Delta.
    • Remove CR
    • Use sources to extract PSF (imUtil.psf.constructPSF ???)
    • Multi itearation PSF fitting with
      • Arguments: Threshold=[100 30 5]
      • In each iteration
        • Measure background, unless remeasureBack=false (not in 1st iteration)
        • Measure Var, or construct Var by adding previous iteration subtracted star to Var map.
        • Find sources using PSF, PSF+delta+extended PSF
        • Measure sources properties
        • Subtract stars
        • go to next iteration

constructPSF - expedite

Can we skip the moment estimation in constructPSF and use values obtained from findMeasureSources? Likely related to multi iteration PSF fitting

New pipeline scheme

  • Master configuration file in LAST.Pipeline..yml
  • Separate pipeline to several steps:
  • Read images and test validity.
    • pipeline.generic.readAndQualify (done but not tested)
    • Log bad images + move to failed dir.
  • Analyze individual images
    • To be implemented in: pipeline.generic.multiRaw2proc
    • dedark, deflat, de-overscan, gain, linearity (using CalibImages)
    • optional: break image to sub images
    • Estimate background and variance - consider use a scalar estimator (see below)
    • Initial star finding (SN>30) and PSF estimation
    • Check if a minimum number of stars were detected, and FWHM is reasonable. If not, do not continue, prep msg, and set Status to false.
    • Multi-iteration PSF fitting, consider modifying the variance image after each iteration by adding the subtracted star to the variance image. In this case, in the next iteration, the variance will take into account the source noise.
    • In each step xcorr with delta fun, PSF + extended PSF.
    • Main parameters of multi-iteration PSF fitting:
      • List of SN steps (e.g., 100, 30, 5)
    • Output: XPEAK, YPEAK, PEAK_VAL, X1, Y1, [X2, Y2, XY], [APER_1/2/3], PSF phot (including FLUX, MAG, ERR, SN), FLAGS
    • Future optional: spikes detection
    • Perform astrometry
    • Add astrometry to catalog
  • Merge catalogs of individual images (consider removing the FITS table).
    • pipeline.generic.mergeCat
    • Generate merged matrices in HDF5 and optional matched table in FITS table.
  • Coadd images.
    • pipeline.generic.procCoadd
    • Require improved coaddition code
  • Image subtraction pipeline.
    • pipeline.generic.transientsDetection
  • Each step will be in its own try/catch and will send failed images to failed dir + log file. treatFailure

Lanczos interpolation function

Requested by : Eran Priority : Low Assigned :

Speed up

exp of cube

Requested by : Eran Priority : Medium Assigned : Dharv

Calculate exp(Cube*Mult).*Norm

functions:

  • Result = tools.array.mex.expCubeNorm_single(Cube, Mult, Norm)
  • Result = tools.array.mex.expCubeNorm_double(Cube, Mult, Norm)
  • Result = tools.array.mex.expCubeNorm(Cube, Mult, Norm)
  • Embed in existing code.

tools.math.ode.rkn1210 / for matrices

Requested by : Eran Priority : Medium Assigned :

tools.math.ode.rkn1210 works on vectors. Required a matrix version (e.g., like ode45).

imUtil.cut.image2cutouts faster version

Requested by : Eran Priority : Medium Assigned :

imUtil.cut.image2cutouts already use a mex code. Check if possible to expedite it using SIMD and multi-threading.

Fast median on array

The fast_median function calculates the median of an entire array faster than the median. We need tools.math.stat.fastMedian that work on the dimension of the array and also can ignore NaNs.

In addition, inspect the fast_median C code and see if it can be expedited.

Faster 2D interpolation

Requested by : Eran Priority : Medium Assigned :

Splinter: https://github.com/bgrimstad/splinter is a fast 2D interpolation code that we can use - we can consider the following improvements:

  • A version for single arrays and double arrays.
  • Extending to nearest interpolation for single, double, uint16, uint32, int16, in32.
  • Check if it is possible to make it faster by using multi-threading, SIMD, BLAS and other tricks.

mfind_bin - faster version

We have a mex version but it is not identical in output to mfind_bin (NaN treatment).

new version of tools.find.mfind_bin is available - x11-x14 faster. Need to test in real pipeline. Default of UseMex is false.

Fast parallel 1D interpolation

Consider using: https://www.mathworks.com/matlabcentral/fileexchange/70015-fastinterpcol-fast-1d-interpolation-along-column-dimension

In existing code.

Fft shift improvement

Requested by : Eran Priority : Medium Assigned :

imUtil.trans.shift_fft is the most important bottleneck for high galactic latitude fields.

Specifically, most of its run time is spent on two lines:

SX = ifft(fft(Image, [], 2).*KernelX, [], 2);
ShiftedImage = real(ifft(fft(SX, [], 1).*KernelY, [], 1);

Note that Image, KernelX, KernelY are cubes!

Since MATLAB already uses FFTW and BLAS it will be hard to improve, but if someone has an idea please suggest...

For additional FFT implementations - check out: https://dsp.stackexchange.com/questions/24375/fastest-implementation-of-fft-in-c

tools.array.conditionalReplace

Function that do: M(A>B)=V.

DONE but not implemented 3-4 times faster.

Need to implement in imUtil.image.moment2

isempty_cell

Requested by : Eran Priority : Low Assigned : Dan

The tools.cell.isempty_cell function is slow. Can we expedite?

DONE

isnan_cell

Requested by : Eran Priority : Low Assigned : Dan

The tools.cell.isnan_cell function is slow. Can we expedite?

DONE

Sphere dist

for small angles Strategy Compare: delta dec Use: sincos Or use Taylor approx.

DONE - celestial.coo.mex.haversine but not implemented.

Aper phot

Chi2 on cubes with flags

Cutouts improvement

New version (much faster) imUtil.cut.mex.mex_cutout2.cpp Test some more, ask for more efficient version, and implement.

imUtil.trans.shift_fft

In all but last iterations can we use a smaller psf with fit radius?

Corr2 on a cube of filters

⚠️ **GitHub.com Fallback** ⚠️