Probably the easiest way to achieve this is by taking the plane to define a rotated and shifted coordinate system. This allows you to construct the matrices for transforming a point in global coordinates into plane coordinates and back. Once you have this, you can simply transform the point into plane coordinates, perform the rounding/projection in a trivial manner, and convert back to world coordinates.
Of course, the problem is underspecified the way you pose the question: the transformation you need has six degrees of freedom, your plane equation only yields three constraints. So you need to add some more information: the location of the origin within the plane, and the rotation of your grid around the plane normal.
Personally, I would start by deriving a plane description in parametric form:
xVec = alpha*direction1 + beta*direction2 + x0
Of course, such a description contains nine variables (three vectors), but you can normalize the two direction vectors, and you can constrain the two direction vectors to be orthogonal, which reduces the amount of freedoms back to six.
The two normalized direction vectors, together with the normalized normal, are the base vectors of the rotated coordinate system, so you can simply construct the rotation matrix by putting these three vectors together. To get the inverse rotation, simply transpose the resulting matrix. Add the translation / inverse translation on the appropriate side of the rotation, and you are done.