Thomas,

Just wanted to say that I find this very interesting. As an astronomer, I
have lots of .fits data from HST observations. I've not worked with them
for a while -- when I did I used IDL routines written by Space Telescope
staff. If I were to go back to them now, I probably would use python, and
especially the tools in the "astropy" distrubition. But it is neat to have
a native tool in J. I've just tried your code & it does display the test
images (though one of my old Hubble images fails). So when I get a bit of
time I'll try to understand what you have done. Anyway, I just wanted to
let you know that there is interest in this project!

Patrick

On Thu, 19 Apr 2018, Thomas Hickey wrote:
This is, of course, great.  Thank you and I'll try to incorporate most of
it.

I noticed that the links in my gmail copy of this post worked, but those
with trailing periods in the archive don't because the link there somehow
incorporates the period.

--Th


On Wed, Apr 18, 2018 at 6:16 PM, Marshall Lochbaum <[email protected]>
wrote:

This is a perfectly good use of the forums in my book. At the very
least, we've had similar posts before and no one's complained.

Overall your program looks pretty good to me. It's not leveraging the
Full Power of J™, but it's a great start. You know what I. does!

There are some ways to use arrays to do more of your work for you. You
write

        'minCut maxCut' =. x
        ...
        minWanted =: (<.minCut*#srtData){srtData
        maxWanted =: (<.maxCut*#srtData){srtData

but this is the same as

        'minWanted maxWanted' =: (<.x*#srtData){srtData

since each of * <. and { has rank zero, that is, maps over arrays, on
the side of x. (Unrelatedly, these should probably use =. rather than =:
to avoid global assignment. It looks like they might have been =: for
debugging, but if so, it's good practice to look over the code again and
change them back when you're done.)

A more sophisticated example is the plane-scaling code in displayFits.
You use different code for each different number of scalings to do:

        r =.  minmax scaleT2D 0{imarray
        select. 3<.{.$imarray   NB. use up to 3 planes
          case. 3 do.
                g =. minmax scaleT2D 1{imarray
                b =. minmax scaleT2D 2{imarray
          case. 2 do.
                b =. minmax scaleT2D 1{imarray
                g =. -: r+b   NB. average of red and blue
          case. 1 do.
                'g b'=. r;r
        end.

Instead, you could do all the scalings in one go:
((#imarray) is like ({.$imarray), but it is a scalar and not a vector.)

        numScales =. 3<.#imarray   NB. use up to 3 planes
        scaled =. minmax scaleT2D"_ _1 numScales {. imarray
        if. numScales < 3 do.
                avg =. (+/ % #) scaled
                scaled =. scaled , (3-numScales)#,:avg
        end.
        'r g b' =. scaled

This example first uses scaleT2D with an explicit rank to operate on
each of the first numScales planes. If you're not familiar with it, my
favorite explanation of the rank conjunction (") is
http://www.jsoftware.com/help/jforc/loopless_code_i_verbs_have_r.htm.
A rank of (_ _1), that is, infinity on the left and minus one on the
right, means to apply the function to the entire left argument and each
subarray of the right argument, effectively mapping over only the right
argument. I could have used a rank of 2 rather than _1 on the right, but
rank _1 is more flexible--it would work even if the planes were instead
three-dimensional.

After this compuation, scaled is a three-dimensional array of the first
numScales planes of imarray. The average of its planes can be found with
the standard averaging function (+/ % #), which happens to be the
canonical example of a fork. If you're interested, these are described
at http://code.jsoftware.com/wiki/Vocabulary/fork.

An exponent like (2^8) unfortunately always evaluates to a
floating-point number. I habitually use (<.@^) when I know the arguments
are integers to avoid accidentally promoting to floats.


I actually gave a talk on a task similar to the parsing you do here. My
task was to parse the header of a simple wave file. A script version
(there's no recording that I'm aware of) is at
http://code.jsoftware.com/wiki/Community/Conference2014/
Talks/UsingDataAsCode
and the current script that I use based on this is
https://github.com/mlochbaum/JSound/blob/master/wav.ijs.
This is certainly difficult material for a beginning J programmer, but
the techniques discussed in the talk might be useful for you as you gain
more skill.

One neat trick similar to the one used there: rather than write

        bitpix =. hdata getHdrVal 'BITPIX'
        naxis  =. hdata getHdrVal 'NAXIS'
        naxis1 =. hdata getHdrVal 'NAXIS1'
        naxis2 =. hdata getHdrVal 'NAXIS2'
        naxis3 =. hdata getHdrVal 'NAXIS3'

you might use

        fields =. ;: 'BITPIX NAXIS NAXIS1 NAXIS2 NAXIS3'
        (fields) =. hdata&getHdrVal&.> fields

and then the fields will be defined in variables BITPIX, etc.

Another very cool trick that I found after the talk would be useful
here:

        select. bitpix
          case. _32 do. adata =. (naxis3,naxis2,naxis1)$ |. _1 fc |. rdata
          case. _64 do. adata =. (naxis3,naxis2,naxis1)$ |. _2 fc |. rdata
          case. 16  do. adata =. (naxis3,naxis2,naxis1)$ |. _1 ic |. rdata
          case. 32  do. adata =. (naxis3,naxis2,naxis1)$ |. _2 ic |. rdata
        end.

The idea is to make an adverb (mine is audioconvert in wav.ijs) which
takes a format specification like bitpix and returns a conversion
function. Something like:

rawconvert =: 1 : 0
        select. u
          case. _32 do. _1&fc
          case. _64 do. _2&fc
          case. 16  do. _1&ic
          case. 32  do. _2&ic
          case. do. 'Invalid format' assert 0
        end.
)

which can then be called with

        cvt =. bitpix rawconvert  NB. A conversion function
        adata =. (naxis3,naxis2,naxis1)$ |. cvt |. rdata

Now the really cool thing is that J knows how to invert each of these
functions. So when you want to store the data again, you just use
(cvt^:_1)! This is a J-only technique: I don't know of any other
language that has the functional capabilities and the automatic
inversion required to do it.

There's always more to say but I think that's enough for one email. Good
luck with your future J projects!

Marshall

On Wed, Apr 18, 2018 at 02:47:02PM -0400, Thomas Hickey wrote:
>  As one of my first J programs, I've written a minimal viewer for FITS
> files, the standard way of sharing astronomical data, especially images.
>
> I'd appreciate it if someone could take a look at it and let me know how
> the code (80 lines of J) looks.  I know I've replaced some really awkward
> code in it, so I imagine there are better ways of doing lots of things.
>
> The code is on GitHub at  https://github.com/ThomasBHickey/JFits, along
> with a sample FITS file in the fitsSample folder.  Pointers on
organizing J
> code on GitHub would be appreciated as well.
>
> Or is there a better way to accomplish this than the programming forum?
I
> know this is longer than most of the questions that are posted here, but
I
> think it is an interesting application of J.
>
> --Th
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm

Reply via email to