001 package util.wavelet;
002
003 import static java.lang.Math.*;
004 import java.io.PrintStream;
005
006 import java.util.List;
007 import util.sdss.Quasar;
008
009 /**
010 * Discrete wavelet transform of a signal using the 1D A-Trous method with a b3-spline smoothing function.
011 * The noise for each position in the signal is also required. The WaveletTransform class computes the
012 * wavelet coefficients, the 1-sigma threshold for wavelet coefficients, the multiresolution support and
013 * it associated fields, structures and object trees for feature extraction and finally a set of signals
014 * corresponding to each individual object detected.
015 * Converted from the Feb 2004 distribution of the
016 * <a href="http://www.eso.org/projects/esomidas/">ESO/MIDAS</a> wavelet code
017 * originally written by Jean-Luc Starck in 25-Feb-2003. In particular version 19_1 of functions
018 * <b>wave_1d_trou()</b> and <b>wave_1d_algo_trou_rec()</b>
019 * located in the file <b>contrib/wavelet/libsrc/wave1d.c</b> available by FTP inside
020 * the gzip file <a href="ftp://ftphost.hq.eso.org/pub/midaspub/04FEB/04FEBpl1.1.tar.gz">
021 * ftp://ftphost.hq.eso.org/pub/midaspub/04FEB/04FEBpl1.1.tar.gz</a>.
022 * For those interested in extending this class to support 2D DWT look into <b>pave_2d_tfo()</b>
023 * and <b>pave_2d_build()</b> located in the file <b>contrib/wavelet/libsrc/pave.c</b>
024 * <b>Reference:</b> Starck,J.-L., Murtagh, F., Bijaoui,A.: 1998, <i>Image Processing and Data Analysis</i>,
025 * ISBN 0521599148 Cambridge Univ. Press.
026 *
027 * @author John Talbot
028 */
029
030 public class WaveletTransform {
031
032 /** Wavelet coefficients */
033 public float[][] w;
034 /** The 1-sigma threshold for wavelet coefficients */
035 public float[][] threshold;
036
037 /** Field labels for contiguous sets of above threshold wavelet coefficients. */
038 public int[][] fieldLabel = null;
039
040 /** A list of Structures in this wavelet transform */
041 public java.util.List<Structure> structures = new java.util.ArrayList<Structure>(300);
042
043 /**
044 * Construct a wavelet transform from a signal with non-uniform noise using
045 * the 1D A-Trous method with a b3-spline smoothing function. <b>Output</b> :
046 * The coefficients of the wavelet transform is in w[0][] to w[scales-2][].
047 * The last row of this array is the smoothed array (w[scales-1][]).
048 * The threshold[][] array is the 1-sigma threshold.
049 * The support[][] array is the multiresolution support.
050 *
051 * @param signal array of float values representing a discrete signal
052 * @param noise the standard deviation of the signal
053 * @param scales number of wavelet scales (if scales==0 then the following equation is solved 2^(scales+1) = 3/4 signal array length)
054 * @param boundaryCondition treatment of boundary conditions
055 */
056 public WaveletTransform(float[] signal, float[] noise, int scales, BoundaryCondition boundaryCondition) {
057 int n = signal.length;
058 if (scales == 0) scales = 1 + (int) (log( 3.0d * ((double)n) / 4.0d) / log(2.0d));
059 if (scales < 1) scales = 1;
060 w = new float[scales][n];
061 threshold = new float[scales][n];
062 fieldLabel= new int[scales][n];
063
064 float[] smooth = new float[n];
065 float[] sigma = getStandardDeviation(n);
066 System.arraycopy(signal, 0, smooth, 0, n);
067 for (int scale = 0; scale < scales - 1; scale++) {
068 System.arraycopy(smooth, 0, w[scale], 0, n);
069 int step = (int)(pow(2.0d, (double) scale) + 0.5d);
070 int step2= step * 2;
071 for (int i = 0; i < n; i++) {
072 int indexMinus2 = boundaryCondition.test(i - step2, n);
073 int indexMinus1 = boundaryCondition.test(i - step, n);
074 int indexPlus1 = boundaryCondition.test(i + step, n);
075 int indexPlus2 = boundaryCondition.test(i + step2, n);
076 // smooth with kernel (1/16, 1/4, 3/8, 1/4, 1/16)
077 smooth[i] = 0.375f * w[scale][i]
078 + 0.25f * ( w[scale][indexMinus1] + w[scale][indexPlus1])
079 + 0.0625f * ( w[scale][indexMinus2] + w[scale][indexPlus2]);
080 // TODO: Create an enum of both methods for computing the noise
081 float maxNoise = 0.0f;
082 for (int j = i - step2; j <= i + step2; j++) // find maximum noise in signal
083 maxNoise = max(maxNoise, noise[boundaryCondition.test(j, n)]);
084
085 threshold[scale][i] = sigma[scale] * maxNoise;
086 }
087 for (int i = 0; i < n; i++) w[scale][i] -= smooth[i];
088 }
089 System.arraycopy(smooth, 0, w[scales - 1], 0, n); // copy smoothed array into last row
090 }
091
092 /**
093 * Reconstruct a signal from a wavelet transform. (see Starck 1998 p.24)
094 * @param nsigma the threshold in units of standard deviation (e.g. 3 is typical, 0 is equivalent to keeping all data)
095 * @param positive if true extract objects that have wavelets coefficients exceeding the threshold, otherwise extract negative features only
096 * @return a reconstructed signal array of float values
097 */
098 public float[] reconstruct(float nsigma, boolean positive) {
099 getFields(nsigma, positive);
100 int n = w[0].length;
101 float[] signal = new float[n]; // automatically filled with 0 in constructor
102 for (int scale = w.length - 1; scale > 0; scale--) // never include the first wavelet scale, its mostly noise
103 for (int i = 0; i < n; i++)
104 if (scale == w.length - 1 || abs(w[scale][i]) > (nsigma * threshold[scale][i]))
105 signal[i] += w[scale][i];
106
107 return signal;
108 }
109
110 // each scale has a List of Structure instances
111
112 /**
113 * Computes the field labels for contiguous sets of above threshold wavelet coefficients.
114 * @param nsigma the significance level in number of sigma units
115 * @param positive restriction to positive significant wavelet coefficients otherwise negative only
116 */
117 protected void getFields(float nsigma, boolean positive) {
118 int scales = w.length;
119 int n = w[0].length;
120 int currentLabel = 1;
121 for (int scale = 0; scale < scales - 1; scale++) {
122 boolean inField = false;
123 boolean fields = false;
124 float extremum = 0.0f;
125 int extremumPosition = 0;
126 for (int i = 0; i < n; i++) {
127 int start = 0;
128 if (w[scale][i] > nsigma * threshold[scale][i]) {
129 if (inField == false) {
130 start = i;
131 inField = true;
132 extremum = w[scale][i];
133 extremumPosition = i;
134 fields = true; // indicates that there are fields at this scale
135 }
136 fieldLabel[scale][i] = currentLabel;
137 float currentValue = w[scale][i];
138 if (currentValue > extremum) {
139 extremum = currentValue;
140 extremumPosition = i;
141 }
142 } else if (!positive && w[scale][i] < -nsigma * threshold[scale][i]) {
143 if (inField == false) {
144 start = i;
145 inField = true;
146 extremum = w[scale][i];
147 extremumPosition = i;
148 fields = true; // indicates that there are fields at this scale
149 }
150 fieldLabel[scale][i] = currentLabel;
151 float currentValue = w[scale][i];
152 if (currentValue < extremum) extremum = currentValue;
153 } else {
154 fieldLabel[scale][i] = 0;
155 if (inField) {
156 int end = i - 1;
157 //structures.add(new Structure(scale, start, end, currentLabel, extremumPosition));
158 currentLabel++;
159 inField = false;
160 }
161 }
162 }
163 // if there were fields at this scale ensure that we start with a fresh label
164 if (fields) currentLabel++;
165 }
166 }
167
168 /**
169 * Get the multiresolution support for this wavelet transform by setting
170 * an array of integers set to 1 for significant wavelet coefficients
171 * (those which exceed the noise) 0 otherwise. (see Starck 1998 p.55).
172 * threshold wavelet coefficients.
173 * @param positive restriction to positive significant wavelet coefficients otherwise negative only
174 */
175 protected float[][] getSupport(float nsigma, boolean positive) {
176 int scales = w.length;
177 int n = w[0].length;
178 float[][] support = new float[scales][n];
179 for (int scale = 0; scale < scales - 1; scale++)
180 for (int i = 0; i < n; i++)
181 if (w[scale][i] > nsigma * threshold[scale][i])
182 support[scale][i] = 1;
183 else if (!positive && w[scale][i] < -nsigma * threshold[scale][i])
184 support[scale][i] = 1;
185 else
186 support[scale][i] = 0;
187 return support;
188 }
189
190
191 /**
192 * Get the multiresolution support visualization for this wavelet transform
193 * as a 1D GIF image. (see Starck 1998 p.56)
194 * Optionally allow negative mr values to produce a reduction in grey intensity.
195 *
196 * @param positive only check for positive significant wavelet coefficients
197 * @return multiresolution support
198 */
199 public static int[] getMultiResolutionSupportVisualization(float[][] support, boolean positive) {
200 int scales = support.length;
201 int n = support[0].length;
202 int[] visual = new int[n]; // automatically set to zero by Java
203 int power;
204 if (positive) power = 1;
205 else power = (int)pow(2.0d, (double)scales);
206
207 for (int scale = 0; scale < scales - 1; scale++) {
208 for (int i = 0; i < n; i++)
209 visual[i] += power * support[scale][i];
210
211 power *= 2;
212 }
213 return visual;
214 }
215
216 /**
217 * Discrete wavelet transform using the 1D A-Trous method with a b3-spline smoothing function.
218 * (see Starck 1998 p.21-26)
219 *
220 * @param signal array of float values representing a discrete signal
221 * @param scales number of wavelet scales (if scales==0 then the following equation is solved 2^(scales+1) = 3/4 signal array length)
222 * @param boundaryCondition treatment of boundary conditions
223 * @return the wavelet transform as a 2-dimensional array of floats with the last row equal to the smoothed array
224 */
225 public static float[][] transform(float[] signal, int scales, BoundaryCondition boundaryCondition) {
226 int n = signal.length;
227 if (scales == 0) scales = 1 + (int) (log( 3.0d * ((double)n) / 4.0d) / log(2.0d));
228 if (scales < 1) scales = 1;
229 float[][] wavelet = new float[scales][n];
230 float[] smooth = new float[n];
231 System.arraycopy(signal, 0, smooth, 0, n);
232 for (int scale = 0; scale < scales - 1; scale++) {
233 System.arraycopy(smooth, 0, wavelet[scale], 0, n);
234 int step = (int)(pow(2.0d, (double) scale) + 0.5d);
235 int step2= step * 2;
236 for (int i = 0; i < n; i++) {
237 int indexMinus2 = boundaryCondition.test(i - step2, n);
238 int indexMinus1 = boundaryCondition.test(i - step, n);
239 int indexPlus1 = boundaryCondition.test(i + step, n);
240 int indexPlus2 = boundaryCondition.test(i + step2, n);
241 // smooth with kernel (1/16, 1/4, 3/8, 1/4, 1/16)
242 smooth[i] = 0.375f * wavelet[scale][i]
243 + 0.25f * ( wavelet[scale][indexMinus1] + wavelet[scale][indexPlus1])
244 + 0.0625f * ( wavelet[scale][indexMinus2] + wavelet[scale][indexPlus2]);
245 }
246 for (int i = 0; i < n; i++) wavelet[scale][i] -= smooth[i];
247 }
248 System.arraycopy(smooth, 0, wavelet[scales - 1], 0, n); // copy smoothed array into last row
249 return wavelet;
250 }
251
252 /**
253 * Get the continuum of a signal. The continuum is computed by iteratively substracting
254 * it from the signal; the iteration is required due to border problems. (see Starck 1998 p.127).
255 *
256 * @param signal array of float values representing a discrete signal
257 * @return the continuum of the signal
258 */
259 public static float[] getContinuum(float[] signal) {
260 int n = signal.length;
261 float[] continuum = new float[n]; // array automatically zeroed by Java VM
262 float[] flattenedSignal = new float[n];
263 for (int iteration = 0; iteration < 5; iteration++) {
264 for (int i = 0; i < n; i++) flattenedSignal[i] = signal[i] - continuum[i];
265 float[][] wavelet = transform(flattenedSignal, 0, BoundaryCondition.Mirror);
266 int scales = wavelet.length;
267 for (int i = 0; i < n; i++) continuum[i] += wavelet[scales - 1][i];
268 }
269 return continuum;
270 }
271
272 public static void main(String[] args) {
273 float[] signal = null;
274 float[] noise = null;
275 Quasar quasar = null;
276 try {
277 List quasars = Quasar.readQuasars();
278 quasar = (Quasar) quasars.get(3);
279 //quasar.showSpectra();
280 signal = quasar.getSpectra();
281 noise = quasar.getNoise();
282 for (int i = 0; i < signal.length; i++) if (signal[i] < 0f) signal[i] = 0f;
283 } catch (java.io.IOException e) {
284 System.out.println(e);
285 System.exit(0);
286 }
287 WaveletTransform wt = new WaveletTransform(signal, noise, 0, BoundaryCondition.Mirror);
288 //outputSignalAndWaveletTransform(System.out, noise, wt.s);
289 float[][] support = wt.getSupport(2f, true);
290 int[] visual = getMultiResolutionSupportVisualization(support, true);
291 float[] reconstruction = wt.reconstruct(2f, true);
292 float[] continuum = getContinuum(signal);
293 int n = visual.length;
294 int signals = wt.w.length;
295 for (int i = 0; i < n; i++) {
296 System.out.print(quasar.getWavelength(i) + " " + visual[i] + " " + signal[i] + " " + reconstruction[i] + " " + continuum[i] + " " + noise[i]);
297 //for (int j = 0; j < signals - 1; j++) System.out.print(mr[j][i] + " ");
298 System.out.println();
299 }
300
301 /*
302 int n = 1200;
303 float[] signal = new float[n];
304 for (int i = 0; i < n; i++) {
305 //Chirp
306 //signal[i] = (float) sin((double)PI * (double) (i*i) / (double) (40*signal.length));
307 //pulse
308 //signal[i] = 0.0f;
309 //signal[n/2] = 1.0f;
310 // emission line
311 int peak = n/2;
312 int x = i - n/2;
313 signal[i] = (float) exp(- (double) (x*x) / (double) (50*50));
314 }
315
316 float[][] wavelet = transform(signal, 0, BoundaryCondition.Mirror);
317 outputSignalAndWaveletTransform(System.out, signal, wavelet);
318 */
319
320
321 }
322
323 //------------------------------------stats-------------------------------------------
324
325 static int iset = 0;
326 static float gset;
327 /**
328 * Get a normally distributed gaussian deviate with zero mean and unit variance.
329 * (see function gasdev(long *idnum) in Press (1992) Numerical Recipes in C, 2nd ed. page 289).
330 *
331 * @return normally distributed floating point number
332 */
333 public static float getGaussianDeviate() {
334 if (iset == 0) {
335 float v1, v2, rsq;
336 do {
337 //TODO: use better random number generator like ran1(idnum) from Num.Rec.C
338 v1 = 2f * (float) random() - 1f;
339 v2 = 2f * (float) random() - 1f;
340 rsq = v1*v1 + v2*v2;
341 } while (rsq >= 1f || rsq == 0f);
342 float fac = (float) sqrt(-2d * log(rsq)/rsq);
343 gset = v1 * fac;
344 iset = 1;
345 return v2 * fac;
346 } else {
347 iset = 0;
348 return gset;
349 }
350 }
351
352 /**
353 * Get the a wavelet transform of a signal comprised of a normally distributed gaussian deviate
354 * @param n the size of the random signal
355 */
356 public static float[][] getGaussianWaveletTransform(int n) {
357 float[] noise = new float[n];
358 for (int i = 0; i < n; i++) {
359 noise[i] = getGaussianDeviate();
360 }
361 float[][] wavelet = transform(noise, 0, BoundaryCondition.Mirror);
362 //outputSignalAndWaveletTransform(System.out, noise, wavelet);
363 return wavelet;
364 }
365
366 /**
367 * Get the standard deviation (sigma) at each scale in the wavelet transform
368 * of a random noise signal comprised of a normally distributed gaussian deviates.
369 * These sigmaNoise values can be used to estimate the sigma of a signal using :
370 * <ul>
371 * <li><b>sigma[scale] = sigma[0] * sigmaNoise[scale] / sigmaNoise[0]</b></li>
372 * </ul>
373 * where sigma[0] is estimated by equating it with the standard deviation of scale 0
374 * of the wavelet transform of a signal (i.e signalWavelet[0][i]).
375 * Then when the absolute value of signalWavelet[scale][i] is greater than
376 * 3*sigma[scale] we have a significant coefficient.
377 * If the signal to noise is non-uniform then another technique is used.
378 *
379 * @param size the size of the random signal
380 * @return a list of standard deviations at each scale of the wavelet transform of the noise signal
381 */
382 public static float[] getStandardDeviation(int size) {
383 float[][] wavelet = getGaussianWaveletTransform(size);
384
385 int scales = wavelet.length;
386 float[] sigma = new float[scales];
387 for (int j = 0; j < scales; j++) {
388 float sum = 0f;
389 float sum2 = 0f;
390 for (int i = 0; i < size; i++) {
391 float f = wavelet[j][i];
392 sum += f;
393 sum2 += f*f;
394 }
395 sigma[j] = (float) sqrt((sum2 - sum * sum / (float)size )/ ((float)size - 1f) );
396 }
397 //System.out.printf("%+10.4f", aFloatNumber); // BUG: this fails in JDK1.5.0 beta2!
398 //for (float s : sigma) System.out.println(" " + s);
399 //when n=1200 sigma[] = {0.750, 0.304, 0.177, 0.112, 0.077, 0.062, 0.045, 0.030, 0.039);
400 return sigma;
401 }
402
403 public static void outputSignalAndWaveletTransform(PrintStream out, float[] signal, float[][]wavelet) {
404 int n = wavelet[0].length;
405 for (int i = 0; i < n; i++) {
406 out.print(signal[i]);
407 for (int j = 0; j < wavelet.length; j++) out.print(" " + wavelet[j][i]);
408 out.println();
409 }
410 }
411
412 public static void outputWaveletTransform(PrintStream out, float[][]wavelet) {
413 int n = wavelet[0].length;
414 for (int i = 0; i < n; i++) {
415 for (int j = 0; j < wavelet.length; j++) out.print(" " + wavelet[j][i]);
416 out.println();
417 }
418 }
419
420
421 }