001 package util.wavelet;
002 import java.util.List;
003 import static java.lang.Math.*;
004
005 /**
006 *A set of wavelet coefficients and a smoothed array computed using the 'atrous' discrete wavelet trasform on a signal.
007 * In more detail a wavelet space contains instances of the following classes :
008 * <ul>
009 * <li>WaveletSpace : contains
010 * <ol>
011 * <li>Original signal (which also contains a noise estimate for each sample)
012 * <li>Ordered list of {@link Scale} instances</li>
013 * <li>Smoothed signal</li>
014 * <li>Set of {@link WObject} instances</li>
015 * </ol>
016 * </li>
017 * <li>{@link Scale} : contains
018 * <ol>
019 * <li>An array of wavelet coeffients</li>
020 * <li>An array of noise threshold coeffients</li>
021 * <li>a set of Structure instances</li>
022 * </ol>
023 * </li>
024 * </ul>
025 *
026 * <p>A WaveletSpace contains zero or more WObject instances which in turn contains
027 * one or more InterscaleRelationship instances which in turn each contain two Structure instances.
028 * A WaveletSpace is derived from a Signal which has an associated noise estimate.
029 * A WaveletSpace contains a certain number of wavelet {@link Scale} instances.</p>
030 *
031 * <p>The discrete wavelet transform of a signal is computed using the 1D A-Trous method with a b3-spline smoothing function.
032 * The constructor computes the wavelet coefficients a a set of scales plus a smoothed signal,
033 * the 1-sigma threshold for wavelet coefficients, structures, interscale relationships and object trees
034 * for feature extraction and makes available a set of signals corresponding to each individual object
035 * detected.</p>
036 *
037 * <b>Reference:</b> Starck,J.-L., Murtagh, F.: 2002, <i>Astronomical Image and Data Analysis</i>,
038 * ISBN 3540428852, Springer.
039 *
040 * @author John Talbot
041 */
042
043 /*
044
045 Operators and variables from Starck (2002) p.102-103 :
046 <table>
047 <tr><td>Oi </td><td>signal corresponding to object i </td></tr>
048 <tr><td>Wi </td><td>wavelet space restricted to Oi Structures, zero elsewhere</td></tr>
049 <tr><td>W </td><td>wavelet transform </td></tr>
050 <tr><td>W^-1 </td><td>inverse transform </td></tr>
051 <tr><td>Pw </td><td>projection onto Oi Structures </td></tr>
052 <tr><td>A=Pw o W</td><td> wavelet transform then projection onto Oi Structures) </td></tr>
053 <tr><td>A~ </td><td>adjoint of A </td></tr>
054 <tr><td>n </td><td>iteration of reconstruction </td></tr>
055 <tr><td>Oi(n) </td><td>object i (a 1D signal) </td></tr>
056 <tr><td>wr{n} </td><td>wavelet residual </td></tr>
057 <tr><td>R(n) </td><td>residual signal </td></tr>
058 <tr><td>alpha(n)</td><td>covergence parameter alpha </td></tr>
059 <tr><td>beta(n) </td><td>covergence parameter beta </td></tr>
060 </table>
061
062 <table>
063 <tr><th colspan="2">Reconstruction algorithm</th></tr>
064 <tr><td>1</td><td>Intialization step</td><td>Oi(0) = W^-1 Wi <br/>
065 wr(0) = Wi - A(Oi(0))<br/>
066 R(0) = A~(wr(0)) <br/></td></tr>
067 <tr><td>2</td><td></td><td></td></tr>
068 <tr><td>3</td><td></td><td></td></tr>
069 <tr><td>4</td><td></td><td></td></tr>
070 <tr><td>5</td><td></td><td></td></tr>
071 <tr><td>6</td><td></td><td></td></tr>
072 <tr><td>7</td><td>Compute convergence parameter beta</td><td>beta(n+1) = mag()^2 / mag()^2</td></tr>
073 <tr><td>8</td><td>Resisdual image computation </td><td>R(n+1) = A~(wr(n+1) + beta(n+1)R(n)</td></tr>
074 <tr><td>9</td><td>Return to step 2</td><td></td></tr>
075
076 ReconstructedSignal
077
078 WaveletResidual
079
080 ResidualSignal
081
082 </li>
083
084 <li>Computation of convergence parameter : alpha = / A(
085
086 </li>
087
088
089 */
090
091 public class WaveletSpace {
092
093 /** An reference to the original signal. */
094 public Signal signal;
095
096 /**
097 * A list of scales containing wavelet coefficients and optionally noise thresholds.
098 * This representation is useful when using decimating wavelet transforms which
099 * have a decreasing number of wavelet coefficients at larger scales.
100 */
101 public java.util.List<Scale> scales;
102
103 /** The last smoothed array as computed by the atrous wavelet transform. */
104 public SimpleSignal smoothedSignal;
105
106 /** The noise computation method. */
107 NoiseComputation noiseComputation = NoiseComputation.UpperLimit;
108
109 /** A list of WObject instances contained in the wavelet space */
110 public java.util.List<WObject> wObjects;
111
112 /**
113 * Construct a discrete wavelet transform using the 1D A-Trous method with a b3-spline
114 * smoothing function. (Starck 2002 p.29). Boundary conditions are handled by the signal.
115 *
116 * @param signal the signal to be wavelet transformed
117 */
118 public WaveletSpace(Signal signal) {
119 this.signal = signal;
120 int maximumScale = signal.estimateNumberOfScales();
121 scales = new java.util.ArrayList<Scale>(maximumScale);
122 smoothedSignal = new SimpleSignal(signal);
123 int n = signal.size();
124 for (int scaleNumber = 0; scaleNumber < maximumScale; scaleNumber++) {
125 Scale scale = new Scale(smoothedSignal, scaleNumber);
126 int step = (int)(pow(2.0d, (double) scaleNumber) + 0.5d);
127 int step2= step * 2;
128 for (int i = 0; i < n; i++) { // smooth with kernel (1/16, 1/4, 3/8, 1/4, 1/16)
129 smoothedSignal.setValue(i, 0.375f * scale.getValue(i)
130 + 0.25f * ( scale.getValue(i-step) + scale.getValue(i+step))
131 + 0.0625f * ( scale.getValue(i-step2) + scale.getValue(i-step2)));
132 }
133 scale.subtract(smoothedSignal);
134 scales.add(scale);
135 }
136 noiseComputation.compute(signal, scales);
137 }
138
139 /**
140 * Construct the WaveletSpace given a signal and a uniform noise value.
141 * @param signal the signal to be wavelet transformed
142 * @param uniformNoise a uniform noise standard deviation
143 */
144 public WaveletSpace(Signal signal, float uniformNoise) {
145 this(signal);
146 // compute noise if (this.noise == null) computeNoise();
147 }
148
149 /**
150 * Get the specified scale as computed by the atrous wavelet transform.
151 * @param scaleNumber the scale number
152 * @return the specified scale as an instance of simple signal
153 */
154 public Scale getScale(int scaleNumber) {
155 return scales.get(scaleNumber);
156 }
157
158 /**
159 * Set the scale at the specified scale number.
160 * @param scaleNumber the scale number
161 * @param scale the specified scale
162 */
163 public void setScale(int scaleNumber, Scale scale) {
164 scales.set(scaleNumber, scale);
165 }
166
167 /** Use the WObject as a mask to eliminate wavelet coefficients not in the object. */
168 public void mask(WObject wObject) {
169 // TODO:
170 }
171
172 /**
173 * Discrete wavelet transform reconstruction using the 1D A-Trous method with a b3-spline smoothing function.
174 * Reconstruction is computed by summing up all the structures belonging to the WObject plus the last smoothed array.
175 * (see Starck 1998 p.21-26).
176 *
177 * @return the reconstructed signal as an instance of simple signal
178 */
179 public SimpleSignal getReconstruction(WObject wObject) {
180 for (int scale = 0; scale < scales.size(); scale++) {
181 //SimpleSignal coefficients = wObject.getStructures(scale);
182 }
183 return getReconstruction();
184 }
185
186 /**
187 * Discrete wavelet transform reconstruction using the 1D A-Trous method with a b3-spline smoothing function.
188 * Reconstruction is computed by summing up all the scales plus the last smoothed array.
189 * (see Starck 1998 p.21-26).
190 *
191 * @return the reconstructed signal as an instance of simple signal
192 */
193 public SimpleSignal getReconstruction() {
194 SimpleSignal reconstruction = new SimpleSignal(signal.size());
195 for (SimpleSignal scale : scales) reconstruction.add(scale);
196 reconstruction.add(smoothedSignal);
197 return reconstruction;
198 }
199
200 //-----------------methods that belong in superclass-----------------------------------
201
202 /** The maximum number of iterations used to obtain the continuum. */
203 static int MAX_ITERATIONS = 5;
204
205 /**
206 * Get the continuum of a signal. The continuum is computed by iteratively substracting
207 * the smoothed signal from the signal; the iteration is required due to border problems.
208 * (see Starck 1998 p.127).
209 * @return the continuum of the signal
210 */
211 public SimpleSignal getContinuum() {
212 SimpleSignal continuum = new SimpleSignal(getSmoothedSignal());
213 SimpleSignal flattenedSignal = new SimpleSignal(signal.size());
214 for (int iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
215 SimpleSignal.subtract(signal, continuum, flattenedSignal);
216 SimpleWaveletSpace simpleWaveletSpace = new SimpleWaveletSpace(flattenedSignal);
217 continuum.add(simpleWaveletSpace.getSmoothedSignal());
218 }
219 return continuum;
220 }
221
222 /**
223 * Get the last smoothed array as computed by the atrous wavelet transform.
224 * @return the last smoothed array as a simple signal
225 */
226 public SimpleSignal getSmoothedSignal() {
227 return smoothedSignal;
228 }
229 }