pmelchior / scarlet2

Scarlet, all new and shiny
MIT License
14 stars 3 forks source link

Source has center #72

Closed pmelchior closed 2 months ago

pmelchior commented 2 months ago

See discussion #71. This PR moves the center attribution from Morphology to Component (which now knows what to produce and where to place it). This means that Morphology and Spectrum now only have shape attributes, not bbox attributes (because they have no placement information).

In addition, the PR brings a simplification for the init procedures. The common init functions now return ArraySpectrum and ArrayMorphology instead of the ndarrays. So, they avoid the second line:

spectrum = init.pixel_spectrum(obs, center)
spectrum = ArraySpectrum(spectrum)

If the old behavior is desired, run the init methods with return_array=True.

Finally, one can convert between different morphology implementation, e.g.

img = <some 2d array>
morph = ArrayMorphology(img)

gaussian_morph = GaussianMorphology.from_image(img)
starlet_morph = StarletMorphology.from_image(morph())
morph_ = ArrayMorphology(starlet_morph())

Together this PR simplifies a number of internal operations, but the important part is for the user interface. A typical call now looks like this:

with Scene(model_frame) as scene:
    try:
        spectrum, morph = init.adaptive_morphology(obs, center, min_corr=0.99)
    except ValueError:
        spectrum = init.pixel_spectrum(obs, center)
        morph = init.compact_morphology()
    Source(center, spectrum, morph)
charlotteaward commented 2 months ago

@pmelchior I'm running into a few issues arising from when I use lists of observations instead of a single one. For example, if I run init_pixel_spectrum on one observation, I get the required ArraySpectrum. But if I run spectrum = init.pixel_spectrum(observations, centerpix) it gives me a list of arrays [Array([103.41209], dtype=float32), Array([74.764824], dtype=float32)] etc that can't be given to Source.

Similarly, since we now get the spectrum and morphology for an extended source in a single call spectrum, morph = init.adaptive_morphology(observations, centerpix, min_corr=0.99), I'm running into issues with giving a list of observations (see error below). (My old workaround for the old version was to run the morphology initialization on a single image but do the spectrum initialization separately on a list).

AttributeError                            Traceback (most recent call last)
Cell In[26], line 98
     81 '''
     82 print('Making other galaxy')
     83 flux = jnp.asarray(np.asarray(initialization.pixel_spectrum(observations_sc2, center))[:,0][inds0])
   (...)
     95     scarlet2.ArrayMorphology(morph_init))
     96 '''
     97 try:
---> 98     spectrum, morph = init.adaptive_morphology(observations_sc2, centerpix, min_corr=0.99)
     99 except IndexError:
    100     spectrum = init.pixel_spectrum(observations_sc2, centerpix)

File ~/scarletcenterless/scarlet2/scarlet2/initialization.py:54, in adaptive_morphology(obs, center, min_size, delta_size, min_snr, min_corr, return_array)
     25 def adaptive_morphology(obs, center, min_size=11, delta_size=3, min_snr=20, min_corr=0.99, return_array=False):
     26     """
     27     Create image of a Gaussian from the centered 2nd moments of the observation.
     28 
   (...)
     52     jnp.ndarrays (for spectrum and morphology) or ArraySpectrum, ArrayMorphology
     53     """
---> 54     assert obs.weights is not None, "Observation weights are required"
     56     if isinstance(center, SkyCoord):
     57         center = obs.frame.get_pixel(center)

AttributeError: 'list' object has no attribute 'weights' 
pmelchior commented 2 months ago

None of these methods were tested for multi-observation inits. The main reason is that it's not clear at all how that should be done. For instance, your example spectrum = init.pixel_spectrum(observations, centerpix) can imply a number of things that are supposed to happen:

  1. observations contain different channels, and the resulting spectrum gets a concatenation in the order of observations. However, this may not be the order of channels in the model frame. In this mode, we should read the model frame channels and associate the observed channel with the model channel accordingly.
  2. observations contain (at least some) identical channels. Should we average the results? Maybe.

So, we could improve this by implementing a channel check as in item 1 and maybe a warning in case 2.

The situation gets more complicated for adaptive_morphology. The method finds the largest bounding box around center to contain light that has the same color as the center pixel. It then measures the moments of the image in that box to create a Gaussian with those moments. If the placement or resolution of the observations is different, we have to make adjustments to how this Gaussian moments are combined. Possible, and a good idea, but not part of this PR.

So, I suggest to check that no observation lists are being passed into the init function until we allow for all init functions to have consistent multi-obs behavior.

charlotteaward commented 2 months ago

That makes sense! I've tested everything and the initializations and sampling over center are working well. My workaround for initializing fluxes for multiple images is now to do something like:

spectrum = jnp.asarray([init.pixel_spectrum(obs, centerpix).data for obs in observations_sc2])[:,0]
spectrum = scarlet2.spectrum.ArraySpectrum(spectrum)
Source(centerpix, spectrum, morph)