Solve for SG's using non-negative least squares
| 82 | |
| 83 | // Solve for SG's using non-negative least squares |
| 84 | static void SolveNNLS(SGSolveParam& params) |
| 85 | { |
| 86 | Assert_(params.XSamples != nullptr); |
| 87 | Assert_(params.YSamples != nullptr); |
| 88 | |
| 89 | // -- Linearly solve for the rgb channels one at a time |
| 90 | Eigen::MatrixXf Ar, Ag, Ab; |
| 91 | Ar.resize(params.NumSamples, int64(params.NumSGs)); |
| 92 | Ag.resize(params.NumSamples, int64(params.NumSGs)); |
| 93 | Ab.resize(params.NumSamples, int64(params.NumSGs)); |
| 94 | Eigen::VectorXf br(params.NumSamples); |
| 95 | Eigen::VectorXf bg(params.NumSamples); |
| 96 | Eigen::VectorXf bb(params.NumSamples); |
| 97 | for(uint32 i = 0; i < params.NumSamples; ++i) |
| 98 | { |
| 99 | // compute difference squared from actual observed data |
| 100 | for(uint32 j = 0; j < params.NumSGs; ++j) |
| 101 | { |
| 102 | float exponent = exp((Float3::Dot(params.XSamples[i], params.OutSGs[j].Axis) - 1.0f) * |
| 103 | params.OutSGs[j].Sharpness); |
| 104 | Ar(i,j) = exponent; |
| 105 | Ag(i,j) = exponent; |
| 106 | Ab(i,j) = exponent; |
| 107 | } |
| 108 | br(i) = params.YSamples[i].x; |
| 109 | bg(i) = params.YSamples[i].y; |
| 110 | bb(i) = params.YSamples[i].z; |
| 111 | } |
| 112 | |
| 113 | Eigen::NNLS<Eigen::MatrixXf> nnlsr(Ar); |
| 114 | Eigen::NNLS<Eigen::MatrixXf> nnlsg(Ag); |
| 115 | Eigen::NNLS<Eigen::MatrixXf> nnlsb(Ab); |
| 116 | nnlsr.solve(br); |
| 117 | nnlsg.solve(bg); |
| 118 | nnlsb.solve(bb); |
| 119 | Eigen::VectorXf rchan = nnlsr.x(); |
| 120 | Eigen::VectorXf gchan = nnlsg.x(); |
| 121 | Eigen::VectorXf bchan = nnlsb.x(); |
| 122 | |
| 123 | for(uint32 j = 0; j < params.NumSGs; ++j) |
| 124 | { |
| 125 | params.OutSGs[j].Amplitude.x = rchan[j]; |
| 126 | params.OutSGs[j].Amplitude.y = gchan[j]; |
| 127 | params.OutSGs[j].Amplitude.z = bchan[j]; |
| 128 | } |
| 129 | } |
| 130 | |
| 131 | // Solve for SG's using singular value decomposition |
| 132 | static void SolveSVD(SGSolveParam& params) |