Metadata-Version: 2.1
Name: spectroflat
Version: 0.5.0
Summary: Flat field and smile correction spectro-polarimetric solar images.
Home-page: https://gitlab.gwdg.de/hoelken/spectroflat
Author: Johannes Hölken
Author-email: hoelken@mps.mpg.de
Project-URL: Bug Tracker, https://gitlab.gwdg.de/hoelken/spectroflat/issues
Classifier: Programming Language :: Python :: 3
Classifier: License :: OSI Approved :: BSD License
Classifier: Operating System :: OS Independent
Requires-Python: >=3.8
Description-Content-Type: text/markdown
License-File: LICENSE

# spectroflat
Python based library to flat field spectro-polarimetric data.

## Theoretical Background 
Generally this is intended to become an extension of the
"[Precise reduction of solar spectra obtained with large CCD arrays](https://www.aanda.org/articles/aa/pdf/2002/42/aa2154.pdf)"
method by Wöhl et al.

## Pre-Flat generation
The pre-flat is a combination of the sensor and slit flat that 
contains most of the "hard" flat field features.
Most prominently the following is corrected:
- Sensor features (e.g. dust on the sensor itself)
- Slit features (e.g. dust on the slit resulting in line features in the spectral direction)
- (Some of) polarimetric fringes.
- Gradient, line-to-line variation in spatial direction and
   col-to-col variation in spectral direction

Please Note: The pre-flat is not a full flat. 
Especially it does not correct for spectral gradients, 
spectrographic curvature nor all fringes.

The algorithm assumes that the spatial direction is vertical (Y) 
and the spectral direction horizontal (X). The pre-flat is 
generated by using averaged flat field frames 
(ideally one per mod state).

### I Smoothed ideal
We fit a polynomial of configurable degree (usually 2 to 7) to 
each column to generate a smoothed ideal of the flat field frame.
We then remove the average spacial gradient from the smoothed ideal.

Finally, if the `gen_column_response_map` is configured, the vertical sensor response is removed from the 
smoothed ideal. For this we create an average column-wise response-map of the sensor and adjust the response 
by adding the offset to one of the two columns.   

### III Pre-flat extraction
To extract the sensor and slit parts of the flat field we divide 
the flat frame by the smoothed ideal. 

## Gain table generation algorithm
This algorithm includes the generation of an `Offset Map` that can (and should) be 
used to perform smile correction on science frames.  

### I. Rotation
From averaged target images (Grid, Random Dot, ...) that show strong horizontal features

-  Detect angle of camera wrt grating (horizontal features)

### II. Smile Detection
Rotate the prepared flat images[1] (one per mod. state) such that the horizontal 
effects really are horizontally.  

For each image do:
- Automatically select lines to fit based on configured window size
- For every line:
    - Every Rth row (configurable) find the line-core (Gaussian or Lorentzian fits)
    - Fit all points with a second order polynomial.
    - Compute horizontal Offsets wrt vertical straight line
- Filter bad lines / fits
    - Not enough points towards the upper or lower border
    - Not enough points at all
    - Strong deviation in leading coefficient (wrt to mean value of leading coefficients)
    - Fitted the wrong line in selection window
    - Poor fit (based on numpy polyfit error value)
- For every row: Fit the horizontal offsets with an Nth order polynomial (configurable, best results 2nd  and 3rd order) and generate an offset map

### III. Desmiling:
Rotate the prepared flat images [1] (one per mod. state) such that the horizontal 
effects really are horizontally.

For each image do:
- Apply offest map from II for each pixel:
    - Use an horizontal (row wise) FFT shift
    - select a discrete Gaussian around the shifted pixel (i. e. 60% central pixel, 30 % the directly adjacent ones and 10 % the next neighbors)
- **Redo II**
    -  If mean correction is above configured values (max distance wrt straight vertical line, value of the leading coefficient, ...):
        - Generate a delta offset-map with these corrections
        - Sum up both offset maps
        - Enforce min(offset_map) = 0
        - **Redo III**
    - If mean correction is not significant or number of max iterations reached:
       - **Jump to IV**

