Реализация специальной функции Matlab fspecial в R

Мне нужно рассчитать индекс структурного сходства (SSIM) в R, и я могу найти только метод, реализованный в Matlab. Переписывание метода в R кажется довольно простым, за исключением двух методов Matlab «fspecial» и «filter2».

fspecial возвращает двумерное распределение Гаусса в матрице 11x11 со стандартным делением 1,5:

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

У кого-нибудь есть предложение, что я сделал неправильно?


r blur matlab gaussian 2d
person Esben Eickhardt    schedule 19.04.2016    source источник
comment
Кстати, вы можете взглянуть на outer и expand.grid. Их можно использовать для облегчения ваших вложенных циклов for.   -  person lmo    schedule 19.04.2016

Ответы (2)


arrow_upward
2
arrow_downward

Вам не хватает корректирующего условия: 1/(2*pi*sigma1*sigma2)

Обратите внимание, что 2*пи*1,5*1,5 = 14,13717.

person lmo    schedule 19.04.2016
comment
Спасибо, мне было интересно, откуда взялась эта цифра! Я пробовал всевозможные манипуляции с пи, но это имеет смысл. - person Esben Eickhardt; 20.04.2016

arrow_upward
1
arrow_downward

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

...похоже, мне подходит!

person rayryeng    schedule 19.04.2016
comment
Спасибо, теперь я выдаю то, что хочу! - person Esben Eickhardt; 20.04.2016