bugGNU Astronomy Utilities - Bugs: bug #50073, ImageWarp 'rotate' does not...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #50073: ImageWarp 'rotate' does not correctly propagate WCS pixel scale information into header

Submitter:  Lee Kelvin <leeskelvin>
Submitted:  Wed 18 Jan 2017 01:20:44 AM UTC
   
 
Category:  All Gnuastro Severity:  4 - Important
Item Group:  Output not reasonable Status:  Fixed
Privacy:  Public Assigned to:  makhlaghi
Open/Closed:  Closed

Jump to the original submission

Fri 20 Jan 2017 06:14:00 PM UTC, comment #7: 

Thank you for all this extra information! What bizarre behaviour! How curious that MJD-OBS should interfere with the plotting of images in DS9...? I do not understand why this should be the case.

Firstly, I followed your suggestion and upgraded my version of DS9. I am now running DS9 version 7.5. This did make a difference (bringing the before and after images into closer agreement), but did not completely fix the issue.

On your suggestion, I upgraded my WCSLIB to version 5.15. I have recompiled GNU Astro and I can now confirm that the CDELTn header keywords print all significant figures and does also print MJD-OBS into the header.

As an aside, task #14315 may help with issues such as this in the future, by propagating all FITS header keywords into any resultant output files.

I attempted to blink between the images and, once more, found a small offset. The offset now is exceptionally small (single-pixel level), and I would not have noticed it had we not started this bug hunt with a much larger offset.

I took a close look at the header to see if anything else was going on. I found a few things:

  • MJD-OBS is missing 1 significant figure when compared with the input file
  • the header keyword RADECSYS has been changed to RADESYS (missing the C)
  • A LATPOLE/LONPOLE have been assumed (rightly so?)


However, on a whim, I tried one more thing... I used DS9 to blink between exactly the same input image! Lo and behold, I found an offset! The same single-pixel level offset I mentioned above. So, I can conclude that the initial bug has now been fixed by upgrading to the latest version of WCSLIB and DS9, and this continuing mis-match must be a DS9 issue.

Thank you for your help with this bug, much appreciated.

Best,
Lee Kelvin

Lee Kelvin <leeskelvin>
Fri 20 Jan 2017 04:59:30 PM UTC, comment #6: 

Sorry for these messages, but this final confirmation idea came to me after posting the previous comment.

I manually added the `MJD-OBS' keyword in your image with this command:


astheader --write=MJD-OBS,55238.2821457 input_rotated_lee.fits


Then, after matching with DS9 it worked (and there was no more shift). So the problem was indeed due to the fact that WCSLIB 4.42 does not print the `MJD-OBS' keyword.

Just for future reference, I am also uploading my output file from the initial command as `input_rotated_mohammad.fits'.

(file #39520)

Mohammad Akhlaghi <makhlaghi>
Group administrator
Fri 20 Jan 2017 04:43:58 PM UTC, comment #5: 

Sorry, I just remembered you had uploaded your version of `input_rotated.fits'. So I downloaded it and tried it out. Indeed the `WCSLIB' keyword wasn't present. But more importantly/strangely after "blink"ing (sorry, in the previous comment, I mistakenly said "flip"!) I also noticed a shift!

So after renaming your output to `input_rotated_lee.fits', I exported both headers to text files to see the difference:


$ astheader input_rotated.fits > mohammad.txt
[mohammad@akhlaghi desktop]$ astheader input_rotated_lee.fits > lee.txt
[mohammad@akhlaghi desktop]$ diff mohammad.txt lee.txt
15,17c15,17
< PC1_1   =             -1.0E+00 / Coordinate transformation matrix element
< CDELT1  =   9.416666295793E-05 / [deg] Coordinate increment at reference point
< CDELT2  =   9.416666295793E-05 / [deg] Coordinate increment at reference point
---
> PC1_1   =                 -1.0 / Coordinate transformation matrix element
> CDELT1  =    9.41666629579E-05 / [deg] Coordinate increment at reference point
> CDELT2  =    9.41666629579E-05 / [deg] Coordinate increment at reference point
28,29d27
< MJD-OBS =        55238.2821457 / [d] MJD of observation
<
44,46c42,43
< DATE    = '2017-01-20T16:37:28' / file creation date (YYYY-MM-DDThh:mm:ss UT)
< CFITSIO = '3.41    '           / CFITSIO version.
< WCSLIB  = '5.15    '           / WCSLIB version.
---
> DATE    = '2017-01-20T14:42:33' / file creation date (YYYY-MM-DDThh:mm:ss UT)
> CFITSIO = '3.37    '           / CFITSIO version.


Besides the missing `3' digit in the `CDELT' keywords, your header file doesn't have `MJD-OBS' while `input.fits' does have it. This must be another feature that has been added to WCSLIB since version 2.24.

Mohammad Akhlaghi <makhlaghi>
Group administrator
Fri 20 Jan 2017 04:28:02 PM UTC, comment #4: 

Thanks for the check.

This cannot be due to the missing `3', the missing decimal is thirteen decimals after the point, so the effect of this change is roughly 10^{-13} of one pixel, not 10^{-1}!

