Stacking and mosaic creation with JSol’Ex 2.0

Here are the main steps of the algorithm I have implemented:

The stitching part works quite differently than with typical stacking. In stacking, we have complete data for each image: we "only" have to align them. With mosaics, there are "missing" parts in the image that we need to fill in. To do this, we have to identify which part of a panel can be blended into the reconstructed image in order to complete it. This means that the alignment process is significanly more complicated than with typical stacking, since we will work on "missing" data. Part of the difficulty is precisely identifying if something is missing or not, that is to say if the signal of a pixel in one of the panels is relevant to the composition of the final image. This is done by comparing it with the estimated background level, but that’s not the only trick.

Despite the fact that our panels are supposedly aligned and that the circles representing the solar disks are supposed to be the same, in practice, depending on the quality of the capture and the ellipse regression success, the disks may be slightly off, with deformations. There can even be slight rotations between panels (because of flexions at capture time, or processing artifacts). As a consequence, a naive approach consisting of trying to minimize the error between 2 panels by moving them a few pixels in each direction like in stacking doesn’t work:

Giving all the issues I described above, I chose to implement an algorithm which would work similarly to stacking, by "warping" a panel into another. This process is iterative, and the idea is to take a "reference" panel, which is the one which has the most "relevant" pixels, and align tiles from the 2d panel into this reference panel.

To do this, we compute a grid of "reference points" which are in the "overlapping" area. These points belong to the reference image, and one difficulty is to filter out points which belong to "incomplete" tiles. Once we have these points, for each of them, we compute an alignment between the reference tile and the tile of the panel we’re trying to integrate. This gives us, roughly, a "model" of how tiles are displaced in the overlapping area. The larger the overlapping area is, the better the model will be, but experience shows that distorsion on one edge of the solar disk can be significanly different at the other edge.

The next step consists of trying to align tiles of the panel we integrate to the reference panel using this model. This is where the iteration process happens. In a nutshell, we have an area where the solar disk is "truncated". Even if we split the image in tiles like with stacking, we cannot really tell whether a tile is "complete" or not, because it depends both on the pixel intensities of the reference panel and the second panel, and the background level. In particular, calcium images may have dark areas within the solar disk which are sometimes as dark as the background.

If you are struggling to understand how difficult it can be to determine if part of the image we consider is relevant or not, let’s illustrate with this image:

Can you see what’s wrong in this image? Let’s increase constrast to make it clearly visible:

Now it should be pretty obvious that below the south edge of the truncated disk, we have an area which has pixels which are above the value of the background, but do not constitute actual signal! This problem took me quite some time to solve, and it’s only recently that I figured out a solution: before mosaicing, I am performing a background neutralization step, by modeling the background and substracting it from the image. While this doesn’t fully solve the problem, it makes it much less relevant for composition.

In addition, we have to compose the image using tiles which are incomplete, and we don’t know the orientation of the panels: they can be assembled north/south, or west/east, and nothing tells us. Potentially, it can even be a combination of these for a large number of panels.

Therefore, the algorithm works by creating a "mask" of the image being reconstructed. This mask tells us "for this particular pixel, I have reconstructed a value, or the value of the reference image is good enough and we won’t touch it". Then, for each tile, we consider tiles for which the mask is incomplete.

In order to determine how to align the truncated disk with data from the other image, we compute an estimate of the distortion of the tile based on the displacements models we have determined earlier. Basically, for a new "tile" to be integrated, we will consider the sample "reference points" which are within a reasonable distance of the tile. For this set of reference points, we know that they are "close enough" to compute an average model of the distorsion, that I call the "local distorsion": we can estimate, based on the distance of each reference point, how much they contribute to the final distorsion model for that particular point.

The key is really to consider enough samples to have a good distorsion model, but not too many because then the "locality" of alignment would become too problematic and we’d face misalignments. Because there are not so many samples for each "incomplete" tile, we are in fact going to reconstruct, naturally, the image from the edges where there’s missing data: when there are no samples, it basically means we cannot compute a model, so we don’t know how to align tiles. If we have enough samples, then we can compute a reliable model of the distorsion, and then we can reconstruct the missing part of each tile, by properly aligning the tiles together. If the number of samples is not sufficient to consider a good model, then we assume that no distorsion happens, which is often the case for "background" tiles.

Most of the difficulty in this algorithm is properly identifying "when" we can stitch tiles together, that is to say when we can tell that the alignment between tiles makes sense and that the alignment is correct. Often, I got good results for one kind of images (e.g, h-alpha images) but horrible results with others (e.g calcium) or the other way around. I cannot really say I took a very scientific approach to this problem, but more an empirical approach, tweaking parameters of my algorithm until it gave good enough results in all cases.

I mentioned that the algorithm is iterative, but didn’t explain why yet: when we compute the tile alignments, we only do so because we have enough local samples for alignment. We do this for all tiles that match this criteria, but we won’t, for example, be able to align a tile which is in the top of the image, if the bottom hasn’t been reconstructed. Therefore, the iteration happens when we have reconstructed all the tiles we could in one step: then we recompute new reference points, and complete the image, not forgetting, of course, to update our mask to tell that some pixels were completed.

Overall the algorithm is fairly fast, and can be stopped once we have completed all tiles, or after a number of fixed iterations in case of difficulties (often due to the background itself).

The algorithm we’ve described provides us with a way to "roughly" reconstruct an image, but it doesn’t work like what you’d intuititvely think of mosaic composition, by "moving" 2 panels until they properly align and blend them toghether. Instead, it will reconstruct an image by assembling tiles together, from what is already reconstructed: it is more fine grained, which will fix a number of the issues we’ve faced before: local distorsions, or images which are not properly aligned because the details at the surface at the sun have moved between the moment the first panel was captured and the second one did.

If we stopped there, we would see an image which looks like this:

We can see a clear horizontal line, which is due to the fact that we’re reconstructing using tiles, and that depending on the alignment of tiles with the "missing" areas, we can have strong or weak artifacts at the borders. Errors are even more visible in this image in calcium K line:

This time it’s very problematic and we are facing several of the issues we attempted to avoid: details have significantly moved between the north and south panel were captured, which leads to "shadowing" artifacts, and there are also tiling artifacts visible.

However, the image we get is good enough to perform one last step: use it as a reference image in the stacking algorithm we described in the first section of this blog post. The reason stacking works well is because we know we have complete images that we can align. Here, we have roughly reconstructed an image that we can use as a "complete reference". The idea is therefore to take each tile of each panel and "blend" it using the reconstructed reference. Of course, there is one big difference between the stacking in the first section and the stacking we have to do now. We’re not really going to use the reconstructed image, except for aligning tiles together and computing a "weight" for each tile, which depends on the relative luminosity between the reference image tile we’re considering and the corresponding panel tile.

This gives us a pretty good result:

The images we got there are not perfect, which is why I’m not fully satisfied yet, but they are however already quite good, given that it’s all done in a few seconds, using the same software that you’d use to reconstruct Sol’Ex images! In other words, the goal of these features is not to get the same level of quality that you’d get by using your favorite post-processing or mosaic composition software, but good enough to get you a reasonable result in a reasonable amount of time.

For example, on my machine, it takes less than one minute to:

It would have taken several minutes, or even more, using external software, especially for the mosaic part.