FIREFLY-116:Evaluate the wavelength readouts from various types of FITS and fix any Errors

src/firefly/java/edu/caltech/ipac/visualize/plot/plotdata/Wavelength.java

@@ -43,7 +43,7 @@ public Wavelength(Header header, BinaryTableHDU tableHDU) {
 
     /**
      * This method will calculate wavelength at any given image point.  The calculation is based on
-     * the paper "Representations of spectral coordinates in FITS", by E. W. Greisen1,M.R.Calabretta2, F.G.Valdes3,
+     * Reference paper: the paper "Representations of spectral coordinates in FITS", by E. W. Greisen1,M.R.Calabretta2, F.G.Valdes3,
      * and S. L. Allen.
      *
      * Current, the Linear, Log and Non-linear algorithms based on the WCS parameters/keywords in the headers
@@ -506,17 +506,31 @@ else if (sArray[1].trim().equalsIgnoreCase("LOG")){
      *
      *  lamda_r : is the reference value,  given by CRVAL3
      *
+     * * Get pixel at given "ImagePt" coordinates
+     *
+     *
+     *   (from the reference paper above) Note that integer pixel numbers refer to the center of the pixel in each axis,
+     *   so that, for example, the first pixel runs from pixel number 0.5 to pixel number 1.5 on every axis.
+     *   Note also that the reference point location need not be integer nor need it even occur within the image.
+     *   The original FITS paper (Wells et al. 1981) defined the pixel numbers to be counted from 1 to NAXIS j ($ \geq $1)
+     *   on each axis in a Fortran-like order as presented in the FITS image[*].
+     *
+     *   This method get the coordinates for given image point (p1, p2, p3) where imagePt=(p1, p2)
+     *   If it is only one plan, the p3=1 since the axis is count staring from 1.
      * @param ipt
      * @return
      * @throws PixelValueException
      */
     public static double[] getPixelCoords(ImagePt ipt, Header header) throws PixelValueException {
 
-        int p0 = (int) Math.round(ipt.getX() - 0.5); //x
-        int p1 = (int) Math.round(ipt.getY() - 0.5); //y
+        //pixel numbers refer to the center of the pixel, so we subtract 0.5, see notes above
+        //As noted above, the pixel is counting from 1 to naxis j where (j=1, 2,... naxis).  Since the p = Math.round(ipt.x - 0.5)
+        //starts from 0.  Thus, 1 is added here.
+        int p0 = (int) Math.round(ipt.getX() - 0.5) +1 ; //x
+        int p1 = (int) Math.round(ipt.getY() - 0.5) +1; //y
+        int p2 = header.getIntValue("SPOT_PL", 1) + 1;  //since SPOT_PL starts from 0
 
-        int p2 = header.getIntValue("SPOT_PL", 0)-1; //z  default is 0  // todo: I think subtracting 1 is wrong, this should be looked at
-        if (p2<0) p2= 0;
+       // if (p2<0) p2= 0;
         int naxis1 = header.getIntValue("NAXIS1");
         int naxis2 = header.getIntValue("NAXIS2");
 

src/firefly/js/visualize/FitsHeaderUtil.js

@@ -86,11 +86,19 @@ export function makeHeaderParse(header, altWcs='') {
             const retAry= [];
             let i= 0;
             for (let headerIdx = startIdx; headerIdx <= endIdx; headerIdx++) {
-                retAry[i]= getNumberHeader(header, `${keyRoot}${headerIdx}${altWcs}`,NaN);
-                if (isNaN(retAry[i])) return def;
+                retAry[i]= getNumberHeader(header, `${keyRoot}${headerIdx}${altWcs}`,def);
                 i++;
             }
             return retAry;
+        },
+        hasKeyStartWith(startKeys){
+          let key;
+          for (key in header) {
+            if (key.startsWith(startKeys)) {
+              return true;
+            }
+          }
+          return false;
         }
     };
 }
