High angular resolution observations at optical wavelengths provide valuable insights into stellar astrophysics1,2, and enable direct measurements of fundamental stellar parameters3,4 and the probing of stellar atmospheres, circumstellar disks5, the elongation of rapidly rotating stars6 and the pulsations of Cepheid variable stars7. The angular size of most stars is of the order of one milliarcsecond or less, and to spatially resolve stellar disks and features at this scale requires an optical interferometer using an array of telescopes with baselines on the order of hundreds of metres. We report on the implementation of a stellar intensity interferometry system developed for the four VERITAS imaging atmospheric Cherenkov telescopes. The system was used to measure the angular diameter of the two sub-milliarcsecond stars β Canis Majoris and ϵ Orionis with a precision of greater than 5%. The system uses an offline approach in which starlight intensity fluctuations that are recorded at each telescope are correlated post observation. The technique can be readily scaled onto tens to hundreds of telescopes, providing a capability that has proven technically challenging to the current generation of optical amplitude interferometry observatories. This work demonstrates the feasibility of performing astrophysical measurements using imaging atmospheric Cherenkov telescope arrays as intensity interferometers and shows the promise for integrating an intensity interferometry system within future observatories such as the Cherenkov Telescope Array.