© ESO, 2018

1. Introduction

The Sun’s magnetic field is concentrated into strong-field elements, or flux tubes, in addition to a weaker, more turbulent background. The strong-field component is thought to be mainly responsible for solar activity. In the lower solar atmosphere, this strong-field flux manifests itself as dark features, such as sunspots and pores, as well as bright structures forming faculae (or plage) and the network. With modern magnetographs and spectropolarimeters both components can be reliably observed. However, before Kitt Peak National Observatory started recording magnetograms in 1974, typically only sunspot records are available. This has negative repercussions for a number of research fields. One of these is the solar influence on climate through its variable irradiance (Haigh 2007; Solanki et al. 2013; Ermolli et al. 2013).

Regular space-based measurements have shown that the total solar irradiance (TSI, the spectrally integrated solar radiative flux at the top of the Earth’s atmosphere at the mean distance of one astronomical unit) varies on different time scales, from minutes to decades (Fröhlich 2013; Kopp 2016). The variability at time scales greater than a day is attributed to changes in the surface distribution of the solar magnetic field. Models have been developed to reconstruct the solar irradiance by accounting for the competing contributions from dark and bright magnetic features on the solar surface. These models have been successful in reproducing the measured TSI (Krivova et al. 2003; Ermolli et al. 2011; Ball et al. 2012; Chapman et al. 2013; Yeo et al. 2014). However, full-disc images of the Sun, typically used by such models, are only available for the last few decades. Reconstructions of past irradiance changes usually rely on observations of sunspots (Lean et al. 1995; Solanki & Fligge 2000; Wang et al. 2005; Krivova et al. 2007, 2010; Dasi-Espuig et al. 2014, 2016), or concentrations of cosmogenic radionuclides, e.g. ¹⁴C and ¹⁰Be (Steinhilber et al. 2009; Delaygue & Bard 2011; Shapiro et al. 2011; Vieira et al. 2011). However, such data provide information about the bright regions only indirectly and on the basis of strong assumptions.

The time series of available bright magnetic feature data can be potentially extended with the help of daily full-disc Ca II K spectroheliograms (SHG, hereafter) that have been stored in photographic archives at various sites since the beginning of the 20th century. Indeed, owing to the relation between the Ca II K brightness and the magnetic field strength averaged over a pixel (e.g. Babcock & Babcock 1955; Skumanich et al. 1975; Schrijver et al. 1989; Harvey & White 1999; Loukitcheva et al. 2009; Kahil et al. 2017), the SHG can be used to describe the temporal and spatial evolution of the solar surface magnetic field, and specifically its bright component as seen in the chromosphere. Since SHG allow direct measurements of the surface coverage and, with a proper calibration, also the photometric properties of bright magnetic features, the historical photographic archives are an important resource for studies of solar magnetism over the last century and the activity and irradiance variations it produced. This has led to an increasing interest towards the recovery and preservation of these early observations, which has recently prompted the digitization of various SHG archives, such as those at the Arcetri (Ermolli et al. 2009a), Coimbra (Garcia et al. 2011), Kodaikanal (Priyal et al. 2013), Meudon (Mein & Ribes 1990), Mitaka (Hanaoka 2013), and Mt. Wilson (Bertello et al. 2010) observatories.

Figure 1 presents typical examples of the images derived from the digitisation of SHG (hereafter referred to as digitised images). These images show non-uniformity across the field-of-view of solar and non-solar origin, as well as (non-solar) spots and scratches on the photographic support.

Analysis of SHG requires accounting for these non-uniformities and for the non-linear relation between the blackening of the photographic support (either glass plate or film, hereafter referred to as plate) and the flux of incident solar radiation during the plate’s exposure. This relation is expressed by the characteristic curve (or Hurter and Driffield curve, CC hereafter), which is specific of each photographic observation (Dainty & Shaw 1974) as it depends on a variety of factors, for example the exposure time, the gelatine of the photographic plate, the composition of the developers (reducing and fixing baths), the duration of plate development and the temperature, and degree of stirring during this step. The CC has in general a sigmoid shape. A typical example of a CC is shown in Fig. 2 (see for example Sect. 3 or Dainty & Shaw 1974, for more details).

Until recent times, the CC was mainly derived using specific exposures (step-wedges, hereafter) acquired with known relative intensity ratios, in addition to the archived observation (e.g. Fredga 1971; Kariyappa & Sivaraman 1994; Kariyappa & Pap 1996; Worden et al. 1998a; Giorgi et al. 2005; Ermolli et al. 2009a). However, in most cases acquisition of these exposures started many years after the start of solar observations. For example, synoptic SHG were taken regularly since 1915 at the Mt. Wilson Observatory, but calibration exposures were imprinted on the plates only from 1961 onwards; at the Kodaikanal site observations started in late 1904 and step-wedges measurements in 1958. Therefore, the bulk of these data lacks specific information on the photometric calibration.

Nevertheless, there have been attempts to calibrate the data lacking wedges. Priyal et al. (2013) calibrated Kodaikanal data by applying an average CC derived from the available step-wedge exposures on all images. Ermolli et al. (2009b) applied the method presented by Mickaelian et al. (2007) for calibration of photographic plates of star surveys, which is based on information stored in the unexposed and darker parts of the plate. Tlatov et al. (2009) suggested to calibrate SHG by using a linear relation between a standardized centre-to-limb variation (CLV, hereafter) profile published by Pierce & Slaughter (1977, derived with a second degree polynomial fit and corresponding to 390.928 nm) and the values on the analysed image at two positions.

Various methods have also been used to process Ca II K images and compensate them for the CLV. Some examples are computations that result in radially symmetric backgrounds (Brandt & Steinegger 1998; Walton et al. 1998; Johannesson et al. 1998; Denker et al. 1999; Zharkova et al. 2003), application of a median filter (Lefebvre et al. 2005; Bertello et al. 2010; Chatterjee et al. 2016), 2D polynomial fittings to the entire image (Worden et al. 1998b; Caccin et al. 1998), or interpolation of the mode pixel values within radial and azimuthal disc sectors (Tlatov et al. 2009). However, all of these methods are unable to thoroughly account for the artefacts affecting historical data. Worden et al. (1998a) presented a method that is able to account for the inhomogeneities in the images by fitting a 1D fifth degree polynomial along image rows and columns to intensity values of the original observation and on the image resulting from its 45° rotation. Regions with values outside the ± 2σ interval of the image, where σ is the standard deviation, were excluded from the analysis and the background is estimated by applying a median, low-pass filter to the average of the fit. Other studies applied variations and combinations of methods previously described in the literature (e.g. Ermolli et al. 2009b; Singh et al. 2012; Priyal et al. 2013). For instance, Priyal et al. (2013) used the method of Denker et al. (1999) and a modified version of the approach suggested by Worden et al. (1998a), applying a third degree polynomial fit instead of the originally proposed fifth degree polynomial.

The long list of references above show that many methods have been employed for processing and calibrating historical Ca II K observations. However, their aptness and accuracy has generally been discussed only qualitatively.

