Более быстрый способ выполнить точечную интерполяцию массива numpy? - PullRequest
Новая
       16

Более быстрый способ выполнить точечную интерполяцию массива numpy?

0 голосов
/ 08 апреля 2011

У меня есть трехмерный куб данных с двумя пространственными измерениями, а третий является многополосным спектром в каждой точке 2D-изображения. 

H[x, y, bands]

Учитывая длину волны (или номер полосы), я хотел бы извлечь 2D-изображение, соответствующее этой длине волны. Это будет просто срез массива, такой как H[:,:,bnd]. Точно так же, учитывая пространственное местоположение (i, j), спектр в этом местоположении равен H[i,j].

Я также хотел бы «сгладить» изображение спектрально, чтобы противостоять шуму при слабом освещении в спектрах. То есть для полосы bnd, я выбираю окно размером wind и подгоняю полином n-степени к спектру в этом окне. С polyfit и polyval я могу найти подходящее спектральное значение в этой точке для полосы bnd.

Теперь, если я хочу получить полное изображение bnd из установленного значения, тогда я должен выполнить это оконное выравнивание для каждого (i,j) изображения. Я также хочу изображение 2-й производной bnd, то есть значение 2-й производной подобранного спектра в каждой точке.

Бегая по точкам, я мог бы polyfit-polyval-polyder каждый из спектров x*y. Хотя это работает, это точечная операция. Есть ли какой-нибудь пито-нумпонский способ сделать это быстрее?

1 Ответ

1 голос
/ 08 апреля 2011

Если вы выполняете полиномиальную наименьших квадратов по точкам (x + dx [i], y [i]) для фиксированного набора dx, а затем вычисляете получившийся полином по x, результат будет (фиксированная) линейная комбинацияиз у [я].То же самое верно для производных многочлена.Так что вам просто нужна линейная комбинация срезов.Посмотрите «Фильтры Савицкого-Голея».

РЕДАКТИРОВАНИЕ, чтобы добавить краткий пример того, как работают фильтры SG.Я не проверил ни одной детали, и поэтому вам не следует полагаться на то, что она верна.

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

Итак, если f (x) ~ = ax ^ 2 + bx + c и f (-2), f (-1), f (0), f (1), f (2) = p, q, r, s, t, тогда мы хотим4a-2b + c ~ = p, a-b + c ~ = q и т. Д. Подгонка по методу наименьших квадратов означает минимизацию (4a-2b + cp) ^ 2 + (a-b + cq) ^ 2 + (cr) ^2 + (a + b + cs) ^ 2 + (4a + 2b + ct) ^ 2, что означает (принимая частные производные по a, b, c):

  • 4 (4a-2b+ cp) + (a-b + cq) + (a + b + cs) +4 (4a + 2b + ct) = 0
  • -2 (4a-2b + cp) - (a-b+ cq) + (a + b + cs) +2 (4a + 2b + ct) = 0
  • (4a-2b + cp) + (a-b + cq) + (cr) + (a+ b + cs) + (4a + 2b + ct) = 0

или, упрощенно,

  • 22a + 10c = 4p + q + s + 4t
  • 10b = -2p-q + s + 2t
  • 10a + 5c = p + q + r + s + t

т. А, b, c = pq/ 2-rs / 2 + t, (2 (tp) + (sq)) / 10, (p +q + r + s + t) / 5- (2p-q-2r-s + 2t).

И, конечно, c - это значение подогнанного многочлена в 0, и, следовательно, это сглаженное значение, которое мыхочу.Таким образом, для каждой пространственной позиции у нас есть вектор входных спектральных данных, из которого мы вычисляем сглаженные спектральные данные путем умножения на матрицу, строки которой (кроме первой и последней пары) выглядят как [0 ... 0 -9 /5 4/5 11/5 4/5 -9/5 0 ... 0], с центральной 11/5 на главной диагонали матрицы.

Таким образом, вы можете выполнить умножение матрицы для каждогопространственное положение;но так как это одинаковая матрица везде, вы можете сделать это одним вызовом tensordot.Поэтому, если S содержит только что описанную мной матрицу (э-э, подождите, нет, транспонирует только что описанной матрицы), а A - это ваш трехмерный куб данных, ваши спектрально-сглаженные данныекуб был бы numpy.tensordot(A,S).

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

...