Looking at PhotonTransferCurveDataset.aMatrix below, it has a higher correlation along with the horizontal direction. It should be in the vertical direction to be matched with the LSST sensor orientation. This finding triggered all the slack discussions attached to this ticket. 

Alex Broughton said

Hi Chris, I was going over the quality of the BF correction, and noticed that it was doing ok, but not great.  So I was trying to cross off some potential sources of error. One route I was considering was that the applied kernel was “flipped” due to difference in the image coordinate system and the numpy ndarray indexing. I know that the real pixels have a lower confining field in vertical, or parallel direction, and so the A matrix component in the vertical direction tends to be larger than it is in the horizontal direction, but if it is not applied in the right orientation, it could be causing the results I’m seeing. I was using fiducial behavior to gauge the difference in image coordinates and array indexing.

Below I’ve plotted the A matrix and Kernel[8:, 8:] matrix, with the indexing (happy to give you the code too). And I noticed that the sensor “parallel” direction looks horizontal on these plots, which seems wrong to me, but might be right with respect to the image coordinate system. I’ve spent a while looking at the indexing here, and I think I might be going crazy, but since you wrote most of the BF code, I thought I’d ask you if this is correct.

Thanks! And sorry about this, I know DM people are very busy right now.

Christopher Waters replied

Can you tell me what config options you've used?  I see a comment on line 412 in the quadratic fit code that the i,j indices are inverted to apply the transposition, as is done in the averaging case.  I can see there's a transpose in the averageCorrelations function (which I'm a bit concerned is not doing the right thing either, as that looks like it needs an axes argument to transpose the correlations individually).

One thing that would be useful to check (generally; this probably should be a unit test) would be to make a PTC dataset that has obvious covariance along one direction only, and confirm that the generated BFK corrects it.  This would test both the generation (is that missing axes actually a problem) and application (is the kernel in numpy coordinates and not afw image coordinates).

As I'm tracing it, makeBrighterFatterKernel takes the PTC covariances, generates a BFK from the transposes (I don't remember why, I'll have to reread the paper) and saves that to disk (either on the AMP or DETECTOR levels).  IsrTask reads that in, using the detKernel directly if it exists, otherwise generating an average from the ampKernel entries.  That averaging again transposes the list of amp-wise kernels (with no axes), so AMP and DETECTOR level products are already different here, meaning there's at least one bug.  This numpy array is then packed into an afwImage's array, which should be transposed relative to the image coordinates (meaning a properly made kernel needs to be transposed once more, if it wasn't already transposed before).

In a separate thread, Christopher Waters added links to the possible suspected places:

• this is for line 412: https://github.com/lsst/cp_pipe/blob/main/python/lsst/cp/pipe/makeBrighterFatterKernel.py#L412
• The first possibly incorrect transpose: https://github.com/lsst/cp_pipe/blob/main/python/lsst/cp/pipe/makeBrighterFatterKernel.py#L375
• The second possibly incorrect transpose: https://github.com/lsst/ip_isr/blob/main/python/lsst/ip/isr/brighterFatterKernel.py#L559
• The different handling of the AMP and DETECTOR kernels: https://github.com/lsst/ip_isr/blob/main/python/lsst/ip/isr/isrTask.py#L1044
• The final pack into afwImage format:
  https://github.com/lsst/ip_isr/blob/fc2fcc5b136c9d21df925c9227bc4a7c33306239/python/lsst/ip/isr/isrFunctions.py#L537