在 R 中实现 Matlab fspecial-function
Implementing the Matlab fspecial-function in R
我需要计算R中的结构相似度（SSIM）指标，只能找到Matlab中实现的方法。除了两个 Matlab 方法 "fspecial" 和 "filter2".
之外，在 R 中重写该方法似乎非常简单
fspecial returns 标准方差为 1.5 的 11x11 矩阵中的二维高斯分布：
h = fspecial('gaussian', 11, 1.5)
所以我在网上找到了一些帮助，实现了一个应该在 R 中做同样事情的函数：
gaussian2D <- function(amplitude) {
# Defining limits of grid
x_min <- 1
x_max <- 11
y_min <- x_min
y_max <- x_max
# Setting parameters of the two-dimensional Gaussian function
# The distribution is centred in [6,6]
x_zero <- 6
y_zero <- 6
# Setting the spread of the filter
sigma_x <- 1.5
sigma_y <- sigma_x
# Running through all x and y combinations applying the 2d-gaussian equation
df <- NULL
for (x_val in c(x_min:x_max)){
for (y_val in c(y_min:y_max)){
cell_value <- amplitude*exp(-( (((x_val-x_zero)^2)/(2*(sigma_x^2))) + (((y_val-y_zero)^2)/(2*(sigma_y^2))) ))
df = rbind(df,data.frame(x_val,y_val, cell_value))
}
}
# Axis labels
x_axis <- c(x_min:x_max)
y_axis <- c(y_min:y_max)
# Populating matrix
gauss_matrix <- matrix(df[,3], nrow = 11, ncol = 11, dimnames = list(x_axis, y_axis))
return(gauss_matrix)
}
h2 = gaussian2D(1)
然而，奇怪的是，当我 运行 这两种方法时，我没有得到相同的结果，而是得到了一个按 14.13 缩放的矩阵：
h2/h
1 2 3 4 5 6 7 8 9 10 11
1 14.13137 14.13185 14.13201 14.13238 14.13187 14.13189 14.13187 14.13238 14.13201 14.13185 14.13137
2 14.13185 14.13186 14.13189 14.13175 14.13164 14.13154 14.13164 14.13175 14.13189 14.13186 14.13185
3 14.13201 14.13189 14.13135 14.13172 14.13176 14.13187 14.13176 14.13172 14.13135 14.13189 14.13201
4 14.13238 14.13175 14.13172 14.13155 14.13209 14.13194 14.13209 14.13155 14.13172 14.13175 14.13238
5 14.13187 14.13164 14.13176 14.13209 14.13194 14.13182 14.13194 14.13209 14.13176 14.13164 14.13187
6 14.13189 14.13154 14.13187 14.13194 14.13182 14.13188 14.13182 14.13194 14.13187 14.13154 14.13189
7 14.13187 14.13164 14.13176 14.13209 14.13194 14.13182 14.13194 14.13209 14.13176 14.13164 14.13187
8 14.13238 14.13175 14.13172 14.13155 14.13209 14.13194 14.13209 14.13155 14.13172 14.13175 14.13238
9 14.13201 14.13189 14.13135 14.13172 14.13176 14.13187 14.13176 14.13172 14.13135 14.13189 14.13201
10 14.13185 14.13186 14.13189 14.13175 14.13164 14.13154 14.13164 14.13175 14.13189 14.13186 14.13185
11 14.13137 14.13185 14.13201 14.13238 14.13187 14.13189 14.13187 14.13238 14.13201 14.13185 14.13137
有没有人对我做错了什么提出建议？
您缺少调整项：1/(2*pi*sigma1*sigma2)
注意 2*pi*1.5*1.5 = 14.13717
在 MATLAB 中，fspecial 还 规范化 内核，使内核中所有元素的总和等于 1。这是因为在执行卷积时这个内核，它避免产生超出与您尝试过滤的信号相关联的数据类型的动态范围的任何输出值。
这也避免了必须使用 @imo 之前声明的任何调整项。很简单，您没有在代码中规范化内核。因此，在你的循环中有一个额外的求和项，它对每个高斯项求和，然后最终矩阵应该将每个条目除以这个数量：
df <- NULL
s <- 0 # Added
for (x_val in c(x_min:x_max)){
for (y_val in c(y_min:y_max)){
cell_value <- amplitude*exp(-( (((x_val-x_zero)^2)/(2*(sigma_x^2))) + (((y_val-y_zero)^2)/(2*(sigma_y^2))) ))
df = rbind(df,data.frame(x_val,y_val, cell_value))
s <- s + cell_value # Added
}
}
最后当你 return 矩阵时：
return(gauss_matrix / s)
仔细检查一下，在 MATLAB 中，这是您调用 fspecial:
时得到的结果
>> format long g
>> h = fspecial('gaussian', 11, 1.5)
h =
Columns 1 through 3
1.05756559815326e-06 7.8144115330536e-06 3.70224770827489e-05
7.8144115330536e-06 5.77411251978637e-05 0.000273561160085806
3.70224770827489e-05 0.000273561160085806 0.0012960555938432
0.000112464355116679 0.000831005429087199 0.00393706926284678
0.000219050652866017 0.00161857756253439 0.00766836382523672
0.000273561160085806 0.00202135875836257 0.00957662749024029
0.000219050652866017 0.00161857756253439 0.00766836382523672
0.000112464355116679 0.000831005429087199 0.00393706926284678
3.70224770827489e-05 0.000273561160085806 0.0012960555938432
7.8144115330536e-06 5.77411251978637e-05 0.000273561160085806
1.05756559815326e-06 7.8144115330536e-06 3.70224770827489e-05
Columns 4 through 6
0.000112464355116679 0.000219050652866017 0.000273561160085806
0.000831005429087199 0.00161857756253439 0.00202135875836257
0.00393706926284678 0.00766836382523672 0.00957662749024029
0.011959760410037 0.0232944324734871 0.0290912256485504
0.0232944324734871 0.0453713590956603 0.0566619704916846
0.0290912256485504 0.0566619704916846 0.070762237763947
0.0232944324734871 0.0453713590956603 0.0566619704916846
0.011959760410037 0.0232944324734871 0.0290912256485504
0.00393706926284678 0.00766836382523672 0.00957662749024029
0.000831005429087199 0.00161857756253439 0.00202135875836257
0.000112464355116679 0.000219050652866017 0.000273561160085806
Columns 7 through 9
0.000219050652866017 0.000112464355116679 3.70224770827489e-05
0.00161857756253439 0.000831005429087199 0.000273561160085806
0.00766836382523672 0.00393706926284678 0.0012960555938432
0.0232944324734871 0.011959760410037 0.00393706926284678
0.0453713590956603 0.0232944324734871 0.00766836382523672
0.0566619704916846 0.0290912256485504 0.00957662749024029
0.0453713590956603 0.0232944324734871 0.00766836382523672
0.0232944324734871 0.011959760410037 0.00393706926284678
0.00766836382523672 0.00393706926284678 0.0012960555938432
0.00161857756253439 0.000831005429087199 0.000273561160085806
0.000219050652866017 0.000112464355116679 3.70224770827489e-05
Columns 10 through 11
7.8144115330536e-06 1.05756559815326e-06
5.77411251978637e-05 7.8144115330536e-06
0.000273561160085806 3.70224770827489e-05
0.000831005429087199 0.000112464355116679
0.00161857756253439 0.000219050652866017
0.00202135875836257 0.000273561160085806
0.00161857756253439 0.000219050652866017
0.000831005429087199 0.000112464355116679
0.000273561160085806 3.70224770827489e-05
5.77411251978637e-05 7.8144115330536e-06
7.8144115330536e-06 1.05756559815326e-06
... 最后在 R 中：
> h2
1 2 3 4 5 6 7
1 1.057566e-06 7.814412e-06 3.702248e-05 0.0001124644 0.0002190507 0.0002735612 0.0002190507
2 7.814412e-06 5.774113e-05 2.735612e-04 0.0008310054 0.0016185776 0.0020213588 0.0016185776
3 3.702248e-05 2.735612e-04 1.296056e-03 0.0039370693 0.0076683638 0.0095766275 0.0076683638
4 1.124644e-04 8.310054e-04 3.937069e-03 0.0119597604 0.0232944325 0.0290912256 0.0232944325
5 2.190507e-04 1.618578e-03 7.668364e-03 0.0232944325 0.0453713591 0.0566619705 0.0453713591
6 2.735612e-04 2.021359e-03 9.576627e-03 0.0290912256 0.0566619705 0.0707622378 0.0566619705
7 2.190507e-04 1.618578e-03 7.668364e-03 0.0232944325 0.0453713591 0.0566619705 0.0453713591
8 1.124644e-04 8.310054e-04 3.937069e-03 0.0119597604 0.0232944325 0.0290912256 0.0232944325
9 3.702248e-05 2.735612e-04 1.296056e-03 0.0039370693 0.0076683638 0.0095766275 0.0076683638
10 7.814412e-06 5.774113e-05 2.735612e-04 0.0008310054 0.0016185776 0.0020213588 0.0016185776
11 1.057566e-06 7.814412e-06 3.702248e-05 0.0001124644 0.0002190507 0.0002735612 0.0002190507
8 9 10 11
1 0.0001124644 3.702248e-05 7.814412e-06 1.057566e-06
2 0.0008310054 2.735612e-04 5.774113e-05 7.814412e-06
3 0.0039370693 1.296056e-03 2.735612e-04 3.702248e-05
4 0.0119597604 3.937069e-03 8.310054e-04 1.124644e-04
5 0.0232944325 7.668364e-03 1.618578e-03 2.190507e-04
6 0.0290912256 9.576627e-03 2.021359e-03 2.735612e-04
7 0.0232944325 7.668364e-03 1.618578e-03 2.190507e-04
8 0.0119597604 3.937069e-03 8.310054e-04 1.124644e-04
9 0.0039370693 1.296056e-03 2.735612e-04 3.702248e-05
10 0.0008310054 2.735612e-04 5.774113e-05 7.814412e-06
11 0.0001124644 3.702248e-05 7.814412e-06 1.057566e-06
...看起来很适合我！
我需要计算R中的结构相似度（SSIM）指标，只能找到Matlab中实现的方法。除了两个 Matlab 方法 "fspecial" 和 "filter2".
之外，在 R 中重写该方法似乎非常简单fspecial returns 标准方差为 1.5 的 11x11 矩阵中的二维高斯分布：
h = fspecial('gaussian', 11, 1.5)
所以我在网上找到了一些帮助，实现了一个应该在 R 中做同样事情的函数：
gaussian2D <- function(amplitude) {
# Defining limits of grid
x_min <- 1
x_max <- 11
y_min <- x_min
y_max <- x_max
# Setting parameters of the two-dimensional Gaussian function
# The distribution is centred in [6,6]
x_zero <- 6
y_zero <- 6
# Setting the spread of the filter
sigma_x <- 1.5
sigma_y <- sigma_x
# Running through all x and y combinations applying the 2d-gaussian equation
df <- NULL
for (x_val in c(x_min:x_max)){
for (y_val in c(y_min:y_max)){
cell_value <- amplitude*exp(-( (((x_val-x_zero)^2)/(2*(sigma_x^2))) + (((y_val-y_zero)^2)/(2*(sigma_y^2))) ))
df = rbind(df,data.frame(x_val,y_val, cell_value))
}
}
# Axis labels
x_axis <- c(x_min:x_max)
y_axis <- c(y_min:y_max)
# Populating matrix
gauss_matrix <- matrix(df[,3], nrow = 11, ncol = 11, dimnames = list(x_axis, y_axis))
return(gauss_matrix)
}
h2 = gaussian2D(1)
然而，奇怪的是，当我 运行 这两种方法时，我没有得到相同的结果，而是得到了一个按 14.13 缩放的矩阵：
h2/h
1 2 3 4 5 6 7 8 9 10 11
1 14.13137 14.13185 14.13201 14.13238 14.13187 14.13189 14.13187 14.13238 14.13201 14.13185 14.13137
2 14.13185 14.13186 14.13189 14.13175 14.13164 14.13154 14.13164 14.13175 14.13189 14.13186 14.13185
3 14.13201 14.13189 14.13135 14.13172 14.13176 14.13187 14.13176 14.13172 14.13135 14.13189 14.13201
4 14.13238 14.13175 14.13172 14.13155 14.13209 14.13194 14.13209 14.13155 14.13172 14.13175 14.13238
5 14.13187 14.13164 14.13176 14.13209 14.13194 14.13182 14.13194 14.13209 14.13176 14.13164 14.13187
6 14.13189 14.13154 14.13187 14.13194 14.13182 14.13188 14.13182 14.13194 14.13187 14.13154 14.13189
7 14.13187 14.13164 14.13176 14.13209 14.13194 14.13182 14.13194 14.13209 14.13176 14.13164 14.13187
8 14.13238 14.13175 14.13172 14.13155 14.13209 14.13194 14.13209 14.13155 14.13172 14.13175 14.13238
9 14.13201 14.13189 14.13135 14.13172 14.13176 14.13187 14.13176 14.13172 14.13135 14.13189 14.13201
10 14.13185 14.13186 14.13189 14.13175 14.13164 14.13154 14.13164 14.13175 14.13189 14.13186 14.13185
11 14.13137 14.13185 14.13201 14.13238 14.13187 14.13189 14.13187 14.13238 14.13201 14.13185 14.13137
有没有人对我做错了什么提出建议？
您缺少调整项：1/(2*pi*sigma1*sigma2)
注意 2*pi*1.5*1.5 = 14.13717
在 MATLAB 中，fspecial 还 规范化 内核，使内核中所有元素的总和等于 1。这是因为在执行卷积时这个内核，它避免产生超出与您尝试过滤的信号相关联的数据类型的动态范围的任何输出值。
这也避免了必须使用 @imo 之前声明的任何调整项。很简单，您没有在代码中规范化内核。因此，在你的循环中有一个额外的求和项，它对每个高斯项求和，然后最终矩阵应该将每个条目除以这个数量：
df <- NULL
s <- 0 # Added
for (x_val in c(x_min:x_max)){
for (y_val in c(y_min:y_max)){
cell_value <- amplitude*exp(-( (((x_val-x_zero)^2)/(2*(sigma_x^2))) + (((y_val-y_zero)^2)/(2*(sigma_y^2))) ))
df = rbind(df,data.frame(x_val,y_val, cell_value))
s <- s + cell_value # Added
}
}
最后当你 return 矩阵时：
return(gauss_matrix / s)
仔细检查一下，在 MATLAB 中，这是您调用 fspecial:
>> format long g
>> h = fspecial('gaussian', 11, 1.5)
h =
Columns 1 through 3
1.05756559815326e-06 7.8144115330536e-06 3.70224770827489e-05
7.8144115330536e-06 5.77411251978637e-05 0.000273561160085806
3.70224770827489e-05 0.000273561160085806 0.0012960555938432
0.000112464355116679 0.000831005429087199 0.00393706926284678
0.000219050652866017 0.00161857756253439 0.00766836382523672
0.000273561160085806 0.00202135875836257 0.00957662749024029
0.000219050652866017 0.00161857756253439 0.00766836382523672
0.000112464355116679 0.000831005429087199 0.00393706926284678
3.70224770827489e-05 0.000273561160085806 0.0012960555938432
7.8144115330536e-06 5.77411251978637e-05 0.000273561160085806
1.05756559815326e-06 7.8144115330536e-06 3.70224770827489e-05
Columns 4 through 6
0.000112464355116679 0.000219050652866017 0.000273561160085806
0.000831005429087199 0.00161857756253439 0.00202135875836257
0.00393706926284678 0.00766836382523672 0.00957662749024029
0.011959760410037 0.0232944324734871 0.0290912256485504
0.0232944324734871 0.0453713590956603 0.0566619704916846
0.0290912256485504 0.0566619704916846 0.070762237763947
0.0232944324734871 0.0453713590956603 0.0566619704916846
0.011959760410037 0.0232944324734871 0.0290912256485504
0.00393706926284678 0.00766836382523672 0.00957662749024029
0.000831005429087199 0.00161857756253439 0.00202135875836257
0.000112464355116679 0.000219050652866017 0.000273561160085806
Columns 7 through 9
0.000219050652866017 0.000112464355116679 3.70224770827489e-05
0.00161857756253439 0.000831005429087199 0.000273561160085806
0.00766836382523672 0.00393706926284678 0.0012960555938432
0.0232944324734871 0.011959760410037 0.00393706926284678
0.0453713590956603 0.0232944324734871 0.00766836382523672
0.0566619704916846 0.0290912256485504 0.00957662749024029
0.0453713590956603 0.0232944324734871 0.00766836382523672
0.0232944324734871 0.011959760410037 0.00393706926284678
0.00766836382523672 0.00393706926284678 0.0012960555938432
0.00161857756253439 0.000831005429087199 0.000273561160085806
0.000219050652866017 0.000112464355116679 3.70224770827489e-05
Columns 10 through 11
7.8144115330536e-06 1.05756559815326e-06
5.77411251978637e-05 7.8144115330536e-06
0.000273561160085806 3.70224770827489e-05
0.000831005429087199 0.000112464355116679
0.00161857756253439 0.000219050652866017
0.00202135875836257 0.000273561160085806
0.00161857756253439 0.000219050652866017
0.000831005429087199 0.000112464355116679
0.000273561160085806 3.70224770827489e-05
5.77411251978637e-05 7.8144115330536e-06
7.8144115330536e-06 1.05756559815326e-06
... 最后在 R 中：
> h2
1 2 3 4 5 6 7
1 1.057566e-06 7.814412e-06 3.702248e-05 0.0001124644 0.0002190507 0.0002735612 0.0002190507
2 7.814412e-06 5.774113e-05 2.735612e-04 0.0008310054 0.0016185776 0.0020213588 0.0016185776
3 3.702248e-05 2.735612e-04 1.296056e-03 0.0039370693 0.0076683638 0.0095766275 0.0076683638
4 1.124644e-04 8.310054e-04 3.937069e-03 0.0119597604 0.0232944325 0.0290912256 0.0232944325
5 2.190507e-04 1.618578e-03 7.668364e-03 0.0232944325 0.0453713591 0.0566619705 0.0453713591
6 2.735612e-04 2.021359e-03 9.576627e-03 0.0290912256 0.0566619705 0.0707622378 0.0566619705
7 2.190507e-04 1.618578e-03 7.668364e-03 0.0232944325 0.0453713591 0.0566619705 0.0453713591
8 1.124644e-04 8.310054e-04 3.937069e-03 0.0119597604 0.0232944325 0.0290912256 0.0232944325
9 3.702248e-05 2.735612e-04 1.296056e-03 0.0039370693 0.0076683638 0.0095766275 0.0076683638
10 7.814412e-06 5.774113e-05 2.735612e-04 0.0008310054 0.0016185776 0.0020213588 0.0016185776
11 1.057566e-06 7.814412e-06 3.702248e-05 0.0001124644 0.0002190507 0.0002735612 0.0002190507
8 9 10 11
1 0.0001124644 3.702248e-05 7.814412e-06 1.057566e-06
2 0.0008310054 2.735612e-04 5.774113e-05 7.814412e-06
3 0.0039370693 1.296056e-03 2.735612e-04 3.702248e-05
4 0.0119597604 3.937069e-03 8.310054e-04 1.124644e-04
5 0.0232944325 7.668364e-03 1.618578e-03 2.190507e-04
6 0.0290912256 9.576627e-03 2.021359e-03 2.735612e-04
7 0.0232944325 7.668364e-03 1.618578e-03 2.190507e-04
8 0.0119597604 3.937069e-03 8.310054e-04 1.124644e-04
9 0.0039370693 1.296056e-03 2.735612e-04 3.702248e-05
10 0.0008310054 2.735612e-04 5.774113e-05 7.814412e-06
11 0.0001124644 3.702248e-05 7.814412e-06 1.057566e-06
...看起来很适合我！