In this paper we present a new method to calibrate the historical photographic Ca II K SHG as well as a novel method to compensate for various patterns on the data. These patterns include the CLV of the solar intensity and patterns of non-solar origin introduced by observational and archival processes. This method, which works without knowledge of the specific CC of the analysed plate and in the absence of step-wedge imprints, is based on the computation of the CLV of density values on quiet Sun (QS, hereafter) regions. Importantly, we rigorously test our technique against modern CCD-based data artificially degraded to correspond to photographically obtained SHG.

The historical and modern observations used in this study are described in Sect. 2. The method is presented in Sect. 3. We assess the accuracy of the method with synthetic data that have known characteristics and artefacts and test the proposed method on a large sample of SHG in Sect. 4. We compare our method with other approaches presented in the literature in Sect. 5 and draw our conclusions in Sect. 6.

2. Observational data

The historical observations, on which our technique is tested, are the digitised SHG from the Arcetri (Ar), Coimbra (Co), Kodaikanal (Ko), Meudon (Me), Mitaka (Mi), McMath-Hulbert (MM), and Mt. Wilson (MW) observatories. These images were taken in the Ca II K line at λ = 393.367 nm, with nominal bandwidths ranging from 0.01 nm to 0.05 nm and scanning time of several minutes needed to cover the solar disc. The data, which were taken between February 1907 and May 2002 with different instrumental set-ups and stored mainly on photographic plates, were digitized by using various devices and methods. This results in solar images with different sizes and characteristics depending on the archive.

The modern observations used to test our technique were taken at the INAF Osservatorio Astronomico di Roma (OAR) with the Precision Solar Photometric Telescope (Rome/PSPT, Ermolli et al. 1998, 2007). These data were acquired from August 2000 to October 2014, by using an interference filter centred at the Ca II K line with a 0.25 nm bandwidth and the exposure time much shorter than for the historical observations, usually about 0.06 s. We analysed images after the usual reduction steps (dark current removal and flat fielding) that greatly reduce instrumental effects and after they were resized to half their original dimension (2 arcsec per pixel after resizing), to roughly match the size of the historical data.

Key information about the datasets used in this study is summarised in Table 1.

The available data include estimates of the solar radius R and centre coordinates in the images. To avoid our results being affected by any uncertainties in the radius estimates, we only analysed the disc pixels within 0.99 R (0.98 R for Ko). This corresponds to μ = cosθ = 0.14 (μ = 0.2 for Ko), where θ is the heliocentric angle.

3. Method description

The digitized SHG reflect the blackening degree of the emulsion, which is proportional to the transparency T of the plate and is related to the density d of the photographic emulsion since: $\mathit{d}\mathrm{=}\log \mathrm{(}\mathrm{1}\mathit{/}{\mathit{T}}^{\mathrm{)}}\mathit{.}$(1)The modern data represent the flux of incident radiation during the observation. The CC is defined as the relation d = f(log E) (shown in Fig. 2), where E = I × t is the plate exposure, I the incident energy per unit area and t the exposure time. Assuming that the plate was exposed evenly, the CC becomes d = f(log I). Therefore the photometric calibration of the data means applying on them the specific and appropriate relation between d and I.

Our method of calibration and processing of the SHG is based on the assumption that, in these observations, the CLV of intensity in the internetwork regions, that is the quietest part of the QS, does not vary with time, in agreement with White & Livingston (1978, 1981), Livingston & Wallace (2003), and Livingston et al. (2007). This assumption is also supported by the results of Bühler et al. (2013) and Lites et al. (2014) that the internetwork magnetic flux, which is the main component of magnetic field populating these quietest parts of the solar surface, remains unchanged over the solar cycle. Before applying the proposed method to the historical images, we convert their values to densities according to Eq. (1). Figure 3a shows the density image corresponding to the Ko observation displayed in Fig. 1b.

The main steps of our method can be summarised as:

calculation of the 2D map of QS regions in each image (includingCLV and inhomogeneities);

extraction of the 1D QS CLV profile in each image;

construction of CC by relating the 1D QS CLV to a reference 1D radial intensity QS CLV from CCD-based Rome/PSPT observations;

calibration of the image using the CC;

compensation for the intensity CLV.

Visual inspection of historical data shows QS density patterns (2D QS background, hereafter) that are in general strongly non-symmetric and inhomogeneous (see, e.g., Fig. 1). This is due to a plethora of problems affecting the data, introduced either during the observation (for example partial coverage by clouds, vignetting, uneven movement of the slit), the development (for example non homogeneous bathing, touching the plate before the process was finished), the storage period (for example dust accumulation, scratches, humidity, ageing burns), or the digitization (for example dust, hair, loss of dynamic range) of the plate. Since results from modern observations indicate the radial uniformity of radiative emission of the internetwork (1D radial QS, Livingston & Sheeley 2008; Peck & Rast 2015), we assume that all departures of the 2D internetwork map in SHG from a radially symmetric map are image artefacts of non-solar origin.

It is worth mentioning that the bandwidths listed in Table 1 are nominal, while the real values can vary even within a single SHG. The varying bandwidths mean that also the heights of the atmosphere that are sampled change. With broader bandwidths, the photospheric contribution is higher and plage regions appear less bright, while sunspots are more visible. As an example, Fig. 4 shows a MW observation where the right part of the image was taken with a broader bandwidth, either intentionally or due to instrument problems. If the bandwidth of the observations was varied or if they were centred at a slightly different wavelength this would have resulted in a different limb darkening profile. However, as information on the real bandwidth is not available, in our analysis we have to assume that the bandwidth of the observations has remained constant over the disc and for all data within an archive.

3.1. Background computation

We first derive the 2D QS background. For this, it is essential to accurately exclude active regions (AR, hereafter), otherwise there is a risk of overestimating the background due to the contribution of remaining AR. The process of calculating the 2D QS background can be outlined as follows:

get a rough estimate of the background;

use the last estimate of the background to identify and exclude AR;

iteratively repeat the process of calculating the background and identifying the AR until sufficient accuracy is achieved.

The first rough estimate of the background is obtained in 4 steps.

Step 1.1. The solar disc is divided into azimuthal slices of 30° that cover the disc in steps of 5°. Within each slice we apply a fifth degree polynomial fit of the form d = f(μ). The best fit values of d are assigned to all pixels within a given slice. Results obtained from contiguous slices are gradually averaged and stitched together (this method will be referred to as rotating slices).

Step 1.2. We compute the density contrast image ${\mathit{C}}_{\mathit{i}}^{\mathit{d}}\mathrm{=}\mathrm{(}{\mathit{d}}_{\mathit{i}}\mathrm{-}{\mathit{d}}_{\mathit{i}}^{\mathrm{QS}}\mathrm{)}\mathit{/}{\mathit{d}}_{\mathit{i}}^{\mathrm{QS}}$, where d_i is the original density image and ${\mathit{d}}_{\mathit{i}}^{\mathrm{QS}}$ the density background resulting from Step 1.1 at the ith pixel.

Step 1.3. We identify the regions in the density image lying outside 1σ of the result of Step 1.2. These are tentatively identified as AR and are excluded from the further analysis (we note that the removal of non-AR pixels in this process does not influence the results; it is more important to discard as many AR pixels as possible).

