Return to Answer

I finally found a relatively short way to do. Not perfect, many special cases should be handled.

With thid code there are the results compared to the routine proposed by @Michael E2 on a uniform grid.

timeAvg = 
  Function[func, 
   Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0, 
     15}], HoldFirst];

\[CapitalOmega] = 
  ImplicitRegion[
   2 x^2 + 3 y^2 + 2 x y - 2 <= 0 \[And] x^2 + y^2 > .1, {x, y}];
\[CapitalOmega]b = RegionBounds[\[CapitalOmega]];

data = Table[Module[{grids, mesh, n},
    grids = 
     N@Map[range \[Function] 
        Range @@ 
         Append[range, -Subtract @@ range/20/2^k], \[CapitalOmega]b];
    (*grids=Union/@Map[g\[Function]{g[[1]],
    Sequence@@(g[[2;;-2]]RandomReal[{1-.02,1+.02},Length[g]-2]),
    g[[-1]]},grids];*)
    n = Length@First@grids;
    mesh = 
     BoundaryDiscretizeRegion[\[CapitalOmega], 
      MaxCellMeasure -> Mean@Flatten[Differences /@ grids]];
    <|
     "Gridlines Count" -> n,
     "AdjustMeshToGrid" -> (AdjustMeshToGrid[mesh, grids, 
         Abs[#1 - #2] <= 10.^-10 &] // timeAvg),
     "snaptogrid" -> (snaptogrid[mesh, Sequence @@ grids] // timeAvg)
     |>
    ], {k, 0, 10}] // Dataset

ListPlot[
 Values@*Normal /@ {data[
    All, {"Gridlines Count", "AdjustMeshToGrid"}], 
   data[All, {"Gridlines Count", "snaptogrid"}]},
 Joined -> True, PlotLegends -> {"AdjustMeshToGrid", "snaptogrid"}, 
 Frame -> True, Mesh -> Full, GridLines -> Automatic, 
 FrameTicks -> Automatic, PlotRange -> All]

I finally found a relatively short way to do. Not perfect, many special cases should be handled. Michael E2 routine permorms better on small uniform grids. I don't know why but on small random grids the rating is reversed.

UPDATE

I sligtly revised the code, mainly to avoid NSolve for obtaining lines intersections and to directly get intersection in the proper order, if possible.

AdjustPolygonToGrid[vertices_, grids_, sameTest_] :=
  Module[{positions, vplist},
   positions = 
    MapThread[
      BinPositions[#1, #2, 
        SameTest -> sameTest] &, {vertices\[Transpose], 
       grids}]\[Transpose];
   vplist = {vertices, positions}\[Transpose];
   Scan[Module[{v1, p1, v2, p2, lines, vlnew},
          {{v1, p1}, {v2, p2}} = #;
          If[AnyTrue[p2, IntegerQ], Sow[v2],
           lines = MapThread[NearestIntegersBetween, {p1, p2}];
           If[lines =!= {{}, {}},
            lines = MapThread[#1[[#2]] &, {grids, lines}];
            
            vlnew = 
             SegmentGridIntersections[Sequence @@ v1, Sequence @@ v2, 
              Sequence @@ lines];
            Sow /@ vlnew;
            ];
           ]
          ] &, Partition[vplist, 2, 1]] // Reap // Last // Last // 
    Append[#, #[[1]]] &
   ];

NearestIntegersBetween = {m, n} \[Function] 
   With[{\[Delta] = n - m, p = IntegerQ[m]}, Which[
     \[Delta] > 2 || \[Delta] == 2 && ! p, {Floor[m] + 1, 
      Ceiling[n - 1]},
     \[Delta] > 1 || \[Delta] == 1 && ! p, {Floor[m] + 1},
     -\[Delta] > 2 || -\[Delta] == 2 && ! p, {Ceiling[m - 1], 
      Floor[n] + 1},
     -\[Delta] > 1 || -\[Delta] == 1 && ! p, {Floor[n] + 1},
     True, {}
     ]];

SegmentGridIntersections =
  {x1, y1, x2, y2, xl, yl} \[Function] 
   Module[{m11 = y2 - y1, m12 = x1 - x2, v1 = x1 y2 - x2 y1},
    Which[
     Length@xl == 0,
     {(v1 - m12*yl)/m11, yl}\[Transpose],
     
     Length@yl == 0,
     {xl, (v1 - m11*xl)/m12}\[Transpose],
     
     True, Join[
        {(v1 - m12*yl)/m11, yl}\[Transpose], {xl, (v1 - m11*xl)/m12}\[Transpose]
        ] // Sort // If[x1 <= x2, #, Reverse@#] &
     ]
    ];

With these updates, for a finer mesh and grid

\[CapitalOmega] = 
  ImplicitRegion[
   2 x^2 + 3 y^2 + 2 x y - 2 <= 0 \[And] x^2 + y^2 > .1, {x, y}];
\[CapitalOmega]g = 
  RegionPlot[\[CapitalOmega], AspectRatio -> Automatic];
\[CapitalOmega]b = RegionBounds[\[CapitalOmega]];
{xg, yg} = 
  N@Map[bound \[Function] 
     Range[bound[[1]], 
      bound[[2]], (bound[[2]] - bound[[1]])/200], \[CapitalOmega]b];
{xg, yg} = 
  Map[g \[Function] {g[[1]], 
     Sequence @@ (g[[2 ;; -2]] RandomReal[{1 - .06, 1 + .06}, 
         Length[g] - 2]), g[[-1]]}, {xg, yg}];

I get an AbsoluteTiming of about 0.3 sec, while the routine of Micheal E2 requires about 5 sec.

MOST RECENT UPDATE

I sligtly revised the code, mainly to avoid NSolve to compute lines intersections and to directly get intersection in the proper order, if possible.

This first helper function, with a binary search algorithm, give the gridline index to wich all coordinates vals given belong, or a half-odd if not on gridline; sameTest allow to configure the accuracy of the check.

BinPositions[vals_, brakes_, sameTest_] :=
 Map[val \[Function] 
   Catch@Module[{lo = 1, mid, hi = Length[brakes], el, res},
     While[lo <= hi, Which[
       TrueQ@sameTest[val, el = brakes[[mid = Floor[(lo + hi)/2]]]], 
       Throw[mid],
       el > val, hi = mid - 1,
       True, lo = mid + 1]];
     lo - 1/2],
  vals]

This second helper function, based on some output of the previous functions, gives the indices of the crossed gridlines. Maybe we can do better.

NearestIntegersBetween = {m, n} \[Function] 
   With[{\[Delta] = n - m, p = IntegerQ[m]}, Which[
     \[Delta] > 1 || \[Delta] == 1 && p, {Ceiling[m], Floor[n]},
     \[Delta] > 0, {Floor[m] + 1},
     \[Delta] == 0 && p, {m, n},
     -\[Delta] > 1 || -\[Delta] == 1 && p, {Floor[m], Ceiling[n]},
     -\[Delta] > 0, {Floor[n] + 1},
     True, {}
     ]];

This third helper function compute the intersections of a line segment with the crossed gridlines, in the order of the segment.

SegmentGridIntersections =
  {x1, y1, x2, y2, xl, yl} \[Function] 
   Module[{m11 = y2 - y1, m12 = x1 - x2, v1 = x1 y2 - x2 y1},
    Which[
     Length@xl == 0,
     {(v1 - m12*yl)/m11, yl}\[Transpose],
     
     Length@yl == 0,
     {xl, (v1 - m11*xl)/m12}\[Transpose],
     
     True, Join[
        {(v1 - m12*yl)/m11, 
          yl}\[Transpose], {xl, (v1 - m11*xl)/m12}\[Transpose]
        ] // Sort // If[x1 <= x2, #, Reverse@#] &
     ]
    ];

This last helper function process a single contour of a MeshRegion:

AdjustPolygonToGrid[vertices_, grids_, sameTest_] :=
  Module[{positions, lines},
   positions = MapThread[BinPositions[#1, #2, SameTest -> sameTest] &, {vertices\[Transpose], grids}]\[Transpose];
   lines = Apply[NearestIntegersBetween,Transpose[Partition[positions, 2, 1], {1, 3, 2}], {2}];
   lines = MapThread[Extract, {grids, Map[List, lines\[Transpose], {2}]}]\[Transpose];
   MapThread[SegmentGridIntersections[Sequence @@ Flatten@#1, Sequence @@ #2] &, {Partition[vertices, 2, 1], lines}] // Flatten[#, 1] & // Append[#, #[[1]]] &
   ];

This main function finally process a whole MeshRegion.

AdjustMeshToGrid[meshRegion_, grids_, sameTest_] :=
 Module[{polygons, vertices, map},
  polygons = meshRegion["BoundaryPolygons"][[All, 1]];
  polygons = AdjustPolygonToGrid[#, grids, sameTest] & /@ polygons;
  
  vertices = DeleteDuplicates[Join @@ polygons];
  map = AssociationThread[vertices, Range@Length@vertices];
  
  (* Restituisce la mesh adattata *)
  BoundaryMeshRegion[vertices, 
   Sequence @@ Line /@ Map[map, polygons, {2}]]
  ]

FIRST ANSWER

UPDATE

I sligtly revised the code, mainly to avoid NSolve for obtaining lines intersections and to directly get intersection in the proper order, if possible.

AdjustPolygonToGrid[vertices_, grids_, sameTest_] :=
  Module[{positions, vplist},
   positions = 
    MapThread[
      BinPositions[#1, #2, 
        SameTest -> sameTest] &, {vertices\[Transpose], 
       grids}]\[Transpose];
   vplist = {vertices, positions}\[Transpose];
   Scan[Module[{v1, p1, v2, p2, lines, vlnew},
          {{v1, p1}, {v2, p2}} = #;
          If[AnyTrue[p2, IntegerQ], Sow[v2],
           lines = MapThread[NearestIntegersBetween, {p1, p2}];
           If[lines =!= {{}, {}},
            lines = MapThread[#1[[#2]] &, {grids, lines}];
            
            vlnew = 
             SegmentGridIntersections[Sequence @@ v1, Sequence @@ v2, 
              Sequence @@ lines];
            Sow /@ vlnew;
            ];
           ]
          ] &, Partition[vplist, 2, 1]] // Reap // Last // Last // 
    Append[#, #[[1]]] &
   ];

NearestIntegersBetween = {m, n} \[Function] 
   With[{\[Delta] = n - m, p = IntegerQ[m]}, Which[
     \[Delta] > 2 || \[Delta] == 2 && ! p, {Floor[m] + 1, 
      Ceiling[n - 1]},
     \[Delta] > 1 || \[Delta] == 1 && ! p, {Floor[m] + 1},
     -\[Delta] > 2 || -\[Delta] == 2 && ! p, {Ceiling[m - 1], 
      Floor[n] + 1},
     -\[Delta] > 1 || -\[Delta] == 1 && ! p, {Floor[n] + 1},
     True, {}
     ]];

SegmentGridIntersections =
  {x1, y1, x2, y2, xl, yl} \[Function] 
   Module[{m11 = y2 - y1, m12 = x1 - x2, v1 = x1 y2 - x2 y1},
    Which[
     Length@xl == 0,
     {(v1 - m12*yl)/m11, yl}\[Transpose],
     
     Length@yl == 0,
     {xl, (v1 - m11*xl)/m12}\[Transpose],
     
     True, Join[
        {(v1 - m12*yl)/m11, yl}\[Transpose], {xl, (v1 - m11*xl)/m12}\[Transpose]
        ] // Sort // If[x1 <= x2, #, Reverse@#] &
     ]
    ];

With these updates, for a finer mesh and grid

\[CapitalOmega] = 
  ImplicitRegion[
   2 x^2 + 3 y^2 + 2 x y - 2 <= 0 \[And] x^2 + y^2 > .1, {x, y}];
\[CapitalOmega]g = 
  RegionPlot[\[CapitalOmega], AspectRatio -> Automatic];
\[CapitalOmega]b = RegionBounds[\[CapitalOmega]];
{xg, yg} = 
  N@Map[bound \[Function] 
     Range[bound[[1]], 
      bound[[2]], (bound[[2]] - bound[[1]])/200], \[CapitalOmega]b];
{xg, yg} = 
  Map[g \[Function] {g[[1]], 
     Sequence @@ (g[[2 ;; -2]] RandomReal[{1 - .06, 1 + .06}, 
         Length[g] - 2]), g[[-1]]}, {xg, yg}];

I get an AbsoluteTiming of about 0.3 sec, while the routine of Micheal E2 requires about 5 sec.