The mextractor function

in the astronomy & astrophysics toolbox for MATLAB

Description:

The mextractor function is a method in the SIM class that enable source finding and measurment in astronomical images. The function uses matched filter to locate sources and perform variery of measurments, including aperture and PSF photometry. The function can automaticaly find the image PSF for matched filter, and perform a careful background subtraction which is critical for many purposes.

This file is accessible through the manual package (i.e. manual.mextractor).

Credit

If you are using this code or products in your scientific publication please give a reference to Ofek (2014; ascl.soft 07005).

License

Unless specified otherwise this code and products are released under the GNU general public license version 3.

Installation

See http://weizmann.ac.il/home/eofek/matlab/doc/install.html for installation instruction and additional documentation.

Reduction steps

mextractor is a flexible code to find, extract and measure sources in astronomical images. mextractor can run on a single or multiple images, and can find sources and or measure user supplied list of sources. mextractor is part of the SIM class and uses the information available within the SIM object, and if some information does not exist, mextractor will generate this information and store it in the SIM object. For example, if the SIM object contains the background and PSF, this information will be used by mextractor, and if it does not exist, mextractor will generate this information.

The user can control the behavior of mextractor via a large number of keywords and full list of the keywords is available by typing:

help mextractor

% or

doc mextractor

mextractor perform (optionally) the following steps:

gain correction - correct the images gain to 1. This is performed using the gain_correct method.

PSF extraction - In order to extract the PSF, mextractor call itself to find the brightest sources and use them to generate the PSF. By default, this is done using the psf_cat_selector and psf_estimator methods.

Mask generation - mextractor propogate and add bit masks to the MASK object associated with the SIM image.

Background and noise estimation - Estimate the background and noise images using the background method, and subtract the background from the image. The background estimation is a critical step for source detection. Any small bias in the background estimation may lead to the detection of false sources, or missing real sources. The default method for background and noise estimation (the 'mode_fit' option in the background function) is to fit a Gaussian to the histogram of counts.

Filter the image - The source detection is done via matched filter. However, mextractor can perform cross-correlation with several filters. For example, the stellar PSF as well as some exdended PSFs for galaxies identification. This allows the user to find both faint stars and galaxies in a single mextractor call.

Measurments - mextractor performs a wide range of measurments on the detected sources. This include, positional information, moments, flag propogation, aperture photometry, PSF photometry, contour photometry, S/N measurments, WCS information, neighbour searches, as well as time-position dependent parameters like, airmass, and parallactic angle, azimuth and altitude. See help for complete list of parameters.

Source flagging - mextractor can remove and/or flag problematic sources. For example, it allows to locate sources on top of diffraction spikes, cosmic rays (via Hypothesis testing), sources in location of background with strong gradients, and more.

Examples

Simple use

mextractor is a method of the SIM class, and therefore it works on SIM objects. The first stage is to read an image(s) into a SIM container. This is done using functions in the FITS static class:

S=FITS.read2sim('test_PTF*100037_c02.fits');

S(1) % display the content of the first element in the SIM object

ans =
  SIM with properties:

               Im: [4096×2048 single]
    ImageFileName: ''
           BackIm: []
            ErrIm: []
         WeightIm: []
              Cat: []
              Col: []
          ColCell: []
         ColUnits: []
         SortedBy: []
      SortedByCol: []
             Name: []
           Source: []
        Reference: []
          Version: []
           Header: {376×3 cell}
              WCS: [1×1 struct]
         UserData: []
             Mask: []
          MaskDic: []
              PSF: []
           ErrPSF: []
           ParPSF: {}
           CooPSF: []

The variable S is a SIM array of one ore more images. In this example, only the Im, Header and WCS fields is populated.

Next, we can run mextractor, with the same output as input:

S=mextractor(S);

  Set gain to 1
mextractor first pass - Constructing optimal PSF
  Set gain to 1
  No filter found - Use DefFilter
  Estimate images background and std

Exiting: Maximum number of iterations has been exceeded
         - increase MaxIter option.
         Current function value: 212933.370859
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)

Exiting: Maximum number of iterations has been exceeded
         - increase MaxIter option.
         Current function value: 922844.488936
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)

Exiting: Maximum number of iterations has been exceeded
         - increase MaxIter option.
         Current function value: 144980.484556
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
  Subtract background
  Filter image number 1 - use 1 filters
    Threshold image and locate local maxima
    Extracting 2599 sources
  Calculate 1st and 2nd moments
Warning: There are some distance smaller than 1 - check why
mextractor second pass - Source detection and measurments
  Subtract background
Warning: CleanByBackGrad was set to false as the user did not requested for the SN, PEAK_GRADBACK, BACK_STD columns
  Filter image number 1 - use 2 filters
    Threshold image and locate local maxima
    Extracting 5551 sources
  Calculate 1st and 2nd moments
    Calculate PSF photometry
    Calculate IL aperture photometry (using aper_phot)

S(1)

ans =
  SIM with properties:

               Im: [4096×2048 single]
    ImageFileName: ''
           BackIm: [4096×2048 double]
            ErrIm: [4096×2048 double]
         WeightIm: []
              Cat: [5551×37 double]
              Col: [1×1 struct]
          ColCell: {1×37 cell}
         ColUnits: []
         SortedBy: {'YWIN_IMAGE'}
      SortedByCol: 2
             Name: []
           Source: []
        Reference: []
          Version: []
           Header: {2×3 cell}
              WCS: []
         UserData: []
             Mask: [4096×2048 uint32]
          MaskDic: []
              PSF: [21×21 double]
           ErrPSF: []
           ParPSF: {}
           CooPSF: []

After running mextractor, several other fields in S are populated. Specifically, the BackIm field contains the background image, the ErrIm field contains the background noise image, the PSF field contains the image point spread function, and Cat contains the catalog of sources detected and measured by mextractor.

For example, the PSF is stored in the PSF field and it can be displayed using:

surface(S(1).PSF)

The following command will display the first 10 lines in the catalog:

S(1).Cat(1:10,:)

ans =
       867.32       2.6424          867            3       1.3096       2.2408     0.026954          NaN          NaN       140.61       1773.4        10.64        3.532       9.8244       14.678       447.02       301.41       15.374      0.73209       125.95       177.23       166.08         1114       2136.6       2223.4       2342.9       3359.8       150.76       265.15        476.1       681.57       885.87       1632.8        39.81       1662.6         45.5            0
       927.97       3.2365          928            3       1.3525       1.8169     0.042438          NaN          NaN       151.21       1778.7        6.744       3.8042       10.932       9.3389       473.35       224.77       15.312      0.51557       35.143       66.033       53.178       593.42       1421.5       2716.5       4906.3       7164.6        136.6       271.46       485.74       694.88       897.44       1627.5       39.748       1642.5       43.297            0
       830.38        3.724          830            4       1.5504       1.3843     0.059015          NaN          NaN       96.036       1730.5       7.8957       2.4117       18.058       11.233       754.76       285.14       14.805      0.41018       36.783       87.997       82.668       760.11       823.24       368.41      -144.98      -406.47       149.75       281.34       502.01       710.51       917.17       1634.5        39.82       1666.5       41.797            0
       927.92       5.5912          927            7       1.3712       2.5147   -0.0059092          NaN          NaN       74.424         1702       3.9299       1.8723       2.4256       10.485       105.56       199.88       16.941       2.0559       39.362       45.291        44.34       517.46       1280.7       3401.7       5170.1       7301.5        147.9       291.66       537.83       756.28       966.38       1627.6       39.749       1642.9       43.519            0
       866.93       5.5929          867            6       1.4138       2.1344     -0.16764          NaN          NaN       184.61       1817.4       9.7667       4.6371       5.9216         17.3       260.62       283.46        15.96       1.1809       126.98       155.53       136.78       987.38         2488       2757.9       3211.3         3386       150.36       295.48       540.51       758.52       969.41       1632.8        39.81       1663.2       44.012            0
       1775.9       6.2621         1776            6      0.90116       1.1351    -0.083165          NaN          NaN       338.83       1968.5       17.469       8.7066        47.86       13.597       1963.6       230.29       13.767      0.12733       35.986       276.78       218.46       2164.2       3474.5       4019.8       3732.2        22539       147.31       294.12        549.5          768       984.09       1629.6       38.917         1628       41.027            0
       1635.1       6.5536         1635            7      0.96377       1.1036     0.026744          NaN          NaN       224.99       1852.5       12.699       5.8274       36.826       10.215       1428.5       279.56       14.113      0.21248       34.408       168.03       140.71       1840.2       2558.6         2891       3432.2       3468.7       156.89       289.66       551.23       771.75       982.31       1627.5       38.609       1626.8       38.791            0
       541.86       7.1964          542            7        0.653      0.74542     0.036833          NaN          NaN       142.31       1760.6       7.5909       3.6684       28.112       4.1977       1117.7       233.74       14.379      0.22705        16.23       67.923        56.42       839.01       1023.3        548.4       553.39      -932.82       142.32       286.22       559.29       785.48       991.35       1618.3       38.793       1617.9        39.76            0
       812.45        7.612          813            8       1.6701       2.1572     -0.15921          NaN          NaN      -6.5472         1627       2.8674     -0.16504      0.82155       8.7736       131.92       219.72       16.699       1.8084       23.679       28.496       27.043       220.76       553.35       2410.9         2609       1502.6       141.58       294.03       576.88       803.24         1016       1633.5       39.669       1651.9       160.57            0
       887.77       7.9181          888            8       1.8026       1.2515    -0.035343          NaN          NaN       96.349       1727.3       5.0589       2.4215       10.498       9.8339       508.98        251.8       15.233      0.53712       21.865       42.775       37.401        169.4      -131.41      -841.68      -551.33      -33.558          142        288.4       572.18       808.85       1023.7         1631       39.789       1666.1       48.484            0

