001 package util.wavelet;
002 import java.util.List;
003 import static java.lang.Math.*;
004
005 /** An enumeration of noise computation methods for thresholding in wavelet space.
006 * Each entry in the enumeration implements a variation of the compute() method
007 * which takes as input the signal and its non-uniform noise map and returns
008 * its output in the noise component of the scales list.
009 *
010 * @author John Talbot
011 */
012 public enum NoiseComputation {
013 UpperLimit() {
014 /**
015 * Compute the upper limit of the noise for each wavelet coefficient using
016 * a box car type range over the noise values affecting a particular wavelet
017 * coefficient and finding the maximim noise value in each range.
018 * The computed noise values are returned in the scales list.
019 * <b>Reference</b> :
020 * <a href="http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1997ApJ...482.1011S">Starck <i>et al.</i> : 1997</a>,
021 * <i>Astrophysics Journal</i> <b>482</b>, 1011.
022 *
023 * @param signal <b>input:</b> the signal with non-uniform noise
024 * @param scales <b>output:</b> noise values returned in noise component of each scale in list
025 */
026 public void compute(Signal signal, java.util.List<Scale> scales) {
027 float[] sigma = getStandardDeviation(signal.size());
028 for (int scaleNumber = 0; scaleNumber < scales.size(); scaleNumber++) {
029 Scale scale = scales.get(scaleNumber);
030 int step2 = (int)(pow(2.0d, (double) scaleNumber) + 0.5d) * 2;
031 for (int i = 0; i < signal.size(); i++) {
032 float maxNoise = 0.0f;
033 for (int j = i - step2; j <= i + step2; j++) maxNoise = max(maxNoise, signal.getNoise(j));
034 scale.setNoise(i, sigma[scaleNumber] * maxNoise);
035 }
036 }
037 }
038
039 },
040 Dirac() {
041 /**
042 * Compute the exact noise threshold values by convolving the square of the wavelet transform of the
043 * dirac delta function with the square of the non-uniform noise of a signal.
044 * (see Stark et al. 2002, Astronomical Image an data Analysis p.43 eq. 2.22)
045 *
046 * @param signal <b>input:</b> the signal with noise
047 * @param scales <b>output:</b> the pre-computed scales of the wavelet transform of the signal
048 */
049 public void compute(Signal signal, java.util.List<Scale> scales) {
050 SimpleWaveletSpace diracWaveletSpace = new SimpleWaveletSpace(SimpleSignal.getDirac(signal.size()));
051 SimpleSignal noise2 = signal.getNoise().square();
052 noise2.setBoundaryCondition(BoundaryCondition.Periodic);
053 for (Scale scale : scales) {
054 SimpleSignal result = new SimpleSignal(noise2);
055 SimpleSignal dirac2 = diracWaveletSpace.getScale(scale.getScaleNumber()).square();
056 scale.setNoise(result.convolve(dirac2));
057 }
058 /* TODO: check to see if the dirac function we are using works by using it to
059 * compute a simple wavelet transform
060 */
061 }
062 };
063
064 public static void main(String[] args) {
065 System.out.println("Testing the dirac method to compute noise thresholds ...");
066 // TODO:
067 }
068
069 /**
070 * Compute the noise assuming unit variance for each point in the signal.
071 *
072 * @param signal <b>input:</b> the signal with noise
073 * @param scales <b>output:</b> the pre-computed scales of the wavelet transform of the signal
074 */
075 public void compute(Signal signal, java.util.List<Scale> scales) {
076 float[] sigma = getStandardDeviation(signal.size());
077 for (int scaleNumber = 0; scaleNumber < scales.size(); scaleNumber++) {
078 Scale scale = scales.get(scaleNumber);
079 int step2 = (int)(pow(2.0d, (double) scaleNumber) + 0.5d) * 2;
080 for (int i = 0; i < signal.size(); i++) {
081 float maxNoise = 0.0f;
082 for (int j = i - step2; j <= i + step2; j++) maxNoise = max(maxNoise, 1f);
083 scale.setNoise(i, sigma[scaleNumber] * maxNoise);
084 }
085 }
086 }
087
088
089
090 /**
091 * Get the standard deviation (sigma) at each scale in the wavelet transform
092 * of a random noise signal comprised of a normally distributed gaussian deviates.
093 * These sigmaNoise values can be used to estimate the sigma of a signal using :
094 * <ul>
095 * <li><b>sigma[scale] = sigma[0] * sigmaNoise[scale] / sigmaNoise[0]</b></li>
096 * </ul>
097 * where sigma[0] is estimated by equating it with the standard deviation of scale 0
098 * of the wavelet transform of a signal (i.e signalWavelet[0][i]).
099 * Then when the absolute value of signalWavelet[scale][i] is greater than
100 * 3*sigma[scale] we have a significant coefficient.
101 * If the signal to noise is non-uniform then another technique is used.
102 * (test data: when n=1200 sigma[] = {0.750, 0.304, 0.177, 0.112, 0.077, 0.062, 0.045, 0.030, 0.039))
103
104 * @param size the size of the random signal
105 * @return a list of standard deviations at each scale of the wavelet transform of the noise signal
106 */
107 public static float[] getStandardDeviation(int size) {
108 SimpleSignal noise = SimpleSignal.getGaussianNoise(size);
109 SimpleWaveletSpace simpleWaveletSpace = new SimpleWaveletSpace(noise);
110 int maximumScale = simpleWaveletSpace.getNumberOfScales();
111 float[] sigma = new float[maximumScale];
112 for (int scaleNumber = 0; scaleNumber < maximumScale; scaleNumber++) {
113 float sum = 0f;
114 float sum2 = 0f;
115 SimpleSignal scale = simpleWaveletSpace.getScale(scaleNumber);
116 for (int i = 0; i < size; i++) {
117 float coefficient = scale.getValue(i);
118 sum += coefficient;
119 sum2 += coefficient*coefficient;
120 }
121 sigma[scaleNumber] = (float) sqrt((sum2 - sum * sum / (float)size )/ ((float)size - 1f) );
122 }
123 return sigma;
124 }
125 }
126
127 /* Wavelet transform of the dirac delta function :
128
129 16*h0 a-trous B3-spline smoothing kernel (h0) 1 4 6 4 1
130 c0 dirac delta funtion 0 0 1 0 0
131 16*c1 smoothed array 1 1 4 6 4 1
132 16*h1 smoothing kernel h1 1 0 4 0 6 0 4 0 1
133 256*c2 smoothed array 2 1 4 10 20 31 40 44 40 31 20 10 4 1
134 16*h2 1 0 0 0 4 0 0 0 6 0 0 0 4 0 0 0 1
135 4096*c3 1 4 10 20 35 56 84 120 161 204 246 284 315 336 344 336 315 284 246 204 161 120 84 56 35 20 10 4 1
136 --------------------------------------------------------------------------------------------------------------------
137 16*w0 (c0-c1) wavelet scale 0 -1 -4 10 -4 -1
138 256*w1 (c1-c2) wavelet scale 1 -1 -4 -10 -20 -15 24 52 24 -15 -20 -10 -4 -1
139 w2 (c2-c3) wavelet scale 2
140
141 256 *w0^2 square of w0 1 16 100 16 1
142 65536*w1^2 square of w1 1 16 100 400 225 576 2704 576 225 400 100 16 1
143 w2^2 square of w2
144 */