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, int resolution, SPDResamplingMethod resamplignMethod)
00036 : SPD() {
00037
00038 int lambdaMin = Ceil2Int(wavelengths[0] / resolution) * resolution;
00039 int lambdaMax = Floor2Int(wavelengths[n-1] / resolution) * resolution;
00040
00041 int sn = (lambdaMax - lambdaMin) / resolution + 1;
00042
00043 float *sam = new float[sn];
00044 float *sd = NULL;
00045
00046 if (resamplignMethod == Linear) {
00047 int k = 0;
00048 for (int i = 0; i < sn; i++) {
00049
00050 float lambda = lambdaMin + i * resolution;
00051
00052 if (lambda < wavelengths[0] || lambda > wavelengths[n-1]) {
00053 sam[i] = 0.0;
00054 continue;
00055 }
00056
00057 for (; k < n; ++k)
00058 if (wavelengths[k] >= lambda)
00059 break;
00060
00061 if (wavelengths[k] == lambda)
00062 sam[i] = samples[k];
00063 else {
00064 float intervalWidth = wavelengths[k] - wavelengths[k - 1];
00065 float u = (lambda - wavelengths[k - 1]) / intervalWidth;
00066 sam[i] = ((1. - u) * samples[k - 1]) + (u * samples[k]);
00067 }
00068 }
00069 }
00070 else {
00071 sd = new float[sn];
00072
00073 calc_spline_data(wavelengths, samples, n, sd);
00074
00075 int k = 0;
00076 for (int i = 0; i < sn; i++) {
00077
00078 float lambda = lambdaMin + i * resolution;
00079
00080 if (lambda < wavelengths[0] || lambda > wavelengths[n-1]) {
00081 sam[i] = 0.0;
00082 continue;
00083 }
00084
00085 while (lambda > wavelengths[k+1])
00086 k++;
00087
00088 float h = wavelengths[k+1] - wavelengths[k];
00089 float a = (wavelengths[k+1] - lambda) / h;
00090 float b = (lambda - wavelengths[k]) / h;
00091
00092 sam[i] = max(a*samples[k] + b*samples[k+1]+
00093 ((a*a*a-a)*sd[k] + (b*b*b-b)*sd[k+1])*(h*h)/6.0, 0.);
00094 }
00095 }
00096
00097 init(lambdaMin, lambdaMax, sam, sn);
00098
00099 delete[] sam;
00100 delete[] sd;
00101 }
00102
00103 void IrregularSPD::init(float lMin, float lMax, const float* const s, int n) {
00104 lambdaMin = lMin;
00105 lambdaMax = lMax;
00106 delta = (lambdaMax - lambdaMin) / (n-1);
00107 invDelta = 1.f / delta;
00108 nSamples = n;
00109
00110 AllocateSamples(n);
00111
00112
00113 for (int i = 0; i < n; i++)
00114 samples[i] = s[i];
00115
00116 Clamp();
00117 }
00118
00119
00120
00121 void IrregularSPD::calc_spline_data(const float* const wavelengths,
00122 const float* const amplitudes, int n, float *spline_data) {
00123 float *u = new float[n-1];
00124
00125
00126 spline_data[0] = u[0] = 0.f;
00127
00128 for (int i = 1; i <= n-2; i++) {
00129 float sig = (wavelengths[i] - wavelengths[i-1]) / (wavelengths[i+1] - wavelengths[i-1]);
00130 float p = sig * spline_data[i-1] + 2.f;
00131 spline_data[i] = (sig-1.0)/p;
00132 u[i] = (amplitudes[i+1] - amplitudes[i]) / (wavelengths[i+1] - wavelengths[i]) -
00133 (amplitudes[i] - amplitudes[i-1]) / (wavelengths[i] - wavelengths[i-1]);
00134 u[i] = (6.0*u[i] / (wavelengths[i+1] - wavelengths[i-1]) - sig*u[i-1]) / p;
00135 }
00136
00137 float qn, un;
00138
00139
00140 qn = un = 0.0;
00141 spline_data[n-1] = (un - qn * u[n-2]) / (qn * spline_data[n-2] + 1.0);
00142 for (int k = n-2; k >= 0; k--) {
00143 spline_data[k] = spline_data[k] * spline_data[k+1] + u[k];
00144 }
00145
00146 delete[] u;
00147 }
00148
00149
00150