The missing `3' is probably due to the precision in printing of floating points. I would guess this is because the version of WCSLIB that you are using (version 4.24) maybe only prints to 12 decimals. This is the result of the operation on my system (with WCSLIB version 5.15):


$ astimgwarp --rotate=360 input.fits

$ astheader input.fits | grep CD
CD1_1   =  -9.416666295793E-05 / Linear projection matrix
CD1_2   =   0.000000000000E+00 / Linear projection matrix
CD2_1   =   0.000000000000E+00 / Linear projection matrix
CD2_2   =   9.416666295793E-05 / Linear projection matrix

$ astheader input_rotated.fits | grep CDELT
CDELT1  =   9.416666295793E-05 / [deg] Coordinate increment at reference point
CDELT2  =   9.416666295793E-05 / [deg] Coordinate increment at reference point


Just to confirm, I manually removed the last `3' digit from the `CDELT' keywords of `input_rotated.fits' and done as you said: opened the two images as different frames in DS9 and matched them with WCS. Then clicking on the "flip" button, I couldn't see any change. Just out of curiosity, what version of DS9 are you using? (mine is version 7.5, maybe this was a bug in DS9 that has been corrected).

I would recommend updating your WCSLIB, it has implemented many new distortions as defined in the WCS paper IV (which will now hopefully be officially published in a journal after 13 years!). As a technical side-note, it also has a function to return the version of WCSLIB. So probably on your system, Gnuastro's output FITS files don't have the `WCSLIB' keyword (which reports the version of WCSLIB).

Mohammad Akhlaghi <makhlaghi>
Group administrator
Fri 20 Jan 2017 03:37:46 PM UTC, comment #3: 

Hello, and thank you for quickly addressing this issue. I can confirm that using the latest build of GNUAstro largely fixes this issue.

However, I still find that there remains a small offset still when matching the input file with an output file by WCS in a FITS viewer such as DS9. I once more make use of the file 'input.fits' attached to the original post. I also once more use the command:


$ astimgwarp --rotate=360 input.fits


Following recent naming convention changes in GNUAstro, this now produces an output file named 'input_rotated.fits', attached to this post.

The input file has a pixel scale of 9.416666295793E-05 degrees per pixel, as encoded by the CD1_1 header keyword.

