# 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 '16 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. – miles Dec 29 '16 at 0:29
• Does the output have to be a decimal number, or could it be exact numbers with constants in them? – JungHwan Min Dec 29 '16 at 3:37
• @JungHwanMin The exact numbers, such as those in Mathematica, are fine. – miles Dec 29 '16 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 '16 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. – Greg Martin Dec 29 '16 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;


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


# 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 '16 at 8:10