Solve for SG's using singular value decomposition
| 130 | |
| 131 | // Solve for SG's using singular value decomposition |
| 132 | static void SolveSVD(SGSolveParam& params) |
| 133 | { |
| 134 | // -- Linearly solve for the rgb channels one at a time |
| 135 | Eigen::MatrixXf Ar, Ag, Ab; |
| 136 | |
| 137 | Ar.resize(params.NumSamples, params.NumSGs); |
| 138 | Ag.resize(params.NumSamples, params.NumSGs); |
| 139 | Ab.resize(params.NumSamples, params.NumSGs); |
| 140 | Eigen::VectorXf br(params.NumSamples); |
| 141 | Eigen::VectorXf bg(params.NumSamples); |
| 142 | Eigen::VectorXf bb(params.NumSamples); |
| 143 | for(uint32 i = 0; i < params.NumSamples; ++i) |
| 144 | { |
| 145 | // compute difference squared from actual observed data |
| 146 | for(uint32 j = 0; j < params.NumSGs; ++j) |
| 147 | { |
| 148 | float exponent = std::exp((Float3::Dot(params.XSamples[i], params.OutSGs[j].Axis) - 1.0f) * |
| 149 | params.OutSGs[j].Sharpness); |
| 150 | Ar(i, j) = exponent; |
| 151 | Ag(i, j) = exponent; |
| 152 | Ab(i, j) = exponent; |
| 153 | } |
| 154 | br(i) = params.YSamples[i].x; |
| 155 | bg(i) = params.YSamples[i].y; |
| 156 | bb(i) = params.YSamples[i].z; |
| 157 | } |
| 158 | |
| 159 | Eigen::VectorXf rchan = Ar.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(br); |
| 160 | Eigen::VectorXf gchan = Ag.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(bg); |
| 161 | Eigen::VectorXf bchan = Ab.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(bb); |
| 162 | |
| 163 | for(uint32 j = 0; j < params.NumSGs; ++j) { |
| 164 | params.OutSGs[j].Amplitude.x = rchan[j]; |
| 165 | params.OutSGs[j].Amplitude.y = gchan[j]; |
| 166 | params.OutSGs[j].Amplitude.z = bchan[j]; |
| 167 | } |
| 168 | } |
| 169 | |
| 170 | // Project sample onto SGs |
| 171 | void ProjectOntoSGs(const Float3& dir, const Float3& color, SG* outSGs, uint64 numSGs) |