Step 1.4. We apply a fifth degree polynomial fit to the density image, excluding AR defined in Step 1.3, along each column and row of the image separately (similarly to Worden et al. 1998a, but without the 45° rotation of the disc). To all pixels of each analysed column/row we assign the density values resulting from the best fit. We get a background map by stitching together the results obtained from the best fit; multiple values derived for the same location on the solar disc are averaged (we will refer to this method as column/row fittings).

The calculations described in Steps 1.1–1.4 provide a first, rough estimate of the background. However, the identification of AR at this stage is rather inaccurate and the obtained values of the background around AR are overestimated. The varying contrast for different μ positions renders the identification near the limb less accurate.

Therefore in order to improve the AR identification and the calculation of the background, we iteratively repeat the following steps:

Step 2.1. We compute the density contrast image ${\mathit{C}}_{\mathit{i}}^{\mathit{d}}\mathrm{=}\mathrm{(}{\mathit{d}}_{\mathit{i}}\mathrm{-}{\mathit{d}}_{\mathit{i}}^{\mathrm{QS}}\mathrm{)}\mathit{/}{\mathit{d}}_{\mathit{i}}^{\mathrm{QS}}$, where d_i is the original density image and ${\mathit{d}}_{\mathit{i}}^{\mathrm{QS}}$ the density background resulting from the previous calculation at the ith pixel. During the first iteration d^QS is the rough background estimate (from Step 1.4), afterwards we use the result of the previous iteration.

Step 2.2. We remove the AR in the original density image, which are identified as the pixels that fulfil the threshold criteria ${\mathit{C}}_{\mathit{i}}^{\mathit{d}}\mathit{>}\mathrm{⟨}{\mathit{C}}^{\mathit{d}}\mathrm{⟩}\mathrm{+}{\mathit{l}}_{\mathrm{1}}\mathit{\sigma }$ or ${\mathit{C}}_{\mathit{i}}^{\mathit{d}}\mathit{<}\mathrm{⟨}{\mathit{C}}^{\mathit{d}}\mathrm{⟩}\mathrm{-}{\mathit{l}}_{\mathrm{2}}\mathit{\sigma }$, where l₁, l₂ are the applied thresholds. For the first 3 iterations l₁ = 1 and l₂ = 4, afterwards l₁ = 0.5 and l₂ = 1. These threshold values were chosen following a series of tests which showed that they effectively removed plage, network and sunspots, leaving just internetwork regions. See discussion in Sect. 4.1.2.

Step 2.3. To fill in the gaps left by the AR over the whole disc in the original density image we apply the column/row fittings on the image obtained in Step 2.2. To avoid artefacts of the fit at the limb due to missing values, the gaps by the AR in the outer 0.1 R are first filled in with the rotating slices method. The column/row fittings method is repeated two more times to the residual of the density image to the calculated background in order to improve the accuracy of the computation (see discussion in Sect. 4.1.2). All three computations are summed together.

Step 2.4. We apply a running window median filter on the image resulting from Step 2.3. The filter window width used is R/ 6 (shown with a rectangle in Fig. 3f; see discussion in Sect. 4.1.2). To avoid inconsistencies, the outer R/ 12 part of the disc is re-sampled outside the disc to fill the space between R and R + R/ 12 and the pixel values of the re-sampled section are adjusted so that the median filter is not affected by the pixel values outside the disc.

Step 2.5. We identify dark and bright linear artefacts affecting many of the SHG (for instance see Fig. 1d) by separately fitting a 1st degree polynomial to every row and column of the residual image between the original density image and the result of Step 2.4, excluding the AR. The fit is repeated 3 times as in Step 2.3.

The sum of the maps derived in Steps 2.4 and 2.5 is the final background map of each iteration. The five-step processing described above is repeated until the difference between two subsequent background maps does not improve the accuracy of the QS background further. Usually 4 iterations allow lowering the relative unsigned mean difference between maps from two consecutive iterations to <0.1%.

The result of the processing described so far is usually an asymmetric map (non-symmetric background, NSB hereafter) that describes the 2D QS background of the original image. The asymmetry is caused by a residual pattern that includes small- and large- scale inhomogeneities due to image problems, for example dark and bright bands and linear artefacts, stray light features, image gradients. Figure 3 shows different steps of the processing on a randomly selected Ko SHG. In particular, the original density image is shown in panel a), the result of the first column/row fittings process in panel b), while panels c)–h) show the results of each step (2.1–2.5) of the iterative process.

Finally, we apply the method of Nesme-Ribes et al. (1996) on the residual image between the original and the NSB to identify and compensate any offset in the average level of the computed NSB and the analysed image. The method by Nesme-Ribes et al. (1996) identifies the level of the QS as the minimum of the average densities of disc regions with density values within ± kσ, where k takes values between 0.1 and 3.0.

3.2. CLV extraction

We compute the radially symmetric background (SB hereafter), which in turn gives the 1D QS CLV, by applying a fifth degree polynomial fit to the density values d = f(μ) of the deduced NSB. All image pixels are considered and the fit is weighted with the local σ defined within 100 concentric and equal area annuli. This is the most important step for the calibration process, because the way the 1D QS CLV is defined, it determines the CC.

Our calculation of NSB includes possible vignetting affecting the original SHG, however this leads to miscalculation of the 1D QS CLV. Under-exposure depends on the intensity level, while vignetting is a purely geometric effect. Unfortunately, the geometry is such that their effects are similar, and vignetting is known to darken solar observations towards the limb, resembling the effects of under-exposure. The effects of vignetting are expected to be restricted mostly at the edges of the image, in which case it might not affect the solar disc at all or it is almost indistinguishable from under-exposure. Notwithstanding the similarities of their effects on the solar observation, some differences in the signatures of vignetting and under-exposure on the QS CLV and CC exist. In fact, vignetting can potentially affect larger parts of the solar observation than under-exposure, by modifying the CLV even near the image centre in contrast to under-exposure that mainly affects observations near the solar limb (see Fig. A.10). This is valid as long as the under-exposure is not extreme. However, in extreme cases the image quality is anyway very poor and such images cannot be processed meaningfully.

We attempt to estimate the vignetting on the analysed image and account for it in the calculation of the 1D QS CLV, by comparing the SB with a rescaled version of the logarithm of the CLV retrieved from the CCD-based Rome/PSPT observations. If there is no vignetting, the SB from the polynomial fit is kept. Otherwise, the SB is computed by rescaling the CLV from the modern observations. Its lower values are adjusted so as to find the minimum density for which the difference between SB and the CLV derived with the fit continuously increases towards the limb.

3.3. Photometric calibration

In our approach, we perform the photometric calibration by seeking for information stored on the solar disc that can be used in a similar manner as the often missing calibration wedges.

For each analysed image, we deduce the CC by relating the density values obtained from its SB at a given μ position to the logarithm of intensity values derived from modern Rome/PSPT observations at the same disc location. The amount of equal area annuli we use to acquire this relation and apply a linear fit to it is equal to the nearest integer of 2 R. This is consistent with the assumption that the QS density values in good-quality observations lie on the linear part of the CC. From the fit we only exclude the last value that corresponds to the regions very close to the limb, due to the higher uncertainties that characterise these regions. We linearly extrapolate the computed relation to the whole range of density values on the image and make use of the fit parameters to photometrically calibrate the original density image. The result of this procedure is illustrated in Fig. 5a. We also photometrically calibrate the estimated NSB. Removal of the QS CLV from the calibrated image and normalization to it then provides the contrast image of the analysed SHG (Fig. 5b).

