Construct a Gaussian Matrix

Question

Gaussian blur is a method used for blurring images smoothly. It involves creating a matrix which will be used by convolving it with the pixels of an image. In this challenge, your task is to construct that matrix used in Gaussian blur. You will take an input r which will be the radius of the blur and an input σ which will be the standard deviation in order to construct a matrix with dimensions (2 r + 1 × 2 r + 1). Each value in that matrix will have an (x, y) value that depends on its absolute distance in each direction from the center and will be used to compute G(x, y) where the formula G is:

$$G(x,y)=\frac{1}{2 \cdot \pi \cdot \sigma ^ 2} e ^ {-\frac{x^2+y^2}{2\cdot \sigma ^2}} $$

For example, if r = 2, we want to generate a 5 x 5 matrix. First, the matrix of (x, y) values is

(2, 2) (1, 2) (0, 2) (1, 2) (2, 2)
(2, 1) (1, 1) (0, 1) (1, 1) (2, 1)
(2, 0) (1, 0) (0, 0) (1, 0) (2, 0)
(2, 1) (1, 1) (0, 1) (1, 1) (2, 1)
(2, 2) (1, 2) (0, 2) (1, 2) (2, 2)

Then, let σ = 1.5 and apply G to each (x, y)

0.0119552 0.0232856 0.0290802 0.0232856 0.0119552
0.0232856 0.0453542 0.0566406 0.0453542 0.0232856
0.0290802 0.0566406 0.0707355 0.0566406 0.0290802
0.0232856 0.0453542 0.0566406 0.0453542 0.0232856
0.0119552 0.0232856 0.0290802 0.0232856 0.0119552

Normally in image blurring, this matrix would be normalized by taking the sum of all the values in that matrix and dividing by it. For this challenge, that is not needed and the raw values calculated by the formula is what the output should be.

Rules

• This is code-golf so the shortest code wins.
• The input r will be a nonnegative integer and σ will be a positive real number.
• The output must represent a matrix. It can be formatted as a 2d array, a string representing a 2d array, or something similar.
• Floating-point inaccuracies will not be counted against you.

Test Cases

(r, σ) = (0, 0.25)
2.54648

(1, 7)
0.00318244 0.00321509 0.00318244
0.00321509 0.00324806 0.00321509
0.00318244 0.00321509 0.00318244

(3, 2.5)
0.00603332 0.00900065 0.0114421  0.012395 0.0114421 0.00900065 0.00603332
0.00900065  0.0134274 0.0170696 0.0184912 0.0170696  0.0134274 0.00900065
 0.0114421  0.0170696 0.0216997  0.023507 0.0216997  0.0170696  0.0114421
  0.012395  0.0184912  0.023507 0.0254648  0.023507  0.0184912   0.012395
 0.0114421  0.0170696 0.0216997  0.023507 0.0216997  0.0170696  0.0114421
0.00900065  0.0134274 0.0170696 0.0184912 0.0170696  0.0134274 0.00900065
0.00603332 0.00900065 0.0114421  0.012395 0.0114421 0.00900065 0.00603332

(4, 3.33)
0.00339074 0.00464913 0.00582484 0.00666854 0.00697611 0.00666854 0.00582484 0.00464913 0.00339074
0.00464913 0.00637454 0.00798657  0.0091434 0.00956511  0.0091434 0.00798657 0.00637454 0.00464913
0.00582484 0.00798657  0.0100063  0.0114556   0.011984  0.0114556  0.0100063 0.00798657 0.00582484
0.00666854  0.0091434  0.0114556   0.013115  0.0137198   0.013115  0.0114556  0.0091434 0.00666854
0.00697611 0.00956511   0.011984  0.0137198  0.0143526  0.0137198   0.011984 0.00956511 0.00697611
0.00666854  0.0091434  0.0114556   0.013115  0.0137198   0.013115  0.0114556  0.0091434 0.00666854
0.00582484 0.00798657  0.0100063  0.0114556   0.011984  0.0114556  0.0100063 0.00798657 0.00582484
0.00464913 0.00637454 0.00798657  0.0091434 0.00956511  0.0091434 0.00798657 0.00637454 0.00464913
0.00339074 0.00464913 0.00582484 0.00666854 0.00697611 0.00666854 0.00582484 0.00464913 0.00339074

Answer 1

Mathematica, 60 54 50 bytes

Thanks @GregMartin for 4 bytes!

