The conventional approach for preoperative parathyroid localization with the dual 99mTc-sestamibi (99mTc-MIBI) and 99mTc-pertechnetate (99mTcO4- ) scintigraphy technique obtains the differential image by subtracting images from the two scans; the former depicts both the thyroid and parathyroid glands while the latter depicts the thyroid only. In this study, we developed a novel method based on Poisson noise modeling and maximum-likelihood estimation to generate the differential image in an iterative fashion using both planar images jointly. We demonstrated improved performance of our joint method as compared with the subtraction method in both phantom and patient studies. In the phantom study, we acquired two noise-free planar datasets using 99mTc on an in-house thyroid phantom and a "lesion" bead (representing a parathyroid gland) with the same attenuation background as the thyroid phantom. These two planar datasets were combined and scaled to simulate noise-free clinical MIBI (four lesion-to-background contrast (LBC) values: 1.2, 1.3, 1.4 and 1.5), and 99mTcO4- datasets. One-hundred Poisson noise realizations were generated for each datasets. The mean and standard deviation (SD) of the lesion contrast in the differential images were computed for both the subtraction and the joint methods. We also applied both the subtraction and the joint methods to one parathyroid patient dataset. The voxel-wise mean-to-SD ratios in four hyperfunctioning parathyroid lesions were calculated. The phantom results showed that the joint method at the 50th iteration yielded a significant SD reduction compared with the subtraction method ranging from 20% to 45% (p <; 0.05). Similarly, the voxel-wise mean-to-SD ratios were substantially improved in the patient study from 0.40-1.60 (subtraction) to 2.68-3.16 (joint).