View Issue Details

IDProjectCategoryView StatusLast Update
0000243Siril[All Projects] Sirilpublic2018-03-21 22:06
ReportersteffensAssigned Tolock42 
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
Product Version0.9.8 
Target Version0.9.9Fixed in Version0.9.9 
Summary0000243: Registration: wrong handling of scale variations in drizzle case
DescriptionDuring star matching in registration, I usually get scaleX / scaleY values slightly different than 1.0, e.g. 0.998 or 1.002. This seems to be correct, as the registered images align perfectly.

Except when I enable the drizzle feature (2x scale up). In this case, stars align perfectly only in the upper left corner of the image, while they form diagonal lines (in direction to the upper left corner) in the lower right corner.

From my experiments, it seems that the scaling is not handled correctly there, so the registered images are slightly too large or too small.
When the star matching found a scale factor of 1.002, then I need to apply a scaling of 1.004 in case of drizzle, i.e. I have to double the deviation from the unit vector.

To be specific, for a homography H found in star matching, I have to apply the following to get perfect results in case of drizzle:

    scaleX = sqrt(H->h00 * H->h00 + H->h01 * H->h01);
    scaleY = sqrt(H->h10 * H->h10 + H->h11 * H->h11);
    H->h00 *= ((1.0 + (scaleX - 1.0) * 2.0) / scaleX);
    H->h01 *= ((1.0 + (scaleX - 1.0) * 2.0) / scaleX);
    H->h10 *= ((1.0 + (scaleY - 1.0) * 2.0) / scaleY);
    H->h11 *= ((1.0 + (scaleY - 1.0) * 2.0) / scaleY);

