R/colorSpec.metrics.R

In colorSpec: Color Calculations with Emphasis on Spectral Data

    
#   x           a colorSpec object with type "responsivity.material" or "responsivity.light"
#
#   value is a list with items
#       generators  a pair of integers, the 1st is the number of original generators,
#                   and the 2nd is the number of 'condensed' generators
#       zeros       vector of wavelengths at which the responsivity is 0 (in all channels)
#       multiples   a list of groups of wavelengths, 
#                   the responsivities in each group are positive multiples of each other
#       salient     a logical where TRUE means that there is some open linear halfspace that contains
#                   all the non-zero generators
#       normal      the inward pointing unit normal for the above halfspace (if salient==TRUE)
#
#   If salient is TRUE then the list also contains:
#       concavities a data.frame with 2 columns: wavelength and extangle
#                   where extangle is the external angle at the wavelength, and is negative.
#                   A negative angle means the chromaticity polygon is concave at that vertex.
#       coneangle   the solid angle of the cone generated by the generators.
#                   This the area of the chromaticity polygon projected onto the unit sphere,
#                   to form a spherical polygon with concavities preserved.
#       cxconeangle the solid angle of the convex cone generated by the generators.
#                   This the area of the convex hull of the chromaticity polygon projected onto the unit sphere,
#                   to form a spherical polygon with no concavities.
#
#   If the type of x is "responsivity.material" then the list also contains:
#       area        the surface are of the color solid
#       volume      the volume of the color solid
    
responsivityMetrics.colorSpec <- function( x )
    {
    theName = deparse(substitute(x))
    
    ok  = (numSpectra(x) == 3L)
    if( ! ok )
        {
        log.string( ERROR, "numSpectra(%s) = %d  != 3.", theName, numSpectra(x) )     
        return(NULL)
        }        

    ok  = type(x)=="responsivity.light"  ||  type(x)=="responsivity.material" 
    if( ! ok )
        {
        log.string( ERROR, "type(%s) = '%s' is invalid.", theName, type(x) )        
        return(NULL)
        }    
        
    x   = linearize(x)
     
    out = list()
    
    wave    = wavelength(x)    
    step    = breakandstep(wave)$stepvec
    W       = step * as.matrix( x )  #;  print(W)    # step   is replicated to all columns of as.matrix(x)

    #   the tolerance here is for collinearity, which can be larger than the face normal differences
    condlist    = condenseGenerators( W, tol=5.e-7 )
    if( is.null(condlist) )    return(NULL)
    
    ncond   = nrow(condlist$Wcond)    
    if( ncond < 3 )
        {
        log.string( ERROR, "Invalid number of condensed generators = %d < 3.", ncond )
        return(NULL)
        }
        
    out$generators  = c( nrow(W), ncond )
        
    out$zeros   = wave[ is.na(condlist$groupidx) ]
    
    #   record the groups of generators that are multiples of each other
    lenvec          = sapply( condlist$group, length )  #; print(lenvec)
    #   print( str(condlist$group[ 1 < lenvec ]) )
    out$multiples   = lapply( condlist$group[ 1 < lenvec ], function(idx) { wave[idx] }  )
        
        
    #   salient
    #   first try the obvious normal=(1,1,1)     
    normal  = NA_real_
    if( all(0 < colSums(W)) )
        {
        normal  = sqrt( c(1,1,1)/3 )
        }
    else
        {
        if( ! requireNamespace( 'quadprog', quietly=TRUE ) )
            {
            log.string( ERROR, "Required package 'quadprog' could not be imported."  )
            return(NULL)
            }     
        
        res = try( quadprog::solve.QP( diag(3), numeric(3), t(x$Wcond), rep(1,ncond) ), silent=TRUE )
        #   print(res)
        
        if( inherits(res,"try-error") )    
            {
            normal  = res$solution
            normal  = normal / sqrt( sum(normal^2) )
            }        
        }
        
    out$salient = all( is.finite(normal) )
    out$normal  = normal
    
    if( ! out$salient ) return(out)
    

        
    #   compute the exterior angles of the spherical polygon
    knext   = c( 2:ncond, 1 )
    kprev   = c( ncond, 1:(ncond-1) )
    
    #   unitize all the rows, to make a spherical polygon
    W2      = .rowSums( condlist$Wcond*condlist$Wcond, nrow(condlist$Wcond), ncol(condlist$Wcond) )  
    Wunit   = condlist$Wcond / sqrt( W2 )   # the denominator is replicated to all columns of Wcond
    
    anglescone  = externalangles( Wunit )
    
    #   compute solid angle of the cone generated by rows of Wcond, and W
    #   this is the same as the area of the spherical polygon
    coneangle   = areaSphericalPolygon( Wunit, normal )   # max( 2*pi - sum(anglescone), 0 )
    theSign     = sign(coneangle)
    
    if( theSign < 0 )
        anglescone  = -anglescone
    
    #   expand these angles from condensed to original generators
    extangle    = rep( NA_real_, nrow(W) )
    for( k in 1:ncond )
        {
        extangle[ condlist$group[[k]] ] = anglescone[k]
        }

    #   make data.frame
    df  = data.frame( wavelength=wave, extangle=extangle )  #; print(df)
    #   only keep negative angles, but not -pi
    mask    = is.finite(extangle)  &  extangle<0  &  -pi<extangle
    out$concavities = df[ mask, ]
    
    #   compute solid angle of the cone generated by rows of Wcond, and W
    #   this is the same as the area of the spherical polygon
    out$coneangle   = theSign * coneangle   # max( 2*pi - sum(anglescone), 0 )
        
    #   compute the convex hull of this spherical polygon
    #   denom is guaranteed to be all positive, by construction
    denom   = as.numeric( condlist$Wcond %*% normal )
    xyz     = condlist$Wcond / denom        # replicated to all columns
    xy      = xyz %*% frame3x2fun(normal)
    
    convexidx       = grDevices::chull( xy[ ,1], xy[ ,2] ) #; print( str(convexidx) )
    #   anglescxcone    = externalangles( Wunit[convexidx, ] )
    out$cxconeangle = abs( areaSphericalPolygon( Wunit[convexidx, ], normal ) )       # max( 2*pi - sum(anglescxcone), 0 )
    
    #   quick check
    ok  = out$coneangle <= out$cxconeangle
    if( ! ok )
        {
        #print( out$coneangle )
        #print( out$cxconeangle )
        #log.string( WARN, "coneangle=%g  >  %g=cxconeangle.  One or both must be incorrect.",
        #                        out$coneangle, out$cxconeangle )
        }
    
    
    #   compute cyclicorder - in the sense of Paul Centore   
    #   out$cyclicorder = cyclicorderfun( condlist, convexidx )

    if( type(x) == "responsivity.material" )
        {
        #   for area and volume, explicitly make the zonohedron.
        #   clustering is not necessary, so set tol2=NULL to skip that.
        zono    = zonohedron( W, tol2=NULL )
        if( is.null(zono) ) return(NULL)
        
        mets    = metrics( zono )
        out$area    = mets$area
        out$volume  = mets$volume
        }

    return( invisible(out) )
    }
    
    
