Mapping the point-spread function over the image

[Darth Takahashi]

Maybe there is, but doesn't he write Perl

That's actually a good point and maybe he's the only one who truly understands it. :D:D

[Darth Takahashi]

They keep it in Fortran because that is the language of choice for scientific/engineering computation. That was true in 1957 and it is even more true now, the latest revision of the standard includes parallelism constructs. And yes, I make and sell Fortran products :-)

Fortan is like black magic to me, you could tell me anything and I'd believe it. :D

I grew up with an ATARI400, then BBC-B before getting my first PC (486 with co-processor). Its only now that I have a MAC that I'm feeling more in love with computer again.

I have a working Sinclair QL and ATARI400 at home. Unfortunately, Fortan was never introduced on home computers to my knowledge? I do remember a promising language called FORTH???

[Darth Takahashi]

I help replace all the obselete **** in nuclear plant processing systems. It might be prehistoric but it is very robust and it just isn't replaceable on a like for like basis so it is huge headache to verify that any new systems are fully compliant and replicate the original systems.

I know, its the verification and validation process that's stopping them from converting this code base. That and the old adage "if it isn't broken don't fix it" right!!!

Some of these nuclear guys look at our products as though they are evil!!!

Anything that has a micro has code which is evil in the nuclear industry! We have to jump through a lot of loops to prove otherwise.

[Photosbykev]

I'm one of those nuclear guys lol, I've been in the industry for 38 years firstly with the CEGB and after numerous name changes now with Magnox Sites.

[michael.h.f.wilkinson]

Fortan is like black magic to me, you could tell me anything and I'd believe it. :D

I grew up with an ATARI400, then BBC-B before getting my first PC (486 with co-processor). Its only now that I have a MAC that I'm feeling more in love with computer again.

I have a working Sinclair QL and ATARI400 at home. Unfortunately, Fortan was never introduced on home computers to my knowledge? I do remember a promising language called FORTH???

Forth is a stack oriented language. "Stack oriented programming destroys more brain cells per hour than any other form of mental self abuse" :D(from Micro Cornucopia, about 1988)

Fortran is available on PCs in the GNU version.

[Macavity]

"Real programmers can write FORTRAN programs in *any* language"?

And something I came to (somewhat) rely on...

Not sure I ever fully grasped the "objectivity" thing - Probably an AGE thing? On the other hand, I suspect programming languages should be hostile enough to keep out the (some) uninitiates. LOL. I still wince that the published plots on my thesis topic were significantly wrong, because I had to defer to some "Old Prof's" BUG-ridden Fortran data analysis.

Hey, that's life? But I kinda hope (some of) these... "Higgs discoverers" know what they're doing.

[michael.h.f.wilkinson]

"Real programmers can write FORTRAN programs in *any* language"?

And something I came to (somewhat) rely on...

Not sure I ever fully grasped the "objectivity" thing - Probably an AGE thing? On the other hand, I suspect programming languages should be hostile enough to keep out the (some) uninitiates. LOL. I still wince that the published plots on my thesis topic were significantly wrong, because I had to defer to some "Old Prof's" BUG-ridden Fortran data analysis.

Hey, that's life? But I kinda hope (some of) these... "Higgs discoverers" know what they're doing.

You can write good, maintainable code in assembler (I have had to write some), just as you can write rotten code in any language. Some languages make life easier. Objects are just part of that. You do not think of data as static items on which you work, but of active entities which can perform certain tasks.

Objects are great for certain parts of a program, less for others. The image processing mainly relies on good old imperative programming, but object oriented (OO) concepts are needed (or at least very, very useful) for some parts. For example, in so called attribute filtering, in which you filter structures in the image based on their properties (or attributes), we first construct a tree structure (a component tree), which contains the structures in the image. OO allows you to view the tree as an object with certain functions you can apply to them. You do not have to know the internal structure of the tree to work with them. If I find a better way to construct the tree, I can insert that into the tree, without having to rewrite or adapt the rest of code.

We have a separate class of object called "AMT_Attribute" which computes some property of each image structure (tree node). A large collection different attributes exist, and they can be inserted into the program at will (two statements added : one to include the header file, one to insert the attribute in the appropriate menu in the GUI). We can do this because the GUI can ask them what their name is, what their minimum and maximum values are, etc.

In the past we built "objects" by making C structs containing pointers to functions (pointers, you have got to love pointers). It works well, but includes a lot of manual administrative stuff, which is error prone. OO takes away that kind of hassle.

[Darth Takahashi]

In the past we built "objects" by making C structs containing pointers to functions (pointers, you have got to love pointers). It works well, but includes a lot of manual administrative stuff, which is error prone. OO takes away that kind of hassle.

