If you really want to get rid of stars, DAOPHOT does a damn good job in my experience.
The work flow is to identify all the objects, then a set of unsaturated and unblended stars. From the latter you build a model of the point spread function. Note that the PSF may be position-dependent. Fit one or more PSF models to each detected object, removing ones which are “good enough” in some sense as you go. Repeat until you no longer have any “good enough” fits.