#   compute cyclicorder - in the sense of Paul Centore    

cyclicorderfun <- function( condlist, convexidx )
    {
    #   first check that the wavelengths in each collinear group are contiguous
    contig  = sapply( condlist$group, contiguous, nrow(condlist$W) )  #;    print( cyclicorder )
    if( ! all(contig) ) return(FALSE)
    
    #   check that every vertex of xy is on the convex hull of xy
    ncond   = nrow(condlist$Wcond)   
    if( length(convexidx) != ncond )    return(FALSE)

    #   check that convexidx can be cyclically rotated to 1:ncond or c(1,ncond:2)
    idx = which( convexidx == 1 )
    if( length(idx) != 1 )  return(FALSE)
    
    if( idx != 1 )
        #   rotate the indexes to put 1 in the first place
        convexidx = convexidx[ c(idx:ncond,1:(idx-1)) ]
        

    ok  = all( convexidx==1:ncond )  ||  all( convexidx==c(1,ncond:2) )
    if( ! ok ) return(FALSE)
    
    return(TRUE)
    }
    
    
    
#   Wunit   Mx3 matrix with each row on the unit sphere
#           no duplicate rows, and no antipodals allowed
#
#   value   M-vector with exterior angle, angle[k] = exterior angle at k
#           the path is from kprev to k to knext    
#           reversing orientation of polygon changes sign of all angles

externalangles <- function( Wunit )
    {
    m       = nrow(Wunit)
    
    knext   = c( 2:m, 1 )
    kprev   = c( m, 1:(m-1) )
    
    #   precompute dot product of each unit vector with the next one
    dotnext = rowSums( Wunit * Wunit[knext, ] )
    
    out     = rep( NA_real_, m )
    for( k in 1:m )
        {
        #   compute 2 tangent vectors at point k
        wk      = Wunit[k, ]
        
        #tannext = -dotnext[k]*wk        + Wunit[knext[k], ]
        #tanprev =  dotnext[kprev[k]]*wk - Wunit[kprev[k], ]

        tanprev = crossproduct( Wunit[kprev[k], ], wk )
        tannext = crossproduct( wk, Wunit[knext[k], ] )
        
        #cat( "---------", k, '\n' )
        #print( tannext )
        #print( tanprev )
        
        out[k]  = angleBetween( tanprev, tannext )
        
        if( is.na(out[k]) ) next
        
        #   now for the sign
        # cross   = crossproduct( tanprev, tannext )
        #theSign = sign( sum( wk * cross ) )         
        #out[k]  = sign( det( rbind( wk, tanprev, tannext ) ) ) * out[k]
        #next
        
        if( det( rbind( wk, tanprev, tannext ) ) < 0 )
            # reverse it
            out[k]  = -out[k]
        }
    
    #   if( sum(out) < 0 )  out = -out      # reverse it
        
    return( out )
    }
    
#   Wunit   3x3 matrix with unit vectors in the rows
#
#   return *signed* area of the triangle    
areaSphericalTriangle  <- function( Wunit )
    {
    anglevec    = externalangles( Wunit )   #; print( anglevec )
    
    if( any( is.na(anglevec) ) )    return(0)   # duplicate vertices
    
    theSum  = sum(anglevec)
    s       = sign(theSum)
    
    return( s * (2*pi - s*theSum) )
    }
    
#   Wunit   3x3 matrix with unit vectors in the rows
#   h       center of hemisphere that contains all points of Wunit
#           h can be any point so that -h is different than any point in Wunit
#
#   return *signed* area of polygon
    
areaSphericalPolygon  <- function( Wunit, h )
    {
    #   accumulate the signed area of all the triangles 
    #   based at h and formed from edges of input polygon
    n       = nrow(Wunit)
    inext   = c(2:n,1)
    
    out = 0
    for( i in 1:n )
        {
        out = out + areaSphericalTriangle( rbind(h,Wunit[ c(i,inext[i]), ]) )
        }

    return( out )
    }
    

#--------       UseMethod() calls           --------------#            
              
responsivityMetrics <- function( x )
    {
    UseMethod("responsivityMetrics")
    }