taskGNU Astronomy Utilities - Tasks: task #15128, TPV and SIP keywords when WCS...

 
 

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

task #15128: TPV and SIP keywords when WCS distortions are necessary

Submitted by:  Mohammad Akhlaghi <makhlaghi>
Submitted on:  Fri 21 Dec 2018 02:27:13 PM UTC  
 
Should Start On:  Fri 21 Dec 2018 12:00:00 AM UTC Should be Finished on:  Fri 21 Dec 2018 12:00:00 AM UTC
Category:  All Gnuastro Priority:  5 - Normal
Item Group:  Enhancement Status:  Postponed
Privacy:  Public Percent Complete:  100%
Assigned to:  Sachin Kumar Singh <sks_15> Open/Closed:  Closed
Effort:  0.00

( Jump to the original submission)

Sat 13 Jun 2020 07:26:22 AM UTC, comment #36: 

Thanks a lot Mohammad for providing great help and support throughout the task. I just went through the commit and I see a lot of work done apart from the library for the interface. Thank you for taking time out of your busy schedule to help.

Also thanks to Raul, Samaeh and Pedram for providing valuable suggestions throughout this task.

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Sat 13 Jun 2020 06:08:24 AM UTC, comment #35: 

I am happy to say that Sachin's great work is now fully merged into Gnuastro master branch (in Commit 7da840d4 and this task is officially closed :-) :-) :-).

The internal library infra-structure we defined is much more general than just the TPV <---> SIP conversion and allows easy addition of other conversions in the future (when the need comes).

For a command-line interface, after playing a little with various scenarios, and considering the ease of implementation (which was really important because I am very busy now!), I thought the best solution is to put this functionality in the Fits program. for now (afterall WCS is a part of the FITS standard!).

So you can use it like this: let's assume you have a FITS file with TPV distortion, called 'tpv.fits'. You can use the command below to convert the distortion to SIP and keep the output data (and new WCS) in 'sip.fits':

astfits tpv.fits --wcsdistortion=SIP --output=sip.fits

I just defined task #15690 to hopefully find better, and more efficient, ways of doing this in the future.

The good thing with this interface is that you don't need to know the type of distortion a priori :-D! What ever it is (as long as its supported!) the output will have a SIP distortion (even if it has a SIP distortion from the start).

I would like to express a BIG THANK YOU for the job that was very well done by Sachin! Digging through that Python implementation and trying to figure out the WCSLIB internal strucure WASN'T EASY!!! As you all see from the changes in the final commit, it was a lot of work! But Sachin finally pulled it off... I just polished it and put the 'astfits' program interface (a very trivial thing!).

Raul and Samaeh, since you already have some datasets with SIP and TPV, can you please play a little with this new feature to see how it works? There may be bugs that Sachin and I haven't found yet and it would be great if we can find and fix them soon.

As always the most recent tarball and its book are available in the links below (currently, version 0.12.8-7da84):

http://akhlaghi.org/gnuastro.pdf
http://akhlaghi.org/gnuastro-latest.tar.lz

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Thu 11 Jun 2020 05:27:29 PM UTC, comment #34: 

Great! Thanks Sachin.

I just went through the changes and made some corrections, you can see them in commit 68f0c071d7e of my fork on a branch with the same name. You can easily fetch this commit over your branch and build over it. Just for reference, I am copying the commit message here (because that commit is on a temporary branch that will be removed in the future after rebasing).

Can you please correct all the remaining points and let me know then we can go into more details ;-).

The corrections can be summarized as follows:

- No TAB should be used anywhere except Makefiles (there was one in
   'THANKS').

- Alignment is generally very important (the names in the copyright
   statement of 'wcsdistortion.h' were not under each other).

- The file summary at the top of 'wcsdistortion.h' was wrong.

- In a '.h' file, only things that are related to OUTSIDE the
   corresponding '.c' file should be included. So the 'max', 'min',
   'gal_wcsdistortion_calc_tpveq' and 'gal_wcsdistortion_calcsip' shouldn't
   be in 'wcsdistortion.h'. Such header-like features that are only for the
   '.c' file should be put at the start of the '.c' file.

- In 'gal_wcs_distortion_identify', we now simply return the corresponding
   macro, instead of saving it to an integer value and returning that, it
   is cleaner like this.

- There is a problem in 'gal_wcs_distortion_identify' (as mentioned in the
   comments, after 'PROBLEM:').

- All functions in 'wcsdistortions' (except the two that go into its
   header) must be static (some weren't), AND, they must not start with
   'gal_'. The 'gal_' prefix is only for functions that are go into the
   headers (used in other files).

- In 'wcsdistortion_intermidate_tpveq', the '+' is put at the end of the
   each line, followed by a '\'. In C when the line isn't a macro we don't
   need a '\'. When you want to break a statement, its more readable if you
   put it in an outer parenthesis and put the statements under that.

- Following the previous point, according to the GNU Coding Standards,
   when you break an expression, its better to put the operator at the
   start of the next line, rather than the end of the previous line. This
   is much more readable and highlights the operator ('+') in this case.

- I done the above two points for 'k[0][2]' and 'l[0][2]', but please
   correct them for the rest also ;-).

- No line should be longer than 75 characters, unless there is absolutely
   no other way! So the 'sprintf' statements in
   'wcsdistortion_add_sipkeywords' were corrected. Please correct the
   similar statements for the other function also.

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Thu 11 Jun 2020 11:34:41 AM UTC, comment #33: 

The new library `wcsdistortion` is tested and pushed on github and gitlab.

Two high-level functions `gal_wcs_distortion_identify` and `gal_wcs_distortion_convert` are added in wcs.c. Final testing is to be done before rebasing and merging it to master branch.

About the interface, I think as the library only has basic functionalities as of now, a common option interface is more apt for it for now. A full-fledged program can be written once the library is mature enough to support more functionalities:-). Any thoughts?

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Tue 09 Jun 2020 10:54:09 AM UTC, comment #32: 

Great Sachin!
I think in the end the best situation is that you can do both parts of the astrometry solution.
But Maybe at this time, the important things that there isn't any good solution for it is distortion.
So I think it would be better to do the first option.

Zahra Sharbaf <samaeh>
Mon 08 Jun 2020 11:12:49 AM UTC, comment #31: 

Great work! I think it would be better to have a program in Gnuastro devoted to WCS-related operations. I think next steps are related with having in Gnuastro a program to obtain the astrometric solution of images. So, I would start having a program for all of this WCS-related things :-)

Raul Infante-Sainz <infantesainz>
Project Member
Sun 07 Jun 2020 10:44:21 PM UTC, comment #30: 

This is great Sachin! This task is almost complete, so I marked it as 80% completed. All that is left is (as you said, sanity checks) and also an interface.

Two options come to my mind (I like both of them, and maybe later we can indeed have both of them! but we'll need to select one for a start).

1. Adding a new common option to all Gnuastro programs, something like '--wcsdistortion', allowing users to set their desired distortion standard for the output of any program.

2. Adding a whole new Gnuastro program devoted to WCS-related operations!

Raul and Samaeh, what do you think?

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Sun 07 Jun 2020 10:29:17 PM UTC, comment #29: 

task #15676 is done and closed now. The final headers are shown there along with the code for the same. Now basic sanitisation of code remains(sanity and error checks and code format) before porting it in Gnuastro.

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Fri 05 Jun 2020 06:48:35 PM UTC, comment #28: 

The TPV->SIP conversions are done and task #15670 is closed. Headers from both the astropy and the test Gnuastro library are shown there along with links to the repository:-)

Next part is task #15676.

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Tue 02 Jun 2020 09:54:26 PM UTC, comment #27: 

Yes, he really helped a lot. So, I'll make the new branch and start exporting changes from the test library to it and also thank him.

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Tue 02 Jun 2020 09:43:55 PM UTC, comment #26: 

This is great! Good work Sachin :-)!

This just occurred to me: Since David helped in this task we should officially thank him in Gnuastro. So after your test projects are complete and you want to put them in Gnuastro, don't forget to add his name in the 'THANKS' file.

Maybe you can even start your Gnuastro branch for this task now and start it by thanking David, for now then committing the future commits of this task over it ;-).

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Tue 02 Jun 2020 09:37:50 PM UTC, comment #25: 

task #15669, which was a subtask of this task is done and closed. I've compared the results of both the sip_tpv and my implementation there and also provided links to the code and scripts.

The next task is task #15670, in which reverse coefficients are to be calculated.

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Sun 31 May 2020 12:33:32 PM UTC, comment #24: 

The reason of not having the equations is that the problem doesn't have an analytical solution? In that case I think the way is to follow the Python implementation (as you said), and let see how it goes. I think the final solution would be avoid using PV, which is only necessary because Swarp uses them. But in that case we will need a program (astwarp?) which where able to take as input the SIP distortion parameters to resample the image. Thanks for the great job!

Raul Infante-Sainz <infantesainz>
Project Member
Sat 30 May 2020 10:37:12 AM UTC, comment #23: 

The forward coefficients were done, but problems arose during the implementation of reverse coefficients. It was also not easy to emulate the python implementation without equations(which were not given in the paper). So I asked the maintainer of sip_tpv module, David Shupe, to give me inner workings of different components to help me with the implementation. He replied back:

"To compute the SIP reverse coefficients, what is done in sip_tpv is to create a grid of points, use the forward coefficients to populate points, and then fit polynomials from that grid of points. The reverse polynomials are an approximation to the forward ones. Some people do not like that; for example, for Hubble Space Telescope data, only the forward coefficients are provided, and an iterative solver must be used for the sky-to-image pixel direction."

Now, WCSLIB does the actual analysis of the polynomials to compute sky <---> pixel translations, and as I don't have the exact equations, I think it's best to use just do the implementation as shown in the python implementation and see how it goes(as now he has explained it, it becomes more understandable to implement it as-is)? Any thoughts?

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Tue 26 May 2020 07:47:29 AM UTC, comment #22: 

After a bit more digging, it seems the WCSLIB doesn't have the functionality for the desired inter-distortion conversions. I tried to convert the SCAMP output containing PV coefficients to TPD. The results are attached in the file. As can be seen, there are no real conversions due to absence of required parameters in the input fits file.

So, it will be convenient and useful to focus on TPV to SIP(and vice-versa) conversions for now. We can, later on, do the calculations to do other distortions conversions as well:-)

(file #49170)

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Mon 18 May 2020 11:34:50 PM UTC, comment #21: 

Good progress! I have had a general look at it and you seem to be on the right path ;-).

Just one small thing: 'wcspih' is the general function that is used to read all WCS keywords from the input FITS header into 'wcsprm'. Its not just limited to the distortion keywords. In Gnuastro, we already use it in the gal_wcs_read_fitsptr function of 'lib/wcs.c'.