### IV. Gaintable generation:
Use smile corrected averaged images (one per mod. State) from III

- Detect, save and subtract global vertical gradient (i. e. average over all columns)
- Detect and mask outliers
- Compute average over all rows and divide by the average
- Re-apply outliers
- Multiply with global vertical gradient

### V. Smoothing (optional):
Use gain tables (one per mod. State) from **IV**

- Detect peaks in horizontal average
- For each row and each peak:
    - Interpolate the value from the left and right side of the peak
    - Replace the peak value by the interpolated on

*[1] Prepared flat images: Average of multiple images in flat field mode with dark and hot pixel corrections
applied.*


## Usage

**NOTE** This library expects the spacial domain on the vertical-axis and 
the spectral domain on the horizontal axis. 

### Detect Rotation
The spectral / vertical rotation will be corrected by the smile correction algorithm.

To detect the horizontal rotation (sensor vs spacial features) of the instrument
you need frames with strong spacial features (e.g. a random dot or pinhole grit target).
It's a good practice to repeat that with every flat field you take.

```python
from spectroflat import RotationAnalysis

detected_rotation = RotationAnalysis.detect_horizontal_rotation(frame_data)
```

If you have that data for multiple mod. states consider using the average
rotation.

### Analyze prepared flat field datacube 

The needed input is a 3D data cube with an averaged flat field 
image per modulation state. The expected shape is `(state, aspacial, spectral)`.
To extract `OffsetMap` and a Gain Table from such a cube, you 
need to adjust the `SmileConfig` and then just run the `Analyzer`

```python 
from spectroflat import SmileConfig, Analyser

# Adjust the default configuration to your needs. 
# You can also load the smile config from a YAML file.
config = SmileConfig(max_iterations=5, rotation=detected_rotation)
an = Analyser(cube, config).run()
```

You can access the gain table via `an.gain_table` as an 3D numpy array with the
same shape as the provided cube with a gain table per mod. state. 
Save it in the format that is suitable for you (e.g. FITS). 

You can access the offset map via `an.offset_map` and dump it to 
a Pickle file with `an.offset_map.dump('/path/to/save/offset_map.pickle')`.

### De-smile science frames
Once you have obtained the offset map you can load it from a pickle
file and correct each frame with

```python 
from spectroflat import OffsetMap, SmileCorrector

om = OffsetMap.from_file('/path/to/offset_map.pickle')
corrected = SmileCorrector(smap: om, img: frame_data, mod_state: 0).run().result
```

## Configuration 
The `SmileConfig` configuration provides a lot of possibilities to 
tweak the algorithm. You can find a description of all values in the 
[API Documentation](https://hoelken.pages.gwdg.de/spectroflat/doc/spectroflat/base/smile_config.html).
However, here are the most important once: 

```python
#: En-/disable post line removal smoothing of the gain-table
smooth: bool = True

#: If True the result will be an average offset map of all mod states
squash: bool = False

#: Image rotation in spacial dimension. If not set, this will be detected 
#: by the algorithm If set, the algorithm will not try to detect the rotation
#: and use this value.
rotation: float = None

#: Integer > 1 to define the minimum distance of two lines 
#: to take into account in pixels.
line_distance: int = 13

#: Degree of horizontal polynomial smile fit
poly_degree: int = 4

#: Max Iterations: Terminate regardless of reaching the goals below
max_iterations: int = 10

#: Maximum of the mean pixel deviation from a straight line after de-smiling.
#: Iteration is aborted when this value is reached with an STD < 08
mean_px_deviation: float = 0.1

#: If the decrease from one iteration to the other of the mean deviation
#: is lower than the given number the iteration is regarded complete.
px_deviation_decrease: float = 0.5
```

## Contact
This code is developed and maintained at the Max Planck Institute for 
Solar System Research (MPS) Göttingen.

### Maintainer 
- Johannes Hoelken ([hoelken@mps.mpg.de](mailto:hoelken@mps.mpg.de))

### Contributions
- Alex Feller 
- Francisco Iglesias

## License
BSD clause-2, see [LICENSE](LICENSE)
