-
Notifications
You must be signed in to change notification settings - Fork 0
/
perv2.py
91 lines (59 loc) · 3.08 KB
/
perv2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
#Project: Human perception of colors [v2]
import numpy as np
import scipy as sp
from numpy import *
# importing daylight spectra [illumination]
daylight_data = np.genfromtxt('/home/aurimas/Amgen/amgen_2013/data/lut.fi/daylight/baso4.asc', delimiter = '') # daylight illumination from reference white
daylight_data_1 = np.genfromtxt('/home/aurimas/Amgen/amgen_2013/data/lut.fi/daylight/tree.asc', delimiter = '') # daylight illumination from spruce tree
daylight_data_2 = np.genfromtxt('/home/aurimas/Amgen/amgen_2013/data/lut.fi/daylight/sky.asc', delimiter = '') # daylight illumination from the sky reflected by the mirror
# importing Munsell chips (Munsell380 glossy) [reflectance]
reflectance_data = np.genfromtxt('/home/aurimas/Amgen/amgen_2013/data/lut.fi/mglossy_all/munsell380_780_1_glossy.asc', delimiter = '') # reflectance spectra of 1600 glossy Munsell color chips,
# importing cone sensitivities [cone sensitivities]
sensitivity_data = np.genfromtxt('/home/aurimas/Amgen/perception/linss2_10e_1.csv', delimiter = ',') # cones sensitivities
baso4 = daylight_data.T;
tree = daylight_data_1.T;
sky = daylight_data_2.T;
print '-' * 10
A = zeros((3,3,1600,)) # defining A matrix
B1 = zeros(1600) # defining B1 matrix
B2 = zeros(1600) # defining B2 matrix
for x in range(1600):
U0 = zeros((3,15)) # defining U0 matrix
U1 = zeros((3,15)) # defining U1 matrix
V0 = zeros((3,15)) # defining V0 matrix
V1 = zeros((3,15)) # defining V1 matrix
for k in range(15):
for i in range(3):
U0[i,k] = np.sum(np.multiply(sensitivity_data[0:388:4,i],tree[0:97,k]))
U1[i,k] = np.sum(np.multiply(sensitivity_data[0:388:4,i],baso4[0:97,k]))
V0[i,k] = np.sum(np.multiply(sensitivity_data[0:388:4,i],reflectance_data[10:398:4,x],tree[0:97,k]))
V1[i,k] = np.sum(np.multiply(sensitivity_data[0:388:4,i],reflectance_data[10:398:4,x],baso4[0:97,k]))
U2 = zeros((3,22)) # defining U2 matrix
V2 = zeros((3,22)) # defining V2 matrix
for l in range(22):
for j in range(3):
U2[j, l] = np.sum(np.multiply(sensitivity_data[0:388:4,j],sky[0:97,l]))
V2[j, l] = np.sum(np.multiply(sensitivity_data[0:388:4,j],reflectance_data[10:398:4,x],sky[0:97,l]))
U = np.concatenate((U0,U1,U2),axis=1) # mixing three U matrixes wrt columns: 3x52
V = np.concatenate((V0,V1,V2),axis=1) # mixing three V matrixes wrt columns: 3x52
U_pinv = np.linalg.pinv(U) # finding pseudoinverse matrix of U
A[:,:,x] = np.dot(V,U_pinv) # finding linear operator [A] that satisfies equation V = AU
D,E = linalg.eig(A[:,:,x]); # eigendecomposition (of linear operator, D - eigenvalues, E - eigenvectors (represent diffent hue compoisitions))
#print D
for x in range (1600):
B1[x] = D[0] / D[1]
#print B1[x]
B2[x] = D[1] / D[2]
b1 = B1[x].max()
# print b1
b2 = B2[x].max()
# print b2
B = zeros(1600) # defining B matrix
B[x] = max(B1[x]/b1,B2[x]/b2) # finding each singularity index
print V.shape
print U_pinv.shape
print A.shape
print D.shape
print E.shape
print sp.diag((D))
print B[x]