@@ -128,10 +136,19 @@ export function makeDoubleHeaderParse(header,zeroHeader,altWcs) {
         },
         getDoubleAry(keyRoot,altWcs, startIdx,endIdx,def=undefined) {
             if (startIdx>endIdx) return def;
-            const retAry= hp.getDoubleAry(keyRoot,altWcs,startIdx,endIdx);
-            if (retAry) return retAry;
+            const retAry= hp.getDoubleAry(keyRoot,altWcs,startIdx,endIdx, def);
+            let validValueArr = retAry.filter( (v)=> v!==def );
+            if (validValueArr.length>0) return retAry;
             return zhp ? zhp.getDoubleAry(keyRoot,altWcs, startIdx,endIdx,def) : def;
         },
+        hasKeyStartWith(startKeys){
+           if (hp.hasKeyStartWith(startKeys)){
+               return true;
+           }
+           else {
+             return zhp ? zhp.hasKeyStartWith(startKeys):false;
+           }
+        },
         isDefinedHeaderList: (list) => hp.isDefinedHeaderList(list) || (zhp && zhp.isDefinedHeaderList(list))
     };
 }

src/firefly/js/visualize/projection/Wavelength.js

@@ -1,7 +1,12 @@
-
-
-
-
+/**
+ *  References:
+ *  1. "Representations of spectral coordinates in FITS", by E. W. Greisen1,M.R.Calabretta2, F.G.Valdes3,
+ *     and S. L. Allen. A&A Volume 446, Number 2, February I 2006
+ *  2. "Representations of world coordinates in FITS" A&A 395, 1061-1075 (2002) DOI: 10.1051/0004-6361:20021326
+ *     E. W. Greisen1 - M. R. Calabretta2
+ *  3. https://www.aanda.org/articles/aa/full/2002/45/aah3859/aah3859.html,
+ *
+ */
 
 export const PLANE  = 'PLANE';
 export const LINEAR = 'LINEAR';
@@ -72,10 +77,11 @@ export function getWavelength(pt, cubeIdx, wlData) {
         return wl;
     }
 }
-
+//TODO check here to see the cubeIdx??
 function getWaveLengthPlane(ipt, cubeIdx, wlData) {
     const {crpix, crval, cdelt} = wlData;
-    const wl = crval + ( cubeIdx - crpix ) * cdelt;
+    //pixel count starts from 1 to naxisn
+    const wl = crval + ( cubeIdx + 1  - crpix ) * cdelt;
     return wl;
 }
 