Pointer -> functions -> array-tables -> to more functions etc...

My head hurts!!!

I have seen some admittedly cleaver "Obfuscated Code" using pointers!!!

Absolute power corrupts absolutely!

Many years ago I started down the Java route because of this and its many other promises (applets). Then I ran into the anonymous inner class!!!

Because everything is a class. Interfaces are neat though.

I believe they are "virtual public functions" in C++

Every language has it's own issues.

[themos]

I do scientific and engineering programming in C/C++ (and MatLab). I worked with Fortran (77/90/95) myself and coming from Pascal/Modula/C background I HATE the lack of proper scope (don't talk to me about common blocks), and the lack of strict typing (still don't like that in C, but at least I MUST declare every variable type in C). The parallel constructs are available in C/C++ with OpenMP in a very transparent way.

I "fondly" remember getting an error message from a NAG routine (with its meaningful name D02BAE) in FORTRAN which read "Impossible error". This was due to a named common block being shared accidentally on a parallel machine, so different instances of the routine were overwriting each other's data. There was no workaround for that (AARGH). They had forgotten to compile the library on our Cray J932 (long since dead) with the --taskcommon switch, which makes private copies of common blocks.

Well, at least you got an error message telling you not to trust the results. We can be proud of that!

COMMON blocks were the source of a lot of headaches and we took a design decision to ban them from the Library.

By the way, we're still debating the wisdom of having descriptive names ( D02PCF vs nagf_ode_ivp_rk_range, say). I am against long names because when you're dealing with a product that has a history of 40 years and a future of as many again (I hope) you can easily "burn" good names when a new interface comes up and we'll end up coming up with ever more baroque names.

[michael.h.f.wilkinson]

Well, at least you got an error message telling you not to trust the results. We can be proud of that!

COMMON blocks were the source of a lot of headaches and we took a design decision to ban them from the Library.

By the way, we're still debating the wisdom of having descriptive names ( D02PCF vs nagf_ode_ivp_rk_range, say). I am against long names because when you're dealing with a product that has a history of 40 years and a future of as many again (I hope) you can easily "burn" good names when a new interface comes up and we'll end up coming up with ever more baroque names.

The NAG library was and is very useful, but I did end up writing my own positive scheme for stiff differential equations. Regarding the names: you get used to that, but in C a set of macros could easily be set up as a translation between old names and new (macros are so useful to enforce in-lining of code to save the overhead of function calls ).

Banning common blocks was a good decision, BTW.

One thing that Fortran supported far to late was recursion. A lot of my processing is tree based, and recursion is the way to handle trees in most cases.

[rfdesigner]

Forth is a stack oriented language. "Stack oriented programming destroys more brain cells per hour than any other form of mental self abuse" (from Micro Cornucopia, about 1988)

How does it comares to Occam? (for transputers)

that's a real head bender....

These days I kind of regard C as a bit high level, preferring to play assembly with my PICs

Derek

[michael.h.f.wilkinson]

How does it comares to Occam? (for transputers)

that's a real head bender....

These days I kind of regard C as a bit high level, preferring to play assembly with my PICs

Derek

Occam is OK, or at least no more difficult than other parallel languages. All parallel programming does your brain in at some point. Forth was worse even when used sequentially.

[Darth Takahashi]

Hum... will any of these parallel languages run on a PlayStation 3 cluster I wonder??? I guess there must be some GPL stuff out there?

[themos]

The NAG library was and is very useful, but I did end up writing my own positive scheme for stiff differential equations.

Didn't we have stiff solvers then? I just noticed that routine D02AGF is still in the Library and it was introduced at Mark 2, early 1970s!

Regarding the names: you get used to that, but in C a set of macros could easily be set up as a translation between old names and new (macros are so useful to enforce in-lining of code to save the overhead of function calls ).

One good thing about HPC library code is that most routines are expected to have so much "meat" in them that the cost of function calls can be ignored.

Banning common blocks was a good decision, BTW.

One thing that Fortran supported far to late was recursion. A lot of my processing is tree based, and recursion is the way to handle trees in most cases.

Yes, I can see that. Ouch.

[themos]

Anyway, enough hijacking. I'd love to see some image processing s/w that helps us analyze our optics and mounts.

[michael.h.f.wilkinson]

The stiff solver was present but would not run in combination with the Cray auto-tasking (similar impossible error ).

I will now stop and get on with some image processing.

[FraserClarke]

get on with some image processing.

I would suggest starting here;

AstrOmatic software – Astromatic.net

You might find a lot of the building blocks are already built...

[themos]

Learning never ends

[michael.h.f.wilkinson]

I would suggest starting here;

