Skip to main content

image processing - Help find a bright object on Mars!


In today's news, scientists found a bright object on one of Curiosity's photos (it's near the bottom of the picture below). It's a bit tricky to find - I actually spent quite some time staring at the picture before I saw it.


Bright object


The question, then, is how one can systematically search for such anomalies. It should be harder than famous How do i find Waldo problem, as we do not necessarily know what we are looking for upfront!


Unfortunately, I know next to nothing about image processing. Playing with different Mathematica functions, I managed to find a transformation which makes the anomaly more visible at the third image after color separation -- but I knew what I was looking for already, so I played with the numerical parameter for Binarize until I found a value (0.55) that separated the bright object from the noise nicely. I'm wondering how can I do such analysis in a more systematic ways.


img = Import["http://www.nasa.gov/images/content/694809main_pia16225-43_946-710.jpg"];
Colorize @ MorphologicalComponents @ Binarize[#, .55] & /@ ColorSeparate[img]

enter image description here



Any pointers would be much appreciated!



Answer



Here's another, slightly more scientific method. One that works for many kinds of anomalies (darker, brighter, different hue, different saturation).


First, I use a part of the image that only contains sand as my training set (I use the high-res image from the NASA site instead of the one linked in the question. The results are similar, but I get much saner probabilities without the JPEG artifacts):


img = Import["http://www.nasa.gov/images/content/694811main_pia16225-43_full.jpg"];
sandSample = ImageTake[img, {0, 200}, {1000, 1200}]

enter image description here


We can visualize the distribution of the R/G channels in this sample:


SmoothHistogram3D[sandPixels[[All, {1, 2}]], Automatic, "PDF",  AxesLabel -> {"R", "G", "PDF"}]


enter image description here


The histogram looks a bit skewed, but it's close enough to treat it as gaussian. So I'll assume for simplicity that the "sand" texture is a gaussian random variable where each pixel is independent. Then I can estimate it's distribution like this:


sandPixels = Flatten[ImageData[sandSample], 1];
dist = MultinormalDistribution[{mR, mG, mB}, {{sRR, sRG, sRB}, {sRG, sGG, SGB}, {sRB, sGB, sBB}}];
edist = EstimatedDistribution[sandPixels, dist];
logPdf = PowerExpand@Log@PDF[edist, {r, g, b}]

Now I can just apply the PDF of this distribution to the complete image (I use the Log PDF to prevent overflows/underflows):


rgb = ImageData /@ ColorSeparate[GaussianFilter[img, 3]];

p = logPdf /. {r -> rgb[[1]], g -> rgb[[2]], b -> rgb[[3]]};

We can visualize the negative log PDF with an appropriate scaling factor:


Image[-p/20]

enter image description here


Here we can see:



  • The sand areas are dark - these pixels fit the estimated distribution from the sand sample

  • Most of the Curiosity area in the image is very bright - it's very unlikely that these pixels are from the same distribution


  • The shadows of the Curiosity probe are gray - they're not from the same distribution as the sand sample, but still closer than the anomaly

  • The anomaly we're looking is very bright - It can be detected easily


To find the sand/non-sand areas, I use MorphologicalBinarize. For the sand pixels, the PDF is > 0 everywhere, for the anomaly pixels, it's < 0, so finding a threshold isn't very hard.


bin = MorphologicalBinarize[Image[-p], {0, 10}]

enter image description here


Here, areas where the Log[PDF] < -10 are selected. PDF < e^-10 is very unlikely, so you won't have to check too many false positives.


Final step: find connected components, ignoring components above 10000 Pixels (that's the rover) and mark them in the image:


components = 

 ComponentMeasurements[bin, {"Area", "Centroid", "CaliperLength"}, 
   10 < #1 < 10000 &][[All, 2]]
Show[Image[img], 
 Graphics[{Red, AbsoluteThickness[5], Circle[#[[2]], 2 #[[3]]] & /@ components}]]

enter image description here


Obviously, the assumption that "sand pixels" are independent gaussian random variables is a gross oversimplification, but the general method would work for other distributions as well. Also, r/g/b values alone are probably not the best features to find alien objects. Normally you'd use more features (e.g. a set of Gabor filters)


Comments

Post a Comment

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more