Figure 6 shows the CC derived from the processing of the Ko observation displayed in Fig. 3a.

4. Performance of the method

4.1. Synthetic data

We tested the method proposed in Sect. 3 on a large number of synthetic images that were obtained from contrast Rome/PSPT images by imposing on them known radially symmetric CLV and a variety of inhomogeneities identified in historical observations (and converted to intensity). We then converted the degraded images from intensities to densities by applying various CC, thus emulating historical observations. We used the 2000 Rome/PSPT images described in Sect. 2, which were processed as described in Sect. 3.1 to derive contrast images. These contrast images are defined as ${\mathit{C}}_{\mathit{i}}^{\mathit{I}}\mathrm{=}\mathrm{(}{\mathit{I}}_{\mathit{i}}\mathrm{-}{\mathit{I}}_{\mathit{i}}^{\mathrm{QS}}\mathrm{)}\mathit{/}{\mathit{I}}_{\mathit{i}}^{\mathrm{QS}}$, where I_i is the intensity, and ${\mathit{I}}_{\mathit{i}}^{\mathrm{QS}}$ the intensity of the QS of the ith pixel. The imposed CLV are in the form of a fifth degree polynomial function of lnμ, as presented by Pierce & Slaughter (1977), with parameters determined from the Rome/PSPT observations. The range of values of the imposed intensity CLV is [0.6, 1.0], while the contrast range of the facular pattern is [0.1, ~0.6], and of the network pattern is [0.026, 0.1]. These ranges match the ones found on images of the Ko, Ar and Mi archives.

We produced 8 subsets of synthetic images with the following features:

Subset 1: Ideal density images with symmetric backgrounds, no linear artefacts or bands, to estimate the precision of the proposed method and the sensitivity of results to the level of solar activity.

Subset 2: Images with imposed CC that are linear functions with varying slopes, to test the sensitivity of the method to different slopes of CC within the range 0.1–4.0.

Subset 3: Images with imposed CC that are non-linear functions with various levels of over- and under-exposure (shown in Fig. A.5), to estimate errors in the retrieved calibration due to exposure problems.

Subset 4: Images with different solar disc sizes, to test the accuracy of the method on the various datasets. The disc diameter was varied between 200 and 1100 pixels.

Subset 5: A vignetting function with different strength levels was added to images, to test the ability of the method to identify possible vignetting effect and to distinguish it from other artefacts. The vignetting used here is axi-symmetric.

Subset 6: Various large scale density patterns of non-solar origin were imposed (shown in Fig. A.11), to evaluate the efficiency of the method to distinguish between the solar and non-solar patterns and the accuracy of accounting for the latter.

Subset 7: Images with a different CLV (shown in Fig. A.13) than the one used during the standardization of the CC, to test the errors of calibrating data obtained with different bandwidths.

Subset 8: A combination of all the artefacts mentioned above was added randomly on each image, to produce an inhomogeneous dataset resembling the historical ones. The CC is defined as a third degree polynomial function with randomly selected parameters within the following ranges: [−4.6,−1.3] for the constant term, [0.5, 4.0] for the linear term, [−0.03, 0.00] for the quadratic term, and [−0.03, 0.00] for the cubic term. The range of intensities within the undegraded images varies and the logarithm of it lies between −1.2 and 0.5 (similar to those shown with the light grey area in Fig. A.5).

Figures 7 and 8 show some examples of the synthetic images derived from the same contrast Rome/PSPT observation (taken on 21/08/2000), illustrating the variety of problems we aim at addressing with the application of the proposed method on these data. In the same figures we show the results obtained for these synthetic data and the pixel-by-pixel errors in the NSB calculation and the calibrated contrast images. The contrast images used to derive the errors are offset so that the mean QS contrast value is 1.

4.1.1. Overall performance of our method on synthetic data

Table 2 summarises key aspects of the various synthetic datasets, and quantitative results obtained by testing the proposed method on them. The results are briefly presented in the following and are described in greater detail in Appendix A. Table 2 summarizes the results derived from the analysis of all data within the various subsets. An exception is made for subsets 3 and 6, for which the presented values correspond to results restricted to images representative of historical data of reasonable quality, that is excluding synthetic images made to suffer from extreme exposure problems affecting the QS regions or images with superposed large-scale patterns whose amplitudes were larger than 0.6 of the CLV.

The results obtained for subsets 1–6 show that our method recovers the QS density CLV with average relative error <3% in the absence of strong non-solar patterns affecting the image, and error <6.5% if the extreme cases of inhomogeneities encompassed in subset 6 are also included. The results derived from subset 7 prove that the above accuracy is maintained as long as the CLV differs by roughly <10% from the one we imposed on the data.

The results obtained from the subset 8 demonstrate that the proposed method retrieves the NSB affecting the observation with a maximum relative error <2% averaged over all the analysed images. We found that, similarly to artefacts on the historical images, the presence of gradients on the synthetic images introduces errors in the calculation of SB, which affects the range of values in the calibrated data.

Figure 9 illustrates the accuracy of the proposed method on the most challenging subset 8. Shown are the maximum and average values of the relative difference as well as the rms difference between the calibrated and processed contrast image and the original undegraded synthetic contrast image. Throughout this analysis the relative differences are given in absolute values (unsigned), but we also provide the signed average difference. These values quantify the final errors of our image processing, which is comprised of the photometric calibration and the removal of the CLV and other image patterns from the analysed image and finally provides the corrected contrast image. We found that the maximum differences are on average <6.5%, while the average differences are <1%. There is, however, a tail consisting of images with higher maximum or average differences. For images representative of low activity periods, maximum errors are approximately 2% lower than those found for images at high activity periods, which illustrates the need for a very careful removal of active areas prior to carrying out any image processing. We discern no significant variation of the mean and rms differences over the solar cycle.

We also tested the accuracy of our proposed method for studies of variations in the fractional coverage of the solar disc by AR and network. Figure 10 shows the relative difference between the disc fraction of various solar features identified on the processed and unprocessed images of subset 8. The disc features were identified by applying a set of constant thresholds on the original (that is prior to degradation) and final (that is after recovery) contrast images. The thresholds used here were defined on the Rome/PSPT data to describe the plage and network regions. The thresholds in the contrast values used are th_p = 0.21 and th_n = 0.03 for the plage and network, respectively. Figure 10 demonstrates that the estimated disc fractions for the processed data typically lie within 2% of the values derived from the original data.

4.1.2. Performance on individual steps of the processing

Our method includes original ideas for processing SHG (for example rotating slices), but it also partly uses ideas from the previously published methods described in Sect. 1. By testing all these methods on synthetic data we identified their drawbacks, which helped us to optimize the steps of the method proposed here.