Now that you have the general outline and want to start digging deeper into implementation details, you can use the 'main()' function I gave before and start playing with the contents of the `wcsprm' to see effect of your changes on the written WCS ;-).

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Mon 18 May 2020 11:12:36 PM UTC, comment #20: 

I read the wcslib's official documentation and now have an understanding of it's working, at least for the current task. I am enlisting the points below:

  • the `dis.h` the header has all the necessary structures and functions for distortion corrections.
  • The main structure, `dispram`, has the first element flag whose value is to be set as 0 when initialising the `dp` dpkey structure. The initialisations are done through an internal function `disinit()` whose value is to be set as -1 for first initialisation only.
  • The keywors for TPD distortion like DPja or DQia are all placed inside the `dpkey` structure.
  • The initialisation in the the `dispram` is done through the function `disset()`. For this particular task, it will automatically be invoked by disp2x() function.
  • The dishdo() function tells wcshdo() to write FITS headers in the form of the TPD translation used internally, but I still have to figure out if it will be of any use.
  • As wcspih() is a more general function to load the TPV and SIP keywords into `dispram`, it makes more sense to use it rather than specific ones(like wcsset() for TPV to TPD conversion of keywords).
  • To specify the particular distortion to be converted into TPD, we can use the `dpram` array in `dispram` to store specific coefficients of distortions in the general TPD coefficients for automatic detection and conversion. I still have to figure out how to go about it.

This is a general overview of how the read and write function's programmatic workflow. If you find any improvements or changes please feel free to comment:-)

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Mon 18 May 2020 11:27:23 AM UTC, comment #19: 

Thanks a lot, Raul! We're trying to fix this conversion and trying to generalise it to accommodate more such distortion conversions into a single distortion coefficient, TPD. Any help and feedback is deeply appreciated:-)

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Mon 18 May 2020 08:33:02 AM UTC, comment #18: 

Great work! I had problems with dealing with SIP and PV for correcting the optical distortion of images. If you need some feedback I can test the programs to check how it works. Thanks a lot guys.

Raul Infante-Sainz <infantesainz>
Project Member
Sun 17 May 2020 10:26:39 PM UTC, comment #17: 

Yes ;-), first focus on the WCS features to possibly do the conversion using WCSLIB itself.

We do later need features to extract the distortion matrix, but that can be better done after the first step above and understanding how to manipulate WCSLIB's internal structure ;-).

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Sun 17 May 2020 10:22:24 PM UTC, comment #16: 

Also, as we now know about the general TPD distortion, I think WCSLIB must already have facilities for conversions of other distortion to TPD, which I think makes focusing on TPV to SIP conversion redundant. So should I proceed with it or should try to make functions that will simply read the distortions and write them as the TPD distortions? What do you suggest?

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Sun 17 May 2020 10:11:20 PM UTC, comment #15: 

I mean with TPD distortion being a more generalized element already present in WCSLIB, I will try to find out how the individual distortion like the TPV to SIP changes can be done so that they can be incorporated in the general TPD distortions.

For example, for TPV to SIP changes we need the forward and reverse coefficient calculation in the current implementation(of which the forward is already implemented). Now I wanted to find out are there any such features/functions already present in `dis.h` to make some portion or calculations easy.

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Sun 17 May 2020 09:28:13 PM UTC, comment #14: 

Okay, so I'll see for all the changes(and additions) in the current implementation(like the usage of reverse coefficients and such) and will discuss the desired changes:-)

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Sun 17 May 2020 08:58:51 PM UTC, comment #13: 

The algorithm implemented in the Python program is a very special case of just converting between SIP to TPV and vice-versa. But if we can find a way to use WCSLIB for this job (possibly through TPD), it is much more robust and general (allowing conversion between even more standards). It also allows us to avoid low-level matrix manipulation.

So I recommend to look in detail into the WCSLIB library and how it behaves on the various datasets (like the PV image you already have).

See how you can access the TPD distortion elements (which is apparently how it stores any of the distortions TPV, SIP, DSS, and WAT), and possibly set WCSLIB to write the headers in a different standard than the input dataset's standard.

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Sun 17 May 2020 08:51:53 PM UTC, comment #12: 

I look into it. But will it still be compatible with the program flow I've written till now? I mean willany major change be required? 

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Sun 17 May 2020 04:20:05 PM UTC, comment #11: 

Great! I am happy the problem with `rval' is clear now ;-).

By the way, while I was looking through the WCSLIB dis.h manual, I encountered this statement (under the TPD title):

"The "Template Polynomial Distortion" (TPD) is a superset of the TPV, SIP, DSS, and WAT (TNX & ZPX) polynomial distortions that also supports 1-D usage and inversions.... Within WCSLIB, TPV, SIP, DSS, and WAT are all implemented as special cases of TPD. ... Moreover, the general Polynomial distortion is translated and implemented internally as TPD whenever possible."

Given its generality, and the fact that WCSLIB can apparently do the translations of all these standards to-and-from TPD internally, a more robust work may be to base all our work on TPD (also for finding the implementing the distortions in task #15637 and task #15636).

You may even find some part of WCSLIB that will automatically do the conversion of TPV to SIP and vice-versa through TPD ;-).

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Sun 17 May 2020 02:51:48 PM UTC, comment #10: 

Also, after a bit more inspection, it seems that the  `rval` value is okay as it states in the SPIE paper that:

"The method also recomputes the PVi_j keywords to match this new CD matrix, constraining PV1_1 = PV2_1 = 1 and PV1_0 = PV1_2 = PV2_0 = PV2_2 = 0, so that only the quadratic and higher-order P V terms are needed."

So it just a way to remove these terms, and so is fixed now.

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Sun 17 May 2020 02:36:42 PM UTC, comment #9: 

