Алгоритм соотнесения графиков с разными шкалами

[aitap]

Добрый вечер!

Имеется два графика, у которых различаются шкалы по оси абсцисс и ординат. Один из них - модельный спектр с известными, правильными длинами волн ($ x_1 = \{ \lambda_i, i = 1 ... N \} $, $ y_1 = \{ y^{mod}(\lambda_i), i = 1 ... N \} $); другой - экспериментальный спектр, снятый с установки, и по оси абсцисс там - номер пикселя на камере ($ x_2 = \{ i, i = 1...N \} $, $ y_2 = \{ y^{exp}_i, i = 1...N \} $). Задача: подобрать такие коэффициенты полинома $ \lambda(i) = a_0 + i a_1 + i^2 a_2 $, чтобы спектры как можно лучше "сходились" (см. ниже). Приблизительно известны начальные значения $ a_0 $ и $ a_1 $, а $ a_2 $ на первых порах можно и вовсе считать нулём, поскольку он отвечает за аберрации спектрографа и является очень маленьким ($ \frac{a_2}{a_1} \approx 10^{-6} $) по сравнению с остальными.

Проблема 1: оба спектра известны в дискретных, в общем случае никогда не совпадающих точках (для модельного спектра это чисто техническая проблема, связанная с тем, что смена сетки выполняется гораздо дольше, чем расчёт спектра). Решение: заранее возьмём модельный спектр с запасом по диапазону длин волн и числу точек, чтобы экспериментальный скорее всего оказался между его границами. Рассчитаем длины волн для экспериментального спектра, исходя из приблизительно известных начальных значений. Объединим обе шкалы и линейно интерполируем оба спектра в точках, в которых значение не известно. Имея общую шкалу, попробуем максимизировать коэффициент корреляции $ \rho_{Pearson}(y_1, y_2) $ по $ a_0, a_1, a_2 $.

Проблема 2: у этой задачи есть несколько очень узких локальных оптимумов, а глобальный оптимум с правильным ответом почти невозможно найти: (картинки больше $ 1000 \times 1000 $ пикселей)

• протабулированный коэффициент корреляции модельного и экспериментального спектров на сетке $ 50 \times 50 $ при $ a_2 = 0 $
• спектры в единственной точке с $ \rho > 0,8 $

"В лоб" это всё решается с трудом: известные мне алгоритмы оптимизации застревают в локальных минимумах и/или норовят сбежать в область, где спектры почти не перекрываются. Табулировать коэффициент корреляции на достаточно мелкой сетке, чтобы потом брать максимум и назначать очень жёсткие ограничения для оптимизатора, долго и ненадёжно.

При склеивании панорам из фотографий сначала находят "особенности", а потом соотносят их друг с другом, но в моих одномерных графиках всего 1 канал и одно измерение, а не двумерная матрица RGB. Плюс, все пики похожи друг на друга.

Как ещё можно найти коэффициенты $ a_0, a_1, a_2 $?

[SVD-d]

Так как тема богатая, на всякий случай сразу уточню два нюанса.
1) Почему нельзя заранее перед моделированием задать такие пределы и такой масштаб, чтобы координаты точно совпадали хотя бы по оси Ox?
2) Почему нельзя зафиксировать камеру так, чтобы на один и тот же пиксель попадала всегда одна и та же часть спектра? Тогда калибровку надо будет провести всего один раз.

[Pavia]

aitap в сообщении #1358544 писал(а):

Плюс, все пики похожи друг на друга.

На изоброжение тоже много похожих писелей. Во вторых 3 канала тут непричём.
Между пиками у вас бесполезный сигнал и его много. Поэтому ваши алгоритмы и стремятся убежать от пиков. Бесполезного сигнала у вас больше чем полезного сигнала, вот они и утаскивают максимум. И фильтеры тут не помогут.
Надо как и с картинками через особые точки.
Вы должны при помощи порога отсечь бесполезный сигнал от пиков. Вычислить координаты пиков.
И вот с этими данными и должны работать. Данных станет меньше надёжность повысится.
Как сравнивать?
Зная координаты пиков. Варьируете Ваши переменные считаете сколько пиков совпало с эталоном. Совпадение проверять с допуском. Где совпадений больше там и оптимум.
А далее если нужно можно их уточнить.

[aitap]

Большое спасибо за ответы!

SVD-d в сообщении #1358758 писал(а):

1) Почему нельзя заранее перед моделированием задать такие пределы и такой масштаб, чтобы координаты точно совпадали хотя бы по оси Ox?
2) Почему нельзя зафиксировать камеру так, чтобы на один и тот же пиксель попадала всегда одна и та же часть спектра? Тогда калибровку надо будет провести всего один раз.

Проблема в том, что дифракционную решётку спектрографа можно поворачивать при помощи специальной ручки и что нам постоянно нужно снимать разные диапазоны длин волн. Конечно, рядом с ручкой есть индикатор, по показаниям которого можно приблизительно оценить, какая длина волны сейчас находится в центре изображения, и теоретические соображения, из которых рассчитывается дисперсия при данном угле поворота решётки, но необходимая точность по длинам волн достигается только путём снятия (при фиксированном положении ручки) "хороших" спектров (с тонкими, не сдвинутыми и во всех отношениях хорошо известными линиями) и соотнесения их координат с табличными длинами волн.

Я мог бы получить начальные значения для центральной длины волны и дисперсии, исходя из информации о данном спектрографе (тогда локальной оптимизации, может быть, и хватило бы), но хочу научиться решать более общую задачу.

Pavia в сообщении #1358815 писал(а):

