00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "irregular.h"
00025
00026 using namespace lux;
00027
00028
00029
00030
00031
00032
00033
00034 IrregularSPD::IrregularSPD(const float* const wavelengths, const float* const samples,
00035 int n, float resolution, SPDResamplingMethod resamplignMethod)
00036 : SPD() {
00037
00038 float lambdaMin = Ceil2Int(wavelengths[0] / resolution) * resolution;
00039 float lambdaMax = Floor2Int(wavelengths[n-1] / resolution) * resolution;
00040
00041 int sn = (lambdaMax - lambdaMin) / resolution + 1;
00042
00043 vector<float> sam(sn);
00044
00045 if (resamplignMethod == Linear) {
00046 int k = 0;
00047 for (int i = 0; i < sn; i++) {
00048
00049 float lambda = lambdaMin + i * resolution;
00050
00051 if (lambda < wavelengths[0] || lambda > wavelengths[n-1]) {
00052 sam[i] = 0.0;
00053 continue;
00054 }
00055
00056 for (; k < n; ++k)
00057 if (wavelengths[k] >= lambda)
00058 break;
00059
00060 if (wavelengths[k] == lambda)
00061 sam[i] = samples[k];
00062 else {
00063 float intervalWidth = wavelengths[k] - wavelengths[k - 1];
00064 float u = (lambda - wavelengths[k - 1]) / intervalWidth;
00065 sam[i] = ((1. - u) * samples[k - 1]) + (u * samples[k]);
00066 }
00067 }
00068 }
00069 else {
00070 vector<float> sd(n);
00071
00072 calc_spline_data(wavelengths, samples, n, &sd[0]);
00073
00074 int k = 0;
00075 for (int i = 0; i < sn; i++) {
00076
00077 float lambda = lambdaMin + i * resolution;
00078
00079 if (lambda < wavelengths[0] || lambda > wavelengths[n-1]) {
00080 sam[i] = 0.0;
00081 continue;
00082 }
00083
00084 while (lambda > wavelengths[k+1])
00085 k++;
00086
00087 float h = wavelengths[k+1] - wavelengths[k];
00088 float a = (wavelengths[k+1] - lambda) / h;
00089 float b = (lambda - wavelengths[k]) / h;
00090
00091 sam[i] = max(a*samples[k] + b*samples[k+1]+
00092 ((a*a*a-a)*sd[k] + (b*b*b-b)*sd[k+1])*(h*h)/6.0, 0.);
00093 }
00094 }
00095
00096 init(lambdaMin, lambdaMax, &sam[0], sn);
00097 }
00098
00099 void IrregularSPD::init(float lMin, float lMax, const float* const s, int n) {
00100 lambdaMin = lMin;
00101 lambdaMax = lMax;
00102 delta = (lambdaMax - lambdaMin) / (n-1);
00103 invDelta = 1.f / delta;
00104 nSamples = n;
00105
00106 AllocateSamples(n);
00107
00108
00109 for (int i = 0; i < n; i++)
00110 samples[i] = s[i];
00111
00112 }
00113
00114
00115
00116 void IrregularSPD::calc_spline_data(const float* const wavelengths,
00117 const float* const amplitudes, int n, float *spline_data) {
00118 vector<float> u(n-1);
00119
00120
00121 spline_data[0] = u[0] = 0.f;
00122
00123 for (int i = 1; i <= n-2; i++) {
00124 float sig = (wavelengths[i] - wavelengths[i-1]) / (wavelengths[i+1] - wavelengths[i-1]);
00125 float p = sig * spline_data[i-1] + 2.f;
00126 spline_data[i] = (sig-1.0)/p;
00127 u[i] = (amplitudes[i+1] - amplitudes[i]) / (wavelengths[i+1] - wavelengths[i]) -
00128 (amplitudes[i] - amplitudes[i-1]) / (wavelengths[i] - wavelengths[i-1]);
00129 u[i] = (6.0*u[i] / (wavelengths[i+1] - wavelengths[i-1]) - sig*u[i-1]) / p;
00130 }
00131
00132 float qn, un;
00133
00134
00135 qn = un = 0.0;
00136 spline_data[n-1] = (un - qn * u[n-2]) / (qn * spline_data[n-2] + 1.0);
00137 for (int k = n-2; k >= 0; k--) {
00138 spline_data[k] = spline_data[k] * spline_data[k+1] + u[k];
00139 }
00140 }
00141
00142
00143