AstrOmatic software – Astromatic.net

You might find a lot of the building blocks are already built...

I know tha stuff (working with the Kapteyn Astronomical Institute on data reduction you cannot avoid it). I am actually trying to replace SExtractor with a hyperconnected morphological counterpart. We had a publication out in IEEE Trans. Pattern Analysis and Machine Intelligence (PAMI) which deals with the theory and a fast algorithm for such filters. The paper can be found here (click on PDF icon, free download, math alert!!). The method is quite general and the best results to date were on historical document processing, but the early astronomy results are encouraging. Most likely we will need a hybrid approach, combining some techniques from SExtractor with hyperconnected filters.

Some initial results are on:

Supplemental material on k-Flat Filtering

The PSFextractor is interesting, but as far as I know it estimates a single PSF, not a spatially variant map. Again, a combination of approaches might be optimal.

[dph1nm]

Michael

The Pan-STARRS project uses spatial variant PSF mapping (but mostly for convolution to a common PSF). They also analyse out of focus stellar images to map the optical defects of the telescope (possibly relying on the presence of the hole in the Cassegrain configuration). I believe it proved impossible to use in-focus images due there being too many potential variables going into constructing the PSF. Of course, the system is more complex than most amateur scopes, having assorted actuators on the primary and secondary mirrors.

I think the weak-lensing people also fit spatially variant PSF - see lensfit

NigelM

[michael.h.f.wilkinson]

Michael

The Pan-STARRS project uses spatial variant PSF mapping (but mostly for convolution to a common PSF). They also analyse out of focus stellar images to map the optical defects of the telescope (possibly relying on the presence of the hole in the Cassegrain configuration). I believe it proved impossible to use in-focus images due there being too many potential variables going into constructing the PSF. Of course, the system is more complex than most amateur scopes, having assorted actuators on the primary and secondary mirrors.

I think the weak-lensing people also fit spatially variant PSF - see lensfit

NigelM

Thanks for the information.

I have heard of several attempts to use spatially-variant PSF estimates (including Pan-STARRS). I think convolution to a common PSF will be essential, to regularity the result. I think my method would focus on estimating the central peak of the PSF, rather than the full PSF. None of the current methods use attribute filters, I am curious to see how they perform.

[NickK]

Hmm ok maybe I was to optimistic about the Belgium part ;-)

And you probably need a supercomputer or at least GPU acceleration.

AIP: deconvolution

Yves.

GPU acceleration is straight forward - OpenCL will give you a good cross platform approach. I worked on GPGPU around 2005-6 there's a few things to appreciate:

1. Precision/Error - the FP precision and errors need to be reviewed as GPUs have only just now started offering IEEE. Note that most GPUs are at 32bit for both mantissa and exponent with some short cuts for handling overflow/underflow.

2. The design of the algorithm needs data-dependency planning (both gathering, i.e. reading a[x,y], and scattering writes a'[x,y]

=value) in true parallel style. GPUs may process fast, however transferring over the PCI-Express bus is the bottleneck so keeping the data put is a core part of the process.

Techniques such as ping-pong etc should be designed in.

3. Memory requirements. I'm lucky as this laptop has 512MB of GPU memory and 8GB of RAM (and I'm considering upping to 16GB).. but most notepads will not have this.

I had to use Fortran 77 and 90 for numerical computation in my software engineering degree (specialised in parallel computation and networks), we also had todo "Formal Specification and Z" for 3 years.. nothing like the mental abuse of translating english into mathematics set theory and then formally proving everything holds true.. .. or that award may go to formally proving concurrency systems as part of Z..

I would advise using C/C++. I would also advise to make the jump away from writing threading in favour of blocks or tasks and let technology do the threading.

[michael.h.f.wilkinson]

I am in collaboration with a group in Sydney, Australia on GPU acceleration of these methods (in OpenCL). One problem is that the processing order is very much data driven (unlike many other image processing methods). This makes GPU acceleration hard. If we manage to do it, we may incorporate it in the project. For the first part I will just rely on multithreaded processing. It really needs to be hand coded, given the nature of this type of algorithm.

[NickK]

Interesting

http://www-group.slac.stanford.edu/kipac/OSU_PDFs_for_lsst_tutorial/peterson_lsst_tutorial.pdf

[michael.h.f.wilkinson]

Interesting

http://www-group.slac.stanford.edu/kipac/OSU_PDFs_for_lsst_tutorial/peterson_lsst_tutorial.pdf

Interesting read. They are really doing the reverse: given an optical system and atmosphere, compute the PSF.What I will try to do is simply measure the shape of the core of the PSF for stars scattered over the FOV. I still have to find a student (or two) for that.