python performance tips for looping calculation

Question

Python 2.7, Windows 7.

I'm looking for tips on how to make a calculation heavy script run faster. First an idea of what I'm doing:

Starting with a given color, I want to generate a list of 30 more colors (rgb values) that are maximally distinctive to the human eye from one another, and with the front of the list more distinctive than the end.

Currently, I estimate that the script will take ~48 hours to complete. I could let it run over the weekend, but I figured I would take the opportunity to learn something about python performance.

An overview what the code does:

gen_colours() contains a loop that runs 30 times. Each time 4 processes run multi(n, l, g) which contains the big loop iterating over each r, g, and b value between 0 and 255 (r value is split between processes so it loops 64 times). The inner most loop contains another loop that checks the rgb value against rgb values already found by calling compute_dist([r, g, b], c).

Anyways, without completely restructuring my code, things to help speed it up would be cool. Also, running all four cpus at max for 48 hours...issues there?

Code:

from math import sqrt, pow, atan2, atan, sin, cos, exp, radians, degrees
from fractions import Fraction
import time
import multiprocessing


def to_xyz(rgb):
    r = rgb[0] / 255.0
    g = rgb[1] / 255.0
    b = rgb[2] / 255.0

    f = Fraction(12, 5)

    if r > 0.04045:
        r = ((r + 0.055) / 1.055) ** f
    else:
        r /= 12.92
    if g > 0.04045:
        g = ((g + 0.055) / 1.055) ** f
    else:
        g /= 12.92
    if b > 0.04045:
        b = ((b + 0.055) / 1.055) ** f
    else:
        b /= 12.92

    r *= 100
    g *= 100
    b *= 100

    # Observer = 2 degrees, Illuminant = D65
    x = r * 0.4124 + g * 0.3576 + b * 0.1805
    y = r * 0.2126 + g * 0.7152 + b * 0.0722
    z = r * 0.0193 + g * 0.1192 + b * 0.9505

    return [x, y, z]


def to_lab(xyz):

    x = xyz[0]
    y = xyz[1]
    z = xyz[2]

    # Observer= 2deg, Illuminant= D65
    x /= 95.047
    y /= 100.0
    z /= 108.883

    f = Fraction(1, 3)

    if x > 0.008856:
        x **= f
    else:
        x = 7.787 * x + 0.13793103448
    if y > 0.008856:
        y **= f
    else:
        y = 7.787 * y + 0.13793103448
    if z > 0.008856:
        z **= f
    else:
        z = 7.787 * z + 0.13793103448

    L = 116 * y - 16
    a = 500 * (x - y)
    b = 200 * (y - z)

    return [L, a, b]


def compute_dist(rgb1, rgb2):
    """ Compute the apparent difference in colours using CIEDE2000 standards """

    xyz1 = to_xyz(rgb1)
    xyz2 = to_xyz(rgb2)

    lab1 = to_lab(xyz1)
    lab2 = to_lab(xyz2)

    a1 = lab1[1]
    a2 = lab2[1]
    b1 = lab1[2]
    b2 = lab2[2]
    L1 = lab1[0]
    L2 = lab2[0]

    c1 = sqrt(a1 * a1 + b1 * b1)
    c2 = sqrt(a2 * a2 + b2 * b2)
    c = (c1 + c2) / 2
    crs = c ** 7
    x = 0.5 - 0.5 * sqrt(crs / (crs + 6103515625))
    temp = (1 + x) * a1
    c1 = sqrt(temp * temp + b1 * b1)
    h1 = hue(temp, b1)
    temp = (1 + x) * a2
    c2 = sqrt(temp * temp + b2 * b2)
    h2 = hue(temp, b2)
    dL = L2 - L1
    dc = c2 - c1
    if c1 * c2 == 0:
        dh = 0
    else:
        temp = round(h2 - h1, 12)
        if abs(temp) <= 180:
            dh = h2 - h1
        else:
            if temp > 180:
                dh = h2 - h1 - 360
            else:
                dh = h2 - h1 + 360
    dh = sqrt(c1 * c2) * sin(radians(dh / 2))
    dh += dh
    lav = (L1 + L2) / 2
    cav = (c1 + c2) / 2
    if c1 * c2 == 0:
       htot = h1 + h2
    else:
        temp = abs(round(h1 - h2, 12))
        if temp > 180:
            if h2 + h1 < 360:
                htot = h1 + h2 + 360
            else:
                htot = h1 + h2 - 360
        else:
            htot = h1 + h2
        htot /= 2

    T = 1 - 0.17 * cos(radians(htot - 30)) + 0.24 * cos(radians(2 * htot)) + 0.32 * cos(radians(3 * htot + 6)) - 0.20 * cos(radians(4 * htot - 63))
    htotdtme = (htot / 25) - 11
    xPH = 30 * exp(-htotdtme * htotdtme)
    cavrs = cav ** 7
    scocp = sqrt(cavrs / (cavrs + 6103515625))
    xRC = scocp + scocp
    lavmf = lav - 50
    lavmfs = lavmf * lavmf
    SL = 1 + 0.015 * lavmfs / sqrt(20 + lavmfs)
    SC = 1 + 0.045 * cav
    SH = 1 + 0.015 * cav * T
    RT = -sin(radians(xPH + xPH)) * xRC
    dL /= SL
    dc /= SC
    dh /= SH
    dE = sqrt(dL * dL + dc * dc + dh * dh + RT * dc * dh)

    return dE


def hue(a, b):          # Function returns CIELAB-Hue value

    c = 0
    if a >= 0 and b == 0:
        return 0
    if a < 0 and b == 0:
        return 180
    if a == 0 and b > 0:
        return 90
    if a == 0 and b < 0:
        return 270
    if a > 0 and b > 0:
        c = 0
    elif a < 0:
        c = 180
    elif b < 0:
        c = 360
    return degrees(atan(b / a)) + c


def multi(p, l, q):
    f = 0
    n = []

    s = p * 64
    e = (p + 1) * 64
    for r in xrange(s, e):
        for g in xrange(256):
            for b in xrange(256):
                s = 1000  # smallest dist

                for c in l:  # compare to existing colours

                    d = compute_dist([r, g, b], c)
                    if d < s:
                        s = d

                if s > f:
                    n = [r, g, b]
                    f = s

    q.put(f)
    q.put(n)


def gen_colours(start_col=[68, 68, 68]):

    out = open('colour_output.txt', 'w')
    l = [start_col]

    if __name__ == '__main__':

        q0 = multiprocessing.Queue()
        q1 = multiprocessing.Queue()
        q2 = multiprocessing.Queue()
        q3 = multiprocessing.Queue()

        for h in xrange(30):  # create 30 more colours

            p0 = multiprocessing.Process(target=multi, args=[0, l, q0])
            p1 = multiprocessing.Process(target=multi, args=[1, l, q1])
            p2 = multiprocessing.Process(target=multi, args=[2, l, q2])
            p3 = multiprocessing.Process(target=multi, args=[3, l, q3])

            p0.start()
            p1.start()
            p2.start()
            p3.start()

            p0.join()
            p1.join()
            p2.join()
            p3.join()

            d0 = q0.get()
            d1 = q1.get()
            d2 = q2.get()
            d3 = q3.get()
            c0 = q0.get()
            c1 = q1.get()
            c2 = q2.get()
            c3 = q3.get()

            d = [d0, d1, d2, d3]
            c = [c0, c1, c2, c3]

            m = max(d)
            i = d.index(m)
            n = c[i]

            l.append(n)
            out.write("[" + str(n[0]) + ", " + str(n[1]) + ", " + str(n[2]) + "]\n")
            print "\nnew colour added: " + str(l)

        out.close()
        print "Done"


gen_colours()

Any tips?

Edit:

An obvious improvement is the fact that I'm calculating Lab values on found rgb colors every time. I added a list to store Lab values for these so that it doesn't need to do this each loop. This reduced time by about 1/4. Not a Python performance improvement that I'm looking for however.

Answer 1

You're doing way too much work.

It appears that you are working in a 24-bit RGB color-space when most monitor/gamut/ambient/eye combinations afford far less discriminability than your CIE calculations produce.

Assuming you are doing this to provide real-world colors for real-eyes, you also have to account for the myriad forms of colorblindness which reduces you to less than 12-bit useful color-space. Even if we are just talking about luminance, as luminance approaches the lower third of the device gamut, noticable differences become ever sparser.

The problem that you have is algorithmic; that is, you are working too hard to get detailed results when the added detail is irrelevant. Between #000 and #fff, there are only 4096 possible colors and the red-green axis can be rejected out of hand.

Answer 2

I'm sure a color that is R:100 G:100 B:101 would not be a "maximally distinctive" solution if color R:100 G:100 B:100 is chosen already.

One quick improvement you could make is to omit checking colors which are similar (ie. R and G values which are the same that have a B value within a given range).