For example, in our calculation of the background, we do not rotate the image by 45° as proposed by Worden et al. (1998a). We found that the rotation does not improve the accuracy of the image processing further, if the outcome of other processing steps has been optimized. In contrast, our iterative fitting improves the accuracy of the AR identification and, in turn, of the QS estimation as compared to both non-iterative computations with a fifth degree polynomial function as suggested by Worden et al. (1998a) and iterative computations with higher degree functions. However, we also noticed that on average more than three computations of the fit per iteration step merely results in an increase of the noise of the final NSB map derived from the processing without improving the accuracy of the result.

In Step 2.2, to identify and exclude AR we apply a thresholding scheme with asymmetric limits, by using the values $\begin{array}{c}\mathrm{+}\mathrm{0.5}\\ -1\end{array}\mathit{\sigma }$ instead of the more widely employed ± 2σ. The asymmetric range allows us to account for both, the potentially inaccurate identification of AR near the solar limb at earlier processing steps and the small disc fraction of dark features in Ca II K observations. For example, using the symmetric limits of ± 1σ on the synthetic image in Fig. 7 S1/1) results in increased errors in the NSB of 3.5% and overestimation of AR intensity values (the errors with our method and the asymmetric limits are 1.8% and are shown in Fig. 7 S1/5). The employed upper limit of +0.5 was found to be a good compromise. Larger values led to the inclusion of significant portions of AR in the QS background, while lower values for the lower limit risk to wrongly exclude large regions of the solar observations near the limb from the 2D QS background calculation.

The window width for the median filter used in Step 2.4 was chosen to scale with R, to achieve consistent results from different data with varying disc size. The adopted width is larger than the typical scale of the network on the analysed observations, in order to avoid effects of small-scale density patterns of solar origin on image processing results, but is small enough to account for rapid changes of the background near the solar limb. We found that window widths in the range R/ 6−R/ 8 perform best on all available data and we adopted the more conservative value of R/ 6. This finding is in contrast to that of Chatterjee et al. (2016) who used a window width of ~R/ 13. Figure 11 shows one example of testing different window widths on an image from subset 1 (shown in Fig. 7 S1/1)). We show widths of R/ 2, R/ 4, R/ 13, and R/ 20 pixels. Large window widths fail very close to the limb, while smaller widths progressively fail on AR and network. Window width of R/ 4 gives comparable results with those of our adopted value, but slightly underestimates the values very close to the limb.

The accuracy of our processing also decreases if we do not identify the density variations of adjacent image lines² and simply derive a smooth background map. Furthermore, if these variations are not properly accounted for in the analysed image, its subsequent analysis aimed at the estimation of for example the photometric properties of AR returns inaccurate results.

Tests on subset 3 that includes non-linear CC, showed that our method is very accurate in recovering the shape of the CC even on observations with strong exposure problems, with relative errors in the computed CC being usually <0.5% under typical conditions as well as for the other synthetic subsets. Even in the extreme cases of over- or under-exposure the relative error in the computed CC lies below 1.7%.

The analysis of the images of the subset 2 suggests that the CC slope only mildly affects the results. Our method also accurately disentangles the vignetting contribution on the CC (see Fig. A.10). Figure 12 shows an example of processing an image from subset 5 without attempting to recover the vignetting (see Fig. 8 S5/4) and S5/6) for the results with the vignetting recovery). The vignetting contribution increases the slope of the CC and reduces the contrast values of the plage regions, with the errors over these regions being 35 times greater than if we account for the vignetting with our method.

4.2. Examples of calibrated SHG

We have applied the proposed method to many SHG randomly selected from the seven available historical (photographic) archives. Figures 13 and 14 show examples of the results obtained from observations taken at periods of high and low solar activity, respectively. From left to right, each panel shows the original observation (density image), the density image saturated to show the background, the background (CLV plus inhomogeneities) deduced from the proposed processing, the calibrated image, the identified inhomogeneities of the analysed image, and the final contrast image after the QS CLV removal. From top to bottom each such set of images is shown for SHG extracted from the Ar, Co, Ko, Me, Mi, MM, and MW archives, respectively. All the calibrated images are shown within the intensity range [0.0, 2.0], while the contrast images (that is images compensated for the CLV) are plotted within [−0.5, 0.5]. For the rest there is a colour bar denoting the range of values. Further examples can be found in Chatzistergos et al. (2016).

The deduced backgrounds describe the different patterns in the analysed images quite accurately. Figures 13 and 14 clearly show that the method works consistently with data extracted from various photographic archives, taken at different activity levels. All calibrated images lie within the same range of values and show a similar CLV pattern. The same is true for the contrast images that return plage regions within the same intensity ranges and no obvious residual large scale artefacts. The method is able to account even for strong inhomogeneities and rather peculiar patterns (for example Fig. 13 Ar, or Fig. 14 Ko) without affecting the plage regions. The inhomogeneities identified here show a CLV that is usually off-centred, in many cases having its highest value towards the limb. Furthermore, images show many dark and bright bands that could occur due to something occluding the Sun for a short period, or not constant exposure over the different rasters.

The tests on historical observations with good quality, resulted in 1D QS CLV very similar to the one from the Rome/PSPT data. This strengthens our argument that the CC of the good quality historical data can be described by a linear relation.

5. Comparison with other methods

5.1. Background calculation methods

We found that the methods presented in the literature that apply radially symmetric computations suffer from their inability to account for the asymmetric patterns affecting the images, while application of median filtering suffers under the presence of AR and its inability to account for density variations along adjacent image lines.

By applying the methods of Caccin et al. (1998) and Worden et al. (1998b) on the degraded synthetic SHG we found that these techniques provide inaccurate QS CLV values towards the limb and do not account for small-scale image artefacts, for example the density variations in adjacent lines. The method by Tlatov et al. (2009) is also unable to account for these linear artefacts.

Figure 15 shows the pixel-by-pixel relative differences between the NSB and the imposed background derived with our method and that by Worden et al. (1998a), on two synthetic images from subset 6. Part of the AR remained undetected by the latter method and so enters the computation of the background. It should be noted that when we applied the method by Worden et al. (1998a) we did not perform the last step of the low-pass filtering, because the information of the window-width or the way they applied it on regions very close to the limb is not described. Applying this step could potentially reduce the errors of isolated pixels, but it would not make a difference in the misidentification of the AR. Thus the method by Worden et al. (1998a) overestimates the actual background in some plage areas and introduces processing errors in recovering the NSB. The maximum relative errors in NSB are lower for observations taken at low solar activity (3.7%), than at periods of high solar activity (26%), but on average the method by Worden et al. (1998a) introduces one order of magnitude higher errors over the disc than obtained from our proposed method.

We applied the method by Priyal et al. (2013) on all data of subset 6 and we found that the method by Priyal et al. (2013) works reasonably well on images with weak anisotropies, but consistently fails to account for the large inhomogeneities affecting the data, by introducing up to 16 (25) times larger maximum (rms) errors than those from our method. Besides, the method by Priyal et al. (2013) applied the fit only to rows of the images and hence does not allow recovering any image patterns that occur in a direction different from the one considered for the fit.

5.2. Calibration methods

