
dotSize[{n1_, n2_}, {N1_, N2_}] :=
  (N1 * n1 + n2 + 2) / (3 (N1^2 + N2));

plotResampling/:
  Options[ plotResampling ] =
    { AspectRatio -> 1,
      Axes -> True };

plotResampling[ resamplingMatrix_, N1_, N2_, opt___ ] :=
  Module[ { allxypairs, cartesianLattice, cartesianPlot, integerTickMarks,
            maxx, maxy, minx, miny, n1, n2, newLattice, newLatticePlot,
            options },

          options =
            Join[ Flatten[ { opt } ],
                  Options[ plotResampling ] ];

          (* Find the extent in x and y that encompasses both plots *)
          allxypairs =
            Join[ Apply[Join,
                        Table[ {n1, n2}, 
                               {n1, 0, N1},
                               {n2, 0, N2} ] ],
                  Apply[Join,
                        Table[ resamplingMatrix . {n1, n2}, 
                               {n1, 0, N1},
                               {n2, 0, N2} ] ] ];
          xvalues = Map[ First, allxypairs ];
          minx = Min[ xvalues ];
          maxx = Max[ xvalues ];
          yvalues = Map[ Part[#, 2]&, allxypairs ];
          miny = Min[ yvalues ];
          maxy = Max[ yvalues ];
          integerTickMarks =
            { Table[ n1, {n1, minx, maxx} ],
              Table[ n2, {n2, miny, maxy} ] };

          AppendTo[ options,
                    PlotRange -> { { minx - 0.5, maxx + 0.5 },
                                   { miny - 0.5, maxy + 0.5 } } ];
          AppendTo[ options, Ticks -> integerTickMarks ];

          (* Create and plot the Cartesian lattice *)
          cartesianLattice =
            Graphics [
              Table [ Disk[ {n1, n2}, dotSize[{n1, n2}, {N1, N2}] ],
                      {n1, 0, N1},
                      {n2, 0, N2} ] ];
          cartesianPlot = Show[ cartesianLattice, options ];

          (* Create and plot lattice corresponding to the resampling matrix *)
          newLattice =
            Graphics [
              Table [ Disk[ resamplingMatrix . {n1, n2},
                            dotSize[{n1, n2}, {N1, N2}] ],
                      {n1, 0, N1},
                      {n2, 0, N2} ] ];
          newLatticePlot = Show[ newLattice, options ];

          { cartesianPlot, newLatticePlot } ];

plotResampling[ {{ 1, 1 }, { 1, -1 }}, 4, 4 ]