Array[s=2#2^2;E^(-{##}.{##}/s)/π/s&,1+2{#,#},-#]&

Takes r and sigma as input, returns the matrix (exact numbers).

Built-in version (58 bytes)

GaussianMatrix[{##},Standardized->1<0,Method->"Gaussian"]&

Of course, Mathematica has a built-in for this, too, but it's too long.

Answer 2

MATL, 20 bytes

_G&:U&+iUE/_Ze5MYP*/

Try it online!

Explanation

_     % Take input r implicitly. Negate
G     % Push r again
&:    % Binary range: [-r -r+1 ... r]
U     % Square, elementwise
&+    % Matrix of all pairwise additions
i     % Take input σ
U     % Square
E     % Multiply by 2. Gives 2σ^2
/     % Divide
_     % Negate
Ze    % Exponential
5M    % Push 2σ^2 again
YP    % Push pi
*     % Multiply
/     % Divide. Display implicitly

Answer 3

Octave, 45 bytes

@(r,s)exp((x=-(-r:r).^2/2/s^2)+x')/2/s^2/pi

Answer 4

Octave, 49 bytes

@(r,q)1./((Q=2*q^2)*pi*e.^(((x=(-r:r).^2)+x')/Q))

Try It Online!

Answer 5

Python, 88 bytes

lambda r,s:[[.5/3.14/s/s/2.72**((x*x+y*y)/2/s/s)for x in range(-r,r+1)]for y in range(-r,r+1)]

Uses the rule where you may hardcode 3.14 and 2.72 at 1 byte cost each.

Answer 6

Perl 6, 71 bytes

->\r,\σ{map ->\y{map ->\x{exp((x*x+y*y)/-2/σ/σ)/2/pi/σ/σ},-r..r},-r..r}

Technically this may be more than 71 bytes if encoded and saved to a file, but I couldn't resist naming the "sigma" input with an actual Greek sigma. It could be renamed to any plain-ASCII letter if one wanted.

Answer 7

SAS Macro Language, 296 bytes

Probably a much more efficient way to do this but it works :)

This code prints out the resulting dataset.

%macro G(r,s);%let l=%eval(2*&r+1);%let pi=%sysfunc(constant(pi));data w;array a[*] t1-t&l;%do i=-&r %to &r;%do j=-&r %to &r;%let t=%sysfunc(sum(&j,&r,1));a[&t]=%sysevalf(1/(2*&pi*&s**2)*%sysfunc(exp(-(%sysfunc(abs(&j))**2+%sysfunc(abs(&i))**2)/(2*&s**2))));%end;output;%end;proc print;run;%mend;

Answer 8

Haskell, 59 bytes

r#s|i<-[-r..r]=[[exp(-(x*x+y*y)/2/s/s)/2/pi/s/s|x<-i]|y<-i]

Usage example:

1#7

[[3.1824449424224664e-3,3.2150851187016326e-3,3.1824449424224664e-3],
 [3.2150851187016326e-3,3.2480600630999047e-3,3.2150851187016326e-3],
 [3.1824449424224664e-3,3.2150851187016326e-3,3.1824449424224664e-3]]

Answer 9

Python 2.7, 167 bytes

A very simple solution:

from __future__ import division;from math import*;r,s=input();s*=2*s;R=range(-r,r+1);print"\n".join("\t".join(str(pow(e,-(x*x+y*y)/s)/(pi*s))[:9]for x in R)for y in R)

Try it here!

Ungolfed:

from __future__ import division
from math import *
r,s = input()                         # Take input
s *= 2*s                              # Set s = 2*s^2; simplifies the expression
R = range(-r,r+1)                     # Range object; used twice

                                   # G(x,y)             # Stripped
print "\n".join("\t".join(str(pow(e,-(x*x + y*y)/s)/(pi*s))[:9] for x in R) for y in R)

Answer 10

JavaScript (Node.js), 260 229 212 208 bytes

(r,σ)=>{t=[];z=[];σ*=σ;b=r*-1;y=r;a=2*r++;while(y!=b-1){x=b;while(r>x){z.push(x,y);x++}y--}for(i=0;i<z.length;i+=2){x=z[i];y=z[i+1];t.push(1/(2*Math.PI*σ)*Math.pow(Math.E,(-1*x**2+y**2)/(2*σ)))}return t}

Attempt This Online!

How?

(r,σ)=>{t=[];z=[];σ*=σ;b=r*-1;y=r;a=2*r++;while(y!=b-1){x=b;while(r>x){z.push(x,y);x++}y--}for(i=0;i<z.length;i+=2){x=z[i];y=z[i+1];t.push(1/(2*Math.PI*σ)*Math.pow(Math.E,(-1*x**2+y**2)/(2*σ)))}return t} : main function
        t=[];z=[]; : Define empty arrays
                  σ*=σ; : set σ to σ squared 
                       a=2*r+1;b=r*-1;y=r; : set variables for looping
                                          while(y!=b-1){ : open while loop for y's
                                                        x=b; : reset x value
                                                            while(r+1>x){ : open while loop for x's
                                                                         z.push(x,y); : push xy values to array z
                                                                                     x++} : increment x and close loop 
                                                                                         y--} : decrement y and close loop
                                                                                             for(i=0;i<z.length;i+=2){ : open for loop based off length of array z
                                                                                                                      x=s(z[i]);y=s(z[i+1]); : retrieve x and y values from z
                                                                                                                                            t.push(....)} : push results to array t +closes for loop 
                                                                                                                                                   1/(2*Math.PI*σ)*Math.pow(Math.E,(-1*x**2+y**2)/(2*σ)) : Gaussian function
                                                                                                                                                                                                          return t} : returns the results and closes the main function

Answer 11

J, 39 bytes

{{(^-a%~1#.*:>{@(,&<)~i:x)%o.a=.+:*:y}}

Attempt This Online!

{{(^-a%~1#.*:>{@(,&<)~i:x)%o.a=.+:*:y}}
{{                                   }} NB. direct definition
                           o.a=.+:*:y   NB. π * (a = 2σ²)
  (^-a%~1#.*:>{@(,&<)~i:x)              NB. computes the eʸ portion
                      i:x               NB. range [-n..n]
              {@(,&<)~                  NB. cartesian product, creates a grid of boxed pairs
             >                          NB. unbox
           *:                           NB. square the pairs
        1#.                             NB. sum¹ the pairs
     a%~                                NB. divide each sum by a, x f~ y -> y f x
    -                                   NB. negate
   ^                                    NB. eʸ
                          %             NB. division

See common use 3 for sum¹. The DD allows a non-tacit definition, which helps implement the function directly, but I am working on refactoring this.