The catalog column names are stored as a cell array in:

S(1).ColCell

ans =
    'XWIN_IMAGE'    'YWIN_IMAGE'    'XPEAK_IMAGE'    'YPEAK_IMAGE'    'X2WIN_IMAGE'    'Y2WIN_IMAGE'    'XYWIN_IMAGE'    'ALPHAWIN_J2000'    'DELTAWIN_J2000'    'PEAKF_VAL'    'PEAKF_VALTOT'    'SN'    'SN_UNF'    'SN_PSF'    'SN_ADD_1'    'FLUX_PSF'    'FLUXERR_PSF'    'MAG_PSF'    'MAGERR_PSF'    'PSF_CHI2'    'PSF_CHI2BACK'    'PSF_CHI2CR'    'FLUX_APER_1'    'FLUX_APER_2'    'FLUX_APER_3'    'FLUX_APER_4'    'FLUX_APER_5'    'FLUXERR_APER_1'    'FLUXERR_APER_2'    'FLUXERR_APER_3'    'FLUXERR_APER_4'    'FLUXERR_APER_5'    'BACK'    'BACK_STD'    'BACK_ANNULUS'    'STD_ANNULUS'    'FLAGS'

or as a structure array in:

S(1).Col

ans =
        XWIN_IMAGE: 1
        YWIN_IMAGE: 2
       XPEAK_IMAGE: 3
       YPEAK_IMAGE: 4
       X2WIN_IMAGE: 5
       Y2WIN_IMAGE: 6
       XYWIN_IMAGE: 7
    ALPHAWIN_J2000: 8
    DELTAWIN_J2000: 9
         PEAKF_VAL: 10
      PEAKF_VALTOT: 11
                SN: 12
            SN_UNF: 13
            SN_PSF: 14
          SN_ADD_1: 15
          FLUX_PSF: 16
       FLUXERR_PSF: 17
           MAG_PSF: 18
        MAGERR_PSF: 19
          PSF_CHI2: 20
      PSF_CHI2BACK: 21
        PSF_CHI2CR: 22
       FLUX_APER_1: 23
       FLUX_APER_2: 24
       FLUX_APER_3: 25
       FLUX_APER_4: 26
       FLUX_APER_5: 27
    FLUXERR_APER_1: 28
    FLUXERR_APER_2: 29
    FLUXERR_APER_3: 30
    FLUXERR_APER_4: 31
    FLUXERR_APER_5: 32
              BACK: 33
          BACK_STD: 34
      BACK_ANNULUS: 35
       STD_ANNULUS: 36
             FLAGS: 37

To display the image with the positions of the extracted sources.

ds9(S(1))

eran     34278  2.9  0.0 180076 77568 pts/39   S    16:06   0:03 ds9
ds9 is already open

ds9.plot(S(1))

We can manipulate and plot columns of the catalog. For example, to display the detection S/N vs. PSF magnitude:

Col = S(1).Col;

plot(S(1).Cat(:,Col.MAG_PSF),S(1).Cat(:,Col.SN),'.')

As expected the S/N is larger for bright stars. However, at the top of the diagram, there are several outliers (with S/N>6000 and Mag>10). You can choose this sources and plot their position on the image:

Col = S(1).Col;

Flag = S(1).Cat(:,Col.MAG_PSF)>10 & S(1).Cat(:,Col.SN)>6000;

ds9(S(1))

eran     34278  3.5  0.0 180076 77568 pts/39   S    16:06   0:05 ds9
ds9 is already open

ds9.plot( S(1).Cat(Flag,[Col.XWIN_IMAGE, Col.YWIN_IMAGE]),'bs')

Now, the nature of these sources is clear - these are saturated objects.

Input parameters

You can change the settings/input parameters of mextractor either uby adding a pairs of keyword,value pairs:

S=mextractor(S,'ZP',27,'MagInLuptitude',false,'Thresh',8);

A complete list of keywords is available in the help of mextractor.

Alternatively you can use the InArg.epar command to edit the arguments and save them as your default parameters:

InArg.epar('mextractor')

Control the output catalog columns

The ColCell input argument controls the mextractor catalog output colums. For example, to produce a catalog with only two columns, the X and Y positions:

S=mextractor(S,'ColCell',{'XWIN_IMAGE','YWIN_IMAGE'});

List of catalog output columns (source attributes)

