Current Shack-Hartmann Wavefront Sensor (SHWFS) reconstruction algorithms for aero-optical research underperform in the presence of flow features with strong density gradients, such as shockwaves in supersonic flow. The large density variations in shockwaves violate the key underlying assumption of SHWFS; that local changes in the wavefront within the lenslet subaperture manifest primarily as tip/tilt. In these cases, the image-plane irradiance pattern of individual lenslets can exhibit high-order aberrations, resulting in non-tip/tilt behaviors. Standard least-squares wavefront reconstruction methods fail to accurately recover the wavefront in the presence of a shock due to the least-squares estimator’s tendency to give too much “influence” to outliers present in the measured SHWFS data, resulting in an underprediction of the optical-path difference (OPD) across the shock. A new algorithm is described to overcome the limitations of the standard least-squares reconstruction method. Two weighting functions are investigated with the aim of using additional intensity information to quantify the degree to which each subaperture is aberrated. The least-squares estimator is replaced with a robust estimator to perform outlier handling and regularizing terms are then used to further constrain the spatial organization of the solution. The problem of wavefront reconstruction is cast as a global functional optimization problem where minimization is achieved iteratively. The algorithm is then evaluated on a sample of Mach 6 SHWFS dot patterns where oblique shocks produce flow discontinuities. The results show that the algorithm is capable of accurately targeting discontinuous flow features as outliers, and subsequently altering those outliers, increasing the OPD as well as the sharpness of the shockwave structure.