These are all great points and I’m glad this is being discussed.
Bob Benjamin of the Spitzer GLIMPSE team recently told me that for the GLIMPSE (right in the MW midplane) data processing it was important for background estimation to take into account and subtract all sources, even 1-2 sigma detections (that were consistent with being point stars). They didn’t put those marginal detections into their final output catalog, but they found it important to subtract them out to find the background otherwise their fluxes were wrong. The expert on this is Brian Babler at UW-Madison. Just thought I’d pass that bit of info along.
For PSF estimation/definition in crowded fields, it is probably best to use an iterative procedure. At the beginning, use a simple analytic PSF model for the brightest stars and only use pixels close to the centroids. Then use that plus the list of detected sources to subtract-off fainter stars, but leave in the bright ones and use them to get a better estimate of the PSF. Iterate until you converge. Since there are faint stars undetected and/or very blended in the original image so you probably need iterative detection/subtraction to find all of those as well. Therefore, you might need two nested iterative loops: one for the PSF determination (the outer loop), and one for detection/subtraction (inner loop).
Having PSF estimation/determination, detection, measurement and deblending more closely linked is probably a good idea. To facilitate that it would be nice if these steps weren’t all so baked-in to processCcd, and you could more easily slap them together as needed.