A numerical model for calculating the electrostatic interaction between two particles of arbitrary shape and topology is described. A key feature of the model is a generalized discretization program, capable of simulating any desired analytical shape as a set of flat, triangular elements. The relative sizes of the elements are adjusted using a density function to better match the desired shape and the spatial variation of the electrical surface properties on each particle. The distribution of either surface potential or surface charge density is then calculated using a boundary element approach to solve the linearized Poisson-Boltzmann equation. Example interaction energy profiles are calculated for three different types of roughness-bumps, pits, and surface waves. It is found that the interaction energy between rough particles remains different from that between two equivalent smooth spheres at all separations, even for gap widths much larger than either the solution Debye length or the characteristic roughness size. This behavior at large gap widths arises from the nature of the decay of the electric potential away from each particle. In addition, the magnitude of the roughness effect is found to depend greatly on the size and shape of the nonuniformity as well as the electrostatic boundary conditions. For example, for a sphere containing asperities of height equal to 0.2 times the particle radius, the interaction energy can be as much as 50% greater than that between two equivalent spheres under the condition of constant surface potential. At constant surface charge density, the ratio of the interaction energies between rough and smooth spheres was found to either diverge or become zero as contact between the two particles is approached, depending on the nature of the roughness. Changes of this magnitude could clearly have a substantial impact on the stability behavior of a dispersion of such particles.