We tested the accuracy of the method by Priyal et al. (2013) by applying an average CC to calibrate a whole dataset with our synthetic data. We derived the average CC from the curves we imposed to all the data of subset 8, and studied the differences between contrast images obtained from the original data and the ones resulting from the calibration with the average CC. We used the imposed background of each image to compensate for the limb darkening, in order to avoid any other uncertainties of our procedure. This error estimate can be considered only as a lower limit, since the CC used to derive the average were the imposed ones, therefore without taking into consideration any errors in the calculation of the individual CC. Figure 16 shows that the errors introduced by the calibration with the average CC as proposed by the method by Priyal et al. (2013) are on average ~50%. These errors reach values as high as 300% for a few cases. We stress, however, that we cannot rule out that in actual historical data sets the CC display a smaller variation than in our subset 8. Therefore, the above test mainly shows the greater versatility of our presented technique for handling a range of CC values.

Figure 17 shows two examples of Ko images processed with our method (left panels) and by Priyal et al. (2013)⁴. The input data employed for this comparison come from different digitizations of the same observation. This limits our analysis of the results to qualitative aspects only. Our images were saturated at the same level to illustrate all AR clearly, however the data by Priyal et al. (2013) were provided in JPG file format and hence we are unable to saturate the images to the same level. Still, Fig. 17 clearly shows that images processed by Priyal et al. (2013) are affected by uncorrected inhomogeneities, to a significantly larger extent than images processed by our method.

In the image processing suggested by Tlatov et al. (2009), all image values are termed intensities and are further scaled linearly to the values of the standardized profile, without conversion to density values as expected by photographic theory, therefore this method merely applies a linear scaling to the image to let it match a desired range of intensity values. Also the standardised CLV that is taken from Pierce & Slaughter (1977) does not accurately represent the Ca II K data since it corresponds to 390.928 nm. Figure 18 shows the relative difference between results derived from the application of the method by Tlatov et al. (2009) and our method to one synthetic image from subset 1 produced from the contrast Rome/PSPT observation shown in Fig. 7 S1. The image calibrated with the method of Tlatov et al. (2009) displays a significant offset ~0.7 and fainter plage and dark regions than in the original image (contrast of ~0.1 and ~0.5 obtained by the method of Tlatov et al. 2009, and ours, respectively).

The method by Ermolli et al. (2009b) cannot be tested on the synthetic data, since it relies on information of the unexposed regions of the plate that we cannot replicate in a meaningful way in the synthetic data.

Finally, we compared the CC computed with our technique and from the calibration wedges stored on Ar data.

Since our method sets the QS near the solar disc centre to be around 1, whereas the wedges describe the response of the whole plate and contain no information as to which range the QS corresponds, there is a scaling factor between the images calibrated with our method and with the wedges. This factor depends mostly on the digitization (that is the range of values of the QS in the digital files), but also on other factors (for example slightly different exposure time should change the location of the QS in the CC). Thus, a direct comparison of the CC derived from the two methods is not straightforward.

To account for the difference in values range, we derived the CC from the wedges by applying a polynomial fit (Ermolli et al. 2009a), used this CC to calibrate the CLV calculated with our method and then rescaled it to match the range of values in the calibrated CLV derived with our method. Since information is lost with the rescaling, this approach does not allow any conclusion on the slope of the CC. Nonetheless, in this way we can test the assumption of using a linear curve to calibrate these data, provided the fitting of the wedge measurements is done accurately enough. However, this may not happen for all Ar data, because of insufficient or inaccurate information stored on the wedges. The wedges of the Ar data usually consist of 3 scans for 7 known exposures, giving 7 points in intensities to fit the sigmoid CC. These values do not necessarily cover the whole range of values on the disc, or even if they do their number may be too low to describe a sigmoid function.

Figure 19 shows an example of a derived CC for one Ar observation including the rescaled wedge measurements and fit. This is one of the good cases where the CC derived with our method matches almost perfectly the one derived from the wedges, missing only a small part of underexposed regions. It is important to note that the scatter of the values derived from the wedges is almost the same as the scatter in the background of the image. We achieve pixel-by-pixel relative differences between the image calibrated with our method and that with the wedge that are <5%, which is consistent with the results presented with the synthetic data from subset 6 (the example tested here was found to have inhomogeneities with a range within 0.6 that of the CLV).

6. Conclusions

We have developed a new method to photometrically calibrate the historical full-disc Ca II K SHG and to correct them for various artefacts of solar and non-solar origin. The method is based on the standardization of the QS CLV intensity pattern to the one resulting from modern observations, under the assumption that it does not vary with time. Modern observations suggest that this holds within the accuracy of the proposed method. We showed that the errors introduced by the above assumptions are relatively small and have minor impact on the CLV estimation unless the analysed observation is of very poor quality. We assume that QS regions store all the information required to construct the CC for the range of brightnesses covered by the QS anywhere on the solar disc. This is not fulfilled for observations with strong over-exposure effects, as these introduce errors into the bright plage regions near the centre of the solar disc, that cannot be calibrated by our method. Therefore very poor quality observations must be rejected. However, they constitute a very small fraction of the available historical data. In addition, it can be that only part of the characteristic curve is accurately represented by the QS density values and the assumed linear relation. This would affect mostly plage regions near disc centre or fainter regions towards the limb and can limit the accuracy of our processing. However, we showed that in spite of these potential shortcomings, our method results in much lower errors than different approaches presented in the literature.

To test the accuracy of the proposed method, we created a large number of synthetic images emulating various problems encountered in historical observations. The maximum error of our method is <6.5% averaged over all the degradations studied here, while the average error is <1%. These errors were derived on synthetic data including extreme cases of imposed artefacts. The maximum errors reduce to <2% if we exclude images with the most extreme artefacts. Application of other methods for the processing of SHG presented in the literature, returns errors that are between 3 and 300 times larger than those derived from our method.

We also estimated the accuracy of processing modern Ca II K data by applying the proposed method to synthetic images unaffected by linear artefacts. The error estimates decreased by almost a factor of 2 with respect to those reported earlier, with maximum relative errors being on average <0.6%.

Finally, we applied for test purposes the proposed method to a sample of images randomly extracted from seven historical SHG archives. We showed that the method allows us to process images from different archives consistently, without having to adjust or tailor the method to observations taken with different instruments at various observatories and at times of different levels of solar activity. The results of the application of the proposed method to various SHG archives will be presented in a forthcoming paper.

It is worth noting that our method to derive the CLV can be applied with minor adjustments to full-disc solar observations taken in other spectral ranges than Ca II K. Examples are archival white-light photographic images used for identifying and measuring sunspot properties (e.g. Ravindra et al. 2013; Willis et al. 2013; Hanaoka 2013), or SHG in the Hα line (e.g. Mein & Ribes 1990; Pötzi 2008; Garcia et al. 2011; Hanaoka 2013).

---

Acknowledgments