Надо как и с картинками через особые точки.
Вы должны при помощи порога отсечь бесполезный сигнал от пиков. Вычислить координаты пиков. <...>
Зная координаты пиков. Варьируете Ваши переменные считаете сколько пиков совпало с эталоном. Совпадение проверять с допуском. Где совпадений больше там и оптимум.

Поясните, пожалуйста, как Вы предлагаете проверять совпадение координат пиков? Взять векторы координат пиков $ \bar{x} = (x_1, x_2, ... x_n ) $ и $ \bar{\lambda} = ( \lambda_1, \lambda_2, ... \lambda_m ) $, а потом найти такое $k$, для которого максимален коэффициент корреляции $r(\bar{x}, \bar{\lambda}_k)$, где $(\bar{\lambda_k})_i = \lambda_{i+k}$? (А то вдруг модельный спектр запросили более широким, чем нужно, и пиков там оказалось больше, чем в экспериментальном?)

[Pavia]

В машинной графике алгоритм сопоставления облака точек называется ICP.
У вас есть эталонный и измеренный спектр. Накладываешь один на другой. Для каждого пика в эталонном спектр находишь наиболее близкий к измеренному. Просто перебором.
$\min\limits_{i,j}(x_i - \lambda_j)$

Цитата:

то вдруг модельный спектр запросили более широким, чем нужно, и пиков там оказалось больше, чем в экспериментальном

Если пики между эталоном и измеренным значениям стоят далеко, то их просто откидываете. А имеющие минимальные подсчитываете.

$count(\min\limits_{i,j}(x_i - \lambda_j))$

Затем найденные расстояния между измеренными и эталонными пиками можете просуммировать - получите некоторую функция для оценки ошибки.

$e=\sum\limits_{p(i)}(\min\limits_{i,j}(x_i - \lambda_j))$

У нас есть масштаб a1 и смещение a0 они так же перебираются с некоторым разумным шагом. И уже ищем минимум для функции ошибки.

$\min\limits_{a0,a1}(\sum\limits_{p(i)} \min\limits_{i,j}(a1*x_i+a0 - \lambda_j))$

Или максимизируете, число совпадений.

$\max\limits_{a0,a1}(count(\min\limits_{i,j}(a1*x_i+a0 - \lambda_j))$

Далее когда масштаб примерно известен совпадающие пары пиков известны. Можно через МНК подобрать оптимум для коэффициентов a1,a0.

[aitap]

Pavia, большое спасибо за наводку на iterative closest point!

Я задал функцию потерь следующим образом:

• $ d_{sum} \leftarrow 0 $
• найдя пики в обоих графиках, считаю все попарные расстояния $ d_{i, j} = | a_0 + x_i a_1 - \lambda_j | $
• в цикле, пока в меньшем из двух множеств найденных пиков не закончатся пики:
  • нахожу $ (i_m, j_m) = \textrm{argmin}_{i,j} d_{i, j} $ и прибавляю $ d_{i_m, j_m} $ к $ d_{sum} $
    
  • всю строку и весь столбец матрицы исключаю из дальнейшего рассмотрения: $ d_{i_m, j} \leftarrow \mathtt{HUGE\_VAL} \, \forall j $; $ d_{i, j_m} \leftarrow \mathtt{HUGE\_VAL}\, \forall i $
    
• $ d_{sum} $ - значение функции потерь

В простых случаях, когда модельный и экспериментальный спектры оба содержат очень характеристичную по длинам волн последовательность пиков, функция очень хорошо минимизируется: 1 2.

Если пиков много, и они все расположены через примерно одинаковые промежутки длин волн (или, хуже, в одном спектре нашлось много больше пиков, чем в другом, но это явно не та ситуация, на которую я должен рассчитывать), ситуация становится хуже: минимум описанной выше функции перестаёт соответствовать тому, что я считаю решением, и соотносит друг с другом совершенно не относящиеся друг к другу пики. Я пробовал вводить штрафы на интенсивность пика:

• $ d_{i, j} = | a_0 + x_i a_1 - \lambda_j | \cdot y_i $ (или $I_j$)
• $ d_{i, j} = | a_0 + x_i a_1 - \lambda_j | \cdot |y_i - I_j| $ ($y_i$ и $I_j$ оба нормированы на $[0; 1]$)
• $ d_{i, j} = | a_0 + x_i a_1 - \lambda_j | \cdot \frac{\max(y_i, I_j)}{\min(y_i, I_j)} $

но особого успеха в этих сложных случаях не добился. Буду продолжать экспериментировать.

[aitap]

Долго не доходили руки попробовать следующий способ учесть интенсивность.

У обоих спектров шкалирую абсциссы и ординаты $ \frac{x_i - \min_i x_i}{\max_i x_i - \min_i x_i} \in [0; 1]$ и запускаю ICP с функцией потерь - суммой Евклидовых расстояний $ d_{ij} = \sqrt{(a_0 + a_1 x_i - \lambda_j)^2 + (b_0 + b_1 y_i - I_j)^2} $. При этом интуитивно получаются очень понятные ограничения на все параметры; кроме того, $ b_0 $ можно считать нулём, т.к. в обоих спектрах обязательно есть участок фона.

При этом спектры начинают соотноситься гораздо лучше, но несколько точек всё равно выпадают. Вероятно, мне придётся после выполнения OLS смотреть на остатки и удалять выбросы. Кроме того, может быть, мне стоит попробовать другой алгоритм поиска пиков (сейчас я делаю maximum overlap discrete wavelet transform и считаю пиками точки, в которых detail coefficient пересекает 0, а интенсивность scale coefficient выше порога).

Смеха ради попробовал запихнуть вообще все точки в ICP, но результата не дождался. И поделом.