20
$\begingroup$

I recently went for a walk with my camera and snapped a few pictures of the sun. The pictures didn't turn out too well, but this was to be expected without a telescope. I thought it might be a nice exercise to segment out the sunspots and see how my results compare to the data of SOHO. Here is what I came up with:

img = Import["http://i.imgur.com/4lDwE33.jpg"];
soho = Import["http://i.imgur.com/XNmroTf.jpg"];

(*determine center and radius of the sun in the the image*)
{center, radius} = 
  ComponentMeasurements[
    MorphologicalComponents[img, 0.2], {"Centroid", 
     "EquivalentDiskRadius"}][[1, 2]];

(*segment the sun from background, polar transformp, post process *)
ImageTransformation[img, 
  center + {Cos[#[[1]]], Sin[#[[1]]]}*#[[2]] &, {2 Pi radius, radius},
   DataRange -> Full, PlotRange -> {{0, 360 °}, {1, radius}}];
polar = ColorConvert[%, "LUV"];
mean = Mean@ImageData[polar][[All, All, 1]][[#]] & /@ 
   Range[Round[radius]];
normalPolar = 
  Image@Rescale[
    ImageData[polar][[All, All, 1]][[#]] - mean[[#]] & /@ 
     Range[Round[radius]]];

(*transform back to cartesien coordinates*)
normalCart = 
  ImageTransformation[
   normalPolar, {ArcTan @@ (radius - #), 
     Norm[# - radius]} &, {2 radius, 2 radius}, 
   DataRange -> {{-180 °, 180 °}, {0, radius}}, 
   PlotRange -> {{0, 2 radius}, {0, 2 radius}}];
(*segmentation of sun spots*)
segmentation = 
  ColorNegate@
   ImageAdd[
    Dilation[
     Binarize[ImageCorrelate[%, GaussianMatrix[6], CosineDistance], 
      0.253], 2.29], Graphics[Disk[center, radius]]];
(*generate spot labels*)
centroids = 
  ComponentMeasurements[
   SelectComponents[MorphologicalComponents[ColorNegate@segmentation],
     "Count", 10], {"Centroid", "EquivalentDiskRadius", "Label"}];
labeled = Show[normalCart,
   Graphics[{Red, Circle @@ # & /@ centroids[[All, 2, 1 ;; 2]],
     MapThread[
      Text, {centroids[[All, 2, 3]] - 1, {0, 6 + #[[2]]} + #[[1]] & /@
         centroids[[All, 2, 1 ;; 2]]}]}]
   ];

(* plot:
"original image", "post processed image"
"segmented sun spots", "labeled post processed image"
"image from Solar and Heliospheric Observatory (SOHO)"
*)
GraphicsGrid[{{ImageCrop[img], normalCart}, {segmentation, 
   labeled}, {soho}}, Frame -> All, ImageSize -> 900]

enter image description here enter image description here

While I'm quite happy with the result, I thought some of you guys might come up with some clever ideas to make the approach more robust and the results even better.

I'm also very thankful for a review of my code. You never stop learning!

$\endgroup$
  • $\begingroup$ Nice question. The part of the code that separates the sun from the background can be replaced by ImageCrop@ImageMultiply[Binarize@img, img] I think. $\endgroup$ – C. E. Sep 2 '14 at 10:49
  • 2
    $\begingroup$ What is the problem with MinDetect[img, .02]? $\endgroup$ – Dr. belisarius Sep 3 '14 at 5:33
14
$\begingroup$

Well I decided to give it a bit of a go...First import the image and convert to grayscale, then crop to focus on the area of interest. Then I used a LaplacianGaussianFilter, which is often used in blob detection.

img = ImageAdjust@ColorConvert[Import["http://i.imgur.com/4lDwE33.jpg"], "Grayscale"];
smallimg = ImageAdjust@ImageTake[img, {200, 500}, {200, 600}];
logimg = ImageAdjust@LaplacianGaussianFilter[smallimg, {6,2}];

(* pretty colours! *)
coloredlogimg = Colorize[ColorNegate@logimg, ColorFunction -> "SunsetColors"];
GraphicsRow[{smallimg, coloredlogimg}, ImageSize -> Full] 

enter image description here

The interesting thing about this is that it starts to reveal the separate sunspots right in the centre of the image, which are distinct in the SOHO comparison image but less clear in your photo. This is what initially drew me to the question!

I then proceed to binarize the image (with an automatic threshold), do a distance transform followed by maxima detection, then finally use WatershedComponents to pick out the sunspots:

binimg = Binarize[logimg, Method -> "MinimumError"];
distimg = ImageAdjust@DistanceTransform@binimg;
maximg = MaxDetect[distimg, 0.2];
components = DeleteBorderComponents[
   WatershedComponents[logimg, maximg, Method -> "Rainfall"]];
measures = ComponentMeasurements[
   components, {"Centroid", "EquivalentDiskRadius", "Label"}];

g1 = Show[coloredlogimg, 
 Graphics[{Red,
   Map[Circle[#, 10] &, measures[[All, 2, 1]]], Blue, 
   MapThread[Text, {measures[[All, 2, 3]], 
     measures[[All, 2, 1]] + ConstantArray[{-5, 15}, Length@measures[[All, 2, 1]]]
      }]
    }]];

End result compared to the SOHO image? Pretty good! It's managed to separate out the central spot into two components (compared to the 3 evident in the SOHO image).

I did try to minimize the amount of user input needed to select appropriate parameters, but sadly it's affected by image noise. I also gave image deconvolution a try, but didn't get very far.

soho = ImageAdjust@Import["http://i.imgur.com/XNmroTf.jpg"]
GraphicsColumn[{g1, ImageAdjust@ImageTake[soho, {175, 375}, {130, 400}]}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.