I'm not exactly sure about the reason and the maths there, maybe there's a better way to fix this...
Steps To ReproduceRegister large images with slight scale differences, once without drizzle (works), and once with drizzle (stars in the lower right corner don't align).
TagsNo tags attached.

Activities

lock42

2018-03-13 09:39

administrator   ~0000490

Last edited: 2018-03-13 09:39

View 2 revisions

I think you're probably right.
There's a bug here.
Let's say that H matrix is equal to:
H00 H01 H02
H10 H11 H12
H20 H21 H22
For now, in the "Drizzle" we just multiply H02 and H12 by 2 (registration.c:732-733) which is probably true when no scale is involved.
I have no images with different scales to really test it. Could you make some test and provide a patch ?

lock42

2018-03-13 10:21

administrator   ~0000493

Could you try this patch please ?
This is the changes you told me if I did well understand.

Drizzle_patch-2 (1,180 bytes)
Index: src/registration/registration.c
===================================================================
--- src/registration/registration.c	(révision 2143)
+++ src/registration/registration.c	(copie de travail)
@@ -728,9 +728,19 @@
 					if (!args->translation_only) {
 						fits_flip_top_to_bottom(&fit);	// this is because in cvTransformImage, rotation center point is at (0, 0)
 						if (args->x2upscale) {
-							cvResizeGaussian(&fit, fit.rx * 2, fit.ry * 2, OPENCV_NEAREST);
-							H.h02 *= 2.0;
-							H.h12 *= 2.0;
+							double upscale = 2.0;
+							double scaleX, scaleY;
+
+							cvResizeGaussian(&fit, fit.rx * upscale, fit.ry * upscale, OPENCV_NEAREST);
+							scaleX = sqrt(H.h00 * H.h00 + H.h01 * H.h01);
+							scaleY = sqrt(H.h10 * H.h10 + H.h11 * H.h11);
+
+							H.h02 *= upscale;
+							H.h12 *= upscale;
+							H.h00 *= ((1.0 + (scaleX - 1.0) * 2.0) / scaleX);
+							H.h01 *= ((1.0 + (scaleX - 1.0) * 2.0) / scaleX);
+							H.h10 *= ((1.0 + (scaleY - 1.0) * 2.0) / scaleY);
+							H.h11 *= ((1.0 + (scaleY - 1.0) * 2.0) / scaleY);
 						}
 						cvTransformImage(&fit, ref, H, args->interpolation);
 						fits_flip_top_to_bottom(&fit);
Drizzle_patch-2 (1,180 bytes)

steffens

2018-03-13 10:24

reporter   ~0000494

Yes, this looks like what I have done. I probably won't have the time to experiment today, but hopefully tomorrow evening...

Also, looking at the code now, I wonder why you first do the
  cvResizeGaussian(&fit, fit.rx * upscale, fit.ry * upscale, OPENCV_NEAREST);
and then the homography transform as a separate step.
Wouldn't it be better, to apply the upscaling to the homography and do just one transformation? Something like

H.h00 *= upscale;
H.h01 *= upscale;
H.h10 *= upscale;
H.h11 *= upscale;
H.h02 *= upscale;
H.h12 *= upscale;
cvTransformImage(&fit, ref, H, args->interpolation);

I'll try if this works when I have some time as well.

lock42

2018-03-13 10:28

administrator   ~0000495

Last edited: 2018-03-13 11:08

View 2 revisions

No.
This is because the interpolation needed for the upscale is different.

steffens

2018-03-13 20:16

reporter   ~0000496

Ah, I understand.

I just quickly compiled again with the attached patch, and it works as expected.
(I changed your prepared patch to use the "upscale" variable instead of the factor 2.0)

Drizzle_patch-steffens (1,245 bytes)
Index: branches/0.9/src/registration/registration.c
===================================================================
--- branches/0.9/src/registration/registration.c	(revision 2143)
+++ branches/0.9/src/registration/registration.c	(working copy)
@@ -728,9 +728,19 @@
 					if (!args->translation_only) {
 						fits_flip_top_to_bottom(&fit);	// this is because in cvTransformImage, rotation center point is at (0, 0)
 						if (args->x2upscale) {
-							cvResizeGaussian(&fit, fit.rx * 2, fit.ry * 2, OPENCV_NEAREST);
-							H.h02 *= 2.0;
-							H.h12 *= 2.0;
+							double upscale = 2.0;
+							double scaleX, scaleY;
+							
+							cvResizeGaussian(&fit, fit.rx * upscale, fit.ry * upscale, OPENCV_NEAREST);
+							scaleX = sqrt(H.h00 * H.h00 + H.h01 * H.h01);
+							scaleY = sqrt(H.h10 * H.h10 + H.h11 * H.h11);
+							
+							H.h02 *= upscale;
+							H.h12 *= upscale;
+							H.h00 *= ((1.0 + (scaleX - 1.0) * upscale) / scaleX);
+							H.h01 *= ((1.0 + (scaleX - 1.0) * upscale) / scaleX);
+							H.h10 *= ((1.0 + (scaleY - 1.0) * upscale) / scaleY);
+							H.h11 *= ((1.0 + (scaleY - 1.0) * upscale) / scaleY);
 						}
 						cvTransformImage(&fit, ref, H, args->interpolation);
 						fits_flip_top_to_bottom(&fit);

Drizzle_patch-steffens (1,245 bytes)

lock42

2018-03-13 22:44

administrator   ~0000497

Many thanks for your report and your patch

lock42

2018-03-13 22:46

administrator   ~0000498

Fix committed to Siril (2144).

steffens

2018-03-18 17:34

reporter   ~0000506

I'm afraid I have to reopen this bug, as my patch doesn't work.
In fact, when I try registering with artificial test images to check alignment when enabling the drizzle mode, then I see that my patch is plain wrong, as the original code produces good alignment, while my patch leads to wrong alignment for frames with scaling differences.

So the best for now is probably to revert that change. Sorry for this.

Right now, I'm trying to reproduce the problem with better test images.
What I found so far is:

- The images have to be rather high resolution, as alignment gets worse, the farther pixels are away from the upper left corner. My DSLR images are 6000x4000, so that's probably around the size it needs
- The error seems to come from applying the two transformation (2x scale up, then the homography transformation) in a row. Probably some small numerical error multiplies there

To check this, I made two test stacks. The first is with the original drizzle registration (without my patch). The second is with a simple change that just multiplies the homography transformations by two and does the upscaling and alignment transformation in one step. In the second case, alignment is perfect - see the attached screenshots.

steffens

2018-03-18 17:35

reporter  

1.siril-drizzle-registration-bug.png (1,209,384 bytes)
2.siril-drizzle-one-transformation.png (1,163,300 bytes)

lock42

2018-03-18 18:33

administrator   ~0000507

Last edited: 2018-03-19 11:23

View 2 revisions

Could you provide the code you use for stack 2 ?

steffens

2018-03-18 18:40

reporter   ~0000508

Just doing the upscaling together with the homography transformation, as I proposed above;

                        if (args->x2upscale) {
                            double upscale = 2.0;
                            H.h02 *= upscale;
                            H.h12 *= upscale;
                            H.h00 *= upscale;
                            H.h01 *= upscale;
                            H.h10 *= upscale;
                            H.h11 *= upscale;
                        }
                        cvTransformImage(&fit, ref, H, args->interpolation);

I know, the upscaling is done with the wrong interpolation method now, but I wanted to check, where the misalignment is coming from.

lock42

2018-03-18 23:05

administrator   ~0000509

Could you share two images with this problem ?

steffens

2018-03-19 10:56

reporter   ~0000510

I'll try to find a minimum set of images that will show the problem. But that may take a while, and the problem is hard to see with few images with lots of noise...

Meanwhile, I managed to create artificial test images that show at least the typical error:
- the images align ok in normal registration
- there are problems in drizzle mode
- the problems are most pronounced in the lower right corner, while upper left corner is good

I created the image pair by doing some transformations (scaling, rotation and translation), and then adding some very minor "warp" errors to the stars so that star matching could not lead to perfect results. Still, it seems to be sufficiently good to produce no visible alignment error in the standard case, while i see an offset of more than 8 pixels in the drizzle case.

Maybe a good first step towards avoiding numerical errors would be to have the center of the transformations in the center of the images. Then, the error would grow towards the borders, but the best resolution would be in the center.

Base.png (318,958 bytes)
Transformed.png (367,427 bytes)

lock42

2018-03-19 10:59

administrator   ~0000511

I think that the best thing to do for now is to replace the code by:

                    if (!args->translation_only) {
                        fits_flip_top_to_bottom(&fit); // this is because in cvTransformImage, rotation center point is at (0, 0)
                        if (args->x2upscale) {
                            double upscale = 2.0;

                            H.h02 *= upscale;
                            H.h12 *= upscale;

                            H.h00 *= upscale;
                            H.h01 *= upscale;
                            H.h10 *= upscale;
                            H.h11 *= upscale;
                            args->interpolation = OPENCV_NEAREST;
                        }
                        cvTransformImage(&fit, ref, H, args->interpolation);
                        fits_flip_top_to_bottom(&fit);
                    }

But be careful. The version currently in the server crashes !!! Yesterday we have merged a lot of features. We have to debug. So do not svn up :).

steffens

2018-03-19 12:34

reporter   ~0000513

> The version currently in the server crashes !!! Yesterday we have merged a lot of features. We have to debug. So do not svn up :).

Off topic, and I'm absolutely not a C guy (Java developer myself), but in stacking.c stack_median method I see a local
    fits fit;
and then later
    copyfits(&fit, &gfit, CP_FORMAT, 0);
    gfit.data = fit.data;

It seems wrong to me to reference fit.data in the global gfit.data after leaving the method.
On the other hand, I didn't look at the copyfits method and know nothing about C ...

vinvin

2018-03-19 14:40

administrator   ~0000514

This is made on purpose to avoid copying data. The data is allocated in the function and never freed after its reference is copied to gfit. CP_FORMAT copies metadata to the image, associated with the data.

lock42

2018-03-21 10:36

administrator  

Drizzle_Patch (2,857 bytes)
Index: src/opencv/opencv.cpp
===================================================================
--- src/opencv/opencv.cpp	(révision 2163)
+++ src/opencv/opencv.cpp	(copie de travail)
@@ -238,6 +238,41 @@
 	return 0;
 }
 