We thank Dr. Tlatov for providing data processed with his method and for providing feedback on his method of calibrating SHG. The authors thank the Arcetri, Coimbra, Kodaikanal, McMath-Hulbert, Meudon, Mitaka, Mt. Wilson, and the Rome Solar Groups. T.C. thanks the INAF Osservatorio astronomico di Roma solar physics group for their hospitality during his lengthy visits. T.C. acknowledges postgraduate fellowship of the International Max Planck Research School on Physical Processes in the Solar System and Beyond. This work was supported by grants PRIN-INAF-2014 and PRIN/MIUR 2012P2HRCR “Il Sole attivo”, COST Action ES1005 “TOSCA”, FP7 SOLID, and by the BK21 plus program through the National Research Foundation (NRF) funded by the Ministry of Education of Korea.

References

Appendix A: Analysis of the synthetic data

In this section we will introduce each subset of synthetic data individually, providing information on their characteristics and how they were created. We will then show all the results for the errors we get when calculating the NSB and the CLV corrected calibrated images by applying our technique. The results obtained are also summarised in Table 2.

Subset 1: Linear relation without inhomogeneities

We selected 500 Rome/PSPT images that sample the period 2000–2014 regularly, we computed the logarithm of each image and applied a linear CC to them to emulate the response of the plate. This set thus represents almost ideal density images without any large scale inhomogeneities or exposure problems, allowing to estimate the precision of the proposed method in the presence of AR.

Figure A.1 shows an example of results obtained on one image from the subset 1, where a) is the density image derived from the Rome/PSPT observation, b) is the imposed CLV, c) is the calculated NSB and d) is the calculated SB.

Figure A.2 shows results derived by comparing the NSB retrieved from the proposed method and the imposed background for all the data of subset 1 (disc positions within 0.98 R). Throughout this analysis the relative differences are given in absolute values (unsigned), but we also provide the signed average difference. The red curve displays the maximum relative difference divided by ten to improve the clarity of the plot. The mean values of the maximum and average relative difference between the retrieved NSB (SB) and the CLV intensity imposed on all images of subset 1 are 0.013 (0.002) and 0.002 (0.0007), respectively; the rms difference is 0.004 (0.001). On average only 0.5% of the pixels show differences exceeding the average value ± 3σ for the NSB. All above values show that the proposed method is able to retrieve radially symmetric backgrounds with great accuracy, even though this assumption is not made during the processing.

Figure A.3 shows results obtained by comparing the contrast image after the calibration performed with the proposed method and the original data; the mean values of the maximum unsigned and average relative difference are 0.024 and 0.004, respectively, and the rms difference is 0.005.

The calculated slopes for the CC derived from the method are in good agreement with the imposed one. The results we get have a slope of on average 0.504 with standard deviation 0.002, compared to the imposed one of 0.5.

To estimate the accuracy of the method when applied on data unaffected by intensity variations between adjacent lines, we analysed subset 1 one more time but now bypassing step 2.5 of the background calculation (Sect. 3.1). We found that the method’s accuracy is on average a factor of 2 higher than obtained by applying the line fittings. In particular, the mean values of the maximum unsigned and average relative difference between the retrieved NSB (contrast) and the CLV intensity (background) imposed on all images of subset 1 are lower than 0.02 (0.013) and 0.0002 (0.0002), respectively; the rms difference is 0.0027 (0.002).

Subset 2: Different linear relations without inhomogeneities

This subset was created to test the efficiency of the proposed method to account for different linear relations of the CC. Due to the changes in the photographic plates or films that were used, the different conditions of observations, development or digitization, the historical data have different CC. Therefore this subset tests the accuracy of our method to restore good quality images with different CC, but with no large-scale inhomogeneities or exposure problems. We selected 10 Rome/PSPT images covering the period 2000–2014 and imposed on them 20 different linear relations with slopes in the range 0.1–4.0. Subset 2 therefore consists of 200 images.

Figure A.4 shows the relative differences obtained by comparing images retrieved using our method and original undegraded data. Panels a) and b) refer to values of NSB, while panels c) and d) refer to values of contrast images. Each box represents a different image of subset 2, the rows (columns) show values derived from synthetic images created with different Rome/PSPT contrast observations (different imposed CC relations). The values on the y-axis represent the year when the used contrast Rome/PSPT observations were taken. The different colours give the errors according to the colour bar next to the plot. The relative differences between the retrieved NSB (SB) and imposed background increase with increasing slope of the relation applied, remaining however relatively small. The maximum unsigned relative difference measured is 0.03 (0.01) for the extreme case, but on average remains <0.0009 and maximum of rms difference is 0.009 (0.004).

The relative difference between the contrast image derived from our method and the original data reaches up to 0.07 for individual pixels in the worst case (on average <0.0009); the rms difference is less than 0.008.

Subset 3: Non-linear relations without inhomogeneities

This subset was created to test the effects of applying a linear calibration on data that were constructed by imposing non-linear CC relations. We used the same 10 images as employed for subset 2, computed the logarithm of each image and applied a series of non-linear relations on them. This way we can evaluate the accuracy of the proposed method for data suffering over- and under-exposure and to quantify the expected deviations introduced when applying a linear calibration on them. Figure A.5 shows the imposed relations. We used 2 separate functions to describe under-exposure and over-exposure. For each exposure problem we used 10 different strength levels and considered all combinations of these functions (in total 100 cases per image). The first 5 cases of over- and under-exposure (labelled 1–5 in Fig. A.5) produce images with intensity values belonging to the QS lying on the linear part of the CC, while for all the rest of the cases only a small part of the QS lies on the linear part of the CC. Subset 3 consists of 1000 images.

Figure A.6 shows the relative differences between the images derived with our method and the original undegraded data of subset 3. Panels a) and b) refer to values of NSB, while panels c) and d) refer to values of contrast images. Each 10 × 10 box is for synthetic images derived from different Rome/PSPT contrast images. In each panel, rows (columns) show values derived from different levels of over-exposure (under-exposure).

Overall the accuracy introduced by the method decreases with increasing exposure problems. Nonetheless the errors remain relatively small (maximum unsigned relative difference is 0.02 for the extreme case and 0.006 maximum unsigned rms difference) and are always comparable with the results from Subset 1. The inaccuracy increases in the cases when over-exposure affects the QS. The differences in SB remain extremely low for no under-exposure problems affecting the QS, maximum unsigned relative difference is 0.007, while maximum unsigned rms difference is 0.005.

The maximum unsigned relative difference obtained by comparing the retrieved and original undegraded contrast images is usually >0.1 and reaches up to 0.95 for the extreme exposure problems (case 10 in Fig. A.5); rms is <0.62. These errors affect mostly the bright features. The maximum unsigned differences remain <0.2 for the maximum over-exposure and under-exposure considered which does not affect the QS (cases 1–5 in Fig. A.5) and decrease with decreasing exposure problems. The errors are usually reasonable for the cases when there is no over-exposure on the QS.

Figure A.7 shows the CC constructed for the extreme exposure problems considered. The green, blue, and red curves display the imposed relation, the measured QS CLV, and the standardized linear relation, respectively. For all the data, our method allows calculating the NSB accurately enough so that the obtained CLV reproduces the imposed relation for almost the whole range of QS values. In particular, the maximum unsigned relative difference between the imposed CC and the one recovered by the method is less than 0.005. The accuracy drops slightly for the last point at the limb and only for the cases of severe under-exposure, rising to 0.017.

