Interpolation of non-rectangular data failed

The problem is that for regular (but not rectangular) meshes the Delaunay mesh is unstable. It's a bug, also mentioned in this question, which was about rectangular grids. The same workaround works here -- just jiggle the points a tiny bit around their perfect lattice positions:

epsilon = 10^-7;
jiggledPlotData = {#1 + RandomReal[epsilon {-1, 1}], #2 + 
      RandomReal[epsilon {-1, 1}], #3} & @@@ plotData;
reg = ConvexHullMesh[jiggledPlotData[[All, 1 ;; 2]]];
f = Interpolation[jiggledPlotData, InterpolationOrder -> 1];
DensityPlot[f[x, y], {x, y} \[Element] reg]