@@ -328,14 +334,29 @@ function searchIndex(indexVec, psi) {
  *
  *  lambda_r : is the reference value,  given by CRVAL3
  *
+ * Get pixel at given "ImagePt" coordinates
+ *   "ImagePt" coordinates have 0,0 lower left corner of lower left pixel
+ *   (from the reference paper above) Note that integer pixel numbers refer to the center of the pixel in each axis,
+ *   so that, for example, the first pixel runs from pixel number 0.5 to pixel number 1.5 on every axis.
+ *   Note also that the reference point location need not be integer nor need it even occur within the image.
+ *   The original FITS paper (Wells et al. 1981) defined the pixel numbers to be counted from 1 to NAXIS j ($ \geq $1)
+ *   on each axis in a Fortran-like order as presented in the FITS image[*].
+ *
+ *   This method get the coordinates for given image point (p1, p2, p3) where imagePt=(p1, p2)
+ *   If it is only one plan, the p3=1 since the axis is count staring from 1.
+ *
  * @param {ImagePt} ipt
  * @param {number} cubeIdx
  * @return {Array.<number>}
  */
 function getPixelCoords(ipt, cubeIdx) {
-    const p0 = Math.round(ipt.x - 0.5); //x
-    const p1 = Math.round(ipt.y - 0.5); //y
-    return [p0, p1, cubeIdx];
+
+    //As noted above, the pixel is counting from 1 to naxis j where (j=1, 2,... naxis).  Since the p = Math.round(ipt.x - 0.5)
+    //starts from 0.  Thus, 1 is added here.
+    //pixel numbers refer to the center of the pixel, so we subtract 0.5, see notes above
+    const p0 = Math.round(ipt.x - 0.5) + 1 ;
+    const p1 = Math.round(ipt.y - 0.5) + 1 ;
+    return [p0, p1, cubeIdx+1]; //since cubeIdx starts from 0
 }
 
 /**
@@ -362,7 +383,7 @@ function getPixelCoords(ipt, cubeIdx) {
  * @return number
  */
 function getOmega(pixCoords, N,  r_j, pc_3j, s_3){
-    if (!pc_3j || !r_j ) return s_3*(pixCoords[0]+pixCoords[1]);
+    //if (!pc_3j || !r_j ) return s_3*(pixCoords[0]+pixCoords[1]);
 
     let omega =0.0;
     for (let i=0; i<N; i++){

src/firefly/js/visualize/projection/WavelengthHeaderParser.js

@@ -7,7 +7,9 @@ import {AWAV, F2W, LINEAR, LOG, PLANE, TAB, V2W, WAVE} from './Wavelength.js';
 
 export function parseWavelengthHeaderInfo(header, altWcs='', zeroHeader, wlTable) {
     const parse= makeDoubleHeaderParse(header, zeroHeader, altWcs);
-    return calculateWavelengthParams(parse,altWcs, wlTable);
+    const which= altWcs?'1':getWCSAXES(parse);
+    const mijMatrixKeyRoot = getPC_ijKey(parse,which);
+    return calculateWavelengthParams(parse,altWcs,which,mijMatrixKeyRoot, wlTable);
 }
 
 
@@ -19,21 +21,170 @@ const allGoodValues= (...numbers)  => numbers.every((n) => !isNaN(n));
 
 const L_10= Math.log(10);
 
-function calculateWavelengthParams(parse, altWcs, wlTable) {
+/**
+ * According to A&A 395, 1061-1075 (2002) DOI: 10.1051/0004-6361:20021326
+ * Representations of world coordinates in FITS E. W. Greisen1 - M. R. Calabretta2
+ * https://www.aanda.org/articles/aa/full/2002/45/aah3859/aah3859.html
+ * the CD_ij is for the older FITs header. The newer FITs header has PC only.
+ * 1. If both PC and CD exist, it is wrong. Instead of throwing exception, no wavelength is displayed.
+ * 2. If PC is not exist, one or more CD is exist, use CD_ij (j=1..N); if any of CD_ij
+ *   is not defined, the default is 0.0;
+ * 3. If any PC_ij is defined or neither PC nor CD is defined, we use PC, the default
+ *    is 0 for j!=i and 1 if j==i
+ *
+ *
+ * @param parser
+ * @param which
+ * @returns {*}
+ */
+function getPC_ijKey(parser,which){
+
+    let hasPC = parser.hasKeyStartWith(`PC${which}_`);
+    let hasCD = parser.hasKeyStartWith(`CD${which}_`);
+
+    if (hasPC && hasCD){
+        return undefined;
+    }
+    if (!hasPC && hasCD){
+        return 'CD';
+    }
+    return 'PC';
+}
+
 
-    const which= altWcs?'1':'3';
+/**
+ * This method will return the value of the WCSAXES.  NOTE: the WCSAXES can be any naxis, it does
+ * not have to be 3.  In general, if the wavelength is depending on two dimensional images, it most likely
+ * is 3.  But the axis 3 can also be other quantity such as Frequency etc.
+ *
+ * If the FITs header has 'WCSAXES', this will be the wavelength axis.
+ * If 'WCSAXES' is not defined, the default will be the larger of naxis and the j where j=1, 2, 3, ..)
+ * @param parse
+ * @returns {*}
+ */
+function getWCSAXES(parse){
+    const  wcsAxis = parse.getValue('WCSAXES');
 
-    const ctype= parse.getValue(`CTYPE${which}${altWcs}`);
-    const crpix= parse.getDoubleValue(`CRPIX${which}${altWcs}`, NaN);
-    const crval= parse.getDoubleValue(`CRVAL${which}${altWcs}`, NaN);
-    const cdelt= parse.getDoubleValue(`CDELT${which}${altWcs}`, NaN);
+    if (wcsAxis) return wcsAxis;
+
+    const nAxis= parse.getIntValue('NAXIS');
+
+    const ctype= parse.getValue(`CTYPE${nAxis+1}`, undefined);
+
+    if(ctype){
+        return (nAxis+1).toString();
+    }
+    else{
+        return nAxis.toString();
+    }
+
+}
+
+/**
+ *
+ * Reference:  A&A 395, 1061-1075 (2002) DOI: 10.1051/0004-6361:20021326
+ * Representations of world coordinates in FITS E. W. Greisen1 - M. R. Calabretta2
+ * paper,  the CD_ij is for the older FITs header. The newer FITs header has PC only.
+ *
+ * This method is checking the pc_ij and r_j array.  If any of the values are undefined,
+ * the default are assigned.  The default values are defined as following based on the reference
+ * paper above:
+ *
+ * PC_ij: an array defined in the FITs header. Fhe size of the dimension is N.  N is defined by
+ *  WCSAXES or naxis
+ *   i: the wcsaxes's value, if the wcsaxes = 2, then i=2
+ *   j: 0, ...N
+ *   If any of the PC_ij is not defined in the header, the following default is assigned:
+ *   1.  PC i_j = 	1.0 when i = j
+ *   2.  PC i_j = 0.0 when i!=j
+ *
+ *  If instead of using PC, CD is used, the following default is assigned:
+ *  CD i_j	0.0  NOTE: CD i_j's default is different from PC i_j
+ *
+ *
+ *  r_j: An array of the CRPIX values, the size of the dimension is N.  N is defined by
+ *  WCSAXES or naxis
+ *  j: 0,...N
+ *
+ *  Default values:
+ *   if r_j  = value of CRPIXj, if it is not defined, the default value is 0.0
+ *
+ * @param inArr
+ * @param keyRoot
+ * @param which
+ * @param N
+ * @returns {Array}
+ */
+function applyDefaultValues(inArr, keyRoot, which, N){
+    let retAry=[];
+    for (let i=0; i<N; i++) {
+        if (inArr && inArr[i] && !isNaN(inArr[i])){
+            //if inArr is defined and it has a valid value
+            retAry[i]=inArr[i];
+            continue;
+        }
+        switch (keyRoot){ //either inArr is undefined, or inArr[i] is undefined
+            case 'PC':
+                retAry[i] = (i+1) === parseInt(which)? 1:0;
+                break;
+            case 'CRPIX' || 'CD':
+                retAry[i] = 0.0;
+                break;
+        };
+    }
+    return retAry;
+}
+function isWaveLength(ctype, pc_3j){
+
+    if (ctype.trim()==='') return false;
+
+    const sArray= ctype.split('-');
+
+    //The header has the axis dependency information, ie. pc_31 (naxis1) or pc_32 (naxis2) are defined.
+    //If no such information, there is no dependency.  The wavelength will not be displayed in the mouse readout.
+    if (sArray[0]==='WAVE'  && ( pc_3j[0]!==0 || pc_3j[1]!==0) ){
+        return true;
+    }
+    return false;
+
+}
+/**
+ * NOTE:
+ *   pc_3j, means the the wavelength axis is 3.  In fact, the wavelength can be in any axis.
+ *   Which means which axis has the wavelength
+ * @param parse
+ * @param altWcs
+ * @param which
+ * @param pc_3j_key
+ * @param wlTable
+ * @returns {*}
+ */
+function calculateWavelengthParams(parse, altWcs, which, pc_3j_key,wlTable) {
+
+    /*
+    * Base on the reference: A&A 395, 1061-1075 (2002) DOI: 10.1051/0004-6361:20021326
+    * Representations of world coordinates in FITS E. W. Greisen.  The default values
+    * defined below:
+    * CDELT i	1.0
+    * CTYPE i	' ' (i.e. a linear undefined axis)
+    * CUNIT i	' ' (i.e. undefined)
+    * NOTE: i is the which variable here.
+    */
+    const ctype= parse.getValue(`CTYPE${which}${altWcs}`, ' ');
+    const crpix= parse.getDoubleValue(`CRPIX${which}${altWcs}`, 0.0);
+    const crval= parse.getDoubleValue(`CRVAL${which}${altWcs}`, 0.0);
+    const cdelt= parse.getDoubleValue(`CDELT${which}${altWcs}`, 1.0);
     const units= parse.getValue(`CUNIT${which}${altWcs}`, '');
     const restWAV= parse.getDoubleValue('RESTWAV'+altWcs, 0);
     const nAxis= parse.getIntValue('NAXIS'+altWcs);
-    // const N= parse.getIntOneOfValue(['WCSAXES'+altWcs, 'WCSAXIS'+altWcs, 'NAXIS'+altWcs], -1);
     const N= parse.getIntOneOfValue(['WCSAXES', 'WCSAXIS', 'NAXIS'], -1);
-    let pc_3j = parse.getDoubleAry('PC3_',altWcs,1,N) || parse.getDoubleAry('CD3_',altWcs,1,N);
-    let r_j = parse.getDoubleAry('CRPIX',altWcs,1,N);
+    let pc_3j = parse.getDoubleAry(`${pc_3j_key}${which}_`,altWcs,1,N, undefined);
+    let r_j = parse.getDoubleAry('CRPIX',altWcs,1,N, undefined);
+
+    //check if any value is not defined, use default
+    pc_3j = applyDefaultValues(pc_3j,pc_3j_key, which, N);
+    r_j   = applyDefaultValues(r_j, 'CRPIX', which, N);
+
     const ps3_0= parse.getValue(`PS${which}_0${altWcs}`, '');
     const ps3_1= parse.getValue(`PS${which}_1${altWcs}`, '');
     const ps3_2= parse.getValue(`PS${which}_2${altWcs}`, '');
@@ -43,12 +194,31 @@ function calculateWavelengthParams(parse, altWcs, wlTable) {
 
     const {algorithm, wlType} = getAlgorithmAndType(ctype,N);
 
-    const canDoPlaneCalc= allGoodValues(crpix,crval,cdelt) && nAxis===3;
-    if (canDoPlaneCalc && (!algorithm || !pc_3j || !r_j)) {
-        return makeSimplePlaneBased(crpix, crval, cdelt, nAxis, wlType, units, 'use PLANE since is cube and parameters missing');
+
+    const canDoPlaneCalc= allGoodValues(crpix,crval,cdelt && nAxis===3 ) ;
+
+    //We only support the standard format in the FITs header.  The standard means the CTYPEka='WAVE-ccc" where
+    //ccc can be 'F2W', 'V2W', 'LOG', 'TAB' or empty that means linear. If the header does not have this kind
+    // CTYPEka defined, we check if it is a spectra cube.  If it is a spectra cube, we display the wavelength
+    // at each plane.
+    const isWL = isWaveLength(ctype, pc_3j);
+
+    //If it is a cube plane FITS, plot and display the wavelength at each plane
+    if (canDoPlaneCalc && !isWL  && nAxis===3 ) {
+        //The FITs has wavelength planes, and each plane has the same wavelength, ie. the third axis
+        //is wavelength
+        return makeSimplePlaneBased(crpix, crval, cdelt,nAxis, wlType, units, 'use PLANE since is cube and parameters missing');
     }
 
-    if (!algorithm || !wlType) return;
+    //Plot and display the wavelength as one of the mouse readout only if the FITs header
+    //contains the required parameters and the wavelength is depending on the image axes.
+    /* We don't show the wavelength in the mouse readout in following four situations:
+     *  1. Algorithm is not defined
+     *  2. wlType is not defined or the type is not supported
+     *  3. The FITs file is not wavelength type (may be plane)
+     *  4. It has both PC_ij and CD_ij, thus, we don't know which parameters to use
+     */
+    if (!algorithm || !wlType  || !isWL || !pc_3j_key) return;
 
 
     if (algorithm===LOG){  //the values in CRPIXk and CDi_j (PC_i_j) are log based on 10 rather than natural log, so a factor is needed.
@@ -90,7 +260,6 @@ function makeSimplePlaneBased(crpix,crval, cdelt, nAxis, wlType, units, reason)
     }
 }
 
-
 /**
  * This method will return the algorithm specified in the FITS header.
  * If the algorithm is TAB, the header has to contain the keyword "EXTNAME".