Astronomical Image and Data Analysis

Jupsat Pro Astronomy Software

Secrets of the Deep Sky

Get Instant Access

Second Edition With 119 Figures

Springer

Jean-Luc Starck Service d'Astrophysique CEA/Saclay Orme des Merisiers, Bat 709 91191 Gif-sur-Yvette Cedex, France

Fionn Murtagh

Dept. Computer Science Royal Holloway University of London Egham, Surrey TW20 0EX, UK

Cover picture: The cover image to this 2nd edition is from the Deep Impact project. It was taken approximately 8 minutes after impact on 4 July 2005 with the CLEAR6 filter and deconvolved using the Richardson-Lucy method. We thank Don Lindler, Ivo Busko, Mike A' Hearn and the Deep Impact team for the processing of this image and for providing it to us.

Library of Congress Control Number: 2006930922

ISSN 0941-7834

ISBN-10 3-540-33024-0 2nd Edition Springer Berlin Heidelberg New York ISBN-13 978-3-540-33024-0 2nd Edition Springer Berlin Heidelberg New York ISBN 3-540-42885-2 1st Edition Springer Berlin Heidelberg New York

This work is subject to copyright. All rights are reserved, whether the whole or part of the material is concerned, specifically the rights of translation, reprinting, reuse of illustrations, recitation, broadcasting, reproduction on microfilm or in any other way, and storage in data banks. Duplication of this publication or parts thereof is permitted only under the provisions of the German Copyright Law of September 9, 1965, in its current version, and permission for use must always be obtained from Springer. Violations are liable to prosecution under the German Copyright Law.

Springer is a part of Springer Science+Business Media springer.com

Springer-Verlag Berlin Heidelberg 2006

The use of general descriptive names, registered names, trademarks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

Typesetting: by the authors

Final layout: Data conversion and production by LE-TEX Jelonek, Schmidt & VocklerGbR, Leipzig, Germany Cover design: design & production GmbH, Heidelberg

Printed on acid-free paper SPIN: 11595496 55/3141 - 5 4 3 2 1 0

Preface to the Second Edition

This book presents material which is more algorithmically oriented than most alternatives. It also deals with topics that are at or beyond the state of the art. Examples include practical and applicable wavelet and other multiresolution transform analysis. New areas are broached like the ridgelet and curvelet transforms. The reader will find in this book an engineering approach to the interpretation of scientific data.

Compared to the 1st Edition, various additions have been made throughout, and the topics covered have been updated. The background or environment of this book's topics include continuing interest in e-science and the virtual observatory, which are based on web based and increasingly web service based science and engineering.

Additional colleagues whom we would like to acknowledge in this 2nd edition include: Bedros Afeyan, Nabila Aghanim, Emmanuel Candes, David Donoho, Jalal Fadili, and Sandrine Pires, We would like to particularly acknowledge Olivier Forni who contributed to the discussion on compression of hyperspectral data, Yassir Moudden on multiwavelength data analysis and Vicent Martinez on the genus function.

The cover image to this 2nd edition is from the Deep Impact project. It was taken approximately 8 minutes after impact on 4 July 2005 with the CLEAR6 filter and deconvolved using the Richardson-Lucy method. We thank Don Lindler, Ivo Busko, Mike A'Hearn and the Deep Impact team for the processing of this image and for providing it to us.

Paris, London June, 2006

Jean-Luc Starck Fionn Murtagh

Preface to the First Edition

When we consider the ever increasing amount of astronomical data available to us, we can well say that the needs of modern astronomy are growing by the day. Ever better observing facilities are in operation. The fusion of information leading to the coordination of observations is of central importance.

The methods described in this book can provide effective and efficient ripostes to many of these issues. Much progress has been made in recent years on the methodology front, in line with the rapid pace of evolution of our technological infrastructures.

The central themes of this book are information and scale. The approach is astronomy-driven, starting with real problems and issues to be addressed. We then proceed to comprehensive theory, and implementations of demonstrated efficacy.

The field is developing rapidly. There is little doubt that further important papers, and books, will follow in the future.

