We present a new technique for obtaining model fittings to VLBI images of astrophysical jets. The method minimizes a performance function proportional to the sum of the squared difference between the model and observed images. The model image is constructed by summing elliptical Gaussian sources characterized by six parameters: two-dimensional peak position, peak intensity, eccentricity, amplitude and orientation angle of the major axis. We present results for the fitting of two main benchmark jets: the first, constructed from three individual Gaussian sources, the second formed by five Gaussian sources. Both jets were analyzed by our cross-entropy technique in finite and infinite signal-to-noise regimes, the background noise chosen to mimic that found in interferometric radio maps. We show that our technique is capable of recovering the parameters of the sources with a similar accuracy to that obtained from the traditional AIPS task IMFIT when the image is relatively simple (e.g., few components). For more complex maps, our method displays superior performance in recovering the parameters of the jet components. Our methodology is also able to show quantitatively the number of individual components present in an image. An additional application of the cross-entropy technique to a real image of a BL Lac object is shown and discussed. Our results indicate that our cross-entropy technique must be used in situations involving the analysis of complex emission regions having more than three sources, even though it is substantially slower than current model fitting tasks (at least 10,000 times slower for a single processor, depending on the number of sources to be optimized). As in the case of any model fitting performed in the image plane, caution is required in analyzing images constructed from a poorly sampled (u,v) plane.