Hi, I have some code dealing with those coordinate transforms which may be useful as an example, though it only covers a subset of the transformations in Greisen&Calabretta:
https://github.com/piero-ranalli/Astro-WCS-PP A (quite old) example using PDL::Transform is here (see lines 34--38 and the pdlwcstransf() and pdlwcstransfinv() functions): https://gist.github.com/piero-ranalli/60fa7b30071de86fd2db7eccbb9e3325 cheers Piero On 31/10/16 17:28, Craig DeForest wrote: > A ton of information on WCS is here: > http://fits.gsfc.nasa.gov/fits_wcs.html > > You probably want to read Paper 1. > > Here’s a brief synopsis of fields you need in an image: > (Note that FITS is old enough to be FORTRAN-based, so all counting > starts at 1 instead of 0). > > The idea for the linear transformations is that you tie the pixel and > scientific coordinate grids together at a single datum, at which you > specify both coordinate systems. Then you specify scale and rotation, > or else a (constant) Jacobian matrix. Nonlinear transformations get > more WCS parameters, but I don’t really use nonlinear WCS yet. > > > > NAXIS - number of dimensions your data points hold > > NAXIS<n> - size of each axis (note that <n> starts at 1 instead of 0) > > CTYPE<n> - name of the <n>th coordinate (e.g., “Elapsed time”, > “Declination”, “Temperature”, or whatever) > > CUNIT<n> - unit used to express the <n>th coordinate (e.g., “Sec”, > “Degree”, “Kelvin”, or whatever) > > CRPIX<n> - <n>th pixel coordinate of the datum (1-based, so the center > of the first pixel has value 1 instead of 0, and the left edge of that > pixel has value 0.5 instead of -0.5). > > CRVAL<n> - <n>th scientific coordinate of the datum (floating point) > > CDELT<n> - if present, this is the scale of the the <n>th scientific > coordinate, in CUNIT<n>s per pixel. But see below > > CD<n>_<m> - the <n,m> value of the Jacobian matrix relating pixel > coordinate to scientific coordinate. Unspecified fields default to 0. > If your data are not rotated or sheared, then this is a diagonal > matrix, and you could get away with using CDELT<n> instead. > > CROTA - you might be tempted to use this legacy keyword. Don’t. It > leads to a world of pain. > > > The PDL::Transform code will see and use these if they are present. > Note that, as currently implemented, header copying takes about as much > time as 10,000 mathematical operations — so PDL disables header copying > by default for many operations. If you want scientific headers to > propagate you have to explicitly turn on header copying in your source PDL. > > Cheers, > Craig > > > >> On Oct 31, 2016, at 10:14 AM, Ingo Schmid <[email protected]> wrote: >> >> Hi Craig, >> >> thanks for explaining. I guess the best approach is to map my >> dicom/nifti/etc. header parameters onto FITS parameters. Do you know >> of a document explaining the individual parameters? So far I've found >> this: >> https://archive.stsci.edu/fits/fits_standard/node40.html#SECTION00942000000000000000 >> >> Specifically, how do I set image orientation, rotation and off-center >> position relative to my WCS in a FITS header hash? >> Just an example, Imagine a stack of images who's orthogonal is tilted >> by 30 degrees off the y-axis and rotated 45 degrees in-plane. >> >> Actual data has a time dimension, is complex and has 2 other >> dimensions, but these are rather simple since rotations don't make >> sense there. >> >> Ingo >> >> On 10/31/2016 04:16 PM, Craig DeForest wrote: >>> Hi, Ingo, >>> >>> PDL::Transform itself deals with coordinate transformations alone. >>> So you are tasked with tracking the transformations between >>> different frames. It provides a collection of parameterized >>> transformation constructors and a basic syntax for composing and >>> inverting those transformations objects. You can apply a >>> transformation to vector data arranged as row vectors (e.g., if your >>> data have four independent variables, they are expected to have a dim >>> 0 of size 4), or map an image whose index variables are related to >>> the values of the independent variables. That functionality doesn’t >>> require tracking anything in particular beyond the coordinate >>> transformation that you are using: all the parameters are >>> encapsulated in the transformation object. >>> >>> But PDL::Transform *does* have the ability to interpret and use the >>> World Coordinate System (WCS), which is a part of the Flexible Image >>> Transport System (FITS) standard. FITS files include metadata stored >>> as keyword/value pairs in ASCII, human-readable punch-card analogs in >>> a fixed-length header record at the beginning of the image. FITS >>> headers are close enough to Perl hashes that you can place one in the >>> hash-ref header field of a PDL. (This is made easier by the >>> Astro::FITS::Header module, which creates a tied hash that behaves >>> exactly like a FITS header, but that isn’t strictly necessary). >>> >>> WCS uses certain keywords in the FITS header to represent parameters >>> of a transformation relating pixel index to actual scientific >>> parameters. PDL::Transform can interpret those keywords >>> automagically, so you can deal with the underlying physical >>> coordinates rather than the raw pixel indices. That part of the code >>> is extremely well tested and mature for images, but less well tested >>> for higher dimensionality (though it should work and has worked in >>> some applications). The FITS header interpretation currently only >>> deals with the linear portion of FITS headers: for about 5-6 years >>> there has been a standard for a family of nonlinear transformations >>> that can be specified (to deal with telescope distortions of various >>> kinds) but these have not been implemented. >>> >>> Cheers, >>> Craig >>> >>> >>>> On Oct 31, 2016, at 3:35 AM, Ingo Schmid <[email protected] >>>> <mailto:[email protected]>> wrote: >>>> >>>> Hi, >>>> >>>> finally I havetoreally face the issue of using real-world >>>> coordinates for my piddles. So far, I've ignored rotations, which I >>>> no longer can do. I guess thereare already a couple of >>>> implementations of at least some of the features I need, so before >>>> Ireinvent the wheelI wanted toask if you could point me to some >>>> solution. >>>> >>>> Basically, I have data in cuboids in space with arbitrary rotation, >>>> translation and pixel scaling. Now I would like toe.g. extract data >>>> from the same spatial volume or add them together, etc. They may >>>> have extra dimensions as well, some without equally spacing, some >>>> with non-numeric coordinates. >>>> >>>> I guess PDL::Transform is good for that. >>>> >>>> How do you effectively store and attach/receive the information? I >>>> guess as translation/rotation/scaling vectors in the header and some >>>> methods to translate the slicing? How to translate this to an >>>> arbitrary number of dimensions? >>>> >>>> >>>> Ingo >>>> >>>> ------------------------------------------------------------------------------ >>>> The Command Line: Reinvented for Modern Developers >>>> Did the resurgence of CLI tooling catch you by surprise? >>>> Reconnect with the command line and become more productive. >>>> Learn the new .NET and ASP.NET <http://asp.net/> CLI. Get your free >>>> copy! >>>> http://sdm.link/telerik_______________________________________________ >>>> pdl-general mailing list >>>> [email protected] >>>> https://lists.sourceforge.net/lists/listinfo/pdl-general >>> >> > > > > ------------------------------------------------------------------------------ > Developer Access Program for Intel Xeon Phi Processors > Access to Intel Xeon Phi processor-based developer platforms. > With one year of Intel Parallel Studio XE. > Training and support from Colfax. > Order your platform today. http://sdm.link/xeonphi > > > > _______________________________________________ > pdl-general mailing list > [email protected] > https://lists.sourceforge.net/lists/listinfo/pdl-general > ------------------------------------------------------------------------------ Developer Access Program for Intel Xeon Phi Processors Access to Intel Xeon Phi processor-based developer platforms. With one year of Intel Parallel Studio XE. Training and support from Colfax. Order your platform today. http://sdm.link/xeonphi _______________________________________________ pdl-general mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/pdl-general
