Message Boards Message Boards

9 Replies
9 Total Likes
View groups...
Share this post:

Computing FWHM from star image data

Posted 10 years ago

I'm trying to use Mathematica to drive an autofocus system for a telescope and ccd camera. Three calculations are commonly used FWHM, HFD (half flux diameter) and sigma. So I need to be able to calculate those from a small (100x100) crop from an image produced by my tethering software. I found a Mathematica program to separate the stars out of the small image that works great so that is all sorted. As I understand it you do a gaussian distribution on the pixel data then calculate the two values to use to figure out the direction and when you are in focus. I'm new to Mathematica discovered it when I bought a Raspberry Pi and since I can buy a copy for about half what commercial software to do autofocus costs I thought it would be the way to go as I can use Mathematica for so much more. But I'm not having any luck so far doing the GD. So any clues tips or links are very much appreciated.


POSTED BY: Dan Pollock
9 Replies
Posted 9 years ago

Ran the first cuda parallel code on the tegra stack last night. Impressive. I wish Mathematica was on those boards. The articles I read for the image processing are these.

You will have to play around with it a bit to get it working. I'll have to do it again myself as I've managed to loose what I had working. I keep forgetting to do a shutdown on the pi's and ruining the SD image. Now I'm trying to do it with OpenCV which is much easier for me to get my head around. It has a function where you click on a star and it will cut out what ever size you want. That is going to be more useful for what I want to do.

POSTED BY: Dan Pollock

Hi Dan,

We mentioned a mathematica program that separate the stars out a small image. Would you mind share it here? I am very interested using mathematica codes in astrophotography.

POSTED BY: Marcelo De Cicco
Posted 9 years ago

Henrik I've taken your idea and ran with it. Put together a cluster of Jetson Tegra K1 dev boards and I'm using the cuFFT library that comes with it to do the FFT's. Its quick. On a single board a 1024x1024 convolution takes 22ms. On a 5in square board that happily runs on a 12v battery. The OpenCV lib provides all the mechanisms to get the fits file into the matrix for the cuFFT routines to operate on. Eventually I'll use one of the Tegras on my scope and one or more on a vision robot I'm working on. Actually giving my mobility scooter the ability to self drive is the goal. I started with a cluster of Pi 2 boards but the Tegras are much better for math. The Raspi cluster is happily running Univ of Hamburgs climate model now but will eventually have a Hadoop cluster running on it. So thanks for that idea its turning out to be an excellent way to get this accomplished.


POSTED BY: Dan Pollock

Hi Dan,

thank you very much for your nice response! I am glad to hear that my little hint was helpful to your project, which really impresses me!

Regards -- Henrik

POSTED BY: Henrik Schachner
Posted 9 years ago

Thanks for the info. I missed these in the deluge of daily e-mail. I'll have to try it out.


POSTED BY: Dan Pollock

Hi Dan,

I do not have any experience in astrophotography and I hope my little remark is not completely off topic: As far as I understand your question, it is about judging sharpness of an image. Here fouriertransform might be an option: The sharper an image, the more spread out are the corresponding fourier values, which can be calculated in terms of standard deviation. This should be an easy, fast and robust method - to be used not just for star images. So I wrote a little function (should be self explaining):

imageSharpness[img_Image] := 
 Module[{imgData, fouData, qRows, qColumns, qData0, qData1},
  imgData = ImageData[ColorConvert[img, "Grayscale"]];
  (* dimensions of one quarter: *)
  fouData = Abs@Fourier[imgData];
  {qRows, qColumns} = Dimensions[fouData]/2;
  (* get the data of the first quarter: *)      
  qData0 = fouData[[;; qRows, ;; qColumns]];
  (* get indizes, remove first line and first column and flatten: *)     
   qData1 = Flatten[MapIndexed[{#2, #1} &, qData0, {2}][[2 ;;, 2 ;;]], 1];
  (* return something like "standard deviation": *)      
  Mean[(#1[[1]]^2 + #1[[2]]^2) #2 & @@@ qData1]

Lets try it out:

enter image description here


imageSharpness /@ imgs
(* Out:   {6.2140,11.9367,25.72055}  *)

Maybe this might somehow helpful.


POSTED BY: Henrik Schachner

Dan, I am not sure if you have seen this, but there is a nice blog post by Tom Sherlock, perhaps interesting to you:

Serial Interface Control of Astronomical Telescopes

POSTED BY: Sam Carrettie
Posted 10 years ago

Thanks that is a huge help.


POSTED BY: Dan Pollock
Posted 10 years ago

Here is one way:

In[1]:= (* make some data *)
data = Table[
   rSquared = x^2 + y^2;
   3. Exp[-rSquared/5^2],
   {x, -10, 10}, {y, -10, 10}

In[2]:= (* generate an {x,y,z} list *)
indexedData = 
  Flatten[Table[{x, y, data[[x, y]]}, {x, Length[data]}, {y, 
     Length[data[[1]]]}], 1];

In[3]:= (* fit to Gaussian to recover parameters *)
 i0 Exp[-((x - x0)^2 + (y - y0)^2)/a^2], {i0, x0, y0, a}, {x, y}]

Out[3]= {i0 -> 3., x0 -> 11., y0 -> 11., a -> 5.}
POSTED BY: David Keith
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract