# Construct a Gaussian Matrix

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

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 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

• How accurate do we need pi and e to be?
– xnor
Dec 29, 2016 at 0:04
• @xnor Good question. If your language allows, you may assume that those values are already stored in a variable or something similar. If not, you may use values to two decimal places making pi = 3.14 and e = 2.72 where you can count each of those values as a single byte. Of course, inaccuracies in the final answer will not be counted against you again. Dec 29, 2016 at 0:29
• Does the output have to be a decimal number, or could it be exact numbers with constants in them? Dec 29, 2016 at 3:37
• @JungHwanMin The exact numbers, such as those in Mathematica, are fine. Dec 29, 2016 at 3:40
• @miles I think it would be easier if you just mandated a certain precision (e.g. 3 decimal places).
– orlp
Dec 29, 2016 at 8:33

# Mathematica, 6054 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.

• You can replace -l by -# at the end (Array will thread that over both dimensions for you); that takes away the need to define l, saving 4 bytes. Dec 29, 2016 at 6:51

# 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


# Octave, 45 bytes

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


# Octave, 49 bytes

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


Try It Online!

# 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.

# 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.

# 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;


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]]


# 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)


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)

• from __future__ import division, really?
– orlp
Dec 29, 2016 at 8:10

# JavaScript (Node.js), 260229212 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


# 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.