;+ ; NAME: ; spec2cie ; PURPOSE: (one line only) ; Convert from reflectance spectrum to CIE Chromaticity coordinates ; DESCRIPTION: ; Conversion from reflectance spectrum to CIE Chromaticity (D65). ; The color conversion is based on a Fortran subroutine ; from A. T. Young, ; SDSU ; See: http://mintaka.sdsu.edu/GF/explain/optics/color/color.html ; ; This routine is designed to let you feed it spectral information that ; may or may not be encoded against other dimensions (most commonly ; position) and convert the spectral dimension to its CIE equivalent. ; CATEGORY: ; Image display ; CALLING SEQUENCE: ; spec2cie,wave,spec,x,y,bigy ; INPUTS: ; wave - scalar or vector, wavelengths, in nanometers for spectrum ; The spectral range of interest is 380-770 nm (visual range) ; This does not need to be regularly spaced but should be monotonical ; increasing. ; spec - Spectrum or spectral image data, this information is expected ; to be in the range of 0-1 where 1 is perfectly reflecting and ; 0 is perfectly absorbing. ; rank(spec) cannot be < rank(wave). Last array index ; for spec matches the index for wave. ; OPTIONAL INPUT PARAMETERS: ; KEYWORD INPUT PARAMETERS: ; SUBSAMPLE - Allows for finer resolution input spectra to be accurately ; transformed. The default value is 1 and should work on anything ; with spectral structure coarser than 10 nm. To properly capture ; the nuances of your input spectra, divide 10 by the sampling ; interval of your spectra. Round this number up to the next ; integer and this is the value you'd give as EXPAND. ; SHOWSUMMARY - Flag, if set will print and plot information for the ; scalar case. Ignored for vector or array data. ; OUTPUTS: ; x,y - CIE chromaticity, rank is one less than spec. ; bigy - intensity ; Note: white is x=0.3127, y=0.3290, bigy=Y=1.0 ; (z=0.3583 because x+y+z=1) ; KEYWORD OUTPUT PARAMETERS: ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; Sorry, this is not an elegant or general coding solution. This could ; be more general but I dont' have the time to figure out the IDL syntax ; to maximize vectorability. This won't be too fast if wave is short and ; spec is really big. ; PROCEDURE: ; The tricky part in using this routine effectively is to understand how ; the wavelength information is imbedded in the input variables. ; Consider the simplest possible case, two wavelenths of information ; and just one points to calculate this for: ; wave=[450,600] ; two wavelengths, don't have to match end points ; spec=[0.5,0.6] ; slight red slope ; will return scalar values for x,y,bigy ; Another case, same wavelength vector, but now two "scans" ; spec -> 100x2 array ; second dimension matches wave vector ; result x,y,bigy are each 100 elements long ; Last case, saveme wavelength vector, but now two images ; spec -> 100x100x2 array ; third dimension matches wave vector ; result x,y,bigy are each 100x100 arrays ; ; MODIFICATION HISTORY: ; Written by Marc W. Buie, Southwest Research Institute, 2009/01/05 ;- pro spec2cie,wave,spec,x,y,bigy,SUBSAMPLE=subsample,SHOWSUMMARY=showsummary self='SPEC2CIE' if badpar(wave,[2,3,4,5],1, caller=self+'(wave) ') then return if badpar(spec,[4,5],[1,2,3],caller=self+'(spec) ', $ dimen=dimen,type=type,rank=rank) then return if badpar(subsample,[0,2,3],0, caller=self+'(SUBSAMPLE) ', $ default=1) then return if badpar(showsummary,[0,1,2,3],0, caller=self+'(SHOWSUMMARY) ', $ default=0) then return refw=findgen(40)*10.0 + 380.0 xbar=[0.004,0.019,0.085,0.329,1.238,2.997,3.975,3.915,3.362,2.272, $ 1.112,0.363,0.052,0.089,0.576,1.523,2.785,4.282,5.880,7.322, $ 8.417,8.984,8.949,8.325,7.070,5.309,3.693,2.349,1.361,0.708, $ 0.369,0.171,0.082,0.039,0.019,0.008,0.004,0.002,0.001,0.001]/100.0 ybar=[0.000,0.000,0.002,0.009,0.037,0.122,0.262,0.443,0.694,1.058, $ 1.618,2.358,3.401,4.833,6.462,7.934,9.149,9.832,9.841,9.147, $ 7.992,6.627,5.316,4.176,3.153,2.190,1.443,0.886,0.504,0.259, $ 0.134,0.062,0.029,0.014,0.006,0.003,0.002,0.001,0.001,0.000]/100.0 zbar=[0.020,0.089,0.404,1.570,5.949,14.628,19.938,20.638,19.299, $ 14.972,9.461,5.274,2.864,1.52,0.712,0.388,0.195,0.086,0.039,0.02, $ 0.016,0.01,0.007,0.002,0.002,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., $ 0.,0.,0.]/100.0 if subsample gt 1 then begin idx=findgen(40*subsample+1)/float(subsample) xbarp = sint(idx,xbar)/float(subsample) ybarp = sint(idx,ybar)/float(subsample) zbarp = sint(idx,zbar)/float(subsample) refwp=findgen(40*subsample+1)*10.0/subsample + 380.0 endif else begin xbarp = xbar ybarp = ybar zbarp = zbar refwp = refw endelse ; wave is vector ; spec is vector ; result is a scalar if rank eq 1 then begin interp,wave,spec,refwp,ispec x = total(ispec*xbarp) ; reduces rank here y = total(ispec*ybarp) ; reduces rank here z = total(ispec*zbarp) ; reduces rank here bigy = y sum = x+y+z x = x/sum y = y/sum setwin,0 allmax=max([xbarp,ybarp,zbarp]) cie2rgb,x,y,bigy,r,g,b color = ishft(long(b),16) + ishft(long(g),8) + long(r) if showsummary then begin print,'CIE x,y,Y ',x,y,bigy print,'RGB values ',r,g,b print,'BGR in hex ',b,g,r,format='(a,1x,3Z2.2)' print,'color hex value ',color,format='(a,1x,Z6.6)' plot,refwp,ispec,yr=[0.0,max([ispec,1.0])],/nodata, $ xtitle='Wavelength (nm)',ytitle='Reflectance' oplot,refwp,ispec,color=color,thick=7 oplot,wave,spec,psym=4 oplot,refw,xbar/allmax/subsample,color='0000ff'xl,psym=4 oplot,refwp,xbarp/allmax,color='0000ff'xl oplot,refw,ybar/allmax/subsample,color='00ff00'xl,psym=4 oplot,refwp,ybarp/allmax,color='00ff00'xl oplot,refw,zbar/allmax/subsample,color='ff0000'xl,psym=4 oplot,refwp,zbarp/allmax,color='ff0000'xl endif endif ; wave is vector ; spec is array ; result is a vector if rank eq 2 then begin bigy=make_array(dimen[0],type=type) x=make_array(dimen[0],type=type) y=make_array(dimen[0],type=type) z=make_array(dimen[0],type=type) for i=0L,dimen[0]-1 do begin interp,wave,reform(spec[i,*]),refwp,ispec x[i] = total(ispec*xbarp) y[i] = total(ispec*ybarp) z[i] = total(ispec*zbarp) bigy[i] = y[i] sum = x[i]+y[i]+z[i] if sum le 0. then begin x[i] = 0.3127 ; set color to white, intensity will be 0. y[i] = 0.3290 endif else begin x[i] = x[i]/sum y[i] = y[i]/sum endelse endfor endif ; wave is vector ; spec is cube ; result is an array if rank eq 3 then begin bigy=make_array(dimen[0:1],type=type) x=make_array(dimen[0:1],type=type) y=make_array(dimen[0:1],type=type) z=make_array(dimen[0:1],type=type) for i=0L,dimen[0]-1 do begin for j=0L,dimen[1]-1 do begin interp,wave,reform(spec[i,j,*]),refwp,ispec x[i,j] = total(ispec*xbarp) y[i,j] = total(ispec*ybarp) z[i,j] = total(ispec*zbarp) bigy[i,j] = y[i,j] sum = x[i,j]+y[i,j]+z[i,j] if sum le 0. then begin x[i,j] = 0.3127 ; set color to white, intensity will be 0. y[i,j] = 0.3290 endif else begin x[i,j] = x[i,j]/sum y[i,j] = y[i,j]/sum endelse endfor endfor endif end