Empirical models of bow shock and magnetic pileup boundary locations are typically based on the identification of individual boundary crossings and their subsequent fitting by properly chosen dependences. Such an approach, however, requires a large set of identified crossings, whose compilation can be easily a source of a significant bias. Moreover, the method is inherently biased by the spacecraft orbit: the more time the spacecraft spends in a given region, the more likely it is for a crossing to be identified in there. We use a different approach based on an automated region identification and Mars Atmosphere and Volatile Evolution (MAVEN) spacecraft data to derive empirical models of both the bow shock and magnetic pileup boundary locations around Mars. We use statistically known parameters in the solar wind, magnetosheath, and induced magnetosphere, along with the observed ratios of measured solar wind parameters, to automatically identify the region where the spacecraft is located at any given time. A simple empirical relation is then assumed for a boundary shape and location, and its free parameters are adapted to optimize the resulting model classification of individual data points. This procedure allows us to model both the bow shock and magnetic pileup boundary locations, reproducing successfully observed variations with the solar wind dynamic pressure, solar ionizing flux, and crustal magnetic fields. However, due to the sparse data coverage, the models are deemed unreliable beyond the terminator. Plain Language Summary Solar wind interaction with Mars results in the formation of two distinct boundaries, bow shock and magnetic pileup boundary. Their exact locations depend on various parameters, and they are typically described by empirical models based on individual boundary crossings identified in spacecraft data. However, this requires many boundary crossings to be identified, and moreover, the locations of these crossings are inherently biased toward the regions where the spacecraft spends more time. We use a different approach based on an automated identification of whether the spacecraft is located inside or outside a given boundary. An empirical dependence of the boundary shape and location is then designed in such a way that the model region classification fits the observations as closely as possible. This removes the necessity to rely on the crossing identification and the related sampling bias. Empirical models of the boundaries derived by our simple and robust method reproduce well observed dependences on solar wind dynamic pressure, solar ionizing flux, and remnant crustal magnetic fields.