FITS image : taking into account distortion coefficients

Discussion in 'Image Processing' started by Alice Pasch, Aug 10, 2016.

  1. Alice Pasch

    Alice Pasch Standard User

    Joined:
    Aug 10, 2016
    Messages:
    8
    Hello all,

    With a telescope, I took sky images in tracking ON. These images were saved in FITS format.
    With MaximDL and PinPoint Atsrometry softwares, we can get the FITS Header of each image, i.e. image data (for instance observation site, CCD camera, focal length, pixels, right ascension of the image center, declination of the image center, etc). I put an example for you (see below at the end of the message).

    What I would like to find is a way to take into account distortion coefficients of an image (we can find them in the FITS Header of the image, see below), in order to convert any pixels of the image into world coordinates.
    I am working with Python.

    Would you know a way to get it ? Or any scripts which could be useful ?

    Thank you in advance for your help.


    Example of the FITS Header of an image :

    SIMPLE = T
    BITPIX = 16 /8 unsigned int, 16 & 32 int, -32 & -64 real
    NAXIS = 2 /number of axes
    NAXIS1 = 2184 /fastest changing axis
    NAXIS2 = 1472 /next to fastest changing axis
    BSCALE = 1.0000000000000000 /physical = BZERO + BSCALE*array_value
    BZERO = 0.00000000000000000 /physical = BZERO + BSCALE*array_value
    DATE-OBS = '2016-08-03T21:53:34' / [ISO 8601] UTC date/time of exposure start


    EXPTIME = 4.00000000000E+000 / [sec] Duration of exposure
    EXPOSURE = 4.00000000000E+000 / [sec] Duration of exposure
    SET-TEMP = -30.000000000000000 /CCD temperature setpoint in C
    CCD-TEMP = -29.968750000000000 /CCD temperature at start of exposure in C
    XPIXSZ = 6.7999999999999998 /Pixel Width in microns (after binning)
    YPIXSZ = 6.7999999999999998 /Pixel Height in microns (after binning)


    XBINNING = 1 / Binning level along the X-axis
    YBINNING = 1 / Binning level along the Y-axis
    XORGSUBF = 0 /Subframe X position in binned pixels
    YORGSUBF = 0 /Subframe Y position in binned pixels
    READOUTM = 'Raw ' / Readout mode of image
    IMAGETYP = 'Light Frame' / Type of image
    FOCALLEN = 1153.0000000000000 /Focal length of telescope in mm
    APTDIA = 254.00000000000000 /Aperture diameter of telescope in mm
    APTAREA = 50670.749319791794 /Aperture area of telescope in mm^2
    EGAIN = 0.80000001192092896 /Electronic gain in e-/ADU
    SBSTDVER = 'SBFITSEXT Version 1.0' /Version of SBFITSEXT standard in effect
    SWCREATE = 'MaxIm DL Version 6.08 160505 1A3Q6' /Name of software
    SWSERIAL = '1A3Q6-HWHSQ-A0990-J2QRF-SUTCT-QC' /Software serial number


    OBJCTRA = '21 28 50.00' / [hms J2000] Target right ascension
    OBJCTDEC = '-06 44 12.0' / [dms +N J2000] Target declination
    OBJCTALT = ' 30.8774' / Nominal altitude of center of image
    OBJCTAZ = '139.4337' / Nominal azimuth of center of image
    OBJCTHA = ' -2.2794' / Nominal hour angle of center of image
    SITELAT = '43 45 11' / Latitude of the imaging location
    SITELONG = '06 55 18' / Longitude of the imaging location
    JD = 2457604.4121990739 /Julian Date at start of exposure
    JD-HELIO = 2457604.4179293886 /Heliocentric Julian Date at exposure midpoint

    AIRMASS = 1.94149884803E+000 / Airmass (multiple of zenithal airmass)
    OBJECT = ' ' / Target object name
    TELESCOP = ' ' / Telescope name
    INSTRUME = 'SBIG STT-3200 3 CCD Camera' / Detector instrument name

    OBSERVER = ' ' / Observer name
    NOTES = ' '
    FLIPSTAT = ' '
    SWOWNER = 'THALES ALENIA SPACE FRANCE' /Licensed owner of software
    INPUTFMT = 'FITS ' / Format of file from which image was read

    HISTORY File was processed by PinPoint 6.0.6 at 2016-08-08T11:19:06
    DATE = '03/08/16' / [old format] UTC date of exposure start
    TIME-OBS = '21:53:34' / [old format] UTC time of exposure start
    UT = '21:53:34' / [old format] UTC time of exposure start
    TIMESYS = 'UTC ' / Default time system
    RADECSYS = 'FK5 ' / Equatorial coordinate system
    RA = '21 28 50.00' / [hms J2000] Target right ascension
    DEC = '-06 44 12.0' / [dms +N J2000] Target declination
    FWHM = 3.31912333965E+000 / [pixels] Mean Full-Width-Half-Max of image star
    ZMAG = 2.22695229625E+001 / Mag zero point for 1 sec exposure
    EQUINOX = 2000.0 / Equatorial coordinates are J2000
    EPOCH = 2000.0 / (incorrect but needed by old programs)
    PA = 3.58967846401E+002 / [deg, 0-360 CCW] Position angle of plate
    CTYPE1 = 'RA---TAN' / X-axis coordinate type
    CRVAL1 = 3.22235640938E+002 / X-axis coordinate value
    CRPIX1 = 1.09200000000E+003 / X-axis reference pixel
    CDELT1 = -3.37671457816E-004 / [deg/pixel] X-axis plate scale
    CROTA1 = 1.03215359918E+000 / [deg] Roll angle wrt X-axis
    CTYPE2 = 'DEC--TAN' / Y-axis coordinate type
    CRVAL2 = -6.73039453583E+000 / Y-axis coordinate value
    CRPIX2 = 7.36000000000E+002 / Y-axis reference pixel
    CDELT2 = -3.38006539857E-004 / [deg/pixel] Y-Axis Plate scale
    CROTA2 = 1.03215359918E+000 / [deg] Roll angle wrt Y-axis
    CD1_1 = -3.37616668483E-004 / Change in RA---TAN along X-Axis
    CD1_2 = 6.08868227955E-006 / Change in RA---TAN along Y-Axis
    CD2_1 = -6.08264627773E-006 / Change in DEC--TAN along X-Axis
    CD2_2 = -3.37951696156E-004 / Change in DEC--TAN along Y-Axis
    TR1_0 = 1.09199958073E+003 / [private] X-axis distortion coefficients
    TR1_1 = 2.18399991151E+003
    TR1_2 = -6.18827693341E-001
    TR1_3 = -1.23992354888E+000
    TR1_4 = -1.21460280616E+000
    TR1_5 = -1.17006143126E+000
    TR1_6 = -2.36278669552E+001
    TR1_7 = 8.13968772630E-001
    TR1_8 = -1.05987616323E+001
    TR1_9 = -2.51664136543E+000
    TR1_10 = 7.75126811018E+000
    TR1_11 = 1.11725347231E+000
    TR1_12 = -1.15863502104E+000
    TR1_13 = 2.78034135924E+000
    TR1_14 = 5.04818239812E+000
    TR2_0 = 7.36000000896E+002 / [private] Y-axis distortion coefficients
    TR2_1 = -1.38948258924E-003
    TR2_2 = 1.47200000121E+003
    TR2_3 = -1.74854087048E+000
    TR2_4 = 3.06867323249E-001
    TR2_5 = -2.74476247742E+000
    TR2_6 = -2.95895994527E-001
    TR2_7 = -1.76582177998E+001
    TR2_8 = 2.41056596524E-001
    TR2_9 = -1.34421356097E+001
    TR2_10 = 6.81670566601E+000
    TR2_11 = -2.72087350085E+000
    TR2_12 = -2.86227884736E+000
    TR2_13 = 7.12777456134E+000
    TR2_14 = 1.02502948504E+001
    HISTORY WCS added by PinPoint 6.0.6 at 2016-08-08T11:19:06
    HISTORY Matched 142 stars from the Gray GSC-ACT Catalog
    HISTORY Average residual was 0.42 arc-seconds
    PLTSOLVD = T / Plate has been solved by PinPoint
    HISTORY Image data was modified, written back to file.
     
  2. Bob Denny

    Bob Denny Cyanogen Customer

    Joined:
    Oct 12, 2014
    Messages:
    439
    Location:
    DC-3 Dreams, SP, Mesa, Arizona +1 480 396 9700
    Python on Windows, with the PyWin package, can use the PinPoint API in Python's native syntax. The API is documented in detail, but you need a full license for the engine in order to use the full capabilities.

    Part of the API is the XYtoSky() method, which converts fractional XY image coordinates to J2000 Equatorial coordinates, including the distortion corrections. There is also a method SkyToXY() to go the other way. The API has an object model of the image's contents and large set of astrometric capabilities, searching for and detecting sources, etc.

    If you have a PinPoint license, you can ask me for help doing specific things, and I will help you approach and solve the problem with API guidance.
     
  3. Alice Pasch

    Alice Pasch Standard User

    Joined:
    Aug 10, 2016
    Messages:
    8
    Thank you Bob.
    I have got a PinPoint license. Should I install the API now ?
     
  4. Alice Pasch

    Alice Pasch Standard User

    Joined:
    Aug 10, 2016
    Messages:
    8
    Sorry the API is included in PinPoint. Therefore I am gonna install PyWin package and see the API documentation !
     
  5. Bob Denny

    Bob Denny Cyanogen Customer

    Joined:
    Oct 12, 2014
    Messages:
    439
    Location:
    DC-3 Dreams, SP, Mesa, Arizona +1 480 396 9700
    Ok, good! The API is already there, as it is used by MaxIm and several other programs. In the Windows Start menu, look for PinPoint and then open the User Guide. The SPI is there, and it is object oriented.

    You will need to install the PyWin package to give Python access to the objects that are registered with Windows, called COM or ActiveX. A little experimenting should result in an AhHa moment.

    If you have a PinPoint license, you have access to our Communication Center.
     
  6. Alice Pasch

    Alice Pasch Standard User

    Joined:
    Aug 10, 2016
    Messages:
    8
    Ok thanks a lot.
    Would you know the mathematical formulae used to convert pixel coordinates into world coordinates by using the distortion coefficients and other parameters of the FITS Header ? The mathematical formulae used behind XyToSky() for instance ?
     
  7. Bob Denny

    Bob Denny Cyanogen Customer

    Joined:
    Oct 12, 2014
    Messages:
    439
    Location:
    DC-3 Dreams, SP, Mesa, Arizona +1 480 396 9700
    Of course I know since I wrote the PinPoint engine, but it would defeat the purpose for our business (our only source of income) if I gave you a way to bypass using the engine, plus it would create a situation where I would not be free to improve/change the internal system for this, having published the semantics of the polynomials. I hope you understand. After all it's a 1-line call with the API :)
     
  8. Alice Pasch

    Alice Pasch Standard User

    Joined:
    Aug 10, 2016
    Messages:
    8
    I understand. I will use the API !
     
  9. Alice Pasch

    Alice Pasch Standard User

    Joined:
    Aug 10, 2016
    Messages:
    8
  10. Bob Denny

    Bob Denny Cyanogen Customer

    Joined:
    Oct 12, 2014
    Messages:
    439
    Location:
    DC-3 Dreams, SP, Mesa, Arizona +1 480 396 9700
    Hello --

    The papers and distortion coefficients you found are the SIP system from 2008. PinPoint's automatic distortion mapping and correction system dates back to 2003. I have not changed it because it works. And as I have explained, I don't want to publish the polynomials. Once again, you will get the full benefit by using the published API to transform between XY and J2000 Equatorial coordinates. Are you having trouble with the API? I might be able to help in that case.
     
  11. Alice Pasch

    Alice Pasch Standard User

    Joined:
    Aug 10, 2016
    Messages:
    8
    Thank you a lot. I have no problem with the API. The function XyToSky() works well ! I just wanted to know if it takes into account the distortion coefficients and I do not think so.
     
  12. Bob Denny

    Bob Denny Cyanogen Customer

    Joined:
    Oct 12, 2014
    Messages:
    439
    Location:
    DC-3 Dreams, SP, Mesa, Arizona +1 480 396 9700
    The functions definitely do include the distortion corrections. Are you finding a problem of some sort?
     
  13. Bob Denny

    Bob Denny Cyanogen Customer

    Joined:
    Oct 12, 2014
    Messages:
    439
    Location:
    DC-3 Dreams, SP, Mesa, Arizona +1 480 396 9700
    Oh, I see.. I just answered this on the DC-3 Dreams Communication Center as well, just a shorter answer. If you ned guidance on using the Pinpoint API please post there as this is for Diffraction/SBIG products.
     

Share This Page