Issues with sinc blot option and compiler warnings/bugs

[mcara]

After cleaning up most warnings on Linux with #208 , there are still 4 warnings remaining that suggest that memory is not accessed correctly and this likely would lead to wrong/random results. For the very least we need to understand the algorithm and the reason for negative indices.

Here are the warnings:

  /project/src/cdrizzleblot.c:673:27: warning: array subscript -1 is below array bounds of ‘float[15]’ [-Warray-bounds=]
    673 |                 px = taper[j - 1] / ax;
        |                      ~~~~~^~~~~~~
  /project/src/cdrizzleblot.c: In function ‘interpolate_sinc’:
  /project/src/cdrizzleblot.c:595:11: note: while referencing ‘taper’
    595 |     float taper[INTERPOLATE_SINC_NCONV];
        |           ^~~~~
  In function ‘interpolate_sinc_’,
      inlined from ‘interpolate_sinc’ at /project/src/cdrizzleblot.c:784:12:
  /project/src/cdrizzleblot.c:681:27: warning: array subscript -1 is below array bounds of ‘float[15]’ [-Warray-bounds=]
    681 |                 py = taper[j - 1] / ay;
        |                      ~~~~~^~~~~~~
  /project/src/cdrizzleblot.c: In function ‘interpolate_sinc’:
  /project/src/cdrizzleblot.c:595:11: note: while referencing ‘taper’
    595 |     float taper[INTERPOLATE_SINC_NCONV];
        |           ^~~~~
  In function ‘interpolate_sinc_’,
      inlined from ‘interpolate_sinc’ at /project/src/cdrizzleblot.c:784:12:
  /project/src/cdrizzleblot.c:685:15: warning: array subscript -1 is below array bounds of ‘float[15]’ [-Warray-bounds=]
    685 |             ac[j - 1] = px;
        |             ~~^~~~~~~
  /project/src/cdrizzleblot.c: In function ‘interpolate_sinc’:
  /project/src/cdrizzleblot.c:596:11: note: while referencing ‘ac’
    596 |     float ac[INTERPOLATE_SINC_NCONV], ar[INTERPOLATE_SINC_NCONV];
        |           ^~
  In function ‘interpolate_sinc_’,
      inlined from ‘interpolate_sinc’ at /project/src/cdrizzleblot.c:784:12:
  /project/src/cdrizzleblot.c:686:15: warning: array subscript -1 is below array bounds of ‘float[15]’ [-Warray-bounds=]
    686 |             ar[j - 1] = py;
        |             ~~^~~~~~~
  /project/src/cdrizzleblot.c: In function ‘interpolate_sinc’:
  /project/src/cdrizzleblot.c:596:39: note: while referencing ‘ar’
    596 |     float ac[INTERPOLATE_SINC_NCONV], ar[INTERPOLATE_SINC_NCONV];
        |                                       ^~

[mcara]

I believe the code for the "sinc" interpolator is wrong. In addition to the possible incorrect memory access suggested by the warnings, fundamentally it appears to me that the kernel is not constructed correctly on several levels.

In my analysis, I looked at the "original" (at least pre-C) Fortran code that was later ported to C used in drizzlepac and drizzle packages: https://github.com/spacetelescope/pydrizzle/blob/09897368f51222d33cff530b1bc9ee89bf3b9d21/src/bieval.f#L345-L549 (SPP version is also available).

Last three loops (labels 210, 250, 290) basically performs convolution with normalized (by sumx * sumy) kernel ac[k]*ar[j]. ac and ak are computed in here basically (hand-wavingly) as taper[j - 1] / ax (or / ay) with ax = N - j + {dx} where {dx} is the fractional part of x and N is NCONV from the code.

Before that, the code computes taper as:

C The taper is a cosbell function which is valid for 0 <= x <= PI / 2
C (Abramowitz and Stegun, 1972, Dover Publications, p 76).

C	   Precompute the taper array. Incorporate the sign change portion
C	   of the sinc interpolator into the taper array.
      if (mod (nsinc, 2) .eq. 0) THEN
         sdx = 1.0
      ELSE
         sdx = -1.0
         do j = -nsinc, nsinc 
            dx2 = sconst * j * j
            taper(j+nsinc+1) = sdx * (1.0 + a2 * dx2 + 
     *      a4 * dx2 * dx2) ** 2
            sdx = -sdx
         ENDDO
      ENDIF

• One issue here is that taper seems to be a Hann window but only when integer part of N / 2 is odd and a rectangular window otherwise. Why?
• It is not immediately clear where the sign change "of the sinc interpolator" comes from (represented by sdx). Why does this sign change affect only Hann window but not the rectangular window?
• The biggest issue (for me) is that I do not see the sinc function used anywhere here. abs(taper) itself is computed as cos(pi * j / (N - 1))**2 which is similar to Hann window as described in https://en.wikipedia.org/wiki/Hann_function or https://en.wikipedia.org/wiki/Window_function#Hann_window except I don't know why it is cos**2 instead of sin**2 - maybe it is a variation of power-of-sine window??? The sinc part of https://en.wikipedia.org/wiki/Whittaker–Shannon_interpolation_formula is missing. taper[j - 1] / ax is not sinc - it is just the window that depends sample index j but on on the interpolation position x (or t as used in https://en.wikipedia.org/wiki/Whittaker–Shannon_interpolation_formula).

So, to me, it seems that this code convolves image data with a weirdly normalized (by ax and ay) Hann window without the sinc function itself. This seems to be simply incorrect. (In addition, because of errors in indexing, in testing, it seems that the resulting image is shifted in comparison to where it should have been.)

Am I missing something? If not, I have two alternative proposals:

1. Simply remove sinc interpolator from blot options: temporarily replace it with Lanczos interpolator while it is being deprecated.
2. Implement correct sinc interpolation (similarly to how Lanczos interpolation is implemented).

CC: @nden @jhunkeler @schlafly

[schlafly]

I probably need to engage more thoroughly with this, but no, I don't get it. I gather the a2 and a4 is an approximation of cos(theta), so I guess I see where you're coming from with the Hann window function, but as you say, that doesn't alternate in sign.

The relationship of sinc interpolation and a rectangular window function may be that they are Fourier transform pairs, but I don't understand why even / odd nsinc should be a good choice for changing the window function like that.

That said, there's a lot of loops in that pydrizzle code that I don't immediately understand so mostly I think I need to invest more effort here.