Ok, I understand now! I'll check it out and will report the progess soon:-)

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Sun 17 May 2020 02:29:05 AM UTC, comment #8: 

I just had a fast look at the WCSLIB manual mentioned before and noticed that you will probably find the dis.h File Reference (section 19.2 of the PDF manual) very useful in this work ;-).

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Sun 17 May 2020 02:14:16 AM UTC, comment #7: 

Sachin, as you dig deeper into WCSLIB and learn its inner structure, please avoid direct manipulation of the raw 'wcsprm' elements in the main conversion function. Rather, define functions to abstract these low-level steps-away.

For example to read the distortion elements, define a function like this:

int
gal_wcs_distortion_read(struct wcsprm *wcs, double *distortion);

This function won't actually touch anything inside the WCS structure, it will just read all the distortion coefficients from the 'wcsprm', and write them into the 'distortion' array (that the caller has already allocated). It will return an integer that specifies a code/macro specifying the distortion (for example SIP, TPV, or other types which you can see from the WCSLIB manual). We can define an 'enum' to assign integers to every distortion.

Later, when you want to write a new distortion into a wcsprm define another function like this:

int
gal_wcs_distortion_write(struct wcsprm *wcs, double *distortion, int type);

You will then use this after converting TPV to SIP here, or writing a whole new distortion in task #15637.

Once you have abstracted these two major operations into such general functions, you can then write a clean and easily-readable function to convert TPV to SIP and vice-versa.

Once the conversion function is written, we can add an even higher-level function to read the distortion as SIP like this (it will internally convert a TPV header to SIP if necessary, and always return a SIP distortion, the return value can be 0 or 1 depending on its success):

int
gal_wcs_distortion_read_to_sip(struct wcsprm *wcs, double *distortion);

For example, you will need this in the final step of your GSoC project (task #15636) where you want to use the distortions to apply a non-linear warping. In that time, you don't want to worry about the format of the input's distortion. So can use this function to always be sure that you have a SIP distortion for example and easily focus on the main job there ;-).

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Fri 15 May 2020 11:39:12 PM UTC, comment #6: 

By the way, just a small technical tip: If the file given to 'gal_fits_img_write' already exists, this function won't over-write it, it will append the given dataset as a new extension/HDU to the existing FITS file.

To simplify your work, before compiling and running your code, delete the output file, something like this:

rm -f test-sip.fits && astbuildprog temp-wcs.c

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Fri 15 May 2020 11:32:01 PM UTC, comment #5: 

This is great progress Sachin, nice work!

You can see my recommendation on the structure of the 'main' function below. In short, you should move all your internal steps/functions into a high-level 'gal_wcs_tpv_to_sip' (which only takes the WCS structure and modifies it internally). After control comes back into 'main', you just read the image data, add the WCS sructure to it and write it.

In short, you need to look closer into the wcsprm structure and see what parts to change, so later when WCSLIB is writing the WCS-related keywords into the FITS image, it writes the newly added SIP structure and there is no more TPV structure.

About the 'rval' issue, what is the most reasonable situation? I don't have time to look into the details of that implementation, but given all the things you have read and also their own coding, you should be able to identify the expected behavior ;-). Ultimately, you can just finish the implementation and see the effect of changing it in practice ;-).

void
gal_wcs_tpv_to_sip(struct wcsprm *wcs)
{
  # This is where you can start the work:
  # Independent of the input/output filenames and HDUs and etc.
  # This task is only inside the input WCS structure ;-).
}

int
main()
{
  /* The only input should be the filename and HDU. */
  char *outfile="test-sip.fits";
  char *infile="test-pv.fits", *inhdu="1";

  /* The outputs. */
  int nwcs;
  gal_data_t *out=NULL;
  struct wcsprm *wcs=NULL;

  /* Read the WCS from the input. */
  wcs=gal_wcs_read(infile, inhdu, 0, 0, &nwcs);

  /* Define this function above to contain all your work, only taking
     the WCS structure as input. */
  gal_wcs_tpv_to_sip(wcs);

  /* Read the data of the input file. */
  out=gal_fits_img_read(infile, inhdu, -1, 1);

  /* Put the modified WCS into the 'out' dataset, then write the
     output dataset. */
  out->wcs=wcs;
  out->nwcs=nwcs;
  gal_fits_img_write(out, outfile, NULL, NULL);

  /* Clean up and return. */
  gal_data_free(out);
  return EXIT_SUCCESS;
}

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Fri 15 May 2020 11:05:16 PM UTC, comment #4: 

The major computation part of the `pv2sip` is complete:). Now only the addition of the keywords in the header remains.