The output file has a pixel scale of 9.41666629579E-05 degrees per pixel, as encoded by the CDELT1 header keyword. (Note that the WCS encoding has been converted from the CD to the PC/CDELT formalism when running ImageWarp as discussed in support thread #109223.)

You will notice that one significant figure (3) is missing from the output pixel scale. All other WCS header keywords appear identical.

When I match the input and output files in DS9 and blink between them there exists an offset of ~0.1 pixels, mostly shifting on the x-axis. I'm unsure if the loss of this significant figure in the pixel scale is the cause for this discrepancy.

Correspondence with Mohammad Akhlaghi has shown that this issue is not reproducible on his system. I therefore copy below the version numbers for the mandatory dependencies currently installed on my machine, for your reference.

Regards,
Lee Kelvin

OS: Ubuntu 15.04 vivid
GNU Astro: 0.2.61-ad122
libgsl0: 1.16+dfsg-2build1
libcfitsio3: 3.340-3
libwcs4: 4.24-1


(file #39518)

Lee Kelvin <leeskelvin>
Wed 18 Jan 2017 02:06:04 AM UTC, comment #2: 

Thanks again for this great description and wonderful suggestion. The issue has been fixed and pushed to the main repository.

The fact that the `CDELTi's were 1 was always bothering me, I just didn't have the motivation to fix it. Thanks a lot for the very clear and motivating description of this bug and also bug #50072 ;-). It would be great if WCSLIB could implement the Singular Value Decomposition also when reading from `CDi_j' inputs instead of simply setting all the `CDELTi' to 1.

You can try out any change that is pushed to the main repo. I have the described in bug #50072 (with links to the book). If you can't do bootstrapping immediately, just let me know and I will send you the tarball for Gnuastro version 0.2.57-9321 privately.

Mohammad Akhlaghi <makhlaghi>
Group administrator
Wed 18 Jan 2017 01:20:44 AM UTC, comment #1: 

This item has been reassigned from the project GNU Astronomy Utilities support tracker to your tracker.

The original report is still available at support #109222

Following are the information included in the original report:

[field #0]                  Item ID: 109222
[field #1]                 Group ID: 11425
[field #2]              Open/Closed: Open
[field #3]                 Severity: 3 - Normal
[field #4]                  Privacy: Public
[field #8]                         : Unknown support Field Display Type
[field #9]                 Category: ImageWarp
[field #10]             Submitted by: leeskelvin
[field #11]              Assigned to: None
[field #12]             Submitted on: Wed 18 Jan 2017 01:36:19 AM JST
[field #13]                  Summary: ImageWarp 'rotate' does not correctly propagate WCS pixel scale information into header
[field #14]      Original Submission: On applying an image rotation using ImageWarp, the World Coordinate System (WCS) pixel scale information is only partially transferred to the output image.

The attached image 'input.fits' is a 3x3 pixel image with WCS information in the header defined using the  CDi_j FITS standard formalism. The pixel scale for this image is 0.339 arcsec/pixel, i.e., 9.4166667e-5 degrees/pixel.

A 360-degree image rotation is performed on this image using the command:

astimgwarp --rotate=360 input.fits

This produces the output image 'image_warped.fits', attached to this report. The pixel flux values have been correctly preserved, and the output image is aligned precisely with the input image, as expected for a 360-degree rotation.

WCS information in the rotated image header has been translated from the CDi_j formalism to the PCi_j formalism. This requires not only PC1_1 and PC2_2 in the header, but also CDELT1 and CDELT2 values. CDELTn values describe the "physical increment at the reference pixel", giving the pixel scale in degrees per pixel.

Whilst the pixel scale is correctly preserved in the output PCi_j header values, both CDELTn values are equal to 1, incorrectly implying that each pixel is 1 degree in size. This creates problems in, e.g., DS9 when viewing the two files attached to this report side by side and matching by WCS. I would expect the pixel scale to be correctly propagated through to the output header CDELTn parameters.

Regards,
Lee Kelvin

[field #16]               Item Group: Output not reasonable
[field #17]                   Status: None
[field #18]        Component Version: None
[field #19]         Operating System: GNU/Linux
[field #20]          Reproducibility: None
[field #21]               Size (loc): None
[field #22]            Fixed Release: None
[field #23]          Planned Release: None
[field #24]                   Effort: 0.00
[field #28]                 Priority: 5 - Normal
[field #31]         Percent Complete: 0%
[field #33]                  Release: None
[field #58]     Custom Select Box #1: None
[field #59]     Custom Select Box #2: None
[field #60]     Custom Select Box #3: None
[field #61]     Custom Select Box #4: None
[field #62]     Custom Select Box #5: None
[field #63]     Custom Select Box #6: None
[field #64]     Custom Select Box #7: None
[field #65]     Custom Select Box #8: None
[field #66]     Custom Select Box #9: None
[field #67]    Custom Select Box #10: None

Mohammad Akhlaghi <makhlaghi>
Group administrator
Wed 18 Jan 2017 01:20:44 AM UTC, original submission:  

On applying an image rotation using ImageWarp, the World Coordinate System (WCS) pixel scale information is only partially transferred to the output image.

The attached image 'input.fits' is a 3x3 pixel image with WCS information in the header defined using the  CDi_j FITS standard formalism. The pixel scale for this image is 0.339 arcsec/pixel, i.e., 9.4166667e-5 degrees/pixel.

A 360-degree image rotation is performed on this image using the command:

astimgwarp --rotate=360 input.fits

This produces the output image 'image_warped.fits', attached to this report. The pixel flux values have been correctly preserved, and the output image is aligned precisely with the input image, as expected for a 360-degree rotation.

WCS information in the rotated image header has been translated from the CDi_j formalism to the PCi_j formalism. This requires not only PC1_1 and PC2_2 in the header, but also CDELT1 and CDELT2 values. CDELTn values describe the "physical increment at the reference pixel", giving the pixel scale in degrees per pixel.

Whilst the pixel scale is correctly preserved in the output PCi_j header values, both CDELTn values are equal to 1, incorrectly implying that each pixel is 1 degree in size. This creates problems in, e.g., DS9 when viewing the two files attached to this report side by side and matching by WCS. I would expect the pixel scale to be correctly propagated through to the output header CDELTn parameters.

Regards,
Lee Kelvin

Lee Kelvin <leeskelvin>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #39518:  input_rotated.fits added by leeskelvin (8KiB - image/fits)
file #39494:  input.fits added by leeskelvin (8KiB - image/fits)
file #39495:  input_warped.fits added by leeskelvin (8KiB - image/fits)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by makhlaghi (Posted a comment)
  • -email is unavailable- added by leeskelvin (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-01-20 makhlaghi Attached File- Added input_rotated_mohammad.fits, #39520
    2017-01-20 leeskelvin Attached File- Added input_rotated.fits, #39518
    2017-01-18 makhlaghi CategoryNone All Gnuastro
        Severity3 - Normal 4 - Important
        Item GroupNone Output not reasonable
        StatusNone Fixed
        Assigned toNone makhlaghi
        Open/ClosedOpen Closed
    2017-01-18 makhlaghi Reassign itemGNU Astronomy Utilities, sr #109222 GNU Astronomy Utilities, bug #50073

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code