+int cvTransformH(Homography *H1, double scale) {
+	Mat H = Mat::eye(3, 3, CV_64FC1);
+	Mat S = Mat::eye(3, 3, CV_64FC1);
+
+	H.at<double>(0, 0) = H1->h00;
+	H.at<double>(0, 1) = H1->h01;
+	H.at<double>(0, 2) = H1->h02;
+	H.at<double>(1, 0) = H1->h10;
+	H.at<double>(1, 1) = H1->h11;
+	H.at<double>(1, 2) = H1->h12;
+	H.at<double>(2, 0) = H1->h20;
+	H.at<double>(2, 1) = H1->h21;
+	H.at<double>(2, 2) = H1->h22;
+
+
+	S.at<double>(0,0) = scale;
+	S.at<double>(1,1) = scale;
+
+	Mat result = S * H * S.inv();
+
+	H1->h00 = result.at<double>(0, 0);
+	H1->h01 = result.at<double>(0, 1);
+	H1->h02 = result.at<double>(0, 2);
+	H1->h10 = result.at<double>(1, 0);
+	H1->h11 = result.at<double>(1, 1);
+	H1->h12 = result.at<double>(1, 2);
+	H1->h20 = result.at<double>(2, 0);
+	H1->h21 = result.at<double>(2, 1);
+    H1->h22 = result.at<double>(2, 2);
+
+	H.release();
+	result.release();
+	return 0;
+}
+
 int cvTransformImage(fits *image, point ref, Homography Hom, int interpolation) {
 	assert(image->data);
 	assert(image->rx);
Index: src/opencv/opencv.h
===================================================================
--- src/opencv/opencv.h	(révision 2163)
+++ src/opencv/opencv.h	(copie de travail)
@@ -21,6 +21,7 @@
 int cvRotateImage(fits *, double, int, int);
 int cvCalculH(s_star *star_array_img,
 		struct s_star *star_array_ref, int n, Homography *H);
+int cvTransformH(Homography *H1, double scale);
 int cvTransformImage(fits *, point, Homography, int);
 int cvUnsharpFilter(fits*, double, double);
 int cvComputeFinestScale(fits *image);
Index: src/registration/registration.c
===================================================================
--- src/registration/registration.c	(révision 2163)
+++ src/registration/registration.c	(copie de travail)
@@ -730,18 +730,10 @@
 						fits_flip_top_to_bottom(&fit);	// this is because in cvTransformImage, rotation center point is at (0, 0)
 						if (args->x2upscale) {
 							double upscale = 2.0;
-							double scaleX, scaleY;
 
 							cvResizeGaussian(&fit, fit.rx * upscale, fit.ry * upscale, OPENCV_NEAREST);
-							scaleX = sqrt(H.h00 * H.h00 + H.h01 * H.h01);
-							scaleY = sqrt(H.h10 * H.h10 + H.h11 * H.h11);
 
-							H.h02 *= upscale;
-							H.h12 *= upscale;
-							H.h00 *= ((1.0 + (scaleX - 1.0) * upscale) / scaleX);
-							H.h01 *= ((1.0 + (scaleX - 1.0) * upscale) / scaleX);
-							H.h10 *= ((1.0 + (scaleY - 1.0) * upscale) / scaleY);
-							H.h11 *= ((1.0 + (scaleY - 1.0) * upscale) / scaleY);
+							cvTransformH(&H, upscale);
 						}
 						cvTransformImage(&fit, ref, H, args->interpolation);
 						fits_flip_top_to_bottom(&fit);
Drizzle_Patch (2,857 bytes)

lock42

2018-03-21 10:36

administrator   ~0000515

I've probably fixed this issue (I hope). Could you please try this patch

steffens

2018-03-21 22:03

reporter   ~0000516

I can confirm that the patch works.
Great, thanks!

lock42

2018-03-21 22:06

administrator   ~0000517

Great,
thanks for the review.

Related Changesets

Siril: 0.9 r2144

2018-03-13 22:46:51

lock42

Details
Fix 0000243: wrong handling of scale variations in drizzle case
Affected Issues
0000243
mod - /branches/0.9/src/registration/registration.c

Siril: 0.9 r2164

2018-03-21 16:30:31

lock42

Details
Trying to fix issue 0000243 Affected Issues
0000243
mod - /branches/0.9/src/opencv/opencv.cpp
mod - /branches/0.9/src/opencv/opencv.h
mod - /branches/0.9/src/registration/registration.c

Issue History

Date Modified Username Field Change
2018-03-13 09:25 steffens New Issue
2018-03-13 09:39 lock42 Assigned To => lock42
2018-03-13 09:39 lock42 Status new => feedback
2018-03-13 09:39 lock42 Note Added: 0000490
2018-03-13 09:39 lock42 Note Edited: 0000490 View Revisions
2018-03-13 09:46 lock42 Severity minor => major
2018-03-13 09:46 lock42 Product Version => 0.9.8
2018-03-13 09:46 lock42 Target Version => 0.9.9
2018-03-13 10:00 lock42 File Added: Drizzle_patch
2018-03-13 10:20 lock42 File Deleted: Drizzle_patch
2018-03-13 10:20 lock42 File Added: Drizzle_patch
2018-03-13 10:21 lock42 File Added: Drizzle_patch-2
2018-03-13 10:21 lock42 Note Added: 0000493
2018-03-13 10:22 lock42 File Deleted: Drizzle_patch
2018-03-13 10:24 steffens Note Added: 0000494
2018-03-13 10:24 steffens Status feedback => assigned
2018-03-13 10:28 lock42 Note Added: 0000495
2018-03-13 11:08 lock42 Note Edited: 0000495 View Revisions
2018-03-13 20:16 steffens File Added: Drizzle_patch-steffens
2018-03-13 20:16 steffens Note Added: 0000496
2018-03-13 22:44 lock42 Note Added: 0000497
2018-03-13 22:46 lock42 Changeset attached => Siril 0.9 r2144
2018-03-13 22:46 lock42 Note Added: 0000498
2018-03-13 22:46 lock42 Status assigned => resolved
2018-03-13 22:46 lock42 Resolution open => fixed
2018-03-13 22:47 lock42 Fixed in Version => 0.9.9
2018-03-18 17:34 steffens Status resolved => feedback
2018-03-18 17:34 steffens Resolution fixed => reopened
2018-03-18 17:34 steffens Note Added: 0000506
2018-03-18 17:35 steffens File Added: 1.siril-drizzle-registration-bug.png
2018-03-18 17:35 steffens File Added: 2.siril-drizzle-one-transformation.png
2018-03-18 18:33 lock42 Note Added: 0000507
2018-03-18 18:40 steffens Note Added: 0000508
2018-03-18 18:40 steffens Status feedback => assigned
2018-03-18 23:05 lock42 Note Added: 0000509
2018-03-19 10:56 steffens File Added: Base.png
2018-03-19 10:56 steffens File Added: Transformed.png
2018-03-19 10:56 steffens Note Added: 0000510
2018-03-19 10:59 lock42 Note Added: 0000511
2018-03-19 11:23 lock42 Note Edited: 0000507 View Revisions
2018-03-19 12:34 steffens Note Added: 0000513
2018-03-19 14:40 vinvin Note Added: 0000514
2018-03-21 10:36 lock42 File Added: Drizzle_Patch
2018-03-21 10:36 lock42 Note Added: 0000515
2018-03-21 16:30 lock42 Changeset attached => Siril 0.9 r2164
2018-03-21 22:03 steffens Note Added: 0000516
2018-03-21 22:06 lock42 Status assigned => resolved
2018-03-21 22:06 lock42 Resolution reopened => fixed
2018-03-21 22:06 lock42 Note Added: 0000517