error in pixel count of a matched polygon region

[kswang1029]

Describe the bug
When there are two images with identical WCS header matched spatially, a polygon region may give inconsistent pixel count as seen with the statistics widget on the matched image. Here is "inconsistent" does not mean they have to be identical on both image as CARTA adopted polygon resampling again in order to perform region reprojection on matched image so it is possible to have few pixel differences. The issue here is the difference appears too large as expected.

To Reproduce
Steps to reproduce the behavior:

1. download the test images at https://drive.google.com/file/d/1uv2f-quEVsf1r9Imt_SwNS6EFv1NvHmO/view?usp=sharing
2. load image01.fits and image02.fits and match them spatially
3. load the region file polygon.crtf
4. load two statistics widget
5. configure one to show image01.fits and the region explicitly and configure the other one to show image02.fits and the region explicitly too
6. check the NumPixels entry and we see quite different pixel counts

Expected behavior
The expected NumPixels from both images should be identical theoretically but it is fine to have just few pixels differences (understandable because of the resampling process).

Screenshots or videos

Platform info (please complete the following information):

• OS [e.g. macOS Monterey]: macOS Monterey
• Browser [e.g. chrome, safari, electron app]: electron app
• Browser version [e.g. 22]:
• Backend branch [e.g. dev, v3b2 release]: v3b3 release
• Frontend branch [e.g. dev, v3b2 release]: v3b3 release

[pford]

@kswang1029 I have been investigating this issue., and here is what I have found out so far. It looks as if the 1000-pt polygon approximation is correct; if I save the points to a .crtf file and import the region it overlays the 11-pt polygon:

If I close image01 so that image02 owns the 11-pt polygon, then export the region as a subimage, it looks as expected and its image number of pixels is 15974, the same as for image01 polygon region:

But if I export the matched polygon as a subimage, it shows NaN pixels within the exported region and its image number of pixels is 14505, the same as for image02 matched polygon region:

I am not sure what to do about this issue.

[kswang1029]

@pford thanks for exploring this. I will try a similar comparison from here too. 🤔️

Do you think it is world1->image1->image2->world2 coordinate conversion chain having numeric errors? In your last image with NaN inside the polygon, we see horizontal border at the bottom. The NaN strip seems to be perfectly horizontal too. 🤔️

[kswang1029]

@pford Out of curiosity I modified

carta-backend/src/Region/Region.h

Line 29 in eae1202

#define DEFAULT_VERTEX_COUNT 1000

to 2000 and 100. With 2000, the issue remains. Surprisingly, with 100, the issue is gone (don't know why though. probably related to numeric error or rounding numbers).

In the following example, we can see the polygon and slightly rotated rectangle have identical pixel counts respectively (differ by a few pixels is fine and understandable). The ellipse, however, has a little error.

It feels like we probably need two sets of DEFAULT_VERTEX_COUNT, one for ellipse/line(? need experiments)/polyline(? need experiments) (such as 1000 as we adopted now) and one for polygon/rectangle (such as 100, need experiments).

[kswang1029]

This might be more complicated as we may need a dynamic VERTEX_COUNT depending on sky curvature (ie sky projection distortion), if the above two-set-DEFAULT_VERTEX_COUNT does not work in general. 🤔

[pford]

@kswang1029 I dug into the casacore code to see how the mask is set for a polygon region. The mask starts as completely false (no pixels included). First it iterates through the control points and if two adjacent points have y-coordinates "near" each other, it flags the point as a horizontal line. Then it iterates through each row of the bounding box (looping over y) and sets the mask to true where the x pixels are in the region. For a row containing a point flagged as horizontal, if the point is close to an integral pixel value then the x-range is masked to the x-value of the next point. This is why in the saved-region image there are a few isolated pixels indicating where the region boundary should be, because the next point is only 1 or 2 pixels away. This is why using fewer points worked, because either the row was not flagged as horizontal or the x-range to the next point was further away.

My solution was first to set the number of vertices in the matched polygon to the pixel length of the original polygon (no need to have more than one point per pixel) or the default 1000, whichever is less. Then I set the number polygon points in a nearly horizontal segment to the delta-y value. So for example instead of 256 points along the bottom of the matched polygon there are only 7 because the difference in y between points 6 and 7 rounded to 7. The non-horizontal portions of the polygon are approximated as usual.

This resulted in a saved-region image from image02 that looks like the one from image01. The number of pixels matches for the region in both images but the other statistics differ; this is also the case when the default number of vertices is changed to 100.

Attached: saved images of the original polygon points for ~horizontal segments, the new points, and the subimage saved from image02 and the matched polygon.

[kswang1029]

@pford Thank you for digging this out! So it is the mask generation algorithm.

Regarding the proposed solution, do you mean we still sample each polygon side with 1000 points first on the reference image, then transfer the new lots-control-point polygon to the matched image, then for nearly horizontal sides (reprojected), we remove unnecessary control points? Or we simply use few control points to sample nearly horizontal sides on the reference image and transfer the re-sampled polygon to the matched image? I wonder if this would introduce a slightly larger field of view area error if the matched image has a very different projection scheme?

[pford]

I applied the horizontal line test for the approximated polygon in the reference image. But I could do this after the points are applied to the matched image and remove points before creating the matched polygon.

[kswang1029]

I applied the horizontal line test for the approximated polygon in the reference image. But I could do this after the points are applied to the matched image and remove points before creating the matched polygon.

I see. I suggest we apply the horizontal line check in the matched image then, as a horizontal line in the reference image may not still be a horizontal line in the matched image (projection scheme dependency), and vice versa.

[pford]

I converted the polygon to the matched image then removed some of the horizontal points. The number of pixels is close but not exact. I only tested with the polygon attached in this issue.

[kswang1029]

I converted the polygon to the matched image then removed some of the horizontal points. The number of pixels is close but not exact. I only tested with the polygon attached in this issue.

do you consider half-pixel when we remove control points in the horizontal line test? IIRC, if the pixel center is inside the polygon or on the polygon side, it counts. Maybe this is why we see slightly different pixel counts if the matched image has exactly the same image dimensions and pixel size. 🤔 I think the goal is to have matched pixel count when the reference and matched images are identical. Then for a more general use cases with different image sizes/pixel sizes, we would get a good approximation using the same horizontal line test.

[pford]

I modified the point selection so the number of pixels now matches.

A point is selected if near a pixel , e.g. y=5048.04 is near 5048. The problem was that with so many points for the polygon, the next points 5048.02 and 5048.00 were also selected, so these additional points triggered the "horizontal line" mask and added 2-3 pixels (delta x between these points) to the x-range for row y=5048. Now only the point y=5048.00 is selected for the matched polygon, which is essentially what LCPolygon does for the mask when it finds where the original line segment with only the two endpoints crosses y=5048.