Segmentation of sunspots

Question

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]

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!

Answer 1

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] 

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}]}]