Improving the speed of a continuous to discrete operator

Question

I have defined an integral operator hm[func_, m_, n_] that discretizes a continuous 2D function into its pixel counterpart where Dm[x_,y_,m_,n_] is the kernel of the operator. This kernel uses Mathematica's 2D UnitBox function.

Definition of the operator:

Dm[x_,y_,m_,n_] := 1/(PixelWidth)^2 UnitBox[(x - m*PixelPitch /2)/PixelWidth, 
                   (y - n*PixelPitch /2)/PixelWidth] /.{PixelWidth-> 1.2, PixelPitch -> 1.5};  

hm[func_,m_,n_] := NIntegrate[Dm[x, y, m, n] func, {x, -Infinity, Infinity},
                              {y, -Infinity, Infinity}, WorkingPrecision -> 10, 
                              PrecisionGoal -> 10, AccuracyGoal -> 10];

I define an example function that becomes the input to hm[func_,m_,n_]:

ExampleFunc[x_, y_] := Exp[-(x^2/2 + y^2/2)]; 

Here I define the pixel coordinates and generate the output, i.e. the discretized function, using hm[func_,m_,n_]:

NumberOfPixelsX = 6;    
NumberOfPixelsY = 10; 

AllPixelCords = Table[{i, j}, {i, NumberOfPixelsY/2, -NumberOfPixelsY/2, -1}, {j, NumberOfPixelsX/2, -NumberOfPixelsX/2, -1}];
FlattenedAllPixelCords = Flatten[AllPixelCords, {1, 2}];

PixelCoordinates = DeleteCases[DeleteCases[FlattenedAllPixelCords, {_, 0}, Infinity], {0, _}, Infinity];  (* Remove pixel coordinates with 0's*)

OperatorResult = hm[ExampleFunc[x, y], PixelCoordinates[[All, 1]], PixelCoordinates[[All, 2]]];

DiscretizedFunction = Partition[OperatorResult, NumberOfPixelsX];
MatrixPlot[DiscretizedFunction]

My problem is that this operation takes a lot of time to complete when I increase the number of pixels, i.e. when I use greater values for NumberOfPixelsX and NumberOfPixelsY. Is there a way to speed this up?.

Answer 1

Skipping Nintegrate, performing $3 \times 3$ Gauss quadrature by hand, and compiling everything into a neat, listable, and parallelized CompiledFunction:

ExampleFunc[x_, y_] := Exp[-0.5 (x^2 + y^2)];

Block[{dx, dy, gausspts, gaussweights, weights},

 {gausspts, gaussweights} = Most[NIntegrate`GaussRuleData[3, $MachinePrecision]];
 weights = Flatten[KroneckerProduct[gaussweightsx, gaussweightsy]];

 cdiscretize = With[{code = Dot[
      Flatten[
       Outer[
        ExampleFunc,
        Compile`GetElement[P, 1] + dx gausspts,
        Compile`GetElement[P, 2] + dy gausspts
        ]
       ],
      Flatten[KroneckerProduct[gaussweights, gaussweights]]
      ]},
   Compile[{{P, _Real, 1}, {dx, _Real}, {dy, _Real}},
    code,
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ];
 ]

Now lets apply this function:

{NumberOfPixelsX, NumberOfPixelsY} = {512, 1024};
{xmin, xmax} = {-1., 1.};
{ymin, ymax} = {-2., 2.};
Δx = (xmax - xmin)/NumberOfPixelsX;
Δy = (ymax - ymin)/NumberOfPixelsY;
lowerx = Most[Subdivide[xmin, xmax, NumberOfPixelsX]];
lowery = Most[Subdivide[ymin, ymax, NumberOfPixelsY]];
lowerleftcorners = Flatten[Outer[List, lowerx, lowery], 1]; // 
  AbsoluteTiming // First

discretizedFunction = Partition[
     cdiscretize[lowerleftcorners, Δx, Δy],
     NumberOfPixelsY
     ]; // AbsoluteTiming // First

ArrayPlot[discretizedFunction]

0.005842

0.02524

As a sanity check, we can compare to the values of the function ExampleFunc sampled on the pixel centers (this is equivalent to using $1 \times 1$ Gauss quadrature per pizel):

sampledFunction = Outer[ExampleFunc, lowerx + 0.5 Δx, lowery + 0.5 Δy];
Max[Abs[discretizedFunction - sampledFunction]]

1.27155*10^-6

Thus, the error introduced by using a $3 \times 3$ Gauss quadrature might be way below machine precision.