16
$\begingroup$

Introduction

I have a greyscale photograph that I want to partition into areas of equal intensity. This means, that the integral over each partition should be (roughly) equal. There is no unique solution to this problem, but although I have some ideas on how to implement it, I am looking for any effective solutions for bitmap input.

1D-example

Let f be my signal, and g it's strictly positive, "padded" version:

f[x_] = 0.2 Cos[5 x] + Sin[x]; const = First[Minimize[f[x], {x}]]; g[x_] := 0.1 + f[x] - const; G[x_] = Integrate[g[s], {s, 0, x}];

My signal is now of total intensity

max = G[4 π];

and I can divide my y-axis into equal segments

k = 20; sols = Table[FindRoot[G[x] - y == 0, {x, y}], {y, Range[0, max, max/(k - 1)]}];

and by intersecting these segments with my integral curve, I can project to the x-axis and get my desired intervals:

Show[ { Plot[{g[x], G[x]}, {x, 0, 4 π}, PlotStyle -> {{Thick, Green}, Red}], Graphics[Point[{#, G[#]}] & /@ (x /. sols)], Graphics[Line[{{0, G[#]}, {#, G[#]}, {#, 0}} & /@ (x /. sols)]] } ]

Which gives me

Mathematica graphics

2D-case

I can imagine that this method generalises to 2D, but it is unclear how to implement it. Also, a big penalty comes from attempts at interpolating the image.

Another idea is to use Sasha's solution for sampling a probability distribution (in my case, the image itself):

RandomVariate from 2-dimensional probability distribution

Thereafter one can generate the Voronoi-diagram for the sample. However, I have been unable to create a feasible solution using this method as well.

It seems to me, that this is a problem better approached by a combination of image processing algorithms. After all, all my input is, is a bitmap. However, after researching this, I could not find a previous description of the problem (or a solution for that matter) in any of the image-processing communities.

Questions

  1. Does this problem have an actual name in the community already?

  2. Can I find an efficient way of clustering a signal, preferably working on the discrete data, into regions of similar intensity?

  3. (Bonus points) Can I generate grids with nice combinatorics? There is a nice Z-combinatorics for the 1D-example of course. Although a Z^2 grid for the 2D-case is too much to ask, perhaps there is a way to get a Quadtree-looking grid?

Improve this question
$\endgroup$
6
  • $\begingroup$ This question seems a better fit for the signal processing site $\endgroup$ Commented Oct 12, 2012 at 12:20
  • $\begingroup$ I checked it out, and though I partly agree, I just named this a "signal" in order to open up the question to more than image processing specialists. In addition, I am specifically looking for a Mathematica solution. $\endgroup$ Commented Oct 12, 2012 at 12:23
  • $\begingroup$ But you need an algorithm first ... I think. Anyway, perhaps someone can come up with something here, mine was only a suggestion $\endgroup$ Commented Oct 12, 2012 at 12:26
  • $\begingroup$ related stats.stackexchange.com/questions/8744/… $\endgroup$ Commented Oct 12, 2012 at 13:26
  • $\begingroup$ Related: Weighted Voronoi Stippling $\endgroup$ Commented Feb 8, 2016 at 2:10

2 Answers 2

Reset to default
7
$\begingroup$

The following seems to work at least for small grayscale images and a few points. My current computing power doesn't allow me to test it for larger examples. As this question was posed more than three years ago I decided to post this, despite the fact that I don't know how it scales.

It makes a Voronoi partition and minimizes the Variance of the region intensities. 

voronoi[pts_] := (* BoundedDiagram fails too often, so I'm using Vitaliy's nice http://mathematica.stackexchange.com/a/39828/193*) Binarize@ImageResize[ImageCrop[ ListDensityPlot[Append[#, 0] & /@ pts, InterpolationOrder -> 0, BoundaryStyle -> Black, Frame -> False, ColorFunction -> (White &)]], id] calcIntensities[{i_Image, voronoi_Image}] := (*probably can be made faster by using the Polygon from Voronoi instead of ComponentMeasurements *) ImageMeasurements[#, "TotalIntensity"] & /@ (ImageMultiply[i, mage@#] & /@ (ComponentMeasurements[voronoi, "Mask", CornerNeighbors -> False][[All, 2]])) obj[i_Image, pts_?(VectorQ[#, VectorQ[#, NumericQ] &] &)] := Variance@calcIntensities[{i, voronoi[pts]}]

Usage:

ed = ExampleData /@ ExampleData["AerialImage"][[1 ;; 5]]; imgs = ColorConvert[ImageResize[#, 100], "Grayscale"] & /@ ed; (*Three segments*) pts = Array[p, {3, 2}]; id = First@ImageDimensions[imgs[[1]]]; nm = NMinimize[{obj[#, pts], And @@ Thread[1 < Flatten@pts < id]}, Flatten@pts, Method -> "NelderMead"] & /@ imgs;

Result:

Mathematica graphics

Result table code

TableForm[{#[[1]], ImageMultiply@@#, Sequence@@{#, Variance@#} &@calcIntensities[#]} &/@ Transpose[{imgs, voronoi /@ (pts /. nm[[All, 2]])}], TableHeadings -> {{}, {"Original", "Segmentation", "Intensity Values", "Variance"}}]

Edit

The following incarnation is 50% faster:

voronoi[pts_] := Module[{g}, g = ListDensityPlot[Append[#, 0] & /@ pts, InterpolationOrder -> 0, Frame -> False, PlotRange -> {{0, id}, {0, id}}]; (Graphics[{FaceForm[White], Polygon@#}, Background -> Black, PlotRange -> {{1, id - 1}, {1, id - 1}}, PlotRangePadding -> 0, ImageSize -> id] & /@ Cases[Cases[g, _GraphicsComplex] // Normal, Polygon[l_, ___] :> l, Infinity])] calcIntensities[{i_Image, voronoi_}] := ImageMeasurements[#, "TotalIntensity"] & /@ (ImageMultiply[i, Image@#] & /@ voronoi) obj[i_Image, pts_?(VectorQ[#, VectorQ[#, NumericQ] &] &)] := Variance@calcIntensities[{i, voronoi[pts]}]

Usage:

i = ColorConvert[ImageResize[ExampleData[{"TestImage", "Mandrill"}], 100], "Grayscale"]; pts = Array[p, {3, 2}]; id = First@ImageDimensions@i; Timing[nm = NMinimize[{obj[i, pts], And @@ Thread[1 < Flatten@pts < id]}, Flatten@pts, Method -> "NelderMead"]] (ImageMultiply[i, #] & /@ voronoi[pts /. nm[[2]]]) ImageMeasurements[#, "TotalIntensity"] & /@ (ImageMultiply[i, #] & /@ voronoi[pts /. nm[[2]]])

Mathematica graphics

Fold[ImageAdd@## &, (ImageMultiply[i, #] & /@ voronoi[pts /. nm[[2]]])]

Mathematica graphics

Just for the record, here is 4-partitioned mandrill image. The performance is horrific:

Mathematica graphics

Improve this answer
$\endgroup$
7
$\begingroup$

The code in this answer implements a 2D version of the approach in the question. It is fairly fast but does not have great precision. The achieved precision might be seen as good enough.

Get image data (as in the set up of the answer by Dr. belisarius):

ed = ExampleData /@ ExampleData["AerialImage"][[1 ;; 5]]; imgs = ColorConvert[ImageResize[#, 500], "Grayscale"] & /@ ed; imgArr = Reverse[imgs[[1]] // ImageData];

The original approach can be re-written to use accumulated sums intestead of integration. The following function does that.

GetSplitCoords[vec_?VectorQ, args___] := GetSplitCoords[{vec}, args]; GetSplitCoords[vecs : {_?VectorQ ..}, qs_: Range[0, 1, 0.1]] := Block[{qvals, accVec = Total /@ Accumulate[Transpose[vecs]]}, qvals = Rescale[qs, {Min[qs], Max[qs]}, {0, Max[accVec]}]; Map[Block[{t = Abs[accVec - #]}, Position[t, Min[t]][[1, 1]]] &, qvals] ];

Here are partitioning coordinates:

qs = Range[0, 1, 0.2];

Find the separation curves in $X$ and $Y$ direction in an alternating manner:

rowCoords = Transpose@ MapIndexed[Thread[{#2[[1]], GetSplitCoords[#1, qs]}] &, Transpose[imgArr]]; colCoords = Transpose@ MapIndexed[Thread[{GetSplitCoords[#1, qs], #2[[1]]}] &, imgArr];

Instead of using the accumulated sum of a single vector to find a splitting point we can use the accumulated sum of several vectors:

pStep = 20; offset = Floor[(pStep - 1)/2]; rowCoords = Transpose@MapIndexed[Thread[{#2[[1]], GetSplitCoords[#1, qs]}] &, Join[Transpose[imgArr][[1 ;; offset]], Partition[Transpose[imgArr], pStep, 1], Transpose[imgArr][[-(pStep - offset - 1) ;; -1]]]]; colCoords = Transpose@MapIndexed[Thread[{GetSplitCoords[#1, qs], #2[[1]]}] &, Join[imgArr[[1 ;; offset]], Partition[imgArr, pStep, 1], imgArr[[-(pStep - offset - 1) ;; -1]]]];

Plot the found parition lines over the image:

Show[{ Graphics[Raster[imgArr]], ListLinePlot[colCoords, PlotStyle -> Yellow], ListLinePlot[rowCoords, PlotStyle -> Yellow]}, Frame -> True]

enter image description here

At this point we want to extract the segments of the image between the grid lines.

Find the $X$ segments:

xSegments = Table[ Thread[Thread[{k, Range[rowCoords[[i, k]][[2]], rowCoords[[i + 1, k]][[2]]]}] -> imgArr[[k, rowCoords[[i, k]][[2]] ;; rowCoords[[i + 1, k]][[2]]]]], {i, 1, Length[rowCoords] - 1}, {k, 1, Length[rowCoords[[i]]]}]; xSegments = SparseArray[Flatten[#], Dimensions[imgArr]] & /@ xSegments;

Find the $Y$ segments:

ySegments = Table[ Thread[Thread[{Range[colCoords[[i, k]][[1]], colCoords[[i + 1, k]][[1]]], k}] -> imgArr[[colCoords[[i, k]][[1]] ;; colCoords[[i + 1, k]][[1]], k]]], {i, 1, Length[colCoords] - 1}, {k, 1, Length[colCoords[[i]]]}]; ySegments = SparseArray[Flatten[#], Dimensions[imgArr]] & /@ ySegments;

Lets plot the $X$ segments:

Graphics[Raster[#]] & /@ xSegments

enter image description here

Precision of the $X$ segments:

Map[Total[#, 2] &, xSegments] (* {16752.8, 18665., 23587., 22058.2, 32367.5} *)

Precision of the $Y$ segments:

Map[Total[#, 2] &, ySegments] (* {31279.3, 24490., 20210.1, 22910.4, 14562.5} *)

Make filtering matrices to extract the grid segments:

xSegments01 = Map[Function[{sm}, SparseArray[Thread[Most[ArrayRules[sm]][[All, 1]] -> 1], Dimensions[sm]]], xSegments]; ySegments01 = Map[Function[{sm}, SparseArray[Thread[Most[ArrayRules[sm]][[All, 1]] -> 1], Dimensions[sm]]], ySegments]; filterMats = Flatten[Table[s1*s2, {s1, xSegments01}, {s2, ySegments01}], 1];

Extract the segments:

segments = Map[imgArr*# &, filterMats];

Precision of the intesity separation over the segments:

t = Map[Total[#, 2] &, segments]; ListPlot[Sort[t/Median[t]], GridLines -> {None, {1}}]

enter image description here

We can see the precision is not that great. With some images I get much better precision with some worse.

Here is what we get with a single vector accumulated sums (the obtained grid can be pretty jagged):

enter image description here

The

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.