eic-opticks/ana/cie.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 cie.py: converts wavelength spectra into XYZ and RGB colorspaces
0023 ===================================================================
0024
0025 Conversion of the binned wavelength spectra into XYZ (using
0026 CIE weighting functions) and then RGB produces a spectrum
0027
0028 [FIXED] Unphysical color repetition
0029 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0030
0031 Uniform scaling by maximal single X,Y,Z or R,G,B
0032 prior to clipping gets rid of the unphysical color repetition
0033 but theres kinda a between the green and the blue, where cyan
0034 should be
0035
0036     #hRGB_raw /= hRGB_raw[0,:,0].max()  # scaling by maximal red, results in muted spectrum
0037     #hRGB_raw /= hRGB_raw[0,:,1].max()  # scaling by maximal green,  OK
0038     #hRGB_raw /= hRGB_raw[0,:,2].max()  # scaling by maximal blue, similar to green by pumps the blues and nice yellow
0039
0040 The entire spectral locus is outside sRGB gamut (the triangle),
0041 so all bins are being clipped.
0042
0043 Not clipping produces a psychedelic mess.
0044
0045
0046 [ISSUE] Blue/Green transition looks unphysical
0047 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0048
0049 Need better way to handle out of gamut ?
0050
0051 Raw numbers show that green ramps up thru 430..480 nm but
0052 its all negative, so that info is clipped.
0053
0054 ::
0055
0056     In [68]: np.set_printoptions(linewidth=150)
0057
0058     In [75]: np.hstack([wd[:-1,None],c.raw[0],c.xyz[0],c.rgb[0]])
0059     Out[75]:
0060     array([[  350.   ,     0.   ,     0.016,     0.102,     0.   ,     0.   ,     0.   ,    -0.   ,     0.   ,     0.   ],
0061            [  370.   ,     0.015,     0.105,     1.922,     0.   ,     0.   ,     0.001,    -0.001,     0.   ,     0.001],
0062            [  390.   ,     1.873,     0.582,    20.444,     0.001,     0.   ,     0.011,    -0.003,     0.   ,     0.012],
0063            [  410.   ,    49.306,     2.691,   205.061,     0.028,     0.002,     0.115,     0.03 ,    -0.019,     0.123],
0064            [  430.   ,   273.393,    10.384,  1386.823,     0.153,     0.006,     0.779,     0.1  ,    -0.105,     0.83 ],
0065            [  450.   ,   343.75 ,    33.415,  1781.385,     0.193,     0.019,     1.   ,     0.098,    -0.11 ,     1.064],
0066            [  470.   ,   191.832,    89.944,  1294.473,     0.108,     0.05 ,     0.727,    -0.091,     0.021,     0.764],
0067            [  490.   ,    32.012,   213.069,   465.525,     0.018,     0.12 ,     0.261,    -0.256,     0.218,     0.253],
0068            [  510.   ,    16.48 ,   500.611,   155.962,     0.009,     0.281,     0.088,    -0.446,     0.522,     0.036],
0069            [  530.   ,   159.607,   869.052,    43.036,     0.09 ,     0.488,     0.024,    -0.472,     0.829,    -0.069],
0070            [  550.   ,   433.715,   994.463,     8.758,     0.243,     0.558,     0.005,    -0.072,     0.812,    -0.095],
0071            [  570.   ,   772.904,   950.107,     1.308,     0.434,     0.533,     0.001,     0.586,     0.58 ,    -0.084],
0072            [  590.   ,  1021.039,   762.587,     0.143,     0.573,     0.428,     0.   ,     1.199,     0.248,    -0.055],
0073            [  610.   ,  1000.205,   500.338,     0.012,     0.561,     0.281,     0.   ,     1.388,    -0.017,    -0.026],
0074            [  630.   ,   656.21 ,   263.667,     0.001,     0.368,     0.148,     0.   ,     0.966,    -0.079,    -0.01 ],
0075            [  650.   ,   283.632,   110.045,     0.   ,     0.159,     0.062,     0.   ,     0.421,    -0.038,    -0.004],
0076            [  670.   ,    80.766,    36.117,     0.   ,     0.045,     0.02 ,     0.   ,     0.116,    -0.006,    -0.002],
0077            [  690.   ,    17.024,    11.172,     0.   ,     0.01 ,     0.006,     0.   ,     0.021,     0.003,    -0.001]])
0078
0079
0080
0081 Chromatic Adaption
0082 ~~~~~~~~~~~~~~~~~~~~
0083
0084 * http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
0085
0086
0087
0088 Refs
0089 ~~~~~
0090
0091 https://github.com/colour-science/colour/issues/191
0092 http://www.scipy-lectures.org/advanced/image_processing/index.html
0093 http://www.scipy-lectures.org/packages/scikit-image/index.html#scikit-image
0094
0095 http://www.scipy.org/scikits.html
0096     separate from scipy, but under the "brand"
0097
0098 """
0099
0100 import os, logging, numpy as np
0101 log = logging.getLogger(__name__)
0102 np.set_printoptions(linewidth=150)
0103
0104 import matplotlib.pyplot as plt
0105 import ciexyz.ciexyz as _cie
0106 from env.graphics.ciexyz.XYZ import Spectrum
0107 from env.graphics.ciexyz.RGB import RGB
0108
0109
0110 class CIE(object):
0111     def __init__(self, colorspace="sRGB/D65", whitepoint=None):
0112         cs = RGB(colorspace)
0113         self.x2r = cs.x2r
0114         self.whitepoint = whitepoint
0115
0116     def hist0d_XYZ(self,w, nb=100):
0117
0118         X = np.sum(_cie.X(w))
0119         Y = np.sum(_cie.Y(w))
0120         Z = np.sum(_cie.Z(w))
0121
0122         hX = np.repeat(X, nb)
0123         hY = np.repeat(Y, nb)
0124         hZ = np.repeat(Z, nb)
0125
0126         raw = np.dstack([hX,hY,hZ])
0127         self.raw = np.copy(raw)
0128         return raw
0129
0130     def hist1d_XYZ(self,w,x,xb):
0131         hX, hXx = np.histogram(x,bins=xb, weights=_cie.X(w))
0132         hY, hYx = np.histogram(x,bins=xb, weights=_cie.Y(w))
0133         hZ, hZx = np.histogram(x,bins=xb, weights=_cie.Z(w))
0134         assert np.all(hXx == xb) & np.all(hYx == xb ) & np.all(hZx == xb)
0135         raw = np.dstack([hX,hY,hZ])
0136         self.raw = np.copy(raw)
0137         return raw
0138
0139     def hist2d_XYZ(self,w,x,y,xb,yb):
0140         bins = [xb,yb]
0141         hX, hXx, hXy = np.histogram2d(x,y,bins=bins, weights=_cie.X(w))
0142         hY, hYx, hYy = np.histogram2d(x,y,bins=bins, weights=_cie.Y(w))
0143         hZ, hZx, hZy = np.histogram2d(x,y,bins=bins, weights=_cie.Z(w))
0144         assert np.all(hXx == xb) & np.all(hYx == xb ) & np.all(hZx == xb)
0145         assert np.all(hXy == yb) & np.all(hYy == yb ) & np.all(hZy == yb)
0146         return np.dstack([hX,hY,hZ])
0147
0148     def norm_XYZ(self, hXYZ, norm=2, scale=1):
0149         """
0150         Trying to find an appropriate way to normalize XYZ values
0151
0152         0,1,2
0153               scale by maximal of X,Y,Z
0154         3
0155               scale by maximal X+Y+Z
0156         4
0157               scale by Yint of an externally determined whitepoint
0158               (problem is that is liable to be with very much more light
0159               than are looking at...)
0160         5
0161               scale by Yint of the spectrum provided, this
0162               is also yielding very small X,Y,Z
0163
0164         >50
0165               scale is used, for normalization with Y value
0166               obtained from the histogram norm identified bin
0167
0168
0169         Hmm, some adhoc exposure factor seems unavoidable given the
0170         range of intensities so perhaps the adhoc techniques are appropriate after all.
0171
0172         Initial thinking was that the out-of-gamut problem was tied up with the
0173         exposure problem, but they are kinda orthogonal: think vectors in XYZ space,
0174         the length of the vector doesnt change the hue.
0175         """
0176         if norm in [0,1,2]:
0177             nscale = hXYZ[:,:,norm].max()
0178         elif norm == 3:
0179             nscale = np.sum(hXYZ, axis=2).max()
0180         elif norm == 4:
0181             assert not self.whitepoint is None
0182             nscale = self.whitepoint[4]
0183         elif norm == 5 or norm > 50:
0184             nscale = scale
0185         else:
0186             nscale = 1
0187         pass
0188
0189         hXYZ /= nscale
0190         self.scale = nscale
0191         self.xyz = np.copy(hXYZ)
0192         return hXYZ
0193
0194     def XYZ_to_RGB(self, hXYZ):
0195         return np.dot( hXYZ, self.x2r.T )
0196
0197     def hist0d(self, w, norm=2, nb=100):
0198         hXYZ_raw = self.hist0d_XYZ(w, nb=nb)
0199         hXYZ = self.norm_XYZ(hXYZ_raw, norm=norm)
0200         hRGB =  self.XYZ_to_RGB(hXYZ)
0201         self.rgb = np.copy(hRGB)
0202         return hRGB,hXYZ,None
0203
0204     def hist1d(self, w, x, xb, norm=2):
0205         hXYZ_raw = self.hist1d_XYZ(w,x,xb)
0206
0207         if norm == 5:
0208             scale = np.sum(_cie.Y(w))
0209         elif norm > 50:
0210             # assume norm is pointing to a bin, the Y value of which is used for scaling
0211             scale = hXYZ_raw[0,norm,1]
0212         else:
0213             scale = 1
0214         pass
0215
0216         hXYZ = self.norm_XYZ(hXYZ_raw, norm=norm, scale=scale)
0217         hRGB =  self.XYZ_to_RGB(hXYZ)
0218         self.rgb = np.copy(hRGB)
0219         return hRGB,hXYZ,xb
0220
0221     def hist2d(self, w, x, y, xb, yb, norm=2):
0222         hXYZ_raw = self.hist2d_XYZ(w,x,y,xb,yb)
0223         self.raw = hXYZ_raw
0224         hXYZ = self.norm_XYZ(hXYZ_raw, norm=norm)
0225         hRGB =  self.XYZ_to_RGB(hXYZ)
0226         extent = [xb[0], xb[-1], yb[0], yb[-1]]
0227         return hRGB,hXYZ,extent
0228
0229     def spectral_plot(self, ax, wd, norm=2):
0230
0231         ndupe = 1000
0232         w = np.tile(wd, ndupe)
0233         x = np.tile(wd, ndupe)
0234         xb = wd
0235
0236         hRGB_raw, hXYZ_raw, bx = self.hist1d(w, x, xb, norm=norm)
0237
0238         hRGB_1d = np.clip(hRGB_raw, 0, 1)
0239         ntile = 100
0240         hRGB = np.tile(hRGB_1d, ntile ).reshape(-1,ntile,3)
0241         extent = [0,2,bx[0],bx[-1]]
0242
0243         #interpolation = 'none'
0244         #interpolation = 'mitchell'
0245         #interpolation = 'hanning'
0246         interpolation = 'gaussian'
0247
0248         ax.imshow(hRGB, origin="lower", extent=extent, alpha=1, vmin=0, vmax=1, aspect='auto', interpolation=interpolation)
0249         ax.yaxis.set_visible(True)
0250         ax.xaxis.set_visible(False)
0251
0252     def swatch_plot(self, wd, norm=2):
0253         ndupe = 1000
0254         w = np.tile(wd, ndupe)
0255         hRGB_raw, hXYZ_raw, bx = self.hist0d(w, norm=norm)
0256
0257         hRGB_1d = np.clip(hRGB_raw, 0, 1)
0258         ntile = 100
0259         hRGB = np.tile(hRGB_1d, ntile ).reshape(-1,ntile,3)
0260         extent = [0,2,0,1]
0261         ax.imshow(hRGB, origin="lower", extent=extent, alpha=1, vmin=0, vmax=1, aspect='auto')
0262         ax.yaxis.set_visible(True)
0263         ax.xaxis.set_visible(False)
0264
0265
0266
0267 def cie_hist1d(w, x, xb, norm=1, colorspace="sRGB/D65", whitepoint=None):
0268     c = CIE(colorspace, whitepoint=whitepoint)
0269     return c.hist1d(w,x,xb,norm)
0270
0271 def cie_hist2d(w, x, y, xb, yb, norm=1, colorspace="sRGB/D65", whitepoint=None):
0272     c = CIE(colorspace, whitepoint=whitepoint)
0273     return c.hist2d(w,x,y,xb,yb,norm)
0274
0275
0276
0277 class Whitepoint(object):
0278     def __init__(self, w):
0279         """
0280         For spectra close to original (think perfect diffuse reflector)
0281         this is expected to yield the characteristic of the illuminant.
0282
0283         XYZ values must be normalized as clearly simulating more photons
0284         will give larger values...
0285
0286         The Yint is hoped to provide a less adhoc way of doing the
0287         normalization.
0288         """
0289         assert w is not None
0290
0291         X = np.sum(_cie.X(w))
0292         Y = np.sum(_cie.Y(w))
0293         Z = np.sum(_cie.Z(w))
0294
0295         Yint = Y
0296
0297         X /= Yint      # normalize such that Y=1
0298         Y /= Yint
0299         Z /= Yint
0300
0301         x = X/(X+Y+Z)  # Chromaticity coordinates
0302         y = Y/(X+Y+Z)
0303
0304         self.wp = np.array([X,Y,Z,Yint,x,y])
0305
0306     def __repr__(self):
0307         return str(self.wp)
0308
0309
0310
0311
0312 def whitepoint(wd):
0313     bb = _cie.BB6K(wd)
0314     bb /= bb.max()
0315
0316     X = np.sum( _cie.X(wd)*bb )
0317     Y = np.sum( _cie.Y(wd)*bb )
0318     Z = np.sum( _cie.Z(wd)*bb )
0319
0320     norm = Y
0321
0322     # Normalize Y to 1
0323     X /= norm
0324     Y /= norm
0325     Z /= norm
0326
0327     return [X,Y,Z], norm
0328
0329
0330
0331 def compare_norms(wdom):
0332     """
0333     norm 0:X,1:Y
0334              look almost identical
0335
0336     norm 2:Z, 3:X+Y+Z
0337               also look the same
0338
0339     0,1 have better yellow and less of a murky gap between green and blue
0340     """
0341     c = CIE()
0342
0343     nplot = 4
0344
0345     for i, norm in enumerate([0,1,2,3]):
0346         ax = fig.add_subplot(1,nplot,i+1)
0347         c.spectral_plot(ax, wdom, norm)
0348
0349         #ax = fig.add_subplot(1,nplot,2*i+2)
0350         #c.swatch_plot(wd, norm)
0351
0352
0353 if __name__ == '__main__':
0354     pass
0355     logging.basicConfig(level=logging.INFO)
0356
0357     plt.close()
0358     plt.ion()
0359
0360     wdom = np.arange(350,720,20)
0361
0362
0363     fig = plt.figure()
0364
0365     compare_norms(wdom)
0366
0367     plt.show()
0368
0369
0370     wp = whitepoint(wdom)
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381