A list of all possible output columns (source attributes) are listed below. Some attributes may have indexed multiple entries (marked by [_#]). Examples, are aperture photometry in multiple apertures.

Measurements in forced positions

mextractor can both find sources, and/or use a user a list of "forced position" supplied by the user.

The following input arguments are relevant for forced position and measurments:

For example, the following command will perform measurment in specified positions:

S=mextractor(S,'ForcePos',[512 512; 600 600],'OnlyForce',true);

Note that the output column name 'ISFORCE' indicate if the source originate from the ForcePos list, or found by mextractor.

Magnitude measurmemnts

mextractor supports several type of flux/magnitude estimate for sources, including, aperture photometry, PSF photometry and isophotal photometry.

Note that the PSF photometry is implemented in SIM/psf_phot, while the aperture photometry is in SIM/aper_phot. The magnitude output can be either im magnitude or luptitude units ('MagInLuptitude' and 'LuptSoft' input arguments), and the magnitude zero point is controlled via the 'ZP' input argument.

Isophotal flux and errors are stored in the FLUX_ISO and FLUXERR_ISO output columns, respectively. The user can specify the isophotal level (in units of the local background noise sigma) using the input argument 'PropThresh' (default is 1.5).

Background estimation

The image background is crucial for source detection and photometry. Small biases in the background can lead to false detection of sources and mis-detection of other sources.

The background image field is populated using the SIM/background method.

The background associated with each source is stored in the output columns. The user can get both the background estimated by the SIM/background method, and the annulys background estimated by SIM/aper_phot.

Filtering

Cosmic rays detection

Bit masks and flags

mextractor use and populate the MASK image. First it can flag specific bits in the MASK image. For example, it looks for diffraction spikes, and mark pixels that are contaminated by diffraction spikes in the MASK image. Second it propagate the MASK image information into the output catalog. It collects all the bits that are found within FlagRadius argument from a source, and propagate them (using or operation) to the catalog FLAGS.

We note that the default mask bit definition of mextractor is in: @MASK.def_bitmask_pipeline. You can use the MASK method val2bitname to convert bit mask value to name.

You can get a logical flag image of pixels that have some bit masks on. The following example return a SIM object in which the Im field is a logical flag indicating if the MASK bits are either 'Bit_Spike 'or 'Bit_ImSaturated':

Im=bitmask_find(S,{'Bit_Spike','Bit_ImSaturated'});

Comparison with SExtractor

To compare mextractor and SExtractor (REF), you can use the matlab sextractor function that calls SExtractor.

Ss=sextractor(S(1),'SexPars',{'DETECT_THRESH','5'});

[1M>
[1A----- SExtractor 2.5.0 started on 2017-01-15 at 09:51:51 with 1 thread

[1M> Setting catalog parameters
[1A[1M> Reading detection filter
[1A[1M> Initializing catalog
[1A[1M> Looking for tp609482ef_26c6_4ccf_bb5f_4e9c0ea41c1d
[1AMeasuring from: "PTF_survey"  / 2048 x 4096 / 32 bits FLOATING POINT data
Detection+Measurement image: [1M> Setting up background maps
[1A[1M> Setting up background map at line:  128
[1A[1M> Setting up background map at line:  256
[1A[1M> Setting up background map at line:  384
[1A[1M> Setting up background map at line:  512
[1A[1M> Setting up background map at line:  640
[1A[1M> Setting up background map at line:  768
[1A[1M> Setting up background map at line:  896
[1A[1M> Setting up background map at line: 1024
[1A[1M> Setting up background map at line: 1152
[1A[1M> Setting up background map at line: 1280
[1A[1M> Setting up background map at line: 1408
[1A[1M> Setting up background map at line: 1536
[1A[1M> Setting up background map at line: 1664
[1A[1M> Setting up background map at line: 1792
[1A[1M> Setting up background map at line: 1920
[1A[1M> Setting up background map at line: 2048
[1A[1M> Setting up background map at line: 2176
[1A[1M> Setting up background map at line: 2304
[1A[1M> Setting up background map at line: 2432
[1A[1M> Setting up background map at line: 2560
[1A[1M> Setting up background map at line: 2688
[1A[1M> Setting up background map at line: 2816
[1A[1M> Setting up background map at line: 2944
[1A[1M> Setting up background map at line: 3072
[1A[1M> Setting up background map at line: 3200
[1A[1M> Setting up background map at line: 3328
[1A[1M> Setting up background map at line: 3456
[1A[1M> Setting up background map at line: 3584
[1A[1M> Setting up background map at line: 3712
[1A[1M> Setting up background map at line: 3840
[1A[1M> Setting up background map at line: 3968
[1A[1M> Filtering background map(s)
[1A[1M> Computing background d-map
[1A[1M> Computing background-noise d-map
[1A(M+D) Background: 802.121    RMS: 28.312     / Threshold: 141.56    
[1M> Scanning image
[1A[1M> Line:   16  Objects:       14 detected /        0 sextracted
[1A[1M> Line:   32  Objects:       33 detected /        0 sextracted
[1A[1M> Line:   48  Objects:       51 detected /        0 sextracted
[1A[1M> Line:   64  Objects:       61 detected /        0 sextracted
[1A[1M> Line:   80  Objects:       77 detected /        0 sextracted
[1A[1M> Line:   96  Objects:       94 detected /        0 sextracted
[1A[1M> Line:  112  Objects:      112 detected /        0 sextracted
[1A[1M> Line:  128  Objects:      125 detected /        0 sextracted
[1A[1M> Line:  144  Objects:      147 detected /        0 sextracted
[1A[1M> Line:  160  Objects:      163 detected /        0 sextracted
[1A[1M> Line:  176  Objects:      176 detected /        0 sextracted
[1A[1M> Line:  192  Objects:      199 detected /        0 sextracted
[1A[1M> Line:  208  Objects:      219 detected /        0 sextracted
[1A[1M> Line:  224  Objects:      243 detected /        0 sextracted
[1A[1M> Line:  240  Objects:      260 detected /        0 sextracted
[1A[1M> Line:  256  Objects:      273 detected /        0 sextracted
[1A[1M> Line:  272  Objects:      290 detected /        0 sextracted
[1A[1M> Line:  288  Objects:      302 detected /        0 sextracted
[1A[1M> Line:  304  Objects:      318 detected /        0 sextracted
[1A[1M> Line:  320  Objects:      339 detected /        0 sextracted
[1A[1M> Line:  336  Objects:      348 detected /        0 sextracted
[1A[1M> Line:  352  Objects:      365 detected /        0 sextracted
[1A[1M> Line:  368  Objects:      384 detected /        0 sextracted
[1A[1M> Line:  384  Objects:      397 detected /        0 sextracted
[1A[1M> Line:  400  Objects:      417 detected /        0 sextracted
[1A[1M> Line:  416  Objects:      432 detected /        0 sextracted
[1A[1M> Line:  432  Objects:      451 detected /        0 sextracted
[1A[1M> Line:  448  Objects:      478 detected /        0 sextracted
[1A[1M> Line:  464  Objects:      491 detected /        0 sextracted
[1A[1M> Line:  480  Objects:      504 detected /        0 sextracted
[1A[1M> Line:  496  Objects:      522 detected /        0 sextracted
[1A[1M> Line:  512  Objects:      534 detected /        0 sextracted
[1A[1M> Line:  528  Objects:      554 detected /        0 sextracted
[1A[1M> Line:  544  Objects:      567 detected /        0 sextracted
[1A[1M> Line:  560  Objects:      580 detected /        0 sextracted
[1A[1M> Line:  576  Objects:      589 detected /        0 sextracted
[1A[1M> Line:  592  Objects:      603 detected /        0 sextracted
[1A[1M> Line:  608  Objects:      615 detected /        0 sextracted
[1A[1M> Line:  624  Objects:      629 detected /        0 sextracted
[1A[1M> Line:  640  Objects:      644 detected /        0 sextracted
[1A[1M> Line:  656  Objects:      661 detected /        0 sextracted
[1A[1M> Line:  672  Objects:      678 detected /        0 sextracted
[1A[1M> Line:  688  Objects:      694 detected /        0 sextracted
[1A[1M> Line:  704  Objects:      713 detected /        0 sextracted
[1A[1M> Line:  720  Objects:      723 detected /        0 sextracted
[1A[1M> Line:  736  Objects:      745 detected /        0 sextracted
[1A[1M> Line:  752  Objects:      767 detected /        0 sextracted
[1A[1M> Line:  768  Objects:      782 detected /        0 sextracted
[1A[1M> Line:  784  Objects:      804 detected /        0 sextracted
[1A[1M> Line:  800  Objects:      820 detected /        0 sextracted
[1A[1M> Line:  816  Objects:      833 detected /        0 sextracted
[1A[1M> Line:  832  Objects:      848 detected /        0 sextracted
[1A[1M> Line:  848  Objects:      856 detected /        0 sextracted
[1A[1M> Line:  864  Objects:      869 detected /        0 sextracted
[1A[1M> Line:  880  Objects:      882 detected /        0 sextracted
[1A[1M> Line:  896  Objects:      901 detected /        0 sextracted
[1A[1M> Line:  912  Objects:      913 detected /        0 sextracted
[1A[1M> Line:  928  Objects:      929 detected /        0 sextracted
[1A[1M> Line:  944  Objects:      950 detected /        0 sextracted
[1A[1M> Line:  960  Objects:      966 detected /        0 sextracted
[1A[1M> Line:  976  Objects:      986 detected /        0 sextracted
[1A[1M> Line:  992  Objects:      994 detected /        0 sextracted
[1A[1M> Line: 1008  Objects:     1010 detected /        0 sextracted
[1A[1M> Line: 1024  Objects:     1028 detected /       14 sextracted
[1A[1M> Line: 1040  Objects:     1042 detected /       32 sextracted
[1A[1M> Line: 1056  Objects:     1061 detected /       51 sextracted
[1A[1M> Line: 1072  Objects:     1076 detected /       61 sextracted
[1A[1M> Line: 1088  Objects:     1095 detected /       70 sextracted
[1A[1M> Line: 1104  Objects:     1114 detected /       86 sextracted
[1A[1M> Line: 1120  Objects:     1133 detected /      102 sextracted
[1A[1M> Line: 1136  Objects:     1148 detected /      112 sextracted
[1A[1M> Line: 1152  Objects:     1160 detected /      131 sextracted
[1A[1M> Line: 1168  Objects:     1177 detected /      145 sextracted
[1A[1M> Line: 1184  Objects:     1190 detected /      163 sextracted
[1A[1M> Line: 1200  Objects:     1205 detected /      178 sextracted
[1A[1M> Line: 1216  Objects:     1230 detected /      199 sextracted
[1A[1M> Line: 1232  Objects:     1240 detected /      213 sextracted
[1A[1M> Line: 1248  Objects:     1251 detected /      228 sextracted
[1A[1M> Line: 1264  Objects:     1271 detected /      241 sextracted
[1A[1M> Line: 1280  Objects:     1286 detected /      251 sextracted
[1A[1M> Line: 1296  Objects:     1307 detected /      269 sextracted
[1A[1M> Line: 1312  Objects:     1320 detected /      278 sextracted
[1A[1M> Line: 1328  Objects:     1339 detected /      293 sextracted
[1A[1M> Line: 1344  Objects:     1361 detected /      308 sextracted
[1A[1M> Line: 1360  Objects:     1379 detected /      322 sextracted
[1A[1M> Line: 1376  Objects:     1394 detected /      336 sextracted
[1A[1M> Line: 1392  Objects:     1411 detected /      349 sextracted
[1A[1M> Line: 1408  Objects:     1428 detected /      363 sextracted
[1A[1M> Line: 1424  Objects:     1440 detected /      380 sextracted
[1A[1M> Line: 1440  Objects:     1452 detected /      396 sextracted
[1A[1M> Line: 1456  Objects:     1475 detected /      414 sextracted
[1A[1M> Line: 1472  Objects:     1493 detected /      424 sextracted
[1A[1M> Line: 1488  Objects:     1511 detected /      441 sextracted
[1A[1M> Line: 1504  Objects:     1527 detected /      455 sextracted
[1A[1M> Line: 1520  Objects:     1541 detected /      461 sextracted
[1A[1M> Line: 1536  Objects:     1558 detected /      477 sextracted
[1A[1M> Line: 1552  Objects:     1578 detected /      490 sextracted
[1A[1M> Line: 1568  Objects:     1591 detected /      504 sextracted
[1A[1M> Line: 1584  Objects:     1607 detected /      512 sextracted
[1A[1M> Line: 1600  Objects:     1627 detected /      520 sextracted
[1A[1M> Line: 1616  Objects:     1638 detected /      528 sextracted
[1A[1M> Line: 1632  Objects:     1663 detected /      543 sextracted
[1A[1M> Line: 1648  Objects:     1674 detected /      556 sextracted
[1A[1M> Line: 1664  Objects:     1685 detected /      571 sextracted
[1A[1M> Line: 1680  Objects:     1697 detected /      590 sextracted
[1A[1M> Line: 1696  Objects:     1713 detected /      599 sextracted
[1A[1M> Line: 1712  Objects:     1732 detected /      616 sextracted
[1A[1M> Line: 1728  Objects:     1745 detected /      627 sextracted
[1A[1M> Line: 1744  Objects:     1761 detected /      651 sextracted
[1A[1M> Line: 1760  Objects:     1780 detected /      669 sextracted

Source detection comparison

We note that in order to have a fair comparison we have to run SExtractor with the same PSF we used for mextractor. The sextractor function takes the PSF field from the SIM object and send it to SExtractor. In this example, we run SExtractor with detection threshold of 5sigma.

Now we can display the results

ds9(S(1),1)

eran     21137  0.7  0.0 217752 116412 pts/39  S    09:36   0:11 ds9
ds9 is already open

ds9.tile([2 1])

ds9.plot(S(1))

ds9(S(1),2)

eran     21137  0.7  0.0 217752 116412 pts/39  S    09:36   0:12 ds9
ds9 is already open

ds9.plot(Ss)

This example demonstrate that SExtractor did not find obvious sources in the image. The main reason for this is a bug in the definition of the detection threshold in SExtractor. This bug can be easily fixed by renormalizing the SExtractor sigma definition:

% Normalization for correcting the SExtractor bug

  Set gain to 1
mextractor first pass - Constructing optimal PSF
  Set gain to 1
  No filter found - Use DefFilter
  Estimate images background and std

Exiting: Maximum number of iterations has been exceeded
         - increase MaxIter option.
         Current function value: 212933.370859
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)

Exiting: Maximum number of iterations has been exceeded
         - increase MaxIter option.
         Current function value: 922844.488936
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)

Exiting: Maximum number of iterations has been exceeded
         - increase MaxIter option.
         Current function value: 144980.484556
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
  Subtract background
  Filter image number 1 - use 1 filters
    Threshold image and locate local maxima
    Extracting 2599 sources
  Calculate 1st and 2nd moments
Warning: There are some distance smaller than 1 - check why
mextractor second pass - Source detection and measurments
  Subtract background
Warning: CleanByBackGrad was set to false as the user did not requested for the SN, PEAK_GRADBACK, BACK_STD columns
  Filter image number 1 - use 2 filters
    Threshold image and locate local maxima
    Extracting 5551 sources
  Calculate 1st and 2nd moments
    Calculate PSF photometry
    Calculate IL aperture photometry (using aper_phot)

Norm = sqrt(sum(S.PSF(:).^2))

Norm =
      0.27729

Ss=sextractor(S,'SexPars',{'DETECT_THRESH',5.*Norm});

[1M>
[1A----- SExtractor 2.5.0 started on 2017-02-20 at 16:59:15 with 1 thread

[1M> Setting catalog parameters
[1A[1M> Reading detection filter
[1A[1M> Initializing catalog
[1A[1M> Looking for tp347a2898_7e21_471b_9433_72cbbd586835.fits
[1AMeasuring from: "Unnamed"  / 2048 x 4096 / 32 bits FLOATING POINT data
Detection+Measurement image: [1M> Setting up background maps
[1A[1M> Setting up background map at line:  128
[1A[1M> Setting up background map at line:  256
[1A[1M> Setting up background map at line:  384
[1A[1M> Setting up background map at line:  512
[1A[1M> Setting up background map at line:  640
[1A[1M> Setting up background map at line:  768
[1A[1M> Setting up background map at line:  896
[1A[1M> Setting up background map at line: 1024
[1A[1M> Setting up background map at line: 1152
[1A[1M> Setting up background map at line: 1280
[1A[1M> Setting up background map at line: 1408
[1A[1M> Setting up background map at line: 1536
[1A[1M> Setting up background map at line: 1664
[1A[1M> Setting up background map at line: 1792
[1A[1M> Setting up background map at line: 1920
[1A[1M> Setting up background map at line: 2048
[1A[1M> Setting up background map at line: 2176
[1A[1M> Setting up background map at line: 2304
[1A[1M> Setting up background map at line: 2432
[1A[1M> Setting up background map at line: 2560
[1A[1M> Setting up background map at line: 2688
[1A[1M> Setting up background map at line: 2816
[1A[1M> Setting up background map at line: 2944
[1A[1M> Setting up background map at line: 3072
[1A[1M> Setting up background map at line: 3200
[1A[1M> Setting up background map at line: 3328
[1A[1M> Setting up background map at line: 3456
[1A[1M> Setting up background map at line: 3584
[1A[1M> Setting up background map at line: 3712
[1A[1M> Setting up background map at line: 3840
[1A[1M> Setting up background map at line: 3968
[1A[1M> Filtering background map(s)
[1A[1M> Computing background d-map
[1A[1M> Computing background-noise d-map
[1A(M+D) Background: 1629.56    RMS: 37.7916    / Threshold: 52.398    
[1M> Scanning image
[1A[1M> Line:   16  Objects:       57 detected /        0 sextracted
[1A[1M> Line:   32  Objects:      106 detected /        0 sextracted
[1A[1M> Line:   48  Objects:      135 detected /        0 sextracted
[1A[1M> Line:   64  Objects:      162 detected /        0 sextracted
[1A[1M> Line:   80  Objects:      180 detected /        0 sextracted
[1A[1M> Line:   96  Objects:      196 detected /        0 sextracted
[1A[1M> Line:  112  Objects:      214 detected /        0 sextracted
[1A[1M> Line:  128  Objects:      233 detected /        0 sextracted
[1A[1M> Line:  144  Objects:      251 detected /        0 sextracted
[1A[1M> Line:  160  Objects:      268 detected /        0 sextracted
[1A[1M> Line:  176  Objects:      296 detected /        0 sextracted
[1A[1M> Line:  192  Objects:      319 detected /        0 sextracted
[1A[1M> Line:  208  Objects:      349 detected /        0 sextracted
[1A[1M> Line:  224  Objects:      371 detected /        0 sextracted
[1A[1M> Line:  240  Objects:      397 detected /        0 sextracted
[1A[1M> Line:  256  Objects:      424 detected /        0 sextracted
[1A[1M> Line:  272  Objects:      441 detected /        0 sextracted
[1A[1M> Line:  288  Objects:      468 detected /        0 sextracted
[1A[1M> Line:  304  Objects:      493 detected /        0 sextracted
[1A[1M> Line:  320  Objects:      516 detected /        0 sextracted
[1A[1M> Line:  336  Objects:      543 detected /        0 sextracted
[1A[1M> Line:  352  Objects:      567 detected /        0 sextracted
[1A[1M> Line:  368  Objects:      581 detected /        0 sextracted
[1A[1M> Line:  384  Objects:      612 detected /        0 sextracted
[1A[1M> Line:  400  Objects:      631 detected /        0 sextracted
[1A[1M> Line:  416  Objects:      653 detected /        0 sextracted
[1A[1M> Line:  432  Objects:      678 detected /        0 sextracted
[1A[1M> Line:  448  Objects:      703 detected /        0 sextracted
[1A[1M> Line:  464  Objects:      718 detected /        0 sextracted
[1A[1M> Line:  480  Objects:      734 detected /        0 sextracted
[1A[1M> Line:  496  Objects:      756 detected /        0 sextracted
[1A[1M> Line:  512  Objects:      784 detected /        0 sextracted
[1A[1M> Line:  528  Objects:      827 detected /        0 sextracted
[1A[1M> Line:  544  Objects:      848 detected /        0 sextracted
[1A[1M> Line:  560  Objects:      875 detected /        0 sextracted
[1A[1M> Line:  576  Objects:      910 detected /        0 sextracted
[1A[1M> Line:  592  Objects:      939 detected /        0 sextracted
[1A[1M> Line:  608  Objects:      969 detected /        0 sextracted
[1A[1M> Line:  624  Objects:     1005 detected /        0 sextracted
[1A[1M> Line:  640  Objects:     1030 detected /        0 sextracted
[1A[1M> Line:  656  Objects:     1061 detected /        0 sextracted
[1A[1M> Line:  672  Objects:     1085 detected /        0 sextracted
[1A[1M> Line:  688  Objects:     1115 detected /        0 sextracted
[1A[1M> Line:  704  Objects:     1135 detected /        0 sextracted
[1A[1M> Line:  720  Objects:     1157 detected /        0 sextracted
[1A[1M> Line:  736  Objects:     1190 detected /        0 sextracted
[1A[1M> Line:  752  Objects:     1217 detected /        0 sextracted
[1A[1M> Line:  768  Objects:     1240 detected /        0 sextracted
[1A[1M> Line:  784  Objects:     1270 detected /        0 sextracted
[1A[1M> Line:  800  Objects:     1299 detected /        0 sextracted
[1A[1M> Line:  816  Objects:     1336 detected /        0 sextracted
[1A[1M> Line:  832  Objects:     1371 detected /        0 sextracted
[1A[1M> Line:  848  Objects:     1406 detected /        0 sextracted
[1A[1M> Line:  864  Objects:     1438 detected /        0 sextracted
[1A[1M> Line:  880  Objects:     1467 detected /        0 sextracted
[1A[1M> Line:  896  Objects:     1508 detected /        0 sextracted
[1A[1M> Line:  912  Objects:     1543 detected /        0 sextracted
[1A[1M> Line:  928  Objects:     1585 detected /        0 sextracted
[1A[1M> Line:  944  Objects:     1630 detected /        0 sextracted
[1A[1M> Line:  960  Objects:     1656 detected /        0 sextracted
[1A[1M> Line:  976  Objects:     1694 detected /        0 sextracted
[1A[1M> Line:  992  Objects:     1733 detected /        0 sextracted
[1A[1M> Line: 1008  Objects:     1768 detected /        0 sextracted
[1A[1M> Line: 1024  Objects:     1800 detected /       30 sextracted
[1A[1M> Line: 1040  Objects:     1866 detected /       59 sextracted
[1A[1M> Line: 1056  Objects:     1960 detected /       85 sextracted
[1A[1M> Line: 1072  Objects:     2060 detected /      107 sextracted
[1A[1M> Line: 1088  Objects:     2152 detected /      127 sextracted
[1A[1M> Line: 1104  Objects:     2193 detected /      145 sextracted
[1A[1M> Line: 1120  Objects:     2240 detected /      159 sextracted
[1A[1M> Line: 1136  Objects:     2277 detected /      175 sextracted
[1A[1M> Line: 1152  Objects:     2320 detected /      190 sextracted
[1A[1M> Line: 1168  Objects:     2365 detected /      211 sextracted
[1A[1M> Line: 1184  Objects:     2401 detected /      228 sextracted
[1A[1M> Line: 1200  Objects:     2437 detected /      251 sextracted
[1A[1M> Line: 1216  Objects:     2471 detected /      273 sextracted
[1A[1M> Line: 1232  Objects:     2510 detected /      288 sextracted
[1A[1M> Line: 1248  Objects:     2557 detected /      320 sextracted
[1A[1M> Line: 1264  Objects:     2594 detected /      344 sextracted
[1A[1M> Line: 1280  Objects:     2638 detected /      363 sextracted
[1A[1M> Line: 1296  Objects:     2698 detected /      376 sextracted
[1A[1M> Line: 1312  Objects:     2727 detected /      404 sextracted
[1A[1M> Line: 1328  Objects:     2769 detected /      432 sextracted
[1A[1M> Line: 1344  Objects:     2821 detected /      447 sextracted
[1A[1M> Line: 1360  Objects:     2878 detected /      467 sextracted
[1A[1M> Line: 1376  Objects:     2939 detected /      485 sextracted
[1A[1M> Line: 1392  Objects:     3014 detected /      511 sextracted
[1A[1M> Line: 1408  Objects:     3129 detected /      530 sextracted
[1A[1M> Line: 1424  Objects:     3262 detected /      548 sextracted
[1A[1M> Line: 1440  Objects:     3374 detected /      570 sextracted
[1A[1M> Line: 1456  Objects:     3442 detected /      589 sextracted
[1A[1M> Line: 1472  Objects:     3482 detected /      604 sextracted
[1A[1M> Line: 1488  Objects:     3523 detected /      619 sextracted
[1A[1M> Line: 1504  Objects:     3551 detected /      636 sextracted
[1A[1M> Line: 1520  Objects:     3586 detected /      660 sextracted
[1A[1M> Line: 1536  Objects:     3616 detected /      682 sextracted
[1A[1M> Line: 1552  Objects:     3641 detected /      705 sextracted
[1A[1M> Line: 1568  Objects:     3663 detected /      728 sextracted
[1A[1M> Line: 1584  Objects:     3688 detected /      756 sextracted
[1A[1M> Line: 1600  Objects:     3711 detected /      780 sextracted
[1A[1M> Line: 1616  Objects:     3730 detected /      798 sextracted
[1A[1M> Line: 1632  Objects:     3753 detected /      816 sextracted
[1A[1M> Line: 1648  Objects:     3777 detected /      835 sextracted
[1A[1M> Line: 1664  Objects:     3798 detected /      855 sextracted
[1A[1M> Line: 1680  Objects:     3815 detected /      873 sextracted
[1A[1M> Line: 1696  Objects:     3837 detected /      900 sextracted
[1A[1M> Line: 1712  Objects:     3854 detected /      914 sextracted
[1A[1M> Line: 1728  Objects:     3871 detected /      931 sextracted
[1A[1M> Line: 1744  Objects:     3884 detected /      951 sextracted
[1A[1M> Line: 1760  Objects:     3905 detected /      978 sextracted

S=mextractor(S(1)); % running mextractor on the same image

ds9(S,1)

eran     34278  0.3  0.0 179416 78660 pts/39   S    16:06   0:11 ds9
ds9 is already open

ds9.tile([2 1])

ds9.plot(S)

ds9(S,2)

ds9.plot(Ss)

ds9.lock_xy

In this example we locked the image coordinates of the two images. There are several important differences:

  1. It is possible to see that SExtractor missed some bright real sources, especially near bright sources. This is likely due to deblending issues - This can be fixed by changing the deblending arguments in SExtractor. However, there is not clear prescription on how to make it work for a wide range of conditions.
  2. SExtractor found some faint sources that are unreal (this can be verified by overlaying the SDSS sources). We speculate that this is due to inaccurate estimation of the background level and noise.
  3. If you will soom out and explore the region around M81, you will see that SExtractor found many fake sources in a ring around M81, and did not identify sources on top of the galaxy. This is likely due to the SExtractor background and noise estimation.
  4. Another important difference is that mextractor found the diffraction spikes of bright stars as multiple sources. In order to mitigate this problem, mextractor locate diffraction spikes using the find_diff_spikes function and flag this sources in two ways: First mextractor populate the MASK image field in the SIM object and flag diffraction spikes. Second, this mask bits are propogated into the catalog so sources on top of diffraction spikes can be removed.

The mextractor (right) and SExtractor (left) sources around M81:

Astrometry comparison

In order to compare the astrometry performences of mextractor with other tools we run it on two PTF images of the same field, matched the stars, register the images (always using the same procedure), and calculated the difference in X/Y position between the two images.

The results indicate that the astrometric precision of mextractor and SExtractor are identical to within the uncertanties.

% retrieve some PTF images using VO.PTF.wget_corrim

%VO.PTF.wget_corrim([100037 2],[],'JD',[celestial.time.julday([20 11 2012]), celestial.time.julday([23 11 2012])],'Filter','R','save','icp');

VO.PTF.wget_corrim([100037 2],[],'JD',[celestial.time.julday([20 11 2014]), celestial.time.julday([23 11 2014])],'Filter','R','save','icp');

PTF_201411205388_i_p_scie_t125555_u023050416_f02_p100037_c02.fits  


PTF_201411204943_i_p_scie_t115144_u023050379_f02_p100037_c02.fits  


PTF_201411225028_i_p_scie_t120406_u022881502_f02_p100037_c02.fits  


PTF_201411225468_i_p_scie_t130720_u022881475_f02_p100037_c02.fits  


PTF_201411205388_i_p_scie_t125555_u023050416_f02_p100037_c02.fits  


PTF_201411204943_i_p_scie_t115144_u023050379_f02_p100037_c02.fits  


PTF_201411225028_i_p_scie_t120406_u022881502_f02_p100037_c02.fits  


PTF_201411225468_i_p_scie_t130720_u022881475_f02_p100037_c02.fits  


PTF_201411205388_c_p_scie_t125555_u023050416_f02_p100037_c02.ctlg  


PTF_201411204943_c_p_scie_t115144_u023050379_f02_p100037_c02.ctlg  


PTF_201411225028_c_p_scie_t120406_u022881502_f02_p100037_c02.ctlg  


PTF_201411225468_c_p_scie_t130720_u022881475_f02_p100037_c02.ctlg  


PTF_201411205388_c_p_scie_t125555_u023050416_f02_p100037_c02.ctlg  


PTF_201411204943_c_p_scie_t115144_u023050379_f02_p100037_c02.ctlg  


PTF_201411225028_c_p_scie_t120406_u022881502_f02_p100037_c02.ctlg  


PTF_201411225468_c_p_scie_t130720_u022881475_f02_p100037_c02.ctlg  


PTF_201411205388_c_d_scie_t125555_u023050416_f02_p100037_c02.ctlg  


PTF_201411204943_c_d_scie_t115144_u023050379_f02_p100037_c02.ctlg  


PTF_201411225028_c_d_scie_t120406_u022881502_f02_p100037_c02.ctlg  


PTF_201411225468_c_d_scie_t130720_u022881475_f02_p100037_c02.ctlg  


PTF_201411205388_c_d_scie_t125555_u023050416_f02_p100037_c02.ctlg  


PTF_201411204943_c_d_scie_t115144_u023050379_f02_p100037_c02.ctlg  


PTF_201411225028_c_d_scie_t120406_u022881502_f02_p100037_c02.ctlg  


PTF_201411225468_c_d_scie_t130720_u022881475_f02_p100037_c02.ctlg  

% Read images to SIM arrays

%S1=FITS.read2sim('PTF_201211213383_i_p_scie_t080710_u014664919_f02_p100037_c02.fits');

%S2=FITS.read2sim('PTF_201211223616_i_p_scie_t084044_u014676172_f02_p100037_c02.fits');

S1=FITS.read2sim('PTF_201411204943_i_p_scie_t115144_u023050379_f02_p100037_c02.fits');

S2=FITS.read2sim('PTF_201411205388_i_p_scie_t125555_u023050416_f02_p100037_c02.fits');

% find stars

S1=mextractor(S1);

  Set gain to 1
mextractor first pass - Constructing optimal PSF
  Set gain to 1
  No filter found - Use DefFilter
  Estimate images background and std
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)

Exiting: Maximum number of iterations has been exceeded
         - increase MaxIter option.
         Current function value: 159696.344677
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
  Subtract background
mextractor image 1 out of 1
  Filter image number 1 - use 1 filters
    Threshold image and locate local maxima
    Extracting 2671 sources
  Calculate 1st and 2nd moments
Warning: There are some distance smaller than 1 - check why
mextractor second pass - Source detection and measurments
  Subtract background
Warning: CleanByBackGrad was set to false as the user did not requested for the SN, PEAK_GRADBACK, BACK_STD columns
mextractor image 1 out of 1
  Filter image number 1 - use 2 filters
    Threshold image and locate local maxima
    Extracting 5505 sources
  Calculate 1st and 2nd moments
    Calculate PSF photometry
    Calculate IL aperture photometry (using aper_phot)

S2=mextractor(S2);

  Set gain to 1
mextractor first pass - Constructing optimal PSF
  Set gain to 1
  No filter found - Use DefFilter
  Estimate images background and std
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)

Exiting: Maximum number of iterations has been exceeded
         - increase MaxIter option.
         Current function value: 148347.162682
Warning: Best fit mode is out of range - set mode to max value and std to sqrt(max)
  Subtract background
mextractor image 1 out of 1
  Filter image number 1 - use 1 filters
    Threshold image and locate local maxima
    Extracting 2547 sources
  Calculate 1st and 2nd moments
Warning: There are some distance smaller than 1 - check why
mextractor second pass - Source detection and measurments
  Subtract background
Warning: CleanByBackGrad was set to false as the user did not requested for the SN, PEAK_GRADBACK, BACK_STD columns
mextractor image 1 out of 1
  Filter image number 1 - use 2 filters
    Threshold image and locate local maxima
    Extracting 5353 sources
  Calculate 1st and 2nd moments
    Calculate PSF photometry
    Calculate IL aperture photometry (using aper_phot)

% align images

TranC = TranClass({@FunOne, [],@FunX,[],@FunY,[],@FunTiltXp,[],@FunTiltXn,[]}, {@FunOne, [],@FunX,[],@FunY,[],@FunTiltYp,[],@FunTiltYn,[] });

%ParDeg = [1 1;2 0;0 2;2 1;1 2];

%TranC = TranClass({@FunOne, [],@FunX,[],@FunY,[],@FunPolyChebyshev1XY,ParDeg}, {@FunOne, [],@FunX,[],@FunY,[],@FunPolyChebyshev1XY,ParDeg});

[S1a,FR] = align(S1,S2,'TranC',TranC); % note WCS of S1a is incorrect

Iep =
     1

Iep =
     2

SIM/transform: WCS was not corrected

% find stars in the aligned image

S1a=mextractor(S1a);

  Set gain to 1
  Subtract background
Warning: CleanByBackGrad was set to false as the user did not requested for the SN, PEAK_GRADBACK, BACK_STD columns
mextractor image 1 out of 1
  Filter image number 1 - use 2 filters
    Threshold image and locate local maxima
    Extracting 5442 sources
  Calculate 1st and 2nd moments
    Calculate PSF photometry
    Calculate IL aperture photometry (using aper_phot)

% match aligned catalogs

S1m = match(S1a,S2,'SkipWCS',true);

% difference in [X,Y] coordinates between the two images

DXY = S1m.Cat(:,1:2) - S2.Cat(:,1:2);

Mag = S2.Cat(:,20); % PSF mag

% do the same using SExtractor

% SExtractor can't work with NaNs

S1a = replace(S1a,[NaN NaN],nanmedian(S1a));

S1s=ImUtil.Im.sextractor(S1a);

[1M>
[1A----- SExtractor 2.5.0 started on 2017-09-05 at 14:16:20 with 1 thread

[1M> Setting catalog parameters
[1A[1M> Reading detection filter
[1A[1M> Initializing catalog
[1A[1M> Looking for tp937023c5_fa33_4f71_bb4f_828fadfd42d3.fits
[1AMeasuring from: "TILU_Spring_2014"  / 2048 x 4096 / 32 bits FLOATING POINT data
Detection+Measurement image: [1M> Setting up background maps
[1A[1M> Setting up background map at line:  128
[1A[1M> Setting up background map at line:  256
[1A[1M> Setting up background map at line:  384
[1A[1M> Setting up background map at line:  512
[1A[1M> Setting up background map at line:  640
[1A[1M> Setting up background map at line:  768
[1A[1M> Setting up background map at line:  896
[1A[1M> Setting up background map at line: 1024
[1A[1M> Setting up background map at line: 1152
[1A[1M> Setting up background map at line: 1280
[1A[1M> Setting up background map at line: 1408
[1A[1M> Setting up background map at line: 1536
[1A[1M> Setting up background map at line: 1664
[1A[1M> Setting up background map at line: 1792
[1A[1M> Setting up background map at line: 1920
[1A[1M> Setting up background map at line: 2048
[1A[1M> Setting up background map at line: 2176
[1A[1M> Setting up background map at line: 2304
[1A[1M> Setting up background map at line: 2432
[1A[1M> Setting up background map at line: 2560
[1A[1M> Setting up background map at line: 2688
[1A[1M> Setting up background map at line: 2816
[1A[1M> Setting up background map at line: 2944
[1A[1M> Setting up background map at line: 3072
[1A[1M> Setting up background map at line: 3200
[1A[1M> Setting up background map at line: 3328
[1A[1M> Setting up background map at line: 3456
[1A[1M> Setting up background map at line: 3584
[1A[1M> Setting up background map at line: 3712
[1A[1M> Setting up background map at line: 3840
[1A[1M> Setting up background map at line: 3968
[1A[1M> Filtering background map(s)
[1A[1M> Computing background d-map
[1A[1M> Computing background-noise d-map
[1A(M+D) Background: 1144.06    RMS: 28.8052    / Threshold: 57.6104    
[1M> Scanning image
[1A[1M> Line:   16  Objects:       26 detected /        0 sextracted
[1A[1M> Line:   32  Objects:       94 detected /        0 sextracted
[1A[1M> Line:   48  Objects:      148 detected /        0 sextracted
[1A[1M> Line:   64  Objects:      190 detected /        0 sextracted
[1A[1M> Line:   80  Objects:      210 detected /        0 sextracted
[1A[1M> Line:   96  Objects:      234 detected /        0 sextracted
[1A[1M> Line:  112  Objects:      248 detected /        0 sextracted
[1A[1M> Line:  128  Objects:      263 detected /        0 sextracted
[1A[1M> Line:  144  Objects:      277 detected /        0 sextracted
[1A[1M> Line:  160  Objects:      291 detected /        0 sextracted
[1A[1M> Line:  176  Objects:      300 detected /        0 sextracted
[1A[1M> Line:  192  Objects:      318 detected /        0 sextracted
[1A[1M> Line:  208  Objects:      337 detected /        0 sextracted
[1A[1M> Line:  224  Objects:      357 detected /        0 sextracted
[1A[1M> Line:  240  Objects:      370 detected /        0 sextracted
[1A[1M> Line:  256  Objects:      385 detected /        0 sextracted
[1A[1M> Line:  272  Objects:      403 detected /        0 sextracted
[1A[1M> Line:  288  Objects:      424 detected /        0 sextracted
[1A[1M> Line:  304  Objects:      434 detected /        0 sextracted
[1A[1M> Line:  320  Objects:      448 detected /        0 sextracted
[1A[1M> Line:  336  Objects:      470 detected /        0 sextracted
[1A[1M> Line:  352  Objects:      489 detected /        0 sextracted
[1A[1M> Line:  368  Objects:      509 detected /        0 sextracted
[1A[1M> Line:  384  Objects:      526 detected /        0 sextracted
[1A[1M> Line:  400  Objects:      534 detected /        0 sextracted
[1A[1M> Line:  416  Objects:      558 detected /        0 sextracted
[1A[1M> Line:  432  Objects:      565 detected /        0 sextracted
[1A[1M> Line:  448  Objects:      577 detected /        0 sextracted
[1A[1M> Line:  464  Objects:      595 detected /        0 sextracted
[1A[1M> Line:  480  Objects:      610 detected /        0 sextracted
[1A[1M> Line:  496  Objects:      622 detected /        0 sextracted
[1A[1M> Line:  512  Objects:      633 detected /        0 sextracted
[1A[1M> Line:  528  Objects:      654 detected /        0 sextracted
[1A[1M> Line:  544  Objects:      688 detected /        0 sextracted
[1A[1M> Line:  560  Objects:      710 detected /        0 sextracted
[1A[1M> Line:  576  Objects:      730 detected /        0 sextracted
[1A[1M> Line:  592  Objects:      751 detected /        0 sextracted
[1A[1M> Line:  608  Objects:      774 detected /        0 sextracted
[1A[1M> Line:  624  Objects:      795 detected /        0 sextracted
[1A[1M> Line:  640  Objects:      815 detected /        0 sextracted
[1A[1M> Line:  656  Objects:      834 detected /        0 sextracted
[1A[1M> Line:  672  Objects:      849 detected /        0 sextracted
[1A[1M> Line:  688  Objects:      872 detected /        0 sextracted
[1A[1M> Line:  704  Objects:      888 detected /        0 sextracted
[1A[1M> Line:  720  Objects:      912 detected /        0 sextracted
[1A[1M> Line:  736  Objects:      929 detected /        0 sextracted
[1A[1M> Line:  752  Objects:      950 detected /        0 sextracted
[1A[1M> Line:  768  Objects:      978 detected /        0 sextracted
[1A[1M> Line:  784  Objects:      995 detected /        0 sextracted
[1A[1M> Line:  800  Objects:     1011 detected /        0 sextracted
[1A[1M> Line:  816  Objects:     1034 detected /        0 sextracted
[1A[1M> Line:  832  Objects:     1059 detected /        0 sextracted
[1A[1M> Line:  848  Objects:     1082 detected /        0 sextracted
[1A[1M> Line:  864  Objects:     1104 detected /        0 sextracted
[1A[1M> Line:  880  Objects:     1128 detected /        0 sextracted
[1A[1M> Line:  896  Objects:     1149 detected /        0 sextracted
[1A[1M> Line:  912  Objects:     1173 detected /        0 sextracted
[1A[1M> Line:  928  Objects:     1207 detected /        0 sextracted
[1A[1M> Line:  944  Objects:     1224 detected /        0 sextracted
[1A[1M> Line:  960  Objects:     1257 detected /        0 sextracted
[1A[1M> Line:  976  Objects:     1282 detected /        0 sextracted
[1A[1M> Line:  992  Objects:     1301 detected /        0 sextracted
[1A[1M> Line: 1008  Objects:     1323 detected /        0 sextracted
[1A[1M> Line: 1024  Objects:     1352 detected /       32 sextracted
[1A[1M> Line: 1040  Objects:     1373 detected /       44 sextracted
[1A[1M> Line: 1056  Objects:     1394 detected /       67 sextracted
[1A[1M> Line: 1072  Objects:     1430 detected /       90 sextracted
[1A[1M> Line: 1088  Objects:     1509 detected /      109 sextracted
[1A[1M> Line: 1104  Objects:     1594 detected /      128 sextracted
[1A[1M> Line: 1120  Objects:     1654 detected /      143 sextracted
[1A[1M> Line: 1136  Objects:     1688 detected /      161 sextracted
[1A[1M> Line: 1152  Objects:     1724 detected /      174 sextracted
[1A[1M> Line: 1168  Objects:     1744 detected /      188 sextracted
[1A[1M> Line: 1184  Objects:     1788 detected /      201 sextracted
[1A[1M> Line: 1200  Objects:     1814 detected /      210 sextracted
[1A[1M> Line: 1216  Objects:     1846 detected /      227 sextracted
[1A[1M> Line: 1232  Objects:     1867 detected /      246 sextracted
[1A[1M> Line: 1248  Objects:     1890 detected /      263 sextracted
[1A[1M> Line: 1264  Objects:     1922 detected /      274 sextracted
[1A[1M> Line: 1280  Objects:     1953 detected /      295 sextracted
[1A[1M> Line: 1296  Objects:     1993 detected /      317 sextracted
[1A[1M> Line: 1312  Objects:     2026 detected /      325 sextracted
[1A[1M> Line: 1328  Objects:     2060 detected /      340 sextracted
[1A[1M> Line: 1344  Objects:     2086 detected /      354 sextracted
[1A[1M> Line: 1360  Objects:     2122 detected /      374 sextracted
[1A[1M> Line: 1376  Objects:     2157 detected /      387 sextracted
[1A[1M> Line: 1392  Objects:     2195 detected /      401 sextracted
[1A[1M> Line: 1408  Objects:     2227 detected /      413 sextracted
[1A[1M> Line: 1424  Objects:     2277 detected /      432 sextracted
[1A[1M> Line: 1440  Objects:     2354 detected /      439 sextracted
[1A[1M> Line: 1456  Objects:     2450 detected /      449 sextracted
[1A[1M> Line: 1472  Objects:     2489 detected /      471 sextracted
[1A[1M> Line: 1488  Objects:     2513 detected /      479 sextracted
[1A[1M> Line: 1504  Objects:     2539 detected /      490 sextracted
[1A[1M> Line: 1520  Objects:     2561 detected /      506 sextracted
[1A[1M> Line: 1536  Objects:     2584 detected /      519 sextracted
[1A[1M> Line: 1552  Objects:     2604 detected /      536 sextracted
[1A[1M> Line: 1568  Objects:     2624 detected /      552 sextracted
[1A[1M> Line: 1584  Objects:     2646 detected /      567 sextracted
[1A[1M> Line: 1600  Objects:     2663 detected /      578 sextracted
[1A[1M> Line: 1616  Objects:     2686 detected /      603 sextracted
[1A[1M> Line: 1632  Objects:     2699 detected /      620 sextracted
[1A[1M> Line: 1648  Objects:     2715 detected /      635 sextracted
[1A[1M> Line: 1664  Objects:     2731 detected /      647 sextracted
[1A[1M> Line: 1680  Objects:     2756 detected /      659 sextracted
[1A[1M> Line: 1696  Objects:     2771 detected /      675 sextracted
[1A[1M> Line: 1712  Objects:     2780 detected /      690 sextracted
[1A[1M> Line: 1728  Objects:     2791 detected /      709 sextracted
[1A[1M> Line: 1744  Objects:     2806 detected /      726 sextracted
[1A[1M> Line: 1760  Objects:     2820 detected /      743 sextracted

S2s=ImUtil.Im.sextractor(S2);

[1M>
[1A----- SExtractor 2.5.0 started on 2017-09-05 at 14:16:25 with 1 thread

[1M> Setting catalog parameters
[1A[1M> Reading detection filter
[1A[1M> Initializing catalog
[1A[1M> Looking for tpf4add99c_99fa_4396_8c67_8db3967d1ca4.fits
[1AMeasuring from: "TILU_Spring_2014"  / 2048 x 4096 / 32 bits FLOATING POINT data
Detection+Measurement image: [1M> Setting up background maps
[1A[1M> Setting up background map at line:  128
[1A[1M> Setting up background map at line:  256
[1A[1M> Setting up background map at line:  384
[1A[1M> Setting up background map at line:  512
[1A[1M> Setting up background map at line:  640
[1A[1M> Setting up background map at line:  768
[1A[1M> Setting up background map at line:  896
[1A[1M> Setting up background map at line: 1024
[1A[1M> Setting up background map at line: 1152
[1A[1M> Setting up background map at line: 1280
[1A[1M> Setting up background map at line: 1408
[1A[1M> Setting up background map at line: 1536
[1A[1M> Setting up background map at line: 1664
[1A[1M> Setting up background map at line: 1792
[1A[1M> Setting up background map at line: 1920
[1A[1M> Setting up background map at line: 2048
[1A[1M> Setting up background map at line: 2176
[1A[1M> Setting up background map at line: 2304
[1A[1M> Setting up background map at line: 2432
[1A[1M> Setting up background map at line: 2560
[1A[1M> Setting up background map at line: 2688
[1A[1M> Setting up background map at line: 2816
[1A[1M> Setting up background map at line: 2944
[1A[1M> Setting up background map at line: 3072
[1A[1M> Setting up background map at line: 3200
[1A[1M> Setting up background map at line: 3328
[1A[1M> Setting up background map at line: 3456
[1A[1M> Setting up background map at line: 3584
[1A[1M> Setting up background map at line: 3712
[1A[1M> Setting up background map at line: 3840
[1A[1M> Setting up background map at line: 3968
[1A[1M> Filtering background map(s)
[1A[1M> Computing background d-map
[1A[1M> Computing background-noise d-map
[1A(M+D) Background: 1253.83    RMS: 35.264     / Threshold: 70.5281    
[1M> Scanning image
[1A[1M> Line:   16  Objects:       28 detected /        0 sextracted
[1A[1M> Line:   32  Objects:       67 detected /        0 sextracted
[1A[1M> Line:   48  Objects:       86 detected /        0 sextracted
[1A[1M> Line:   64  Objects:      109 detected /        0 sextracted
[1A[1M> Line:   80  Objects:      120 detected /        0 sextracted
[1A[1M> Line:   96  Objects:      137 detected /        0 sextracted
[1A[1M> Line:  112  Objects:      149 detected /        0 sextracted
[1A[1M> Line:  128  Objects:      164 detected /        0 sextracted
[1A[1M> Line:  144  Objects:      177 detected /        0 sextracted
[1A[1M> Line:  160  Objects:      190 detected /        0 sextracted
[1A[1M> Line:  176  Objects:      199 detected /        0 sextracted
[1A[1M> Line:  192  Objects:      214 detected /        0 sextracted
[1A[1M> Line:  208  Objects:      220 detected /        0 sextracted
[1A[1M> Line:  224  Objects:      237 detected /        0 sextracted
[1A[1M> Line:  240  Objects:      250 detected /        0 sextracted
[1A[1M> Line:  256  Objects:      265 detected /        0 sextracted
[1A[1M> Line:  272  Objects:      282 detected /        0 sextracted
[1A[1M> Line:  288  Objects:      299 detected /        0 sextracted
[1A[1M> Line:  304  Objects:      311 detected /        0 sextracted
[1A[1M> Line:  320  Objects:      320 detected /        0 sextracted
[1A[1M> Line:  336  Objects:      340 detected /        0 sextracted
[1A[1M> Line:  352  Objects:      356 detected /        0 sextracted
[1A[1M> Line:  368  Objects:      371 detected /        0 sextracted
[1A[1M> Line:  384  Objects:      386 detected /        0 sextracted
[1A[1M> Line:  400  Objects:      392 detected /        0 sextracted
[1A[1M> Line:  416  Objects:      414 detected /        0 sextracted
[1A[1M> Line:  432  Objects:      420 detected /        0 sextracted
[1A[1M> Line:  448  Objects:      428 detected /        0 sextracted
[1A[1M> Line:  464  Objects:      445 detected /        0 sextracted
[1A[1M> Line:  480  Objects:      455 detected /        0 sextracted
[1A[1M> Line:  496  Objects:      467 detected /        0 sextracted
[1A[1M> Line:  512  Objects:      476 detected /        0 sextracted
[1A[1M> Line:  528  Objects:      497 detected /        0 sextracted
[1A[1M> Line:  544  Objects:      531 detected /        0 sextracted
[1A[1M> Line:  560  Objects:      560 detected /        0 sextracted
[1A[1M> Line:  576  Objects:      578 detected /        0 sextracted
[1A[1M> Line:  592  Objects:      599 detected /        0 sextracted
[1A[1M> Line:  608  Objects:      624 detected /        0 sextracted
[1A[1M> Line:  624  Objects:      649 detected /        0 sextracted
[1A[1M> Line:  640  Objects:      672 detected /        0 sextracted
[1A[1M> Line:  656  Objects:      700 detected /        0 sextracted
[1A[1M> Line:  672  Objects:      719 detected /        0 sextracted
[1A[1M> Line:  688  Objects:      741 detected /        0 sextracted
[1A[1M> Line:  704  Objects:      759 detected /        0 sextracted
[1A[1M> Line:  720  Objects:      780 detected /        0 sextracted
[1A[1M> Line:  736  Objects:      793 detected /        0 sextracted
[1A[1M> Line:  752  Objects:      810 detected /        0 sextracted
[1A[1M> Line:  768  Objects:      829 detected /        0 sextracted
[1A[1M> Line:  784  Objects:      845 detected /        0 sextracted
[1A[1M> Line:  800  Objects:      857 detected /        0 sextracted
[1A[1M> Line:  816  Objects:      875 detected /        0 sextracted
[1A[1M> Line:  832  Objects:      897 detected /        0 sextracted
[1A[1M> Line:  848  Objects:      920 detected /        0 sextracted
[1A[1M> Line:  864  Objects:      938 detected /        0 sextracted
[1A[1M> Line:  880  Objects:      956 detected /        0 sextracted
[1A[1M> Line:  896  Objects:      972 detected /        0 sextracted
[1A[1M> Line:  912  Objects:      990 detected /        0 sextracted
[1A[1M> Line:  928  Objects:     1018 detected /        0 sextracted
[1A[1M> Line:  944  Objects:     1035 detected /        0 sextracted
[1A[1M> Line:  960  Objects:     1055 detected /        0 sextracted
[1A[1M> Line:  976  Objects:     1069 detected /        0 sextracted
[1A[1M> Line:  992  Objects:     1082 detected /        0 sextracted
[1A[1M> Line: 1008  Objects:     1101 detected /        0 sextracted
[1A[1M> Line: 1024  Objects:     1125 detected /       26 sextracted
[1A[1M> Line: 1040  Objects:     1137 detected /       33 sextracted
[1A[1M> Line: 1056  Objects:     1149 detected /       47 sextracted
[1A[1M> Line: 1072  Objects:     1164 detected /       66 sextracted
[1A[1M> Line: 1088  Objects:     1200 detected /       81 sextracted
[1A[1M> Line: 1104  Objects:     1279 detected /       98 sextracted
[1A[1M> Line: 1120  Objects:     1357 detected /      108 sextracted
[1A[1M> Line: 1136  Objects:     1405 detected /      121 sextracted
[1A[1M> Line: 1152  Objects:     1434 detected /      134 sextracted
[1A[1M> Line: 1168  Objects:     1458 detected /      144 sextracted
[1A[1M> Line: 1184  Objects:     1487 detected /      159 sextracted
[1A[1M> Line: 1200  Objects:     1508 detected /      166 sextracted
[1A[1M> Line: 1216  Objects:     1534 detected /      177 sextracted
[1A[1M> Line: 1232  Objects:     1553 detected /      191 sextracted
[1A[1M> Line: 1248  Objects:     1573 detected /      206 sextracted
[1A[1M> Line: 1264  Objects:     1594 detected /      215 sextracted
[1A[1M> Line: 1280  Objects:     1625 detected /      232 sextracted
[1A[1M> Line: 1296  Objects:     1657 detected /      251 sextracted
[1A[1M> Line: 1312  Objects:     1675 detected /      263 sextracted
[1A[1M> Line: 1328  Objects:     1697 detected /      274 sextracted
[1A[1M> Line: 1344  Objects:     1717 detected /      293 sextracted
[1A[1M> Line: 1360  Objects:     1741 detected /      308 sextracted
[1A[1M> Line: 1376  Objects:     1776 detected /      321 sextracted
[1A[1M> Line: 1392  Objects:     1802 detected /      334 sextracted
[1A[1M> Line: 1408  Objects:     1833 detected /      343 sextracted
[1A[1M> Line: 1424  Objects:     1894 detected /      360 sextracted
[1A[1M> Line: 1440  Objects:     1972 detected /      366 sextracted
[1A[1M> Line: 1456  Objects:     2023 detected /      374 sextracted
[1A[1M> Line: 1472  Objects:     2042 detected /      395 sextracted
[1A[1M> Line: 1488  Objects:     2059 detected /      400 sextracted
[1A[1M> Line: 1504  Objects:     2076 detected /      412 sextracted
[1A[1M> Line: 1520  Objects:     2089 detected /      423 sextracted
[1A[1M> Line: 1536  Objects:     2103 detected /      436 sextracted
[1A[1M> Line: 1552  Objects:     2117 detected /      452 sextracted
[1A[1M> Line: 1568  Objects:     2136 detected /      465 sextracted
[1A[1M> Line: 1584  Objects:     2154 detected /      478 sextracted
[1A[1M> Line: 1600  Objects:     2168 detected /      498 sextracted
[1A[1M> Line: 1616  Objects:     2184 detected /      520 sextracted
[1A[1M> Line: 1632  Objects:     2197 detected /      540 sextracted
[1A[1M> Line: 1648  Objects:     2208 detected /      554 sextracted
[1A[1M> Line: 1664  Objects:     2221 detected /      566 sextracted
[1A[1M> Line: 1680  Objects:     2240 detected /      579 sextracted
[1A[1M> Line: 1696  Objects:     2257 detected /      591 sextracted
[1A[1M> Line: 1712  Objects:     2266 detected /      604 sextracted
[1A[1M> Line: 1728  Objects:     2276 detected /      620 sextracted
[1A[1M> Line: 1744  Objects:     2287 detected /      632 sextracted
[1A[1M> Line: 1760  Objects:     2301 detected /      645 sextracted

S1sm=match(S1s,S2s,'SkipWCS',true);

DXYs=S1sm.Cat(:,1:2) - S2s.Cat(:,1:2);

Mags = S2s.Cat(:,5); % aper mag.

Nm = size(S2.Cat,1);

Ns = size(S2s.Cat,1);

CC = zeros(0,5);

K = 0;

for Im=1:1:Nm

D = Util.Geom.plane_dist(S2.Cat(Im,1),S2.Cat(Im,2),S2s.Cat(:,1),S2s.Cat(:,2));

Id = find(D<2);

if (numel(Id)==1)

K = K + 1;

CC(K,:) = [DXY(Im,:), DXYs(Id,:), Mag(Im)];

end

end

Dm=sqrt(sum(CC(:,1:2).^2,2));

Ds=sqrt(sum(CC(:,3:4).^2,2));

nanmedian(Dm-Ds) % mextractor is a bit better: (19+/-5) mpix

ans =
      0.14207

nanstd(Dm-Ds)./sqrt(numel(Ds))

ans =
     0.005009

semilogy(CC(:,5)+6.5,Dm,'k.')

axis([13 22 1e-3 1e0]);

%print MextractorAstrom.eps -depsc2

D1 = FITS.read_table('PTF_201411204943_c_d_scie_t115144_u023050379_f02_p100037_c02.ctlg');

D2 = FITS.read_table('PTF_201411205388_c_d_scie_t125555_u023050416_f02_p100037_c02.ctlg');

D1.ColCell{2} = 'XWIN_IMAGE';

D1.ColCell{3} = 'YWIN_IMAGE';

D1.ColCell{4} = 'ALPHAWIN_J2000';

D1.ColCell{5} = 'DELTAWIN_J2000';

D2.ColCell{2} = 'XWIN_IMAGE';

D2.ColCell{3} = 'YWIN_IMAGE';

D2.ColCell{4} = 'ALPHAWIN_J2000';

D2.ColCell{5} = 'DELTAWIN_J2000';

D1 = colcell2col(D1);

D2 = colcell2col(D2);

D1m = match(D1,D2);

Photometry comparison