Colleagues we would like to acknowledge include: Alexandre Aussem, Albert Bijaoui, Francois Bonnarel, Jonathan G. Campbell, Ghada Jammal, Rene Gastaud, Pierre-Francois Honore, Bruno Lopez, Mireille Louys, Clive Page, Eric Pantin, Philippe Querre, Victor Racine, Jerome Rodriguez, and Ivan Valtchanov.

The cover image is from Jean-Charles Cuillandre. It shows a five minute exposure (5 60-s dithered and stacked images), R filter, taken with CFH12K wide field camera (100 million pixels) at the primary focus of the CFHT in July 2000. The image is from an extremely rich zone of our Galaxy, containing star formation regions, dark nebulae (molecular clouds and dust regions), emission nebulae (Ha), and evolved stars which are scattered throughout the field in their two-dimensional projection effect. This zone is in the constellation of Saggitarius.

Paris, Belfast June, 2002

Jean-Luc Starck Fionn Murtagh

Table of Contents

1. Introduction to Applications and Methods 1

1.1 Introduction 1

1.2 Transformation and Data Representation 3

1.2.1 Fourier Analysis 5

1.2.2 Time-Frequency Representation 6

1.2.3 Time-Scale Representation: The Wavelet Transform . . 9

1.2.4 The Radon Transform 12

1.2.5 The Ridgelet Transform 12

1.2.6 The Curvelet Transform 14

1.3 Mathematical Morphology 15

1.4 Edge Detection 18

1.4.1 First Order Derivative Edge Detection 18

1.4.2 Second Order Derivative Edge Detection 20

1.5 Segmentation 23

1.6 Pattern Recognition 24

1.7 Chapter Summary 27

2. Filtering 29

2.1 Introduction 29

2.2 Multiscale Transforms 31

2.2.1 The A Trous Isotropic Wavelet Transform 31

2.2.2 Multiscale Transforms Compared to Other Data Transforms 33

2.2.3 Choice of Multiscale Transform 36

2.2.4 The Multiresolution Support 37

2.3 Significant Wavelet Coefficients 38

2.3.1 Definition 38

2.3.2 Noise Modeling 39

2.3.3 Automatic Estimation of Gaussian Noise 40

2.3.4 Detection Level Using the FDR 48

2.4 Filtering and Wavelet Coefficient Thresholding 50

2.4.1 Thresholding 50

2.4.2 Iterative Filtering 51

2.4.3 Other Wavelet Denoising Methods 52

2.4.4 Experiments 54

2.4.5 Iterative Filtering with a Smoothness Constraint 56

2.5 Filtering from the Curvelet Transform 57

2.5.1 Contrast Enhancement 57

2.5.2 Curvelet Denoising 59

2.5.3 The Combined Filtering Method 61

2.6 Haar Wavelet Transform and Poisson Noise 63

2.6.1 Haar Wavelet Transform 63

2.6.2 Poisson Noise and Haar Wavelet Coefficients 64

2.6.3 Experiments 67

2.7 Chapter Summary 70

3. Deconvolution 71

3.1 Introduction 71

3.2 The Deconvolution Problem 74

3.3 Linear Regularized Methods 75

3.3.1 Least Squares Solution 75

3.3.2 Tikhonov Regularization 75

3.3.3 Generalization 76

3.5 Bayesian Methodology 78

3.5.1 Definition 78

3.5.2 Maximum Likelihood with Gaussian Noise 79

3.5.3 Gaussian Bayes Model 79

3.5.4 Maximum Likelihood with Poisson Noise 80

3.5.5 Poisson Bayes Model 81

3.5.6 Maximum Entropy Method 81

3.5.7 Other Regularization Models 82

3.6 Iterative Regularized Methods 84

3.6.1 Constraints 84

3.6.2 Jansson-Van Cittert Method 85

3.6.3 Other Iterative Methods 85

3.7 Wavelet-Based Deconvolution 86

3.7.1 Introduction 86

3.7.2 Wavelet-Vaguelette Decomposition 87

3.7.3 Regularization from the Multiresolution Support 90

3.7.4 Wavelet CLEAN 93

3.7.5 The Wavelet Constraint 98

3.8 Deconvolution and Resolution 104

3.9 Super-Resolution 105

3.9.1 Definition 105

3.9.2 Gerchberg-Saxon Papoulis Method 106

3.9.3 Deconvolution with Interpolation 107

