[Math] Algorithm for unique “additive” color mixing

MATLAB

my question is
if i have only 3 basic colors (each made of rgb):

color1 : R:150, B:zero, G:255
color2 : R:255, B:150, G:zero
color3 : R:zero, B:255, G:150

 
They can be mixed using the formula :

new_color = floor(X*0.9)+floor(Y*0.1)

X and Y can be a basic color or a new color already created by using the formula.
 
for example, if i want to mix color1 as main with color3 :

new_color(R,B,G) = (floor(0.9*150)+floor(0.1*0) , floor(0.9*0)+floor(0.1*255) , floor(0.9*255)+floor(0.1*150) ) = (135, 25, 244).

 
I need to find a way to mix many colors in order to get a desired color, for example : R:187 B:135 G:201
so far i wrote a "brute force" program which go all over the combinations of basic colors (runing for 7 days now got up to 16 mixing steps) and a bit smarter AStar algorithm (got good reults for small sequnces, letting it run for the week end…).
Since i know some colors are made of more then 100 mixing steps, it will take me forever to go over all of the results…
 
Help me find a smarter and faster way to solve this.

Thanks.

Best Answer

(EDIT: I've completely rewritten this first algorithm for clarity)

If an approximate answer is good enough, I think there is a simple algorithm. First, consider the exact version of your problem: ExactMix(C1, C2) = 0.9 * C1 + 0.1 * C2.

We can represent any color as a unique linear combination of the basic colors:

target = a * color1 + b * color2 + c * color3

A slight modification that will simplify things is if we change from (R,G,B) coordinates to (x,y,z) coordinates by

(R,G,B) = x * color1 + y * color2 + (1 + z - x - y) * color3

In (x,y,z) coordinates, we have

color1 = (1,0,0)
color2 = (0,1,0)
color3 = (0,0,0)

and ExactMix has the same form in (x,y,z) coordinates as it does in (R,G,B)-coordinates.

The thing to note is that given any three colors in the z=0 plane, you can use ExactMix to form arbitrarily good approximations to any color within the triangle they define. The algorithm is essentially a binary search: if the color T lies in the triangle PQR, then let S=ExactMix(P,Q). Then either T lies in SQR or it lies in PSR.

By iterating this algorithm (and rotating which edge we split with ExactMix), we get an infinite sequence of triangles whose corners converge to T.

If we are so inclined, we might consider splitting the segment PQ into three segments: P-ExactMix(P,Q), ExactMix(P,Q)-ExactMix(Q,P), and ExactMix(Q,P)-Q rather than just two segments, and choose which of the three triangles T lies within.

Now, Mix introduces round-off error. We can still use the same binary search algorithm to get the (x,y) coordinates approximately right, but it will also spontaneously decrement the z coordinate as well.

If our target z coordinate is very negative, we can use the same binary search algorithm to get the (x,y) coordinates right, then keep doing calculations until round-off error gives us the correct z coordinate. (this might take some cleverness to avoid getting stuck in a situation where there is no round-off error)

If our target z coordinate is nearly zero, we can still run the algorithm and hope for the best, getting something that is at least close to our target color.

When implementing the inexact version, it is probably best to redo the coordinate change every time you pick a new triangle. That is, rather than always using the (x,y,z) coordinates defined by the initial three colors, you instead compute a new set of (x,y,z) coordinates based off of your current set of three colors.

For an exact answer, I expect exhaustive brute force to be feasible if done efficiently: a one-time calculation to compute every color you can possibly produce.

The natural data structure to hold the results would be to keep a structure that stores, for every color, one of the following three factoids:

  • I do not know how to make this color
  • This is a basic color
  • This color is produced by mix(color1, color2)

Then keep a second list of every color ever produced. Initialize it with your three basic colors, then execute the following algorithm:

  • i = 0
  • while i < length(array of generated colors)
    • for j = 0 to i
      • c = mix(color[i], color[j])
      • if c is marked as "I don't know how to make":
        • Append c to the list of generated colors
        • Mark c as produced by mix(color[i], color[j])
    • for j = 0 to i-1
      • c = mix(color[j], color[i])
      • if c is marked as "I don't know how to make":
        • Append c to the list of generated colors
        • Mark c as produced by mix(color[j], color[i])

Efficient computation of the mix function is essential. You can probably get some advantage from the fact that your inner loop is over j with i fixed (and the compiler might recognize that for you). Also, you will need an efficient data structure for the main data. For example, you could store the data as follows:

  • A 256x256x256 array of pair<color,color>
  • A 256x256x256 bit array to represent which entries are "do not know how to make"
  • A list of basic colors

The point of the bit array is that it will fit in cache (or nearly so), so that accessing it should be fast. Note that in the above algorithm you don't need the list of basic colors, so you could completely omit that part for your calculation.

It may or may not be slightly better to replace the first array with an array of 'pointers', and keep a second, flat array around of all of the pairs ever stored. This should result in a more compact data structure, since a large fraction of colors will never be produced, and may be faster to create. (less cache pollution)

If every color was reachable, then this algorithm would require 2^48 steps. Even if you could get the time of the inner loop down to a dozen ticks on a 1 GHz computer, it would take 36 days.

However, a large fraction of the colors are not reachable from the basic colors, and the algorithm designed above is meant to calculate only the entries that are reachable. I don't have the tools handy to estimate the fraction easily, but I fully expect the calculation to take only a few days with a very well optimized inner loop. A C/C++ implementation (or Fortran, I suppose) is probably a requirement for this calculation (unless it happens to be even smaller than I expect it to be).

Also, you could parallelize the algorithm, so as to use multiple cores on your CPU or multiple computers. It will be a little tricky to do so, as there are two obstacles:

  • You have to ensure that you iterate over all pairs of generated colors: so you have to combine the lists produced by different threads intelligently
  • You have to devise a scheme of having thread work independently for a while then synchronize results -- if multiple threads tried to update a shared data structure, you will get killed by the synchronization costs. (or get corrupt data if you forego synchronization)
Related Question