001 package ui.model;
002
003 import ui.recognizer.DataCoordinate;
004 import ui.recognizer.Calibration;
005 import ui.recognizer.Units;
006
007 import java.io.*;
008 import java.util.List;
009 import java.util.ArrayList;
010
011 /** A Spectra is a kind of data model which contains a list of DataCoordinates
012 * which could represent intensity versus wavelength. This list can represents
013 * a continuous spectrum or a discrete line list such as atomic transitions.
014 * A spectra can be contrusted from a pre-existing list of datacoordinates or
015 * from a file. There are several methods which can return the minimum and
016 * maximum values of wavelength and intensity. The wavelength unit is a number
017 * by which the DataCoordinate wavelength value is multiplied; it can be chosen
018 * among a small list of predifined units or can be chosen by the user.
019 *
020 * @author John Talbot
021 */
022 public class Spectra extends ArrayList<DataCoordinate> {
023
024 /** The spectra attribute for the object's name (should use FITS header conventions for these) */
025 public static final String ATTRIBUTE_NAME = "NAME";
026
027 /** The spectra attribute for the object's name (should use FITS header conventions for these) */
028 public static final String ATTRIBUTE_WAVELENGTH_UNIT = "WAVELENGTH_UNIT";
029
030 static final int INITIAL_CAPACITY = 1000; // The initial capacity of the ArrayList
031
032 // special characters used in the spectra data file
033 protected char separator = ' ';
034 protected char commentChar = '#';
035
036 // the data ranges (will be different for most plots)
037 protected double minX = Double.MAX_VALUE;
038 protected double maxX = -Double.MAX_VALUE;
039 protected double minY = Double.MAX_VALUE;
040 protected double maxY = -Double.MAX_VALUE;
041
042 protected Calibration calibration = new Calibration();
043
044 // the wavelenght unit for the data coordinates
045 protected double wavelengthUnit = Units.ANGSTROMS;
046
047 // true if there is an equal difference between successive wavelength samples
048 protected boolean uniformSampling = true;
049
050 protected java.util.HashMap attributes = new java.util.HashMap();
051
052 /** Create a spectra by loading file with the specified file name.
053 * @param aFileName a file name pointing to file containing the spectra data
054 */
055 public Spectra(String aFileName) {
056 this(new File(aFileName));
057 }
058
059 /** Create a spectra by loading a file.
060 * @param aFile a file containing the spectra data
061 */
062 public Spectra(File aFile) {
063 this(aFile, Units.ANGSTROMS);
064 }
065
066 /** Create a spectra by loading file with the specified wavelength unit.
067 * @param aFile a file containing the spectra data
068 * @param aWavelengthUnit the wavelength unit, all x coordinates in the data will be multipled by this number.
069 */
070 public Spectra(File aFile, double aWavelengthUnit) {
071 this(aWavelengthUnit);
072 setAttribute(ATTRIBUTE_NAME, aFile.getName());
073 load(aFile);
074 }
075
076 /** Create a spectra using an array of DataCoordinate instances of the specified wavelength unit.
077 *
078 * @param aData an array of data coordinates representing spectra points
079 * @param aWavelengthUnit the wavelength unit, all x coordinates in the data will be multipled by this number.
080 */
081 public Spectra(DataCoordinate[] aData, double aWavelengthUnit) {
082 this(java.util.Arrays.asList(aData), aWavelengthUnit);
083 }
084
085 /** Create a spectra using a list of DataCoordinate instances with the specified wavelength unit.
086 * @param aData an ordered collection of data coordinates representing spectra points
087 * @param aWavelengthUnit the wavelength unit, all x coordinates in the data will be multipled by this number.
088 */
089 public Spectra(List<DataCoordinate> aData, double aWavelengthUnit) {
090 this(aWavelengthUnit);
091 addAll(aData);
092 }
093
094 /** Create a blank spectra with the specified wavelength unit.
095 * DataCoordinates can be added later by using the add() method.
096 * @param aWavelengthUnit the wavelength unit, all x coordinates in the data will be multipled by this number.
097 */
098 public Spectra(double aWavelengthUnit) {
099 super(INITIAL_CAPACITY);
100 getCalibration().setUnitX(aWavelengthUnit);
101 }
102
103 /** Create a spectra using another spectra (copy constructor).
104 * @param aSpectra a spectra to be cloned
105 */
106 public Spectra(Spectra aSpectra) {
107 // Should use deep clone() because only references to DataCoordinates are used
108 // if we invoke scale(), the parent Spectra will also be affected.
109 super(INITIAL_CAPACITY);
110 setWavelengthUnit(aSpectra.getWavelengthUnit());
111 addAll(aSpectra);
112 Object[] keys = aSpectra.getAttributeKeys();
113 for (int i = 0; i < keys.length; i++) {
114 setAttribute(keys[i], aSpectra.getAttribute(keys[i]) );
115 }
116 }
117
118 /** Load the spectra from a file and interpret the wavelength using the wavelengthUnit.
119 * @param aFile a file containing xy pairs, one per line, x and y separatated by a space
120 */
121 public void load(File aFile) {
122 try {
123 LineNumberReader in = new LineNumberReader(new FileReader(aFile));
124 String line = in.readLine();
125 double dataX, dataY;
126 clear();
127
128 // variables required for determining uniform sampling (uniform is the default default)
129 boolean firstPoint = true;
130 double lastX = 0.0;
131 double dx = -1.0;
132 // how much dx must deviate from average before sampling is considered non-uniform
133 double dxFudgeFactor = 5.0e-3;
134
135 while (line != null) {
136 line = line.trim();
137 //String[] numbers = line.split("\\s+"); // one or more spaces between numbers
138 int whiteSpacePosition = line.indexOf(getSeparator());
139 int commentPosition = line.indexOf(getCommentChar());
140 if (commentPosition != -1) { // extract attributes from comment lines
141 String attributeLine = line.substring(commentPosition+1).trim();
142 String attributeKey;
143 String attributeValue;
144 int attributeSeparator;
145 if ((attributeSeparator = attributeLine.indexOf("=")) > -1) { // FITS attribute
146 attributeKey = attributeLine.substring(0, attributeSeparator).trim();
147 attributeValue = attributeLine.substring(attributeSeparator+1).trim();
148 } else {
149 attributeSeparator = attributeLine.indexOf(getSeparator());
150 attributeKey = attributeLine.substring(0, attributeSeparator).trim();
151 attributeValue = attributeLine.substring(attributeSeparator+1).trim();
152 }
153 if (attributeKey != null && attributeKey.equals(ATTRIBUTE_WAVELENGTH_UNIT)) {
154 setWavelengthUnit(Double.parseDouble(attributeValue));
155 } else {
156 setAttribute(attributeKey, attributeValue);
157 }
158 } else if (whiteSpacePosition != -1) { // extract x y pair from line
159 try {
160 // TODO: make the parsing more robust
161 // 1. ignore junk after the second number
162 // 2. skip multiple spaces between numbers or at beginning of line
163
164 dataX = Double.parseDouble(line.substring(0, whiteSpacePosition));
165
166 // is sampling uniform ?
167 if (firstPoint)
168 firstPoint = false;
169 else if (isUniformSampling()) {
170 double newdx = dataX - lastX;
171 if (dx < 0.0)
172 dx = newdx;
173 else if ( Math.abs((newdx - dx) / dx) > dxFudgeFactor)
174 setUniformSampling(false);
175 }
176 String restOfLine = line.substring(whiteSpacePosition+1).trim();
177 int secondWhiteSpacePosition = restOfLine.indexOf(getSeparator());
178 String secondNumber = restOfLine;
179 if (secondWhiteSpacePosition != -1) {
180 secondNumber = restOfLine.substring(0, secondWhiteSpacePosition);
181 //System.out.println(secondNumber + ":" + secondWhiteSpacePosition);
182 }
183 dataY = Double.parseDouble(secondNumber);
184 minX = Math.min(minX, dataX);
185 maxX = Math.max(maxX, dataX);
186 //System.out.println("x=" + dataX + " maxX=" + maxX + " minX=" + minX);
187 minY = Math.min(minY, dataY);
188 maxY = Math.max(maxY, dataY);
189 DataCoordinate entry = new DataCoordinate(dataX, dataY);
190 add(entry);
191 lastX = dataX;
192 } catch(NumberFormatException e) {
193 System.out.println("Could not parse number on line :" + in.getLineNumber());
194 System.out.println(line);
195 }
196 }
197 line = in.readLine();
198 }
199 // can invoke map() on the calibration object to convert from wavelength to pixel and vice-versa
200 // the rectangle represents this spectra if it were one pixel high and size() pixels wide
201 if (maxY == minY) minY = 0; // kludge to prevent divide by zero errors
202 if (maxX == minX) minX = 0; // kludge to prevent divide by zero errors
203 if (maxY == 0) maxY = 1.0d; // kludge to prevent divide by zero errors
204 setCalibration(new Calibration(new java.awt.Rectangle(0, 0, size(), 1),
205 minX, minY, maxX, maxY, getWavelengthUnit(), 1.0d/maxY));
206 in.close();
207 } catch (IOException ioe) {
208 ioe.printStackTrace();
209 }
210 }
211
212 /** Save the spectra as a file with properties encoded as comments.
213 * @param aFile a file in which xy pairs will be added, one per line, x and y separatated by a space
214 */
215 public void save(File aFile) {
216 try {
217 PrintStream out = new PrintStream((OutputStream) new FileOutputStream(aFile));
218 Object[] keys = getAttributeKeys();
219 for (int i = 0; i < keys.length; i++) {
220 out.println(getCommentChar() + " " + keys[i] + " " + getAttribute(keys[i]) );
221 }
222 for (java.util.Iterator iterator = this.iterator(); iterator.hasNext(); ) {
223 out.println(iterator.next());
224 }
225 out.close();
226 } catch (IOException ioe) {
227 ioe.printStackTrace();
228 }
229 }
230
231 /** Append the contents of given spectra to this spectra.
232 * If regions overlap then the new spectra will replace the old.
233 * The added spectra does not have to be sequential, it come before this spectra.
234 * A future version of this method could optionally scale the source spectra to more
235 * smoothly 'sew' the two spectra together in the overlap region.
236 * A future version of this method could contain a scaling value if the other spectra differs in vertical scale.
237 * @param aSpectra a spectra which will be appended to this spectra
238 */
239 public void append(Spectra aSpectra) {
240 if (getCalibration().getMaxX() < aSpectra.getCalibration().getMinX()) { // new spectra is placed after this spectra
241 addAll(aSpectra);
242 } else if(getCalibration().getMinX() > aSpectra.getCalibration().getMaxX()) { // new spectra is placed before this spectra
243 ArrayList<DataCoordinate> temp = new ArrayList<DataCoordinate>(this);
244 clear();
245 addAll(aSpectra);
246 addAll(temp);
247 } else { //new spectra overlaps this spectra
248 // TODO: make new spectra replace the old in the overlap region.
249 // Optionally scale the source spectra to more smoothly 'sew' the two spectra together.
250 }
251 }
252
253 /** Append the contents of given spectra to this spectra.
254 * If regions overlap then the new spectra will replace the old.
255 * The added spectra does not have to be sequential, it come before this spectra.
256 * A future version of this method could optionally scale the source spectra to more
257 * smoothly 'sew' the two spectra together in the overlap region.
258 * @param aSpectra a spectra which will be appended to this spectra
259 * @param aScale a scaling value to be used if the other spectra differs in vertical scale.
260 */
261 public void append(Spectra aSpectra, double aScale) {
262 Spectra tempSpectra = new Spectra(aSpectra);
263 tempSpectra.scaleIntensity(aScale);
264 append(tempSpectra);
265 }
266
267 /** Scale the spectra intensity by the given factor.
268 * @param aScale a scale factor to multiply the spectra intensity values.
269 */
270 public void scaleIntensity(double aScale) {
271 scale(1.0d, aScale);
272 }
273
274 /** Scale the spectra wavelength by the given factor.
275 * @param aScale a scale factor to multiply the spectra wavelength values.
276 */
277 public void scaleWavelength(double aScale) {
278 scale(aScale, 1.0d);
279 }
280
281 /** Scale the spectra wavelength and intensity by the given factors.
282 * @param wavelengthScale a scale factor to multiply the spectra wavelength values.
283 * @param intensityScale a scale factor to multiply the spectra intensity values.
284 */
285 public void scale(double wavelengthScale, double intensityScale) {
286 for (java.util.Iterator iterator = this.iterator(); iterator.hasNext(); ) {
287 DataCoordinate data = (DataCoordinate) iterator.next();
288 data.setLocation(data.getX() * wavelengthScale, data.getY() * intensityScale);
289 }
290 }
291
292 /** Get the spectra wavelength at the given integer position.
293 * @param position the integer position at which to obtain the wavelength
294 */
295 public double getWavelength(int position) {
296 DataCoordinate data = (DataCoordinate) get(position);
297 return data.getX();
298 }
299
300 /** Get the spectra intensity at the given integer position.
301 * @param position the integer position at which to obtain the intensity
302 */
303 public double getIntensity(int position) {
304 DataCoordinate data = (DataCoordinate) get(position);
305 return data.getY();
306 }
307
308 /** Get the spectra intensity given an arbitrary wavelength, interpolate if necessary.
309 *
310 * @param aWavelength the wavelength position at which the intensity is obtained
311 */
312 public double getIntensity(double aWavelength) {
313 double intensity = 0.0;
314 // Linear interpolation
315 if (isUniformSampling()) {
316 double dx = (maxX - minX) / (size() - 1); // sampling interval
317 double n = (aWavelength - minX) / dx;
318 int n1 = (int) Math.floor(n);
319 int n2 = (int) Math.ceil(n);
320 if (n1 >= 0 && n2 < size()) {
321 //System.out.println("size()=" + size() + " maxX=" + maxX + " minX=" + minX + "dx=" + dx + " n=" + n + " n1=" + n1 + " n2=" + n2 + " wavelength=" + aWavelength);
322 intensity = interpolate(n1, n2, aWavelength);
323 }
324 } else { // non-uniformly sampled spectra
325 // TODO: FIXME binary search to locate position in spectra
326 if (aWavelength > minX && aWavelength < maxX) {
327 int guess = size() / 2; // initial guess for wavelength position
328 boolean located = false; // flag indicating that we have found the bin
329 int interval = guess;
330 //System.out.println("Wavelength=" + aWavelength);
331 while (!located) {
332 interval = interval / 2;
333 double x = getWavelength(guess);
334 //System.out.println(" Guess=" + guess + " x=" + x + " dx=" + (x-aWavelength));
335 if (x == aWavelength || interval == 0)
336 located = true;
337 else if (x > aWavelength)
338 guess = guess - interval;
339 else if (x < aWavelength)
340 guess = guess + interval;
341 }
342 // find the points between which the wavelength is located
343 double x = getWavelength(guess);
344
345 if (x == aWavelength) {
346 intensity = getIntensity(guess);
347 } else {
348 int point1 = guess;
349 int point2 = guess;
350 if (x < aWavelength)
351 point2 = guess + 1;
352 else if (x > aWavelength)
353 point1 = guess - 1;
354
355 //KLUDGE: shift bracketting points slightly to make aWavelength fall within interval
356 double x1 = getWavelength(point1);
357 double x2 = getWavelength(point2);
358 boolean outside = aWavelength < x1 || aWavelength > x2;
359 int increment = 1;
360 if (aWavelength < x1)
361 increment = -1;
362 else if(aWavelength > x2)
363 increment = 1;
364
365 while(outside) {
366 point1 += increment;
367 point2 += increment;
368 x1 = getWavelength(point1);
369 x2 = getWavelength(point2);
370 outside = aWavelength < x1 || aWavelength > x2;
371 }
372 //END_KLUDGE
373
374 intensity = interpolate(point1, point2, aWavelength);
375 }
376
377 } else if (aWavelength == minX) {
378 intensity = getIntensity(0);
379 } else if (aWavelength == maxX) {
380 intensity = getIntensity(size()-1);
381 }
382
383 }
384 return intensity;
385 }
386
387 /** Interpolate linearly the intensity at a specified wavelength
388 * position between two data points.
389 * Preconditions:
390 * (1) The separation between points 1 and 2 must be greater than or equal to 1.
391 * (2) The wavelength must lie within the interval
392 * @param point1 integer representing the first point
393 * @param point2 integer representing the second point
394 * @param aWavelength
395 */
396 protected double interpolate(int point1, int point2, double aWavelength) {
397 double x1 = getWavelength(point1);
398 double x2 = getWavelength(point2);
399 double y1 = getIntensity(point1);
400 double y2 = getIntensity(point2);
401 double intensity;
402 if (x1 == x2)
403 intensity = y1; // pre-condition failure (should not happen) TODO: throw exception here
404 else {
405 double slope = (y2 - y1)/(x2 - x1);
406 intensity = slope * (aWavelength - x1) + y1;
407 }
408 return intensity;
409 }
410
411 /** Get the intensity minus the continuum. The spectra is obtained by substracting
412 * a polynomical which approximate the average continuum.
413 * This creates a flattened spectrum in wich emission lines dominate.
414 * @param aWavelength
415 */
416 public double getIntensityMinusContinuum(double aWavelength) {
417 return getIntensity(aWavelength) - getContinuum(aWavelength);
418 }
419
420 /** The order of the polynomial used to fit the continuum */
421 static final int ORDER = 6;
422
423 /** An internal flag which indicates that the coefficients to
424 * continuum polynomial need to be computed.
425 */
426 protected boolean continuumFitted = false;
427
428 /** The polynomial fit coefficients to the continuum */
429 protected double[] continuum = new double[ORDER];
430
431 /** Get the continuum of the spectra, obtained by fitting a
432 * polynomial to the regions of the spectra without emission lines.
433 * @param aWavelength
434 */
435 public double getContinuum(double aWavelength) {
436 if (!continuumFitted) {
437 // compute the coefficients to the continuum polynomial
438 //
439 // continuumFitted = true;
440 }
441 // evaluate the 6th order polynomial continuum fit at aWavelength
442 return 0.0;
443 }
444
445 protected double[] fitPolynomial() {
446
447 return new double[6];
448 }
449
450 //------Accessor methods-----------------------------------------------------------------
451
452 /** Get the data coordinate for a given index in the spectra array.
453 *
454 * @param anIndex the position in the spectral point array
455 */
456 public DataCoordinate getDataCoordinate(int anIndex) {
457 return get(anIndex);
458 }
459
460 /** Get the data coordinate scaled to unit intensity for a given index in the spectra array and with the wavelength in meters.
461 * The units are obtained from the calibration object in this spectra.
462 * @param anIndex the position in the spectral point array
463 */
464 public DataCoordinate getScaledDataCoordinate(int anIndex) {
465 DataCoordinate data = new DataCoordinate(getDataCoordinate(anIndex)); // make a copy
466 data.scale(getCalibration());
467 return data;
468 }
469
470 /** Set the wavelength unit, all x coordinates in the data will be multipled by this number.
471 * @param aWavelengthUnit the wavelength unit
472 */
473 private void setWavelengthUnit(double aWavelengthUnit) {
474 wavelengthUnit = aWavelengthUnit;
475 setAttribute(ATTRIBUTE_WAVELENGTH_UNIT, new Double(aWavelengthUnit));
476 if (getCalibration() != null) {
477 getCalibration().setUnitX(aWavelengthUnit);
478 }
479 }
480
481 /** Get the wavelength unit, all x coordinates in the data will be multipled by this number.
482 * @return the wavelength unit in double precision
483 */
484 private double getWavelengthUnit() {
485 return wavelengthUnit;
486 }
487
488 /** Get the calibration for this spectra, required when wavelength to pixel/point conversion is required
489 */
490 public Calibration getCalibration() {
491 return calibration;
492 }
493
494 /** Set the calibration for this spectra, required when wavelength to pixel/point conversion is required
495 * @param aCalibration a calibration object
496 */
497 public void setCalibration(Calibration aCalibration) {
498 calibration = aCalibration;
499 }
500
501 /** Get the separator character used between the wavelength and the intensity in the the spectra file.
502 * @return a character used to separate the wavelength and intensity
503 */
504 public char getSeparator() {
505 return separator;
506 }
507
508 /** Set the separator character used between the wavelength and the intensity in the the spectra file.
509 * @param aSeparator a character used to separate the wavelength and intensity
510 */
511 public void setSeparator(char aSeparator) {
512 separator = aSeparator;
513 }
514
515 /** Get the comment character which indicates a comment line to be ignored.
516 * A future option could be to parse comment lines for FITS headers.
517 * @return a character used to mark a comment line
518 */
519 public char getCommentChar() {
520 return commentChar;
521 }
522
523 /** Set the comment character which indicates a comment line to be ignored.
524 * @param aCommentChar a character used to mark a comment line
525 */
526 public void setCommentChar(char aCommentChar) {
527 commentChar = aCommentChar;
528 }
529
530 /** Get all the attribute keys.
531 * @return an array or keys one for each attribute
532 */
533 public Object[] getAttributeKeys() {
534 return attributes.keySet().toArray();
535 }
536
537 /** Get a spectra attribute using a key.
538 * @param aKey an object representing a key
539 * @return an object which represents the value corresponding to the key
540 */
541 public Object getAttribute(Object aKey) {
542 return attributes.get(aKey);
543 }
544
545 /** Set an attribute using a key and its corresponding value.
546 * @param aKey an object representing a key
547 * @param aValue an object which represents the value corresponding to the key
548 */
549 public void setAttribute(Object aKey, Object aValue) {
550 attributes.put(aKey, aValue);
551 }
552
553 /** Is the wavelength sampling interval uniform. This should
554 * be automatically detected from the data file by substracting
555 * subsequent wavelength samples and determining that their separation
556 * is the same throughout the spectra.
557 */
558 public boolean isUniformSampling() {
559 return uniformSampling;
560 }
561
562 public void setUniformSampling(boolean aUniformSampling) {
563 uniformSampling = aUniformSampling;
564 }
565
566
567 }