For that, I searched the documentation and found out about it but in the fits.c libraries. How to do the writing of keywords using the `struct wcspram` structure?

Also, the doubtful value is `rval` for which I raised an issue is yet not confirmed, so what should I do with that calculation? I've attached the independent file, which I tested with the given fits file. I'll modify the arguments and style further when I import it to the library:)

Also, after this review, I'll version control it as well.

(file #49086)

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Wed 13 May 2020 10:38:18 PM UTC, comment #3: 

Thanks for the nice progress report Sachin.

About your question on the need for TPV, the whole problem is SCAMP! It is currently one of the only programs that does the distortion calculations properly. But unfortunately it uses TPV (like the example FITS file I attached in my previous comment). Through this task, we want to be able to internally convert such TPV headers into SIP, so even when the input to any Gnuastro program has TPV headers, it will be converted to SIP on the fly and programs can always assume that distortions (when necessary) are in SIP. In the inverse scenario, we want to let people choose to write the distortions in TPV if they want, probably as an extra command-line option.

Here are some points:

  • Instead of using 'gal_fits_key_read_from_ptr', I recommend to use WCSLIB's wcsprm structure. Reading all the keywords with the various histories of the WCS standard is not easy! Its best to leave the low-level WCS keyword reading to WCSLIB ;-). Just set your function to accept 'struct wcsprm *wcs' as input, similar to all the other functions in Gnuastro's WCS-related functions. If your test implementation isn't yet in Gnuastro, you can use 'gal_wcs_read_fitsptr' or 'gal_wcs_read' to read the 'wcsprm' structure from a given FITS pointer of file. Fortunately WCSLIB's 'wcs.h' (for example in '/usr/include/wcslib/wcs.h') has a very nice explanation on the contents of this structure as comments.
  • In the end, you should also write your results into this 'wcsprm' structure. This way, the writing of the keywords will also be left to WCSLIB and we don't have to worry about it ;-).
  • Essentially, your function should take a `wcsprm' structure and correct it internally, so it's declaration looks something like this:

void
gal_wcs_tpv_to_sip(struct wcsprm *wcs);

  • About the if-else case in 'calcsip' function, I guess its best to contact is author, maybe you can post an issue on Github for them to answer ;-).
  • About what to write in the output headers, as mentioned above, I think its best to leave the final writing to WCSLIB (which is written by the author of the WCS standard and is effectively the WCS standard!). We shouldn't make the same mistake as the author of SCAMP and make our own custom formats/conventions ;-).
Mohammad Akhlaghi <makhlaghi>
Project Administrator
Wed 13 May 2020 10:14:43 PM UTC, comment #2: 

To help play with a real-world example, I just upload this example file: http://akhlaghi.org/data/astrometry/2923938133-pv.fits

It has PV headers and created by SCAMP. Ideally, we want to convert these to SIP. To test if the conversion has been done correctly we can use DS9:

1. Load both images (as separate "frame"s), first the original one in PV, then the converted one with SIP.

2. Load the following "region" file (contents shown at the end of this message, put these in a plain-text `.reg' file) into both "frames" (Select "Region" menu --> "Load Regions ..." --> Select file containing the contents below).

3. By flipping through the two "frames" the regions of both images should be identical.

Here is the contents of the region file:

# Region file format: DS9 version 4.1
global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1
fk5
circle(113.7174541,31.8758297,40")
circle(114.9751658,30.7796376,40")
circle(115.1066211,33.1931181,40")
circle(112.9084759,33.2238771,40")
circle(112.8294555,31.5815660,40")
circle(114.3327922,32.3422483,40")
circle(112.7558973,30.7745099,40")
circle(113.9851253,31.4624990,40")

Mohammad Akhlaghi <makhlaghi>
Project Administrator
Wed 13 May 2020 08:47:43 PM UTC, comment #1: 

I found these ELF 64-bit LSB binary executable files for sip2pv and pv2sip whose link for python program is given(the python programs are its prototypes).

Here are some points:

  • This file really made it easy for me to understand the SPIE paper better.
  • I went through the python code and for the first part, I've made a raw C implementation for `pv2sip`(after this `sip2pv` can be easily done). I'll use `gal_fits_key_read_from_ptr` to read the keywords. In place of `sympy expressions` I think using a 2x2 matrix according to the indices of the u-v expression. The resulting matrix will be an upper triangular matrix so we can default the remaining values to be null(or something). For calculating the expression for TPV distortion, I'll just multiply the matrix with the calculated inverse of a 2x2 CD matrix(But doesn't the paper suggests to use only SIP and PV keys so that there should be no use of TPV keywords, so why do we need this TPV expression anyway?).
  • This if-else condition in`calcsip` function seems wrong as it gives the same value of `rval` either way. Please confirm.
  • Should I calculate the reverse coefficients too?
  • And finally should I add options to remove `SIP` or `PV` keywords from the final header before merging it in the fits data. I think having both keywords produces uniformity so there shouldn't be any removal.

What are your suggestions:-)

Sachin Kumar Singh <sks_15>
Project MemberIn charge of this item.
Fri 21 Dec 2018 02:27:13 PM UTC, original submission:  

Currently, Gnuastro just uses WCSLIB's raw output when it wants to write WCS keywords that are distorted. In particular, it will use the SIP distortion coefficients (which has become the commonly used standard recently).

But, as discussed in Proc. of SPIE Vol. 8451, 84511M,  there are many software packages (Astromatic software in particular) that use TPV.

In that SPIE paper, they suggest to include both conventions in the header (with `CTYPE' referring to SIP, but TPV also present as `PV*' keywords.

It would be good if we apply this standard in Gnuastro's outputs also (to allow easy integration with Astromatic).

WCSLIB does appear to have a flag for writing TPV keywords (`WCSHDO_PVn_ma', when calling `wcshdo'), but I couldn't see any change when I changed the current value (`WCSHDO_safe') in the call to `wcshdo' of `lib/fits.c' (using WCSLIB 6.2).

So we might have to do this manually. That paper does give some of the basic equations, but I couldn't find the `sip2pv' and `pv2sip' C programs it talks about! However, they have apparently we-written it in Python and are distributing that here.

So if there is no better solution, we can translate this Python code into C and include it in Gnuastro's WCS library (for general use).

Afterwards, it is trivial to also add a common option to all programs (maybe called `--tpvinwcs'). When it is called, Gnuastro's programs can follow the SPIE paper's suggestion and have both headers present.

For the time being, I am busy and won't have the chance to work on this, so if you are interested please post a comment (for encouragement!), or just take up the job (its just a relatively simple translation from Python).

Mohammad Akhlaghi <makhlaghi>
Project Administrator

 

Attached Files
file #49170:  header.txt added by sks_15 (4KiB - text/plain)
file #49086:  temp-wcs.c added by sks_15 (7KiB - text/x-csrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by samaeh (Posted a comment)
  • -email is unavailable- added by makhlaghi
  • -email is unavailable- added by infantesainz (Posted a comment)
  • -email is unavailable- added by makhlaghi
  • -email is unavailable- added by pedram
  • -email is unavailable- added by sks_15 (Posted a comment)
  • -email is unavailable- added by makhlaghi (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only logged-in users can vote.

     

     

     

    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-06-13 makhlaghi Item GroupNone => Enhancement
        Percent Complete80% => 100%
        Assigned toNone => sks_15
        Open/ClosedOpen => Closed
    2020-06-07 makhlaghi Percent Complete0% => 80%
        Carbon-Copy- => Added zahra.sharbaf2@gmail.com
    2020-05-26 sks_15 Attached File- => Added header.txt, #49170
    2020-05-15 makhlaghi Carbon-Copy- => Added infantesainz@gmail.com
    2020-05-15 sks_15 Attached File- => Added temp-wcs.c, #49086
    2020-05-13 pedram Carbon-Copy- => Added pedram

    Back to the top


    Powered by Savane 3.5