Make an Lch histogram from an image. Try `lch_histo.py --wheel` and then `lch_histo.py --histofancy image.png wheel.png`.

lch_histo.py

#!/usr/bin/env python3

'''
Author: Yotam Gingold <yotam (strudel) yotamgingold.com>
License: Public Domain [CC0](http://creativecommons.org/publicdomain/zero/1.0/)
On GitHub as a gist: https://gist.github.com/yig/2de8f832b96f6ddffd788bb9c3e72ccf
'''

## pip install scikit-image drawsvg

import numpy as np
from numpy import *
import skimage.io
from skimage import color
import os

## For reproducibility
random.seed(1)

LUMINANCE = 60.
CHROMA = 33.

VERIFY = True

def generate_one_random_lch():
    '''
    Lch has Luminance in the range [0,100],
    chroma in the range [0,100],
    and hue in radians [0,2*pi].
    '''
    
    return ( LUMINANCE, CHROMA, 2*pi*random.random() )

def lch2rgb( lch ):
    '''
    Given a color in Lch color space as a tuple of three numbers,
    returns a color in sRGB color space as a tuple of three numbers.
    '''
    assert len( lch ) == 3
    assert float( lch[0] ) == lch[0]
    assert float( lch[1] ) == lch[1]
    assert float( lch[2] ) == lch[2]
    
    if VERIFY and not lch_color_is_in_sRGB_gamut( lch ):
        print( "Color out of sRGB gamut!" )
    
    return color.lab2rgb( color.lch2lab( asfarray(lch)[None,None,:] ) )[0,0]

def lch_color_is_in_sRGB_gamut( lch ):
    '''
    Given a color in Lch color space as a tuple of three numbers,
    returns True if the color can be converted to sRGB and False otherwise.
    '''
    
    c_lab = color.lch2lab( asfarray(lch)[None,None,:] )
    c_lab2 = color.rgb2lab( color.lab2rgb( c_lab ) )
    return abs( c_lab - c_lab2 ).max() < 1e-5

# def hue_too_close( hue1, hue2 ):
#    return lch_too_close( ( LUMINANCE, CHROMA, hue1 ), ( LUMINANCE, CHROMA, hue2 ) )
def lch_too_close( col1, col2 ):
    ## 2.3 is JND in Lab.
    return linalg.norm( lch2lab(col1) - lch2lab(col2) ) <= 2.3*10
def lch2lab( lch ): return color.lch2lab( asfarray(lch)[None,None,:] )[0,0]

def sort_lch_colors_relative_to_first( seq ):
    seq = asfarray( seq )
    hue0 = seq[0][2]
    for i in range(1,len(seq)):
        while seq[i][2] < hue0: seq[i][2] += 2*pi
    
    seq = seq.tolist()
    seq_sorted = seq[:1] + sorted( seq[1:], key = lambda x: x[2] )
    ## mod the numbers back
    seq_sorted = [ ( c[0], c[1], fmod( c[2], 2*pi ) ) for c in seq_sorted ]
    return seq_sorted

def ij2xy( i, j, RESOLUTION ):
    '''
    Given:
        i,j: Row, column integer pixel coordinates in a square image.
        RESOLUTION: The width and height of the image.
    Returns:
        x, y: Floating point coordinates in [-1,1]^2.
    '''
    center = (RESOLUTION-1)/2
    # return ( i - center )/center, ( j - center )/center
    return j/center - 1, i/center - 1

def xy2ij( x, y, RESOLUTION ):
    '''
    Given:
        x,y: Floating point coordinates in [-1,1]^2.
        RESOLUTION: The width and height of the image.
    Returns:
        i,j: Row, column integer coordinates in the range [0,RESOLUTION]^2.
    '''
    center = (RESOLUTION-1)/2
    # return x * center + center, y * center + center
    return int(round( ( y + 1 ) * center )), int(round( ( x + 1 ) * center ))

def generate_wheel_as_circle():
    '''
    Returns an sRGB image of all hues around the Lch color wheel.
    '''
    
    ## Resolution should be even.
    ## If resolution is odd, there will be a weird color at the center.
    RESOLUTION = 1000
    img = zeros( (RESOLUTION,RESOLUTION,4), dtype = float )
    for i in range(RESOLUTION):
        for j in range(RESOLUTION):
            ## Convert to coordinates with the origin in the center of the image
            x, y = ij2xy( i,j, RESOLUTION )
            ## Get the hue as an angle in radians
            hue_radians = ( arctan2( y, x ) + 2*pi ) % (2*pi)
            ## Create an lch color
            lch = ( LUMINANCE, CHROMA, hue_radians )
            ## Skip this pixel if out of gamut
            # if not lch_color_is_in_sRGB_gamut( lch ): continue
            ## Otherwise, get the RGB color and write it to the image
            img[i,j,:3] = lch2rgb( lch )
            ## Make the pixel opaque
            img[i,j,3] = 1
    
    return img

