Color Deconvolution at Wolfram Mathematica

This article was inspired by a recent article on monkey guts . Since the Chukchi is not a reader, the Chukchi is a writer, I decided to try to do this myself. Moreover, the task does not seem complicated and a lot of code is not required.

image

The simplest algorithm that comes to mind looks like this:
  • We define several basic colors of the picture. The RGB components of these colors will be used as basis vectors.
  • The color of each pixel is decomposed into a linear combination of basic ones.
  • Display an image for each base color.
  • Self-esteem is automatically increased.

Further, in more detail for each item.

Base colors


Basic colors are just the middle colors of the image segments. Segmentation algorithms - the sea is poured, choose to taste, some are implemented in Mathematica . I must say right away that little depends on the choice of segmentation algorithm, the main thing is that the number of clusters coincides with the desired number of components:

original = 
  Import["https://hsto.org/files/f42/962/1c5/f429621c55cd46a1beb4d4eb5aefed53.jpg"];
sizeX = First@ImageDimensions[original];
clusters =
  ClusteringComponents[
   original,
   3,
   Method -> "KMeans",
   DistanceFunction -> CosineDistance
   ];

Another simple function that we need minIndexis to return the position of the minimum element in the list, i.e., a compiled replacement Composition[First, Ordering]:
not for truth, but for truth
minIndex =
  Compile[
   {
    {list, _Real, 1}
    },
   Module[
    {
     i = 0,
     min = list[[1]],
     minPos = 1
     },
    Do[
     If[list[[i]] < min, minPos = i; min = list[[i]]],
     {i, 1, Length@list}
     ];
    minPos
    ],
   CompilationTarget -> "C"
   ];


Now a function that returns the basis vectors. removeDarkestremoves the darkest color - this is the background, we do not need it:

removeDarkest[list_] := Delete[list, minIndex[Norm /@ list]]
getBasisColors[image_, clusters_] :=
 Module[
  {
   data = ImageData[image],
   components = Union[Flatten @ clusters]
   },
  Table[
   Median[Join @@ Pick[data, clusters, component]],
   {component, components}
   ]
  ]

Launching ...

B1 = removeDarkest@getBasisColors[ColorNegate@original, clusters]

and voila:
image

And no, hey, not voila. It makes no sense to decompose a picture into such basic vectors: the resulting components will more or less repeat the clusters. Here is the moment of truth, i.e. the only inevitable heuristic step. In order to make the components of the picture more complete, you need to slightly expand the basis vectors.

If we want from the vector and subtract vector b , at the same time want to keep the components of the positive, then the maximum that we can make c  =  a  - min ( a i / b i )  b . Moreover, one of the components cwill go to zero. Naturally, we do not want to push so much apart. Just subtract a little the first from the second:

alpha = 0.5;
basis = {B1[[1]], B1[[2]] - alpha Min[B1[[2]]/B1[[1]]] B1[[1]]};
metric = Outer[Dot, basis, basis, 1];

You can, for example, subtract the average vector from all, it does not matter, the main thing is to somehow increase the angle between the vectors. Freedom in this step is the inevitable payment for the uncertainty of the task. The coefficient 0 <<  alpha 1 is the only fitting parameter of the algorithm. For example, in order to get a picture at the beginning of the post, I had to put alpha = 0.95.

The matrix metricis the metric of the linear span of the basis vectors; it is also the Gram matrix, which we will need in the next section.

Baseline decomposition


Baseline decomposition is a standard procedure. The only subtlety: from the statement of the problem it is clear that the expansion coefficients must be positive. The linear span of a set of vectors with positive coefficients is an infinite simplex (a pyramid whose vertex rests at the origin, and the base is carried to infinity). The faces of this simplex are simplexes of lower dimension. All we need is to project the given vector onto the simplex. If the vector falls inside the simplex (all expansion coefficients along the basis are positive), then the problem is solved, if not, then it is necessary to project on the face.

Here, in fact, the whole function:

It looks irony only because of my style of writing code.
reduceMetric =
  Compile[
   {
    {metric, _Real, 2},
    {index, _Integer}
    },
   Transpose@Delete[Transpose@Delete[metric, index], index],
   CompilationTarget -> "C"
   ];
getComponents =
 Compile[
  {
   {vec, _Real, 1},
   {basis, _Real, 2},
   {metric, _Real, 2}
   },
  Module[
   {
    covariant,
    contravariant = Table[0., {Length[basis]}],
    flag = True,
    subspace
    },
   covariant = basis.vec;
   flag =
    If[
     Det[metric] != 0,
     contravariant = Inverse[metric].covariant;
     FreeQ[Sign[contravariant], -1],
     False
     ];
   Chop@If[
     flag,
     contravariant,
     subspace =
      Table[
       Insert[
        getComponents[vec, Delete[basis, i], reduceMetric[metric, i]], 0., i],
       {i, 1, Length@basis}
       ];
     subspace[[minIndex[-(subspace.covariant)]]]
     ]
   ],
  {
   {_minIndex, _Integer},
   {_reduceMetric, _Real, 2},
   {_getComponents, _Real, 1}
   },
  CompilationTarget -> "C",
  RuntimeAttributes -> {Listable},
  Parallelization -> True
  ]


Check that it works correctly:

getComponents[{0.1, 0.9}.basis, basis, metric]
{0.1, 0.9}

Below are comments for those who want to understand the code. The expansion coefficients of the vector v along the basis e (i) are called the contravariant coordinates v  =  v i e (i) . First of all, we compute the covariant coordinates v i  = ( e (i) ,  v ). If the metric g ij  = ( e (i) ,  e (j) ) is invertible, then the contravariant coordinates are calculated by the inverse metric g ij = ( g −1 ) ij , i.e. v i =  g ij v j . That's all. The flag is flag = Falsegenerated in two cases: the metric is irreversible (a simplex of a smaller dimension than we thought), or at least one of the contravariant coordinates is negative (did not get inside the simplex). In these cases, we stupidly recursively sort through all the faces, and select the projection onto which is maximum. The projection onto the subspace of the basis e (i) is the convolution of covariant and contravariant coordinates v i  v i . Now for sure.

results


Although all functions are compiled in C, Mathematica stubbornly does not want to run getComponents on multiple cores. To understand laziness, let it remain on the conscience of the developers. So just grind all 2,736,000 pixels:

pixels = Join @@ ImageData[ColorNegate@original];
components =
   Table[
    getComponents[pixel, basis, metric],
    {pixel, pixels}
    ];

The first component:
data1 = (({1, 0} #).basis) & /@ components;
ColorNegate[Image@Partition[data1, sizeX]]

image

The second component:
data2 = (({0, 1} #).basis) & /@ components;
ColorNegate@Image@Partition[data2, sizeX]

image

Can I have more?


Can. Shove in getComponentsany, even degenerate / redundant basis, and you will be happy:
Module[
 {
  basis = {{1, 0}, {0, 1}, {1, 1}},
  metric
  },
 metric = Outer[Dot, basis, basis, 1];
 getComponents[{0.1, 0., 0.9}.basis, basis, metric]
 ]
{0.1, 0., 0.9}

You can download the Mathematica file here .

Also popular now: