Hi @Joel_David_Berkson, I meant to get to this earlier and got caught up in the holidays, so my apologies. The Southwell integration method can have some serious issues particularly for non-symmetric or ‘funny’ shaped apertures. This is strongly apparent in cases were there are a) not a number (NaNs) inside of your aperture, where the aperture has an edge (polygons, rectangles, etc), and c) where there is ‘bad’ data in your slope data. This is because the zone is what is being considered in the integration approach, and this has ‘bleed out’ effects, where one point can impact surrounding data.
As you said, Chebyshev polynomials are particularly convenient for rectangular data, and Zernike polynomials are useful for circular aperture data. The historical limit to these approaches are that is can be really resource intensive to generate enough polynomials to fit high frequency data AND that there can be a risk of overfitting. As an extra ‘and’, you also should be cautious of how the polynomials are being generated. If they are hard-coded for example, did the author check up to whatever order they coded to verify there are no errors? If they are recursively generated, is there any error that may stack on higher orders?
Maham Aftab, as you are familiar with, came up with the recursive Chebyshev Gradient polynomial set which does a great job of generating large numbers of polynomials, a great solution for rectangular apertures. You can also apply these to polygon apertures, but try to consider how this is impacted by the rectangular orthonormality of the base set (http://www.loft.optics.arizona.edu/documents/journal_articles/Aftab2019_Article_ModalDataProcessingForHighReso.pdf). A few comments on this approach, check in the slope domain to see if when you code it there is a scaling issue. The weighting factor referenced in the paper can be dropped when you code it. Recursive generation of these sets using a matrix approach is just stupidly resource intensive when coded in Matlab, I would recommend a vector approach or coding it in some other language.
You can also now find various fast Zernike polynomial set generators for fitting data. See:
Specifically, I think this recent article may be particularly useful to you:
Finally, @bdube has made an excellent open source optics engine called Prysm (prysm — prysm 0.20 documentation) which offers many general optics tools, including a fast and stable wavefront fitting tool. I leave the technicalities of how this could be usable to you, but the tool offers Zernike, Chebyshev, and Lagrange fitting, all relating back to the Jacobi bases these share ( How Zernike, Legendre, and Chebychev Polynomials relate to Jacobi Polynomials). This is also something we are focused on providing in our toolset as we move forward (we also have suffered through slope fitting 
).
TL;DR: Slope fitting is complex, check out Maham’s method and particularly check out this recent paper for a Zernike approach: Average gradient of Zernike polynomials over polygons
This is a great question that I hope other members will chime in on. I will plan to write up a guide though in the near future on slope fitting methods!