def generate_lch_slice():
    '''
    Returns an sRGB image of all hues and chromas.
    '''
    
    RESOLUTION = 1000
    CLIP_TO_CIRCLE = False
    
    img = zeros( (RESOLUTION,RESOLUTION,4), dtype = float )
    for i in range(RESOLUTION):
        for j in range(RESOLUTION):
            ## Convert to coordinates with the origin in the center of the image
            x, y = ij2xy( i,j, RESOLUTION )
            ## Get the hue as an angle in radians
            hue_radians = ( arctan2( y, x ) + 2*pi ) % (2*pi)
            ## Chroma is the magnitude scaled (and clipped?) to [0,100]
            chroma = sqrt( x**2 + y**2 )*100
            if CLIP_TO_CIRCLE and chroma > 100: continue
            ## Create an lch color
            lch = ( LUMINANCE, chroma, hue_radians )
            ## Skip this pixel if out of gamut
            if not lch_color_is_in_sRGB_gamut( lch ): continue
            ## Otherwise, get the RGB color and write it to the image
            img[i,j,:3] = lch2rgb( lch )
            ## Make the pixel opaque
            img[i,j,3] = 1
            
    return img

def generate_wheel_as_strip():
    '''
    Returns an sRGB image of all hues around the Lch color wheel.
    '''
    
    out_of_gamut_count = 0
    
    img = zeros( (50,360,3), dtype = float )
    c = asfarray(( LUMINANCE, CHROMA, 0.0 ))
    for col, h in enumerate( linspace( 0, 2*pi, 360 ) ):
        c[2] = h
        img[:,col,:] = asfarray(lch2rgb(c))[None,:]
        
        ## verify
        if not lch_color_is_in_sRGB_gamut( c ):
            print( "Color out of sRGB gamut!" )
            out_of_gamut_count += 1
    
    print( "Wheel: %s colors out of gamut (%.2f%%)" % ( out_of_gamut_count, (100.*out_of_gamut_count) / 360 ) )
    
    return img

def save_image_to_file( img, filename, clobber = False ):
    if os.path.exists( filename ) and not clobber:
        raise RuntimeError( "File exists, will not clobber: %s" % filename )
    
    skimage.io.imsave( filename, skimage.img_as_ubyte( img ) )
    print( "Saved:", filename )

def circle_histogram_for_imagepath( inpath, clobber = False ):
    mode = 'histogram'
    outpath = os.path.splitext( inpath )[0] + f'-{mode}.svg'
    if os.path.exists( outpath ) and not clobber:
        raise RuntimeError( "File exists, will not clobber: %s" % outpath )
    
    # Load the image
    print( "Loading:", inpath )
    img = skimage.io.imread( inpath )
    
    polyline = circle_histogram_for_image( img, clobber = False )
    
    import drawsvg
    
    # Save the polyline as an SVG
    d = drawsvg.Drawing(100, 100)
    d.append(drawsvg.Lines(*np.array(polyline).ravel(), close='true', stroke='black', stroke_width=2, fill='none'))
    
    d.save_svg( outpath )
    print( "Saved:", outpath )

def circle_histogram_for_imagepath_masking_background( inpath, backgroundpath, clobber = False ):
    mode = 'histofancy'
    outpath = os.path.splitext( inpath )[0] + f'-{mode}.svg'
    if os.path.exists( outpath ) and not clobber:
        raise RuntimeError( "File exists, will not clobber: %s" % outpath )
    
    # Load the image
    print( "Loading:", inpath )
    img = skimage.io.imread( inpath )
    
    polyline = circle_histogram_for_image( img, clobber = False )
    
    import drawsvg
    
    # Save the polyline as an SVG
    PADDING = 10
    STROKE_WIDTH = 0.5
    d = drawsvg.Drawing( 100 + 2*PADDING + STROKE_WIDTH, 100 + 2*PADDING + STROKE_WIDTH, origin = ( -PADDING - 0.5*STROKE_WIDTH, -PADDING - 0.5*STROKE_WIDTH ) )
    
    clip_path = drawsvg.Path()
    ## The inner histogram
    polyline.reverse()
    clip_path.M( *polyline[0] )
    for pt in polyline[1:]: clip_path.L( *pt )
    clip_path.Z()
    ## The outer circle
    clip_path.M( 100 + PADDING, 50 )
    clip_path.A( 50 + PADDING/2, 50 + PADDING/2, 0, 1, 1, -PADDING, 50 )
    clip_path.A( 50 + PADDING/2, 50 + PADDING/2, 0, 0, 1, 100 + PADDING, 50 )
    
    clip = drawsvg.ClipPath()
    clip.append( clip_path )
    
    d.append( drawsvg.Image( -PADDING, -PADDING, 100 + 2*PADDING, 100 + 2*PADDING, path = backgroundpath, embed = True, clip_path = clip ) )
    
    ## The same data on top as a line
    d.append( drawsvg.Lines( *np.array(polyline).ravel(), close='true', stroke='black', stroke_width=STROKE_WIDTH, fill='none' ) )
    d.append( drawsvg.Circle( 50, 50, 60, stroke='black', stroke_width=STROKE_WIDTH, fill='none' ) )
    
    d.save_svg( outpath )
    print( "Saved:", outpath )

