Bezdepbonusu7Программирование: статьи и ответы
Меню

Вопрос

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

Автор вопроса: Esben EickhardtОтветов: 2

Мне нужно рассчитать индекс структурного сходства (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

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

  • Кстати, вы можете взглянуть на outer и expand.grid. Их можно использовать для облегчения ваших вложенных циклов for. lmo 19.04.2016

rblurmatlabgaussian2d

Ответы

Ответ 1 · рейтинг 2 · lmo

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

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

Ответ 2 · рейтинг 1 · rayryeng

В 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

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

Перевод вопроса и ответов сообщества Stack Overflow, лицензия CC BY-SA. Оригинал на Stack Overflow.

Похожие материалы