s-image.c: fix offsets when inputting an image.

[kglotfelty]

This is for #50

There are a couple of fixes here

• The binary_search routine is finding the last pixel below the random value. However, what is needed is the 1st pixel that includes (>=) the random value. This is causing a 1 pixel shift in X (or a missed wrap if at the edge of the image).
• We want to use floor(x/len) when computing the X coordinate. Even though this is integer divided by integer, it's still producing a floating point value -- but we don't want the fractional part of division. This will also cause a fractional shift in the X coordinate.

The other +1 changes are there to help clarify the change from C-index (0:n-1) to pixel index (1:n).

This looks to fix the case in #50, though it's easier to see when working with an image with odd axis lengths

dmcopy delta.img"[bin #1=1:101:1,#2=1:101:1]" delta101.img
pset marx_delta S-ImageFile=delta101.img

Will definitely need to be run through marx's test suite.

[hamogu]

Will definitely need to be run through marx's test suite.

If I only had such a thing. As far as I can tell, previous development model was to never make errors and thus no tests are needed. (Or more seriously, many things are hard to test, because it's (a) Monte-Carlo, and (b) the truth is often not known, since we can't asssume CIAO to be correct and static either. So, tests are mostly "does it run without a crash?" and "Does it look fine by eye?".

In other words, I'll have to look at this is detail because I can't rely on tests telling me if this is the correct way to change things.

[kglotfelty]

That's unfortunate 😞 ... but I'm also familiar with that modus operandi. Let me know if you need more info.

FWIW: to debug I cut the input image down to a 5x5 pixel image so that it was easier to see exactly which pixel was being used.

dmcopy "delta.img[sky=box(4097,4097,5,5)]" 5x5.img

[DougBurke]

Do we know if this will make if doe the December 2023 release?

[hamogu]

Yes, it will.

[hamogu]

I usually prepare the next marx release when I have the new version of the CALDB, because the contamination files get reformated and included into the marx release.

[hamogu]

So, with this fix, it moves the center of the simulated distribution from (4096, 4097) to (4097.5, 4097.5). I thought that in Chandra, the pixel center is an integer number, so it should have been moved to (4097.0, 4097.0), right?

So, maybe

	/* Convert from C-array index to image pixel index */
	y += 1;
	x += 1;

should be changed to

	y += 0.5;
	x += 0.5;

These questions on "is the corner (0,0) or (1,1) or (0.5, 0.5)" always confuse me...

[hamogu]

I should say that I get "the center" as

r = pycrates.read_file('delta_out.fits[cols x,y]')
 r.get_column('x').values.mean()
 r.get_column('y').values.mean()

because I'm afraid I'm shooting myself in the foot when I do:

dmcopy delta_out.fits"[bin x=4046.5:4146.5:1,y=4046.5:4146.5:1][opt type=i4]" delta_out.img clob+
dmstat "delta_out.img" cen+

Or is looking at the X,Y in the event file not valid and it only becomes OK when I look at the binned image?

[kglotfelty]

I don't like comparing centroids computed from images w/ those computed from events for this kind of reason ; with image the weight is all at the center of the pixel, with events it's {whatever} (to within the pixel boundaries). So they're really not the same.

It's been a while since I've looked, but the +1 there is, I think, just going from C-arrays starting at 0, and the WCS defining the low-left pixel pixel as (1,1).

The bin edge really depends on the bin syntax, so in this example the 4046.5:4146.5:1, the 1st pixel includes x values from 4046.5 up to but not including 4047.5, so 4047 would be the middle of the pixel. But if the bin was 4046:4146:1, then the xcenter of the 1st pixel would be 4047.5.

[hamogu]

Thanks. I'll figure it out ;-)

Note: I renamed branches in the repro for unrelated reasons after I used a bunch of cherry-picking to back out (but keep in a separate branch) some commits I don't want to release yet. When I reset this PR to go onto the branch that's currently called "main" (but had a different name when this PR was opened), all of those cherry-picked commits get added to this PR. Thus, I've manually taken this one commit (by using another cherry-pick!) and put int locally onto my main; I'll either update this PR later or merge locally and push directly to main.

[hamogu]

Suggested change

	y = (double) floor(ofs / X_Image_Size); // floor() because we want the integer part only here
	y = (double) (ofs / X_Image_Size);

Non need for "floor" - / is integer division in C, the result will always be the integer part, it's not going to do any hamr iether though: https://stackoverflow.com/questions/12240228/c-integer-division-and-floor

[hamogu]

Suggested change

	// binary_search is finding the last pixel less than JDMrandom().
	// We want the 1st pixel that includes JDMrandom(), so +1.
	ofs = JDMbinary_search_f (JDMrandom (), Image, Image_Size) +1;
	ofs = JDMbinary_search_f (JDMrandom (), Image, Image_Size);

The comment "last pixel less than" I think is based on a misunderstanding of the indexing scheme. It's easier to think of this a zero-indexed array - and it returns the index into a zero-indexed array.

[hamogu]

Suggested change

	if (ofs >= Image_Size) {
	ofs = Image_Size-1;
	}

[hamogu]

It turns out that this is not a problem from the image source, but instead there was an underlying bug in the binary search in JDmath, that's been there for at least two decades and also effects other parts of marx - for example when selecting with energy bin to put a photon in a CCD level ARF. About 50% of the time, the result of the binary search would be off by 1. Fortunately, our ARF binning typically oversamples so it's not a big impact in practice, but it's annoying nevertheless and shows why one should use well-tested standard libraries instead of implementing your own math functions!

I'm also adding tests based on the example in #50 to the marx test suite which is slowly growing now that I've found a way to use a single framwork (pytest) to drive both tests of the functions C code itself (through CFFI) and end-to-end tests that involve CIAO tools (through ciao_contrib.runtool).

Thank you so much for @kglotfelty for your original investigation and I apologize that it's taken me so long to take the time to get to the bottom of this and fix it.