Implementation of MATLAB's "numgrid" function

Question

I am looking for an implementation of MATLAB's numgrid function, in particular the "B" mode.

For example, I want to get the matrix corresponding to the following matrix in MATLAB:

numgrid('B', 80)

Answer 1

The source code of some of the MATLAB functions is available to its users. In this example, if you type edit numgrid, its file will be open, and to my surprise, it wasn't very complicated, so I decided to rewrite it with Mathematica's built-in functions (I'm sure there is room for improvement).

Update

Thanks to @chyanog's suggestion to remove Thread in building SparseArray, it performs almost 3 times faster.

ClearAll[NumGridB];
NumGridB[n_Integer] :=
 Block[{x = ConstantArray[Subdivide[-1., 1, n - 1], n], y, t, r, temp,
    Indeterminate = 0.},
y = Transpose[x[[All, -1 ;; 1 ;; -1]]];
t = Quiet@ArcTan[x, y];
  r = Sqrt[x^2 + y^2];
temp = SparseArray[
   1 - Unitize[
     UnitStep[
       r - (Sin[2 t] + 0.2*Sin[8 t])] + (1 -
        UnitStep[Abs[x] - 1]) + (1 - UnitStep[Abs[y] - 1]) - 3]];
SparseArray[
 SortBy[temp["ExplicitPositions"], Last] ->
  Range[temp["ExplicitLength"]], {n, n}, 0]
  ]

Notes:

• I used the same variable names as MATLAB, except temp, so you can follow its logic.
• In Block scope, I used Indeterminate=0. with Quiet because in Mathematica ArcTan[0, 0] returns Indeterminate with a Message, but in MATLAB, atan2 simply return 0 - @chyanog suggests using Arg[x+I*y], but since it's a little slower, I decided not to use it.

Example:

Mathematica:

ArrayPlot[NumGridB[7], Frame -> False]

MATLAB:

imagesc(numgrid('B',7))

Higher number (n=500):

Speed

n	MATLAB (second - tic;toc;)	Mathematica (second - RepeatedTiming)
500	0.006	0.18 0.044
5000	0.55	17.9 4.75

Without touching *Compile, for not large n, its timing is reasonable. If you did improve the speed, feel free to post it here.

Update 2

Other answers show clearly compiling will further boost the performance. I decide to use the wolfram-library-link-rs crate to use Rust for building a LibraryLink library.

You can get the code and compiled library from this GitHub page.

At first glance, it wasn't much faster than others (clearly slower than @chyanog solution!), but if you open Task Manager while running MATLAB numgrid, you'll notice is not a single-threaded command, So this became an opportunity to leverage Rust's power in multithreading to set a new record.

After downloading the .dll file from the Release Section of the GitHub page, execute the following command and you're good to go:

(* Single-threaded *)
NumGridBCompiled =
 LibraryFunctionLoad[
  "C:\\numgrid.dll", "numgrid_b", {Integer},
  LibraryDataType["NumericArray", "UnsignedInteger32", 2]]
(* Multi-threaded *)
NumGridBParallelCompiled =
 LibraryFunctionLoad[
  "C:\numgrid.dll",
  "numgrid_b_parallel", {Integer},
  LibraryDataType["NumericArray", "UnsignedInteger32", 2]]

Benchmark

Speed

• For "Ben Izd", AbsoluteTiming was used.

Table (all the timings are in seconds):

n	MATLAB	NumGridBParallelCompiled	NumGridBCompiled	chyanog	xzczd	Ben Izd
1,000	0.023	0.0039	0.033	0.016	0.061	0.16
5,000	0.61	0.25	0.84	0.63	1.49	4.72
10,000	2.4	1.47	3.39	2.06	5.98	46.1

Memory

Although I'm not fully aware of Mathematica accuracy on LibraryLink calls, here is the plot:

Table (all the sizes are in bytes):

n	MATLAB	NumGridBParallelCompiled	NumGridBCompiled	chyanog	xzczd	Ben Izd
1,000	7,864,000	4,000,256	4,000,256	24,020,064	16,000,432	95,388,352
5,000	195,764,000	100,000,256	100,000,256	600,085,664	400,000,432	2,390,195,048
10,000	782,712,000	400,000,256	400,000,256	2,400,163,872	1,600,000,432	9,563,573,936

Notes:

• MATLAB timings were based on multiple tic;toc; and I try a slightly different size and the timings were around the same number.
• MATLAB memory usage was based on its profiler with profile('-memory','on'), and its Peak Memory (values were in Kb, so it was converted accordingly)
• MATLAB version is 9.12.0 (2022b update 2), Mathematica 13.1.0, and system is AMD Ryzen 3950X

Answer 2

After searching for a while I found the source code of numgrid here. So, here's my trial with FunctionCompile:

bfnumgrid =
   FunctionCompile@
    Function[Typed[n, "Real64"], 
     With[{atan2 = {y, x} |-> If[x == 0 && y == 0, 0., ArcTan[x, y]], 
           sqrt = Sqrt, sin = Sin}, 
       Module[{t, r, i = 1},
         Table[t = atan2[y, x];
               r = sqrt@(x^2 + y^2); 
               If[r >= (sin@(2*t) + .2*sin@(8*t)) && -1 < x < 1 && -1 < y < 1, 
                  i++, 0], 
               {x, -1, 1, 2/(n - 1)}, 
               {y, 1, -1, -2/(n - 1)}]\[Transpose]]]]; // AbsoluteTiming
(* {9.24076, Null} *)
bfshape[5000]; // AbsoluteTiming
(* {1.57289, Null} *)

Be careful with the subtle differences between MATLAB and Mathematica e.g. atan2 v.s. ArcTan, row-major v.s. column-major, etc. I don't choose the Arg trick shown by chyanog because it turns out to be a little slower compared with the If[…] work-around in my current implementation.

For comparison, Ben lzd's solution takes about 6.55 s on my laptop for n == 5000.

Visualization:

ArrayPlot[bfnumgrid[500], ColorFunction -> "AvocadoColors"]

Answer 3

Clear[cf, numGridB];
cf = Compile[{{n, _Integer}, x, {range, _Real, 1}},
   Table[With[{r = Sqrt[x^2 + y^2], t = Arg[x + I y]},
     Boole[r >= Sin[2 t] + 0.2*Sin[8 t] && Abs@x < 1 && Abs@y < 1]], {y, range}
    ], RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed", 
       CompilationTarget -> "C"
   ];

numGridB[n_] :=
  With[{range = Range[-1., 1, 2/(n - 1)]},
    Transpose@ArrayReshape[# Accumulate@#, {n, n}] &@Flatten@cf[n, -range, range]
   ];

numGridB[10] // MatrixForm
numGridB[5000]; // RepeatedTiming

$Version

13.1.0 for Microsoft Windows (64-bit) (June 16, 2022)

For GCC (mingw-w64) compiler, there is no significant difference between -O2, -O3, -Ofast flag

{0.840892, Null}

For VC compiler, /O2 flag

{0.596376, Null}

For Clang compiler, same flag as GCC

{0.561526, Null}

Compared with MATLAB

Another implementation:

Clear[numGridB2];
numGridB2 = Compile[{{n, _Integer}},
   Module[{x, y, t, r, bin, tol = 2.*^-16},
    y = ConstantArray[Range[-1., 1, 2/(n - 1)], n];
    x = -Transpose[y];
    t = ArcTan[x + tol, y];
    r = Sqrt[x^2 + y^2];
    bin = Flatten@BitAnd[UnitStep[r - (Sin[2 t] + 0.2*Sin[8 t])], 
       UnitStep[1 - tol - Abs@x], UnitStep[1 - tol - Abs@y]];
    Transpose@Partition[bin Accumulate@bin, n]
    ]];
numGridB2[10] // Grid
numGridB2[5000]; // RepeatedTiming

{0.980508, Null}

Answer 4

This is a little slower than Ben's version, but I thought I'd share this anyway since the "ExplicitLength" and "ExplicitPositions" properties are not known to earlier versions, and I wanted to show off the use of CoordinateBoundsArray[]:

numGridB[n_Integer?Positive] := Module[{pos, t, temp, x, y},
   {x, y} = Transpose[CoordinateBoundsArray[{{-1., 1.}, {-1., 1.}},
                                            {Into[n - 1], Into[n - 1]}],
                      {3, 2, 1}];
   y = Reverse[y]; t = ArcTan[x + $MachineEpsilon, y];
   temp = SparseArray[1 - Unitize[UnitStep[Sqrt[x^2 + y^2] -
                                           (Sin[2 t] + 0.2*Sin[8 t])] -
                      UnitStep[Abs[x] - 1] - UnitStep[Abs[y] - 1] - 1]];
   pos = temp["NonzeroPositions"];
   SparseArray[SortBy[pos, Last] -> Range[Length[pos]], {n, n}, 0]]

For visualization purposes, let me include this as well:

(* https://tpfto.wordpress.com/2019/04/17/faking-a-birds-colors-and-miscellanea/ *)
birdie[x_?NumericQ] /; 0 <= x <= 1 := 
Blend[{RGBColor[0.242, 0.15, 0.66], RGBColor[0.27, 0.2135, 0.83], 
       RGBColor[0.281, 0.298, 0.939], RGBColor[0.271, 0.385, 0.99], 
       RGBColor[0.202, 0.479, 0.99], RGBColor[0.172, 0.564, 0.94], 
       RGBColor[0.13, 0.639, 0.892], RGBColor[0.048, 0.7025, 0.827], 
       RGBColor[0.07, 0.746, 0.726], RGBColor[0.203, 0.78, 0.61], 
       RGBColor[0.36, 0.8, 0.463], RGBColor[0.58, 0.79, 0.29], 
       RGBColor[0.786, 0.757, 0.16], RGBColor[0.946, 0.7285, 0.217], 
       RGBColor[0.995, 0.79, 0.2], RGBColor[0.961, 0.891, 0.153], 
       RGBColor[0.977, 0.984, 0.08]}, x]

and then, we have

ArrayPlot[numGridB[7], ColorFunction -> birdie, Frame -> False]

ArrayPlot[numGridB[500], ColorFunction -> birdie, Frame -> False]

---

Edit 7/10/2022

The stuff above was written a little before I went to bed. Now with the benefit of getting sleep and seeing Ben's comment and chyanog's new answer, I came up with the following:

numGridB[n_Integer?Positive] := Module[{pos, sgn2, temp, zz},
   zz = {1, I} . Transpose[CoordinateBoundsArray[{{-1., 1.}, {1., -1.}},
                                                 {Into[n - 1], Into[n - 1]}], {3, 2, 1}];
   sgn2 = Sign[zz]^2;
   temp = SparseArray[BitOr[1 - UnitStep[Abs[zz] - Im[sgn2 (1 + 0.2 sgn2^3)]],
                            UnitStep[Abs[Re[zz]] - 1], UnitStep[Abs[Im[zz]] - 1]],
                      Automatic, 1];
   pos = temp["NonzeroPositions"];
   SparseArray[SortBy[pos, Last] -> Range[Length[pos]], {n, n}, 0]]

How this version works is an exercise I will leave for the interested.