Subset 4: Different disc sizes with linear CC without inhomogeneities in the image

The photographic plates from different observatories have different sizes. Even among the images of the same archive there are variations due to changes of the spectroheliograph and of the photographic plates. The digitization affects the number of pixels covered by the disc as well. The disc size is also affected by the seasonal variation of the Sun–Earth distance. Therefore, this subset was created to test the efficiency of the proposed method when applied to datasets that have a different disc size. We selected 10 Rome/PSPT images covering the period 2000–2014 and re-sampled them, so that each was represented in the form of solar images with 10 different radii between 100 and 550 pixels. The range of radii was defined to include the sampling present in most available datasets, however the upper limit was dictated by the number of pixels covered by the disc on Rome/PSPT observations. Subset 4 consists of 100 images.

Figure A.8 shows the relative differences between the images derived with our method and the original Rome/PSPT data. Panels a) and b) refer to values of NSB, while panels c) and d) refer to values of contrast images, each (coloured) box for a different image of subset 4. In each panel, rows (columns) show values derived from different observations (imposed radii). The relative difference between the retrieved NSB (SB) and imposed background decreases with increasing radius of the disc; the maximum unsigned difference measured is 0.03 (0.007) for the images with the smallest radius, average difference is <0.0008, and maximum of rms difference is 0.004 (0.001).

The relative difference between the contrast image derived with our method and the original data reaches up to 0.056; the average differences are <0.007; the rms differences are <0.009.

Subset 5: Vignetting with linear CC without inhomogeneities

This subset was created to evaluate the effects of vignetting on the performance of the proposed method. We took the same 10 Rome/PSPT observations as employed above and imposed on them the same radially symmetric CLV as in subset 1, along with the vignetting function c × cos⁴s (Gardner 1947; Goldman & Chen 2005) in intensity, where s is the distance of each pixel to the centre of the disc and c is a constant that takes 10 discrete values between 0 and 0.5. This vignetting function causes maximum decrease in normalised intensities between 0 near the disc centre and 0.3 (which corresponds to 50%) near the limb in the density images. Subset 5 consists of 100 images.

Figure A.9 shows the relative differences between the images derived with our method and the original data. Panels a) and b) refer to values of NSB, while panels c) and d) refer to values of contrast images, each box corresponds to a different image of subset 5. In each panel, rows (columns) show values derived from different observations (imposed vignetting). The relative difference between the retrieved NSB (SB) and imposed background increases with increasing strength of vignetting (value of c); the maximum unsigned difference measured is 0.018 (0.014) for the extreme case, the average differences are <0.0004, and maximum of rms difference is 0.005 (0.016).

The relative difference between the contrast image derived with our method and the original data reaches up to 0.065; the average differences are <0.0003; the rms differences are <0.008.

Figure A.10 shows the effect of the vignetting on the CC. As intended, vignetting decreases non-linearly the density towards the limb. This results in significantly wrong slopes for the computed CC if the vignetting is not accounted for. Calibrating the data with the wrong slope of the CC results in lower contrasts than the ones from the undegraded data and the contrast tends to decrease with increasing magnitude of vignetting.

Subset 6: Varying magnitude of inhomogeneities with linear CC

This subset was created to evaluate the effects of the varying levels of inhomogeneities on the results derived with the proposed method. The magnitude of the inhomogeneities is compared to the range of values of the CLV and we construct for each background 10 cases with inhomogeneities whose amplitude lies in the range [0.1, 1.0] of the CLV. The backgrounds used here are shown in Fig. A.11. The facular patterns overlain on these backgrounds were derived from 10 Rome/PSPT observations in the period 2000–2014. Subset 6 consists of 2000 images, resulting from 20 backgrounds, 10 amplitudes for each background, each applied to the 10 chosen Rome/PSPT images.

The relative differences of the NSB derived with our method and the original data of subset 6 can be seen in Fig. A.12. We show the results for only two Rome/PSPT images, one taken at high and one at low solar activity periods, in order to make it easier to discern the differences between the results obtained from the various inhomogeneities and levels. In particular, in Fig. A.12 the x axis represents the different backgrounds, the circles show results for the high activity image, while the squares are for the low activity image and the colours signify the varying levels of inhomogeneities (black implying low amplitudes, while red is for large amplitudes).

Unsurprisingly, the errors introduced when retrieving the intensity image using the proposed method increase with magnitude of the inhomogeneities on the analysed image. In many images the errors remain relatively low, for example background No. 6 in Fig. A.11, where the maximum difference is 0.007. The maximum unsigned difference is obtained for the images with very bright/dark small artefacts, while larger-scale inhomogeneities and small artefacts with mild brightness are reliably processed by our method. The errors are highest when density variations in adjacent image lines are significantly brighter than their surroundings, with a maximum relative difference of 0.5. The maximum differences for NSB remain <0.1 for most of the backgrounds even in the extreme cases that were considered.

The results for the rest of the data in subset 6 are very similar to the ones presented here. We did not find any significant dependence of the results on the level of solar activity.

Subset 7: Different CLV with linear CC

Subset 7 was created to study the effects of applying our method on data that have a different CLV than the one we use to standardize them. The data could have a different CLV because the observation was centred at a different wavelength than the core of the Ca II K line, or possibly due to a different bandwidth (for example due to a different spectral resolution or slit width). We used again 10 Rome/PSPT images to derive the facular pattern and on each one of those imposed 10 different functions for the quiet Sun (seen in Fig. A.13), describing deviations from the one we measure in modern Rome/PSPT observations. Subset 7 consists of 100 images.

Figure A.14 shows the relative differences between the corrected and calibrated images derived with our method and the original data. Panels a) and b) refer to values of NSB, while panels c) and d) refer to values of contrast images, each box for a different image of subset 7. In each panel, rows (columns) show values derived from different observations (imposed CLV). The relative difference between the retrieved NSB (SB) and imposed background has a maximum value of 0.019 (0.007), the average differences are <0.0004, and the maximum rms difference is 0.004 (0.003).

We found that the different CLV affect more the calculation of the slope of the CC. The errors for the contrast images increase with the rms differences getting up to 0.038, while the maximum differences remain <10% for 3 CLV surrounding the average curve we use.

Subset 8: Random problems

Subset 8 was created to study the performance of our method on data where every image suffers from a different random set of problems. We used 2000 Rome/PSPT images to derive the facular pattern and on each one of those imposed a random function for the quiet Sun (shown in Fig. A.13), describing deviations from the one we measure in modern Rome/PSPT observations. We also added a random pattern of inhomogeneities that were used for subset 6 (Fig. A.11), with a random level of strength for each image. We convert the images to density by applying a CC in the form of a third degree polynomial with randomly selected parameters within the following ranges: [−4.6,−1.3] for the constant term, [0.5, 4.0] for the linear term, [−0.03, 0.00] for the quadratic term, and [−0.03, 0.00] for the cubic term. A random level of vignetting was also added. Subset 8 consists of 2000 images.

The results for this subset are presented in the main text in Sect. 4.1.1, here in Fig. A.15 we also show the plot with the relative errors for the NSB calculation.

All Tables

All Figures