SourceXtractorPlusPlus  0.11
Please provide a description of the project.
measurement_images.py
Go to the documentation of this file.
1 # -*- coding: utf-8 -*-
2 
3 # Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université
4 #
5 # This library is free software; you can redistribute it and/or modify it under
6 # the terms of the GNU Lesser General Public License as published by the Free
7 # Software Foundation; either version 3.0 of the License, or (at your option)
8 # any later version.
9 #
10 # This library is distributed in the hope that it will be useful, but WITHOUT
11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 # FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
13 # details.
14 #
15 # You should have received a copy of the GNU Lesser General Public License
16 # along with this library; if not, write to the Free Software Foundation, Inc.,
17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 from __future__ import division, print_function
19 
20 import os
21 import re
22 import sys
23 
24 from astropy.io import fits
25 
26 import _SourceXtractorPy as cpp
27 
28 if sys.version_info.major < 3:
29  from StringIO import StringIO
30 else:
31  from io import StringIO
32 
33 
34 measurement_images = {}
35 
36 
37 class MeasurementImage(cpp.MeasurementImage):
38  """
39  A MeasurementImage is the processing unit for SourceXtractor++. Measurements and model fitting can be done
40  over one, or many, of them. It models the image, plus its associated weight file, PSF, etc.
41 
42  Parameters
43  ----------
44  fits_file : str or astropy.fits.HDUList
45  The path to a FITS image, or an instance of astropy.fits.HDUList
46  psf_file : str
47  The path to a PSF. It can be either a FITS image, or a PSFEx model.
48  weight_file : str or astropy.fits.HDUList
49  The path to a FITS image with the pixel weights, or an instance of astropy.fits.HDUList
50  gain : float
51  Image gain. If None, `gain_keyword` will be used instead.
52  gain_keyword : str
53  Keyword for the header containing the gain.
54  saturation : float
55  Saturation value. If None, `saturation_keyword` will be used instead.
56  saturation_keyword : str
57  Keyword for the header containing the saturation value.
58  flux_scale : float
59  Flux scaling. Each pixel value will be multiplied by this. If None, `flux_scale_keyword` will be used
60  instead.
61  flux_scale_keyword : str
62  Keyword for the header containing the flux scaling.
63  weight_type : str
64  The type of the weight image. It must be one of:
65 
66  - none
67  The image itself is used to compute internally a constant variance (default)
68  - background
69  The image itself is used to compute internally a variance map
70  - rms
71  The weight image must contain a weight-map in units of absolute standard deviations
72  (in ADUs per pixel).
73  - variance
74  The weight image must contain a weight-map in units of relative variance.
75  - weight
76  The weight image must contain a weight-map in units of relative weights. The data are converted
77  to variance units.
78  weight_absolute : bool
79  If False, the weight map will be scaled according to an absolute variance map built from the image itself.
80  weight_scaling : float
81  Apply an scaling to the weight map.
82  weight_threshold : float
83  Pixels with weights beyond this value are treated just like pixels discarded by the masking process.
84  constant_background : float
85  If set a constant background of that value is assumed for the image instead of using automatic detection
86  image_hdu : int
87  For multi-extension FITS file specifies the HDU number for the image. Default 1 (primary HDU)
88  psf_hdu : int
89  For multi-extension FITS file specifies the HDU number for the psf. Defaults to the same value as image_hdu
90  weight_hdu : int
91  For multi-extension FITS file specifies the HDU number for the weight. Defaults to the same value as image_hdu
92  """
93 
94  def __init__(self, fits_file, psf_file=None, weight_file=None, gain=None,
95  gain_keyword='GAIN', saturation=None, saturation_keyword='SATURATE',
96  flux_scale=None, flux_scale_keyword='FLXSCALE',
97  weight_type='none', weight_absolute=False, weight_scaling=1.,
98  weight_threshold=None, constant_background=None,
99  image_hdu=1, psf_hdu=None, weight_hdu=None
100  ):
101  """
102  Constructor.
103  """
104  if isinstance(fits_file, fits.HDUList):
105  hdu_list = fits_file
106  file_path = fits_file.filename()
107  else:
108  hdu_list = fits.open(fits_file)
109  file_path = fits_file
110 
111  if isinstance(weight_file, fits.HDUList):
112  weight_file = weight_file.filename()
113 
114  super(MeasurementImage, self).__init__(os.path.abspath(file_path),
115  os.path.abspath(psf_file) if psf_file else '',
116  os.path.abspath(weight_file) if weight_file else '')
117 
118  if image_hdu <= 0 or (weight_hdu is not None and weight_hdu <= 0) or (psf_hdu is not None and psf_hdu <= 0):
119  raise ValueError('HDU indexes start at 1')
120 
121  self.meta = {
122  'IMAGE_FILENAME': self.file,
123  'PSF_FILENAME': self.psf_file,
124  'WEIGHT_FILENAME': self.weight_file
125  }
126 
127  self.meta.update(hdu_list[image_hdu-1].header)
128  self._load_header_file(file_path, image_hdu)
129 
130  if gain is not None:
131  self.gain = gain
132  elif gain_keyword in self.meta:
133  self.gain = float(self.meta[gain_keyword])
134  else:
135  self.gain = 0.
136 
137  if saturation is not None:
138  self.saturation = saturation
139  elif saturation_keyword in self.meta:
140  self.saturation = float(self.meta[saturation_keyword])
141  else:
142  self.saturation = 0.
143 
144  if flux_scale is not None:
145  self.flux_scale = flux_scale
146  elif flux_scale_keyword in self.meta:
147  self.flux_scale = float(self.meta[flux_scale_keyword])
148  else:
149  self.flux_scale = 1.
150 
151  self.weight_type = weight_type
152  self.weight_absolute = weight_absolute
153  self.weight_scaling = weight_scaling
154  if weight_threshold is None:
155  self.has_weight_threshold = False
156  else:
157  self.has_weight_threshold = True
158  self.weight_threshold = weight_threshold
159 
160  if constant_background is not None:
162  self.constant_background_value = constant_background
163  else:
164  self.is_background_constant = False
165  self.constant_background_value = -1
166 
167  self.image_hdu = image_hdu
168 
169  if psf_hdu is None:
170  self.psf_hdu = image_hdu
171  else:
172  self.psf_hdu = psf_hdu
173 
174  if weight_hdu is None:
175  self.weight_hdu = image_hdu
176  else:
177  self.weight_hdu = weight_hdu
178 
179  global measurement_images
180  measurement_images[self.id] = self
181 
182  # overrides the FITS headers using an ascii .head file if it can be found
183  def _load_header_file(self, filename, hdu):
184  pre, ext = os.path.splitext(filename)
185  header_file = pre + ".head"
186  current_hdu = 1
187 
188  if os.path.exists(header_file):
189  print("processing ascii header file: " + header_file)
190 
191  with open(header_file) as f:
192  for line in f:
193  line = re.sub("\\s*$", "", line)
194  if line == "":
195  continue
196 
197  if line.upper() == "END":
198  current_hdu += 1
199  elif current_hdu == hdu:
200  m = re.match("([^=]{1,8})=([^\\/]*)(.*)", line)
201  if m:
202  keyword = m.group(1).strip().upper()
203  value = m.group(2).strip()
204  self.meta[keyword] = value
205 
206  def __str__(self):
207  """
208  Returns
209  -------
210  str
211  Human readable representation for the object
212  """
213  return 'Image {}: {} / {}, PSF: {} / {}, Weight: {} / {}'.format(
214  self.id, self.meta['IMAGE_FILENAME'], self.image_hdu, self.meta['PSF_FILENAME'], self.psf_hdu,
215  self.meta['WEIGHT_FILENAME'], self.weight_hdu)
216 
217 
218 def print_measurement_images(file=sys.stderr):
219  """
220  Print a human-readable representation of the configured measurement images.
221 
222  Parameters
223  ----------
224  file : file object
225  Where to print the representation. Defaults to sys.stderr
226  """
227  print('Measurement images:', file=file)
228  for i in measurement_images:
229  im = measurement_images[i]
230  print('Image {}'.format(i), file=file)
231  print(' File: {}'.format(im.file), file=file)
232  print(' PSF: {}'.format(im.psf_file), file=file)
233  print(' Weight: {}'.format(im.weight_file), file=file)
234 
235 
236 class ImageGroup(object):
237  """
238  Models the grouping of images. Measurement can *not* be made directly on instances of this type.
239  The configuration must be "frozen" before creating a MeasurementGroup
240 
241  See Also
242  --------
243  MeasurementGroup
244  """
245 
246  def __init__(self, **kwargs):
247  """
248  Constructor. It is not recommended to be used directly. Use instead load_fits_image or load_fits_images.
249  """
250  self.__images = []
251  self.__subgroups = None
252  self.__subgroup_names = set()
253  if len(kwargs) != 1 or ('images' not in kwargs and 'subgroups' not in kwargs):
254  raise ValueError('ImageGroup only takes as parameter one of "images" or "subgroups"')
255  key = list(kwargs.keys())[0]
256  if key == 'images':
257  if isinstance(kwargs[key], list):
258  self.__images = kwargs[key]
259  else:
260  self.__images = [kwargs[key]]
261  if key == 'subgroups':
262  self.__subgroups = kwargs[key]
263  for name, _ in self.__subgroups:
264  self.__subgroup_names.add(name)
265 
266  def __len__(self):
267  """
268  See Also
269  --------
270  is_leaf
271 
272  Returns
273  -------
274  int
275  How may subgroups or images are there in this group
276  """
277  if self.__subgroups:
278  return len(self.__subgroups)
279  else:
280  return len(self.__images)
281 
282  def __iter__(self):
283  """
284  Allows to iterate on the contained subgroups or images
285 
286  See Also
287  --------
288  is_leaf
289 
290  Returns
291  -------
292  iterator
293  """
294  if self.__subgroups:
295  return self.__subgroups.__iter__()
296  else:
297  return self.__images.__iter__()
298 
299  def split(self, grouping_method):
300  """
301  Splits the group in various subgroups, applying a filter on the contained images. If the group has
302  already been split, applies the split to each subgroup.
303 
304  Parameters
305  ----------
306  grouping_method : callable
307  A callable that receives as a parameter the list of contained images, and returns
308  a list of tuples, with the grouping key value, and the list of grouped images belonging to the given key.
309 
310  See Also
311  --------
312  ByKeyword
313  ByPattern
314 
315  Raises
316  -------
317  ValueError
318  If some images have not been grouped by the callable.
319  """
320  if self.__subgroups:
321  #if we are already subgrouped, apply the split to the subgroups
322  for _, sub_group in self.__subgroups:
323  sub_group.split(grouping_method)
324  else:
325  subgrouped_images = grouping_method(self.__images)
326  if sum(len(p[1]) for p in subgrouped_images) != len(self.__images):
327  self.__subgroups = None
328  raise ValueError('Some images were not grouped')
329  self.__subgroups = []
330  for k, im_list in subgrouped_images:
331  assert k not in self.__subgroup_names
332  self.__subgroup_names.add(k)
333  self.__subgroups.append((k, ImageGroup(images=im_list)))
334  self.__images = []
335 
336  def add_images(self, images):
337  """
338  Add new images to the group.
339 
340  Parameters
341  ----------
342  images : list of, or a single, MeasurementImage
343 
344  Raises
345  ------
346  ValueError
347  If the group has been split, no new images can be added.
348  """
349  if self.__subgroups is not None:
350  raise ValueError('ImageGroup is already subgrouped')
351  if isinstance(images, MeasurementImage):
352  self.__images.append(images)
353  else:
354  self.__images.extend(images)
355 
356  def add_subgroup(self, name, group):
357  """
358  Add a subgroup to a group.
359 
360  Parameters
361  ----------
362  name : str
363  The new of the new group
364 
365  group : ImageGroup
366  """
367  if self.__subgroups is None:
368  raise Exception('ImageGroup is not subgrouped yet')
369  if name in self.__subgroup_names:
370  raise Exception('Subgroup {} alread exists'.format(name))
371  self.__subgroup_names.add(name)
372  self.__subgroups.append((name, group))
373 
374  def is_leaf(self):
375  """
376  Returns
377  -------
378  bool
379  True if the group is a leaf group
380  """
381  return self.__subgroups is None
382 
383  def __getitem__(self, name):
384  """
385  Get a subgroup.
386 
387  Parameters
388  ----------
389  name : str
390  The name of the subgroup.
391 
392  Returns
393  -------
394  ImageGroup
395  The matching group.
396 
397  Raises
398  ------
399  ValueError
400  If the group has not been split.
401  KeyError
402  If the group has not been found.
403  """
404  if self.__subgroups is None:
405  raise ValueError('ImageGroup is not subgrouped yet')
406  try:
407  return next(x for x in self.__subgroups if x[0] == name)[1]
408  except StopIteration:
409  raise KeyError('Group {} not found'.format(name))
410 
411  def print(self, prefix='', show_images=False, file=sys.stderr):
412  """
413  Print a human-readable representation of the group.
414 
415  Parameters
416  ----------
417  prefix : str
418  Print each line with this prefix. Used internally for indentation.
419  show_images : bool
420  Show the images belonging to a leaf group.
421  file : file object
422  Where to print the representation. Defaults to sys.stderr
423  """
424  if self.__subgroups is None:
425  print('{}Image List ({})'.format(prefix, len(self.__images)), file=file)
426  if show_images:
427  for im in self.__images:
428  print('{} {}'.format(prefix, im), file=file)
429  else:
430  print('{}Image sub-groups: {}'.format(prefix, ','.join(str(x) for x, _ in self.__subgroups)), file=file)
431  for name, group in self.__subgroups:
432  print('{} {}:'.format(prefix, name), file=file)
433  group.print(prefix + ' ', show_images, file)
434 
435  def __str__(self):
436  """
437  Returns
438  -------
439  str
440  A human-readable representation of the group
441  """
442  string = StringIO()
443  self.print(show_images=True, file=string)
444  return string.getvalue()
445 
446 
447 class ImageCacheEntry(object):
448  def __init__(self, image, kwargs):
449  self.image = image
450  self.kwargs = kwargs
451 
452  def match_kwargs(self, kwargs):
453  mismatch = []
454  for key, value in kwargs.items():
455  if key not in self.kwargs:
456  mismatch.append('{} {} != undefined'.format(key, value))
457  elif self.kwargs[key] != value:
458  mismatch.append('{} {} != {}'.format(key, value, self.kwargs[key]))
459  return mismatch
460 
461 _image_cache = {}
462 
463 def load_fits_image(image, psf=None, weight=None, **kwargs):
464  """Creates an image group with the images of a (possibly multi-HDU) single FITS file.
465 
466  If image is multi-hdu, psf and weight can either be multi hdu or lists of individual files.
467 
468  In any case, they are matched in order and HDUs not containing images (two dimensional arrays) are ignored.
469 
470  :param image: The FITS file containing the image(s)
471  :param psf: psf file or list of psf files
472  :param weight: FITS file for the weight image or a list of such files
473 
474  :return: A ImageGroup representing the images
475  """
476 
477  image_hdu_list = fits.open(image)
478  image_hdu_idx = [i for i, hdu in enumerate(image_hdu_list) if hdu.is_image and hdu.header['NAXIS'] == 2]
479 
480  # handles the PSFs
481  if isinstance(psf, list):
482  if len(psf) != len(image_hdu_idx):
483  raise ValueError("The number of psf files must match the number of images!")
484  psf_list = psf
485  psf_hdu_idx = [0] * len(psf_list)
486  else:
487  psf_list = [psf] * len(image_hdu_idx)
488  psf_hdu_idx = range(len(image_hdu_idx))
489 
490  # handles the weight maps
491  if isinstance(weight, list):
492  if len(weight) != len(image_hdu_idx):
493  raise ValueError("The number of weight files must match the number of images!")
494  weight_list = weight
495  weight_hdu_idx = [0] * len(weight_list)
496  elif weight is None:
497  weight_list = [None] * len(image_hdu_idx)
498  weight_hdu_idx = [0] * len(weight_list)
499  else:
500  weight_hdu_list = fits.open(weight)
501  weight_hdu_idx = [i for i, hdu in enumerate(weight_hdu_list) if hdu.is_image and hdu.header['NAXIS'] == 2]
502  weight_list = [weight_hdu_list] * len(image_hdu_idx)
503 
504  image_list = []
505  for hdu, psf_file, psf_hdu, weight_file, weight_hdu in zip(
506  image_hdu_idx, psf_list, psf_hdu_idx, weight_list, weight_hdu_idx):
507  image_list.append(MeasurementImage(image_hdu_list, psf_file, weight_file,
508  image_hdu=hdu+1, psf_hdu=psf_hdu+1, weight_hdu=weight_hdu+1, **kwargs))
509 
510  return ImageGroup(images=image_list)
511 
512 def load_fits_images(images, psfs=None, weights=None, **kwargs):
513  """Creates an image group for the given images.
514 
515  Parameters
516  ----------
517  images : list of str
518  A list of relative paths to the images FITS files. Can also be single string in which case,
519  this function acts like load_fits_image
520  psfs : list of str
521  A list of relative paths to the PSF FITS files (optional). It must match the length of image_list or be None.
522  weights : list of str
523  A list of relative paths to the weight files (optional). It must match the length of image_list or be None.
524 
525  Returns
526  -------
527  ImageGroup
528  A ImageGroup representing the images
529 
530  Raises
531  ------
532  ValueError
533  In case of mismatched list of files
534  """
535 
536  if isinstance(images, list):
537  if len(images) == 0:
538  raise ValueError("An empty list passed to load_fits_images")
539 
540  psfs = psfs or [None] * len(images)
541  weights = weights or [None] * len(images)
542 
543  if not isinstance(psfs, list) or len(psfs) != len(images):
544  raise ValueError("The number of image files and psf files must match!")
545 
546  if not isinstance(weights, list) or len(weights) != len(images):
547  raise ValueError("The number of image files and weight files must match!")
548 
549  groups = []
550  for f, p, w in zip(images, psfs, weights):
551  groups.append(load_fits_image(f, p, w, **kwargs))
552 
553  image_list = []
554  for g in groups:
555  image_list += g
556 
557  return ImageGroup(images=image_list)
558  else:
559  load_fits_image(images, psfs, weights, **kwargs)
560 
561 class ByKeyword(object):
562  """
563  Callable that can be used to split an ImageGroup by a keyword value (i.e. FILTER).
564 
565  Parameters
566  ----------
567  key : str
568  FITS header keyword (i.e. FILTER)
569 
570  See Also
571  --------
572  ImageGroup.split
573  """
574 
575  def __init__(self, key):
576  """
577  Constructor.
578  """
579  self.__key = key
580 
581  def __call__(self, images):
582  """
583  Parameters
584  ----------
585  images : list of MeasurementImage
586  List of images to group
587 
588  Returns
589  -------
590  list of tuples of str and list of MeasurementImage
591  i.e. [
592  (R, [frame_r_01.fits, frame_r_02.fits]),
593  (G, [frame_g_01.fits, frame_g_02.fits])
594  ]
595  """
596  result = {}
597  for im in images:
598  if self.__key not in im.meta:
599  raise KeyError('The image {}[{}] does not contain the key {}'.format(
600  im.meta['IMAGE_FILENAME'], im.image_hdu, self.__key
601  ))
602  if im.meta[self.__key] not in result:
603  result[im.meta[self.__key]] = []
604  result[im.meta[self.__key]].append(im)
605  return [(k, result[k]) for k in result]
606 
607 
608 class ByPattern(object):
609  """
610  Callable that can be used to split an ImageGroup by a keyword value (i.e. FILTER), applying a regular
611  expression and using the first matching group as key.
612 
613  Parameters
614  ----------
615  key : str
616  FITS header keyword
617  pattern : str
618  Regular expression. The first matching group will be used as grouping key.
619 
620  See Also
621  --------
622  ImageGroup.split
623  """
624 
625  def __init__(self, key, pattern):
626  """
627  Constructor.
628  """
629  self.__key = key
630  self.__pattern = pattern
631 
632  def __call__(self, images):
633  """
634  Parameters
635  ----------
636  images : list of MeasurementImage
637  List of images to group
638 
639  Returns
640  -------
641  list of tuples of str and list of MeasurementImage
642  """
643  result = {}
644  for im in images:
645  if self.__key not in im.meta:
646  raise KeyError('The image {}[{}] does not contain the key {}'.format(
647  im.meta['IMAGE_FILENAME'], im.image_hdu, self.__key
648  ))
649  group = re.match(self.__pattern, im.meta[self.__key]).group(1)
650  if group not in result:
651  result[group] = []
652  result[group].append(im)
653  return [(k, result[k]) for k in result]
654 
655 
656 class MeasurementGroup(object):
657  """
658  Once an instance of this class is created from an ImageGroup, its configuration is "frozen". i.e.
659  no new images can be added, or no new grouping applied.
660 
661  Parameters
662  ----------
663  image_group : ImageGroup
664  """
665 
666  _all_groups = list()
667 
668  def __init__(self, image_group, is_subgroup=False):
669  """
670  Constructor.
671  """
672  self.__images = None
673  self.__subgroups = None
674  if image_group.is_leaf():
675  self.__images = [im for im in image_group]
676  else:
677  self.__subgroups = [(n, MeasurementGroup(g, is_subgroup=True)) for n,g in image_group]
678  if not is_subgroup:
679  MeasurementGroup._all_groups.append(self)
680 
681  def __iter__(self):
682  """
683  Returns
684  -------
685  iterator
686  """
687  if self.__subgroups:
688  return self.__subgroups.__iter__()
689  else:
690  return self.__images.__iter__()
691 
692  def __getitem__(self, index):
693  """
694  The subgroup with the given name or image with the given index depending on whether this is a leaf group.
695 
696  Parameters
697  ----------
698  index : str or int
699  Subgroup name or image index
700 
701  Returns
702  -------
703  MeasurementGroup or MeasurementImage
704 
705  Raises
706  ------
707  KeyError
708  If we can't find what we want
709  """
710 
711  if self.__subgroups:
712  try:
713  return next(x for x in self.__subgroups if x[0] == index)[1]
714  except StopIteration:
715  raise KeyError('Group {} not found'.format(index))
716  else:
717  try:
718  return self.__images[index]
719  except:
720  raise KeyError('Image #{} not found'.format(index))
721 
722  def __len__(self):
723  """
724  Returns
725  -------
726  int
727  Number of subgroups, or images contained within the group
728  """
729  if self.__subgroups:
730  return len(self.__subgroups)
731  else:
732  return len(self.__images)
733 
734  def is_leaf(self):
735  """
736  Returns
737  -------
738  bool
739  True if the group is a leaf group
740  """
741  return self.__subgroups is None
742 
743  def print(self, prefix='', show_images=False, file=sys.stderr):
744  """
745  Print a human-readable representation of the group.
746 
747  Parameters
748  ----------
749  prefix : str
750  Print each line with this prefix. Used internally for indentation.
751  show_images : bool
752  Show the images belonging to a leaf group.
753  file : file object
754  Where to print the representation. Defaults to sys.stderr
755  """
756  if self.__images:
757  print('{}Image List ({})'.format(prefix, len(self.__images)), file=file)
758  if show_images:
759  for im in self.__images:
760  print('{} {}'.format(prefix, im), file=file)
761  if self.__subgroups:
762  print('{}Measurement sub-groups: {}'.format(prefix, ','.join(
763  x for x, _ in self.__subgroups)), file=file)
764  for name, group in self.__subgroups:
765  print('{} {}:'.format(prefix, name), file=file)
766  group.print(prefix + ' ', show_images, file=file)
767 
768  def __str__(self):
769  """
770  Returns
771  -------
772  str
773  A human-readable representation of the group
774  """
775  string = StringIO()
776  self.print(show_images=True, file=string)
777  return string.getvalue()
def __init__(self, fits_file, psf_file=None, weight_file=None, gain=None, gain_keyword='GAIN', saturation=None, saturation_keyword='SATURATE', flux_scale=None, flux_scale_keyword='FLXSCALE', weight_type='none', weight_absolute=False, weight_scaling=1., weight_threshold=None, constant_background=None, image_hdu=1, psf_hdu=None, weight_hdu=None)
def __init__(self, image_group, is_subgroup=False)
def load_fits_image(image, psf=None, weight=None, kwargs)
ELEMENTS_API auto join(Args &&... args) -> decltype(joinPath(std::forward< Args >(args)...))
def load_fits_images(images, psfs=None, weights=None, kwargs)
def print(self, prefix='', show_images=False, file=sys.stderr)
def print(self, prefix='', show_images=False, file=sys.stderr)