def circle_histogram_for_image( img, filter = None, num_bins = None, clobber = False ):
    '''
    Given:
        img: An sRGB image as a width-by-height-by-3 numpy.array with floating point values in [0,1].
        filter: An optional filter parameter. Only colors with an ||ab||/128 greater than this value are counted. The default is 0.1.
        num_bins: The number of histogram bins. Too many bins can lead to a noisy appearance. The default is 180.
    '''
    
    if filter is None: filter = 0.1
    if num_bins is None: num_bins = 180
    
    # Convert the image to LAB
    if img.shape[2] == 4:
        assert ( skimage.util.img_as_float( img )[:,:,3] == 1.0 ).all()
        img = img[:,:,:3]
    image_lch = color.lab2lch( color.rgb2lab( img ) )
    
    # Flatten to a 1D sequence of pixels.
    pixels_lch = image_lch.reshape(-1,3)
    
    # Filter out pixels with chroma magnitude < 0.1.
    # Chroma range is [0,100]: https://scikit-image.org/docs/stable/api/skimage.color.html#skimage.color.lch2lab
    if filter is not None: pixels_lch = pixels_lch[ pixels_lch[:,1]/100 > filter ]
    
    # Make a degree histogram
    hist, bins = np.histogram( pixels_lch[:,2], bins = num_bins, range = ( 0, 2*pi ), density = False )
    # Get the maximum bin value for normalization purposes
    max_hist = hist.max()
    # Get the distance from a bin edge to the midpoint
    bin_edge_to_midpoint = .5 * ( bins[0] + bins[1] )
    
    # Create a polyline connecting bin tips
    polyline = []
    for value, bin_edge in zip( hist, bins ):
        radians = bin_edge + bin_edge_to_midpoint
        # We want a radius based on the histogram value.
        # A value of 0 should have r = 1.
        # The maximum value should have r = 0.
        r = 1 - value/max_hist
        x, y = r*np.cos( radians ), r*np.sin( radians )
        # A little transformation to put points in the unit square [0,100]^2
        polyline.append( ( 50*x + 50, 50*y + 50 ) )
    
    return polyline

def main():
    import argparse
    parser = argparse.ArgumentParser( description = "Save an lch color wheel." )
    parser.add_argument( "--clobber", action = 'store_true', help="If set, this will overwrite existing files." )
    parser.add_argument( "-L", type = float, help="Override the default luminance." )
    parser.add_argument( "-c", type = float, help="Override the default chroma." )
    parser.add_argument( "--seed", type = int, help="The random seed. Default 1." )
    parser.add_argument( "--wheel", action = 'store_true', help="Generate a wheel saved to 'wheel.png'." )
    parser.add_argument( "--slice", action = 'store_true', help="Generate a slice saved to 'slice.png'." )
    parser.add_argument( "--strip", action = 'store_true', help="Generate a strip saved to 'strip.png'." )
    parser.add_argument( "--histogram", type = str, help="Path to image to create a histogram for." )
    parser.add_argument( "--histofancy", type = str, nargs = 2, help="Path to image to create a histogram for and the background image. Generates output similar to Color Harmonization [Cohen-Or et al. 2006]." )
    args = parser.parse_args()
    
    if args.seed is not None:
        random.seed( args.seed )
    
    global LUMINANCE, CHROMA
    if args.L is not None:
        LUMINANCE = args.L
    if args.c is not None:
        CHROMA = args.c
    
    print( "Luminance:", LUMINANCE )
    print( "Chroma:", CHROMA )
    
    if args.strip:
        W = generate_wheel_as_strip()
        save_image_to_file( W, 'strip.png', clobber = args.clobber )
    
    if args.wheel:
        W = generate_wheel_as_circle()
        save_image_to_file( W, 'wheel.png', clobber = args.clobber )
    
    if args.slice:
        W = generate_lch_slice()
        save_image_to_file( W, 'slice.png', clobber = args.clobber )
    
    if args.histogram:
        circle_histogram_for_imagepath( args.histogram, clobber = args.clobber )
    
    if args.histofancy:
        circle_histogram_for_imagepath_masking_background( args.histofancy[0], args.histofancy[1], clobber = args.clobber )

if __name__ == '__main__':
    main()