3.9.4 Undersampled Point Spread Function 107

3.10 Conclusions and Chapter Summary 109

4. Detection 111

4.1 Introduction 111

4.2 From Images to Catalogs 112

4.3 Multiscale Vision Model 116

4.3.1 Introduction 116

4.3.2 Multiscale Vision Model Definition 117

4.3.3 From Wavelet Coefficients to Object Identification 117

4.3.4 Partial Reconstruction 120

4.3.5 Examples 122

4.3.6 Application to ISOCAM Data Calibration 122

4.4 Detection and Deconvolution 126

4.5 Detection in the Cosmological Microwave Background 130

4.5.1 Introduction 130

4.5.2 Point Sources on a Gaussian Background 132

4.5.3 Non-Gaussianity 132

4.6 Conclusion 135

4.7 Chapter Summary 135

5. Image Compression 137

5.1 Introduction 137

5.2 Lossy Image Compression Methods 139

5.2.1 The Principle 139

5.2.2 Compression with Pyramidal Median Transform 140

5.2.3 PMT and Image Compression 142

5.2.4 Compression Packages 145

5.2.5 Remarks on these Methods 146

5.2.6 Other Lossy Compression Methods 148

5.3 Comparison 149

5.3.1 Quality Assessment 149

5.3.2 Visual Quality 150

5.3.3 First Aladin Project Study 151

5.3.4 Second Aladin Project Study 155

5.3.5 Computation Time 159

5.3.6 Conclusion 160

5.4 Lossless Image Compression 161

5.4.1 Introduction 161

5.4.2 The Lifting Scheme 161

5.4.3 Comparison 166

5.5 Large Images: Compression and Visualization 167

5.5.1 Large Image Visualization Environment: LIVE 167

5.5.2 Decompression by Scale and by Region 168

5.5.3 The SAO-DS9 LIVE Implementation 169

5.6 Hyperspectral Compression for Planetary Space Missions 170

5.7 Chapter Summary 173

6. Multichannel Data 175

6.1 Introduction 175

6.2 The Wavelet-Karhunen-Loeve Transform 176

6.2.1 Definition 176

6.2.2 Correlation Matrix and Noise Modeling 178

6.2.3 Scale and Karhunen-Loeve Transform 179

6.2.4 The WT-KLT Transform 179

6.2.5 The WT-KLT Reconstruction Algorithm 180

6.3 Noise Modeling in the WT-KLT Space 180

6.4 Multichannel Data Filtering 181

6.4.1 Introduction 181

6.4.2 Reconstruction from a Subset of Eigenvectors 181

6.4.3 WT-KLT Coefficient Thresholding 183

6.4.4 Example: Astronomical Source Detection 183

6.5 The Haar-Multichannel Transform 183

6.6 Independent Component Analysis 184

6.6.1 Definition 184

6.6.3 FastICA 186

6.7 CMB Data and the SMICA ICA Method 189

6.7.1 The CMB Mixture Problem 189

6.8 ICA and Wavelets 193

6.8.2 Covariance Matching in Wavelet Space: WSMICA ... 194

6.8.3 Numerical Experiments 195

6.9 Chapter Summary 198

7. An Entropic Tour of Astronomical Data Analysis 201

7.1 Introduction 201

7.2 The Concept of Entropy 204

7.3 Multiscale Entropy 210

7.3.1 Definition 210

7.3.2 Signal and Noise Information 212

7.4 Multiscale Entropy Filtering 215

7.4.1 Filtering 215

7.4.2 The Regularization Parameter 215

7.4.3 Use of a Model 217

7.4.4 The Multiscale Entropy Filtering Algorithm 218

7.4.5 Optimization 219

7.4.6 Examples 220

7.5 Deconvolution 220

7.5.1 The Principle 220

7.5.2 The Parameters 224

7.5.3 Examples 225

7.6 Multichannel Data Filtering 225

7.7 Relevant Information in an Image 228

7.8 Multiscale Entropy and Optimal Compressibility 230

7.9 Conclusions and Chapter Summary 231

8. Astronomical Catalog Analysis 233

8.1 Introduction 233

8.2 Two-Point Correlation Function 234

8.2.1 Introduction 234

8.2.2 Determining the 2-Point Correlation Function 235

8.2.3 Error Analysis 236

8.2.4 Correlation Length Determination 237

8.2.5 Creation of Random Catalogs 237

8.2.6 Examples 238

8.2.7 Limitation of the Two-Point Correlation Function: Toward Higher Moments 242

8.3 The Genus Curve 245

8.4 Minkowski Functionals 247

8.5 Fractal Analysis 249

8.5.1 Introduction 249

8.5.2 The Hausdorff and Minkowski Measures 250

8.5.3 The Hausdorff and Minkowski Dimensions 251

8.5.4 Multifractality 251

8.5.5 Generalized Fractal Dimension 253

8.5.6 Wavelets and Multifractality 253

8.6 Spanning Trees and Graph Clustering 257

8.7 Voronoi Tessellation and Percolation 259

8.8 Model-Based Clustering 260

8.8.1 Modeling of Signal and Noise 260

8.8.2 Application to Thresholding 262

8.9 Wavelet Analysis 263

8.10 Nearest Neighbor Clutter Removal 265

8.11 Chapter Summary 266

9. Multiple Resolution in Data Storage and Retrieval 267

9.1 Introduction 267

9.2 Wavelets in Database Management 267

9.3 Fast Cluster Analysis 269

9.4 Nearest Neighbor Finding on Graphs 271

9.5 Cluster-Based User Interfaces 272

9.6 Images from Data 273

9.6.1 Matrix Sequencing 273

9.6.2 Filtering Hypertext 277

9.6.3 Clustering Document-Term Data 278

9.7 Chapter Summary 282

10. Towards the Virtual Observatory 285

10.1 Data and Information 285

10.2 The Information Handling Challenges Facing Us 287

Appendix

A. A Trous Wavelet Transform 291

B. Picard Iteration 297

C. Wavelet Transform Using the Fourier Transform 299

D. Derivative Needed for the Minimization 303

E. Generalization of the Derivative

Needed for the Minimization 307

F. Software and Related Developments 309

Bibliography 311

Index 331

1. Introduction to Applications and Methods

1.1 Introduction

"May you live in interesting times!" ran the old Chinese wish. The early years of the third millennium are interesting times for astronomy, as a result of the tremendous advances in our computing and information processing environment and infrastructure. The advances in signal and image processing methods described in this book are of great topicality as a consequence. Let us look at some of the overriding needs of contemporary observational astronomical.

Unlike in Earth observation or meteorology, astronomers do not want to interpret data and, having done so, delete it. Variable objects (supernovae, comets, etc.) bear witness to the need for astronomical data to be available indefinitely. The unavoidable problem is the sheer overwhelming quantity of data which is now collected. The only basis for selective choice for what must be kept long-term is to associate more closely the data capture with the information extraction and knowledge discovery processes. We have got to understand our scientific knowledge discovery mechanisms better in order to make the correct selection of data to keep long-term, including the appropriate resolution and refinement levels.

The vast quantities of visual data collected now and in the future present us with new problems and opportunities. Critical needs in our software systems include compression and progressive transmission, support for differential detail and user navigation in data spaces, and "thinwire" transmission and visualization. The technological infrastructure is one side of the picture.

Another side of this same picture, however, is that our human ability to interpret vast quantities of data is limited. A study by D. Williams, CERN, has quantified the maximum possible volume of data which can conceivably be interpreted at CERN. This points to another more fundamental justification for addressing the critical technical needs indicated above. This is that selective and prioritized transmission, which we will term intelligent streaming, is increasingly becoming a key factor in human understanding of the real world, as mediated through our computing and networking base. We need to receive condensed, summarized data first, and we can be aided in our understanding of the data by having more detail added progressively. A hyperlinked and networked world makes this need for summarization more and more acute. We need to take resolution scale into account in our information and knowledge spaces. This is a key aspect of an intelligent streaming system.

A further area of importance for scientific data interpretation is that of storage and display. Long-term storage of astronomical data, we have already noted, is part and parcel of our society's memory (a formulation due to Michael Kurtz, Center for Astrophysics, Smithsonian Institute). With the rapid obsolescence of storage devices, considerable efforts must be undertaken to combat social amnesia. The positive implication is the ever-increasing complementarity of professional observational astronomy with education and public outreach.

Astronomy's data centers and image and catalog archives play an important role in our society's collective memory. For example, the SIMBAD database of astronomical objects at Strasbourg Observatory contains data on 3 million objects, based on 7.5 million object identifiers. Constant updating of SIMBAD is a collective cross-institutional effort. The MegaCam camera at the Canada-France-Hawaii Telescope (CFHT), Hawaii, is producing images of dimensions 16000 x 16000, 32-bits per pixel. The European Southern Observatory's VLT (Very Large Telescope) is beginning to produce vast quantities of very large images. Increasingly, images of size 1 GB or 2 GB, for a single image, are not exceptional. CCD detectors on other telescopes, or automatic plate scanning machines digitizing photographic sky surveys, produce lots more data. Resolution and scale are of key importance, and so also is region of interest. In multiwavelength astronomy, the fusion of information and data is aimed at, and this can be helped by the use of resolution similar to our human cognitive processes. Processing (calibration, storage and transmission formats and approaches) and access have not been coupled as closely as they could be. Knowledge discovery is the ultimate driver.

Many ongoing initiatives and projects are very relevant to the work described in later chapters.

Image and Signal Processing. The major areas of application of image and signal processing include the following.

— Visualization: Seeing our data and signals in a different light is very often a revealing and fruitful thing to do. Examples of this will be presented throughout this book.

— Filtering: A signal in the physical sciences rarely exists independently of noise, and noise removal is therefore a useful preliminary to data interpretation. More generally, data cleaning is needed, to bypass instrumental measurement artifacts, and even the inherent complexity of the data. Image and signal filtering will be presented in Chapter 2.

— Deconvolution: Signal "deblurring" is used for reasons similar to filtering, as a preliminary to signal interpretation. Motion deblurring is rarely important in astronomy, but removing the effects of atmospheric blurring, or quality of seeing, certainly is of importance. There will be a wide-ranging discussion of the state of the art in deconvolution in astronomy in Chapter 3.

— Compression: Consider three different facts. Long-term storage of astronomical data is important. A current trend is towards detectors accommodating ever-larger image sizes. Research in astronomy is a cohesive but geographically distributed activity. All three facts point to the importance of effective and efficient compression technology. In Chapter 5, the state of the art in astronomical image compression will be surveyed.

— Mathematical morphology: Combinations of dilation and erosion operators, giving rise to opening and closing operations, in boolean images and in greyscale images, allow for a truly very esthetic and immediately practical processing framework. The median function plays its role too in the context of these order and rank functions. Multiple scale mathematical morphology is an immediate generalization. There is further discussion on mathematical morphology below in this chapter.

— Edge detection: Gradient information is not often of central importance in astronomical image analysis. There are always exceptions of course.

— Segmentation and pattern recognition: These are discussed in Chapter 4, dealing with object detection. In areas outside astronomy, the term feature selection is more normal than object detection.

— Multidimensional pattern recognition: General multidimensional spaces are analyzed by clustering methods, and by dimensionality mapping methods. Multiband images can be taken as a particular case. Such methods are pivotal in Chapter 6 on multichannel data, 8 on catalog analysis, and 9 on data storage and retrieval.

— Hough and Radon transforms, leading to 3D tomography and other applications: Detection of alignments and curves is necessary for many classes of segmentation and feature analysis, and for the building of 3D representations of data. Gravitational lensing presents one area of potential application in astronomy imaging, although the problem of faint signal and strong noise is usually the most critical one. Ridgelet and curvelet transforms (discussed below in this chapter) offer powerful generalizations of current state of the art ways of addressing problems in these fields.

A number of outstanding general texts on image and signal processing are available. These include Gonzalez and Woods (1992), Jain (1990), Pratt (1991), Parker (1996), Castleman (1995), Petrou and Bosdogianni (l999), Bovik (2000). A text of ours on image processing and pattern recognition is available on-line (Campbell and Murtagh, 2001). Data analysis texts of importance include Bishop (1995), and Ripley (1995).

1.2 Transformation and Data Representation

Many different transforms are used in data processing, - Haar, Radon, Hadamard, etc. The Fourier transform is perhaps the most widely used. The goal of these transformations is to obtain a sparse representation of the data, and to pack most information into a small number of samples. For example, a sine signal f (t) = sin(2^vt), defined on N pixels, requires only two samples (at frequencies —v and v) in the Fourier domain for an exact representation. Wavelets and related multiscale representations pervade all areas of signal processing. The recent inclusion of wavelet algorithms in JPEG 2000 - the new still-picture compression standard - testifies to this lasting and significant impact. The reason for the success of wavelets is due to the fact that wavelet bases represent well a large class of signals. Therefore this allows us to detect roughly isotropic elements occurring at all spatial scales and locations. Since noise in the physical sciences is often not Gaussian, modeling in wavelet space of many kind of noise - Poisson noise, combination of Gaussian and Poisson noise components, non-stationary noise, and so on - has been a key motivation for the use of wavelets in scientific, medical, or industrial applications. The wavelet transform has also been extensively used in astronomical data analysis during the last ten years. A quick search with ADS (NASA Astrophysics Data System, adswww.harvard.edu) shows that around 500 papers contain the keyword "wavelet" in their abstract, and this holds for all astrophysical domains, from study of the sun through to CMB (Cosmic Microwave Background) analysis:

— Sun: active region oscillations (Ireland et al., 1999; Blanco et al., 1999), determination of solar cycle length variations (Fligge et al., 1999), feature extraction from solar images (Irbah et al., 1999), velocity fluctuations (Lawrence et al., 1999).

— Solar system: asteroidal resonant motion (Michtchenko and Nesvorny, 1996), classification of asteroids (Bendjoya, 1993), Saturn and Uranus ring analysis (Bendjoya et al., 1993; Petit and Bendjoya, 1996).

— Star studies: Ca II feature detection in magnetically active stars (Soon et al., 1999), variable star research (Szatmary et al., 1996).

— Interstellar medium: large-scale extinction maps of giant molecular clouds using optical star counts (Cambresy, 1999), fractal structure analysis in molecular clouds (Andersson and Andersson, 1993).

— Planetary nebula detection: confirmation of the detection of a faint planetary nebula around IN Com (Brosch and Hoffman, 1999), evidence for extended high energy gamma-ray emission from the Rosette/Monoceros Region (Jaffe et al., 1997).

— Galaxy: evidence for a Galactic gamma-ray halo (Dixon et al., 1998).

— QSO: QSO brightness fluctuations (Schild, 1999), detecting the non-Gaussian spectrum of QSO Lya absorption line distribution (Pando and Fang, 1998).

— Gamma-ray burst: GRB detection (Kolaczyk, 1997; Norris et al., 1994) and GRB analysis (Greene et al., 1997; Walker et al., 2000).

— Black hole: periodic oscillation detection (Steiman-Cameron et al., 1997; Scargle, 1997)

— Galaxies: starburst detection (Hecquet et al., 1995), galaxy counts (Aus-sel et al., 1999; Damiani et al., 1998), morphology of galaxies (Weistrop et al., 1996; Kriessler et al., 1998), multifractal character of the galaxy distribution (Martinez et al., 1993a).

— Galaxy cluster: sub-structure detection (Pierre and Starck, 1998; Krywult et al., 1999; Arnaud et al., 2000), hierarchical clustering (Pando et al., 1998a), distribution of superclusters of galaxies (Kalinkov et al., 1998).

— Cosmic Microwave Background: evidence for scale-scale correlations in the Cosmic Microwave Background radiation in COBE data (Pando et al., 1998b), large-scale CMB non-Gaussian statistics (Popa, 1998; Aghanim et al., 2001), massive CMB data set analysis (Gorski, 1998).

Cosmology: comparing simulated cosmological scenarios with observations (Lega et al., 1996), cosmic velocity field analysis (Rauzy et al., 1993).

This broad success of the wavelet transform is due to the fact that astronomical data generally gives rise to complex hierarchical structures, often described as fractals. Using multiscale approaches such as the wavelet transform, an image can be decomposed into components at different scales, and the wavelet transform is therefore well-adapted to the study of astronomical data.

This section reviews briefly some of the existing transforms. 1.2.1 Fourier Analysis

The Fast Fourier Transform. The Fourier transform of a continuous function f (t) is defined by:

and the inverse Fourier transform is:

The discrete Fourier transform is given by:

and the inverse discrete Fourier transform is:

In the case of images (two variables), this is:

Since f (u,v) is generally complex, this can be written using its real and imaginary parts:

MN 1

MN uk vl

It can also be written using its modulus and argument: f (u,v) = I f (u,v) I e4 arg

I f (u, v) |2 is called the power spectrum, and 0(u, v) = arg f (u, v) the phase.

Two other related transforms are the cosine and the sine transforms. The discrete cosine transform is defined by:

2N 2N

1.2.2 Time-Frequency Representation

The Wigner-Ville Transform. The Wigner-Ville distribution (Wigner, 1932; Ville, 1948) of a signal s(t) is

ir 2nv dr

where s* is the conjugate of s. The Wigner-Ville transform is always real (even for a complex signal). In practice, its use is limited by the existence of interference terms, even if they can be attenuated using specific averaging approaches. More details can be found in (Cohen, 1995; Mallat, 1998).

The Short-Term Fourier Transform. The Short-Term Fourier Transform of a 1D signal f is defined by:

If g is the Gaussian window, this corresponds to the Gabor transform. The energy density function, called the spectrogram, is given by:

Fig. 1.1 shows a quadratic chirp s(t) = sin(3N2), N being the number of pixels and t G {1,.., N}, and its spectrogram.

50 40

Fig. 1.1. Left: a quadratic chirp and right: its spectrogram. The y-axis in the spectrogram represents the frequency axis, and the x-axis the time. In this example, the instantaneous frequency of the signal increases with the time.

The inverse transform is obtained by:

Example: QPO Analysis. Fig. 1.2, top, shows an X-ray light curve from a galactic binary system, formed from two stars of which one has collapsed to a compact object, very probably a black hole of a few solar masses. Gas from the companion star is attracted to the black hole and forms an accretion disk around it. Turbulence occurs in this disk, which causes the gas to accrete

Wigner Ville Gif

Fig. 1.1. Left: a quadratic chirp and right: its spectrogram. The y-axis in the spectrogram represents the frequency axis, and the x-axis the time. In this example, the instantaneous frequency of the signal increases with the time.

1. Introduction to Applications and Methods 250

1. Introduction to Applications and Methods 250

63.9

57.5

25.6

1000 2000 Time Spectrogram

3000

1000 2000 Time Spectrogram

3000

2400

1200 1800 Time (sec)

Fig. 1.2. Top: QPO X-ray light curve, and bottom: its spectrogram.

3000

2400

1200 1800 Time (sec)

Fig. 1.2. Top: QPO X-ray light curve, and bottom: its spectrogram.

3000

slowly to the black hole. The X-rays we see come from the disk and its corona, heated by the energy released as the gas falls deeper into the potential well of the black hole. The data were obtained by RXTE, an X-ray satellite dedicated to the observation of this kind of source, and in particular their fast variability which gives us information on the processes in the disk. In particular they show sometimes a QPO (quasi-periodic oscillation) at a varying frequency of the order of 1 to 10 Hz (see Fig. 1.2, bottom), which probably corresponds to a standing feature rotating in the disk.

1.2.3 Time-Scale Representation: The Wavelet Transform

The Morlet-Grossmann definition (Grossmann et al., 1989) of the continuous wavelet transform for a 1-dimensional signal f (x) G L2 (R), the space of all square integrable functions, is:

where:

— W(a, b) is the wavelet coefficient of the function f (x)

— b is the position parameter

The inverse transform is obtained by:

1 f+™ f+™ 1 fx - b\ dadb f (x) = ^ / ^W a b)H — —T"

Reconstruction is only possible if C^ is defined (admissibility condition) which implies that ^(0) = 0, i.e. the mean of the wavelet function is 0.

Fig. 1.3. Mexican hat function.

Fig. 1.3. Mexican hat function.

Fig. 1.3 shows the Mexican hat wavelet function, which is defined by:

This is the second derivative of a Gaussian. Fig. 1.4 shows the continuous wavelet transform of a 1D signal computed with the Mexican Hat wavelet. This diagram is called a scalogram. The y-axis represents the scale.

Continuous Wavelet Transform Haar
Fig. 1.4. Continuous wavelet transform of a 1D signal computed with the Mexican Hat wavelet.

The Orthogonal Wavelet Transform. Many discrete wavelet transform algorithms have been developed (Mallat, 1998; Starck et al., 1998a). The most widely-known one is certainly the orthogonal transform, proposed by Mallat (1989) and its bi-orthogonal version (Daubechies, 1992). Using the orthogonal wavelet transform, a signal s can be decomposed as follows:

k k j = l with $j,i(x) = 2-j$(2-jx — l) and ^j,l(x) = 2-j^(2-jx — l), where $ and ^ are respectively the scaling function and the wavelet function. J is the number of resolutions used in the decomposition, Wj the wavelet (or detail) coefficients at scale j, and cJ is a coarse or smooth version of the original signal s. Thus, the algorithm outputs J +1 subband arrays. The indexing is such that, here, j = 1 corresponds to the finest scale (high frequencies). Coefficients Cj,k and Wjk are obtained by means of the filters h and g:

k where h and g verify:

k and the reconstruction of the signal is performed with:

Cjii = 2^[h(k + 2l)cj+ik + g(k + 2l)Wj+1,k ] (1.20)

k where the filters h and g must verify the conditions of dealiasing and exact reconstruction:

The two-dimensional algorithm is based on separate variables leading to prioritizing of horizontal, vertical and diagonal directions. The scaling function is defined by 4>(x, y) = 4>(x)4>(y), and the passage from one resolution to the next is achieved by:

cj + l(kK ,ky E E h(lx - 2kx)h(ly - 2 ky )fj (l x ,l y ) (1.22)

The detail signal is obtained from three wavelets:

which leads to three wavelet subimages at each resolution level. For three dimensional data, seven wavelet subcubes are created at each resolution level, corresponding to an analysis in seven directions. Other discrete wavelet transforms exist. The a trous wavelet transform which is very well-suited for astronomical data is discussed in the next chapter, and described in detail in Appendix A.

1.2.4 The Radon Transform

The Radon transform of an object f is the collection of line integrals indexed by (0, t) e [0, 2n) x R given by where 6 is the Dirac distribution. The two-dimensional Radon transform maps the spatial domain (x,y) to the Radon domain (0,t), and each point in the Radon domain corresponds to a line in the spatial domain. The transformed image is called a sinogram (Liang and Lauterbur, 2000).

A fundamental fact about the Radon transform is the projection-slice formula (Deans, 1983):

This says that the Radon transform can be obtained by applying the one-dimensional inverse Fourier transform to the two-dimensional Fourier transform restricted to radial lines going through the origin.

This of course suggests that approximate Radon transforms for digital data can be based on discrete fast Fourier transforms. This is a widely used approach, in the literature of medical imaging and synthetic aperture radar imaging, for which the key approximation errors and artifacts have been widely discussed. See (Toft, 1996; Averbuch et al., 2001) for more details on the different Radon transform and inverse transform algorithms. Fig. 1.5 shows an image containing two lines and its Radon transform. In astronomy, the Radon transform has been proposed for the reconstruction of images obtained with a rotating Slit Aperture Telescope (Touma, 2000), for the BATSE experiment of the Compton Gamma Ray Observatory (Zhang et al., 1993), and for robust detection of satellite tracks (Vandame, 2001). The Hough transform, which is closely related to the Radon transform, has been used by Ballester (1994) for automated arc line identification, by Llebaria (1999) for analyzing the temporal evolution of radial structures on the solar corona, and by Ragazzoni and Barbieri (1994) for the study of astronomical light curve time series.

1.2.5 The Ridgelet Transform

The two-dimensional continuous ridgelet transform in R2 can be defined as follows (Candes and Donoho, 1999). We pick a smooth univariate function ^ : R ^ R with sufficient decay and satisfying the admissibility condition

Rf (0,t) = f (x1,x2)S(x1 cos 0 + x2 sin 0 — t) dx1dx2,

Was this article helpful?

0 0
Learn Photoshop Now

Learn Photoshop Now

This first volume will guide you through the basics of Photoshop. Well start at the beginning and slowly be working our way through to the more advanced stuff but dont worry its all aimed at the total newbie.

Get My Free Ebook


Post a comment