Message Boards

WOLFRAM COMMUNITY

20156 Views

9 Replies

9 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Astronomy Image Processing Mathematica Graphics and Visualization

Computing FWHM from star image data

Dan Pollock

Dan Pollock, Retired

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. Dan

POSTED BY: Dan Pollock

9 Replies

Sort By:

David Keith

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 *) FindFit[indexedData, 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.}

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 *)
FindFit[indexedData, 
 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

Dan Pollock

Dan Pollock, Retired

Posted 10 years ago

Thanks that is a huge help. Dan

POSTED BY: Dan Pollock

Sam Carrettie

Sam Carrettie, Freelancer

Posted 10 years ago

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

Henrik Schachner

Henrik Schachner, Radiation Therapy Center, Weilheim, Germany

Posted 10 years ago

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: then imageSharpness /@ imgs ( Out: {6.2140,11.9367,25.72055} *) Maybe this might somehow helpful. Henrik

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

then

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

Maybe this might somehow helpful.

Henrik

POSTED BY: Henrik Schachner

Dan Pollock

Dan Pollock, Retired

Posted 10 years ago

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

POSTED BY: Dan Pollock

Dan Pollock

Dan Pollock, Retired

Posted 10 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. Dan

POSTED BY: Dan Pollock

Henrik Schachner

Henrik Schachner, Radiation Therapy Center, Weilheim, Germany

Posted 10 years ago

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

Marcelo De Cicco

Marcelo De Cicco, INMETRO

Posted 10 years ago

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

Dan Pollock

Dan Pollock, Retired

Posted 10 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. http://blog.wolfram.com/2014/08/14/fixing-bad-astrophotography-ii-imaging-mars-with-mathematica/ http://blog.wolfram.com/2010/12/27/fixing-bad-astrophotography-using-mathematica-8-and-advanced-image-deconvolution/ 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

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback