Как извлечь только 3 собственных вектора изображения nxn в opencv?

Question

Я пытаюсь преобразовать RGB-изображение в оттенки серого, используя следующую статью.

Основной алгоритм, используемый в статье: алгоритм на основе нового PCA для преобразования изображений в оттенки серого

Однако, когда я пытаюсь извлечь собственные векторы из изображения, я получаю 500 собственных значений, а не 3, как требуется. Насколько я знаю, матрица NxN обычно дает N собственных векторов, но я не совсем уверен, что мне следует делать здесь, чтобы получить только 3 собственных вектора.

Любая помощь относительно того, что я должен сделать? Вот мой код до сих пор:

import numpy as np
import cv2

def pca_rgb2gray(img):
    """
    NOVEL PCA-BASED COLOR-TO-GRAY IMAGE CONVERSION
    Authors:
    -Ja-Won Seo
    -Seong Dae Kim
    2013 IEEE International Conference on Image Processing
    """
    I_re = cv2.resize(img, (500,500))
    Iycc = cv2.cvtColor(I_re, cv2.COLOR_BGR2YCrCb)
    Izycc = Iycc - Iycc.mean()
    eigvals = []
    eigvecs = []
    final_im = []
    for i in range(3):
        res = np.linalg.eig(Izycc[:,:,i])
        eigvals.append(res[0])
        eigvecs.append(res[1])
    eignorm = np.linalg.norm(eigvals)
    for i in range(3):
        eigvals[i]/=eignorm
        eigvecs[i]/=np.linalg.norm(eigvecs[i])
        temp = eigvals[i] * np.dot(eigvecs[i], Izycc[:,:,i])
        final_im.append(temp)
    final_im = final_im[0] + final_im[1] + final_im[2]
    return final_im
if __name__ == '__main__':
    img = cv2.imread('image.png')
    gray = pca_rgb2gray(img)

Answer 1

Фон

Когда Сео и Ким просят lambda_i, v_i <- PCA(Iycc), за i = 1, 2, 3, они хотят:

from numpy.linalg import eig
lambdas, vs = eig(np.dot(Izycc, Izycc.T))

для массива 3×N Izycc, То есть им нужны три собственных значения и собственные векторы ковариационной матрицы 3×3 Izycc массив 3×N (для вас N = 500*500).

Однако вы почти никогда не хотите вычислять ковариационную матрицу, а затем находите ее собственное разложение из-за численной нестабильности. Есть намного лучший способ получить то же самое lambdas, vs с использованием разложения по сингулярным числам (SVD) Izycc напрямую (см. этот ответ). Код ниже показывает, как это сделать.

Просто покажи мне код

Сначала загрузите http://cadik.posvete.cz/color_to_gray_evaluation/img/155_5572_jpg/155_5572_jpg.jpg и сохраните его как peppers.jpg,

Затем выполните следующее:

import numpy as np
import cv2
from numpy.linalg import svd, norm

# Read input image
Ibgr = cv2.imread('peppers.jpg')

# Convert to YCrCb
Iycc = cv2.cvtColor(Ibgr, cv2.COLOR_BGR2YCR_CB)

# Reshape the H by W by 3 array to a 3 by N array (N = W * H)
Izycc = Iycc.reshape([-1, 3]).T

# Remove mean along Y, Cr, and Cb *separately*!
Izycc = Izycc - Izycc.mean(1)[:, np.newaxis]

# Make sure we're dealing with zero-mean data here: the mean for Y, Cr, and Cb
# should separately be zero. Recall: Izycc is 3 by N array.
assert(np.allclose(np.mean(Izycc, 1), 0.0))

# Compute data array's SVD. Ignore the 3rd return value: unimportant.
(U, S) = svd(Izycc, full_matrices=False)[:2]

# Square the data's singular vectors to get the eigenvalues. Then, normalize
# the three eigenvalues to unit norm and finally, make a diagonal matrix out of
# them. N.B.: the scaling factor of `norm(S**2)` is, I believe, arbitrary: the
# rest of the algorithm doesn't really care if/how the eigenvalues are scaled,
# since we will rescale the grayscale values to [0, 255] anyway.
eigvals = np.diag(S**2 / norm(S**2))

# Eigenvectors are just the left-singular vectors.
eigvecs = U;

# Project the YCrCb data onto the principal components and reshape to W by H
# array.
Igray = np.dot(eigvecs.T, np.dot(eigvals, Izycc)).sum(0).reshape(Iycc.shape[:2])

# Rescale Igray to [0, 255]. This is a fancy way to do this.
from scipy.interpolate import interp1d
Igray = np.floor((interp1d([Igray.min(), Igray.max()],
                            [0.0, 256.0 - 1e-4]))(Igray))

# Make sure we don't accidentally produce a photographic negative (flip image
# intensities). N.B.: `norm` is often expensive; in real life, try to see if
# there's a more efficient way to do this.
if norm(Iycc[:,:,0] - Igray) > norm(Iycc[:,:,0] - (255.0 - Igray)):
    Igray = 255 - Igray

# Display result
if True:
    import pylab
    pylab.ion()
    pylab.imshow(Igray, cmap='gray')

# Save result
cv2.imwrite('peppers-gray.png', Igray.astype(np.uint8))

Это приводит к следующему изображению в градациях серого, которое, кажется, соответствует результату на рисунке 4 документа (хотя см. Предостережение внизу этого ответа!):

Ошибки в вашей реализации

Izycc = Iycc - Iycc.mean() НЕПРАВИЛЬНО. Iycc.mean() выравнивает изображение и вычисляет среднее. Ты хочешь Izycc так что канал Y, канал Cr и канал Cb имеют нулевое среднее значение. Вы могли бы сделать это в for dim in range(3) -Нет, но я делал это выше, с массивом вещания. У меня также есть утверждение выше, чтобы убедиться, что это условие выполняется. Трюк, в котором вы получаете собственное разложение ковариационной матрицы из SVD массива данных, требует каналов Y/Cr/Cb с нулевым средним.

np.linalg.eig(Izycc[:,:,i]) НЕПРАВИЛЬНО. Вклад этой статьи заключается в использовании основных компонентов для преобразования цвета в оттенки серого. Это означает, что вы должны комбинировать цвета. Обработка, которую вы выполняли выше, проходила по каналам - без сочетания цветов. Более того, было совершенно неправильно разлагать массив 500×500: ширина / высота массива не имеют значения, только пиксели. По этой причине я изменяю три канала входа в 3 раза и работаю с этой матрицей. Убедитесь, что вы понимаете, что происходит после конвертации BGR в YCrCb и до SVD.

Не столько ошибка, сколько предостережение: при звонке numpy.linalg.svd, full_matrices=False ключевое слово важно: это делает SVD "экономичного размера", вычисляя только три левых / правых сингулярных вектора и только три сингулярных значения. Полноразмерный SVD попытается создать массив N×N правых сингулярных векторов: с N = 114270 пикселей (изображение 293 на 390), массив N×N float64 будет N ** 2 * 8 / 1024 ** 3 или 97 гигабайт.

Конечная нота

Волшебство этого алгоритма действительно в одной строке из моего кода:

Igray = np.dot(eigvecs.T, np.dot(eigvals, Izycc)).sum(0) # .reshape...

Здесь математика самая толстая, поэтому давайте разберемся с ней.

• Izycc массив 3×N, строки которого имеют среднее значение нуля;
• eigvals является диагональным массивом 3×3, содержащим собственные значения ковариационной матрицы dot(Izycc, Izycc.T) (как упоминалось выше, вычисляется с помощью ярлыка, используя SVD Izycc),
• eigvecs является ортонормированной матрицей 3×3, столбцы которой являются собственными векторами, соответствующими этим собственным значениям этой ковариации.

Потому что это Numpy array а не matrix ES, мы должны использовать dot(x,y) для матрицы-матрицы-умножения, а затем мы используем sum и оба они затеняют линейную алгебру. Вы можете проверить сами, но приведенный выше расчет (до .reshape() вызов) эквивалентно

np.ones([1, 3]) · eigvecs.T · eigvals · Izycc = dot([[-0.79463857, -0.18382267,  0.11589724]], Izycc)

где · истинно матрица-матрица-умножение, и sum заменяется предварительным умножением на вектор-строку из единиц. Эти три числа,

1. -0,79463857 умножение Y-канала каждого пикселя (яркость),
2. -0.18382267 умножение Cr (красная разница), и
3. 0,11589724 умножение Cb (синяя разница),

задайте "идеальное" средневзвешенное значение для этого конкретного изображения: каналы Y/Cr/Cb каждого пикселя выровнены с ковариационной матрицей изображения и суммированы. Говоря численно, значение Y каждого пикселя слегка ослаблено, его значение Cr значительно ослаблено, а значение Cb еще более ослаблено, но с противоположным знаком - это имеет смысл, мы ожидаем, что яркость будет наиболее информативной для шкалы яркости. поэтому его вклад самый высокий.

Незначительное предостережение

Я не совсем уверен, откуда происходит преобразование OpenCV в RGB в YCrCb. Документация для cvtColor, в частности, раздел о RGB ↔︎ YCrCb JPEG, похоже, не соответствует ни одному из преобразований, указанных в Википедии. Когда я использую, скажем, пакет Matlasb Transformations Colorspace для преобразования RGB в YCrCb (который цитирует статью в Википедии), я получаю более хорошее изображение в градациях серого, которое кажется более похожим на рисунок 4:

Я совершенно не в себе, когда речь заходит об этих преобразованиях цвета - если кто-то может объяснить, как получить эквиваленты преобразования цветов в Википедии или Matlab в Python/OpenCV, это было бы очень любезно. Тем не менее, это предостережение о подготовке данных. После того, как вы делаете Izycc массив данных с нулевым средним 3×N, приведенный выше код полностью определяет оставшуюся обработку.

Answer 2

Принятый ответ Ахмеда, к сожалению, неверен в математике PCA, что приводит к совершенно другому результату, чем рукопись. Вот изображения экрана, снятые из рукописи.

Среднее центрирование и SVD должны быть выполнены по другому измерению, причем каналы рассматриваются как разные выборки. Среднее центрирование нацелено на получение среднего отклика пикселя, равного нулю, а не среднего отклика канала, равного нулю.

Связанный алгоритм также четко утверждает, что проекция модели PCA предполагает умножение изображения сначала на оценки, а на это произведение - на собственные значения, а не наоборот, как в другом ответе.

Для получения дополнительной информации о математике см. Мой ответ по математике PCA здесь.

Разницу в коде можно увидеть в выходных данных. Поскольку в рукописи не содержится пример выходных данных (что я обнаружил), между результатами могут быть тонкие различия, поскольку рукописные снимки являются снимками экрана.

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

Сначала результат из кода Ахмеда:

Тогда результат из обновленного кода:

Исправленный код (основанный на Ахмеде для удобства сравнения)

import numpy as np
import cv2
from numpy.linalg import svd, norm

# Read input image
Ibgr = cv2.imread('path/peppers.jpg')
#Convert to YCrCb
Iycc = cv2.cvtColor(Ibgr, cv2.COLOR_BGR2YCR_CB)

# Reshape the H by W by 3 array to a 3 by N array (N = W * H)
Izycc = Iycc.reshape([-1, 3]).T

# Remove mean along Y, Cr, and Cb *separately*!
Izycc = Izycc - Izycc.mean(0)  #(1)[:, np.newaxis]
# Mean across channels is required (separate means for each channel is not a 
# mathematically sensible idea) - each pixel's variation should centre around 0

# Make sure we're dealing with zero-mean data here: the mean for Y, Cr, and Cb
# should separately be zero. Recall: Izycc is 3 by N array. 
# Original assertion was based on a false presmise. Mean value for each pixel should be 0
assert(np.allclose(np.mean(Izycc, 0), 0.0))

# Compute data array's SVD. Ignore the 3rd return value: unimportant in this context. 
(U, S, L) = svd(Izycc, full_matrices=False)

# Square the data's singular vectors to get the eigenvalues. Then, normalize
# the three eigenvalues to unit norm and finally, make a diagonal matrix out of
# them. 
eigvals = np.diag(S**2 / norm(S**2))

# Eigenvectors are just the right-singular vectors.
eigvecs = U;

# Project the YCrCb data onto the principal components and reshape to W by H
# array.
# This was performed incorrectly, the published algorithm shows that the eigenvectors 
# are multiplied by the flattened image then scaled by eigenvalues
Igray = np.dot(eigvecs.T, np.dot(eigvals, Izycc)).sum(0).reshape(Iycc.shape[:2])
Igray2 = np.dot(eigvals, np.dot(eigvecs, Izycc)).sum(0).reshape(Iycc.shape[:2])
eigvals3 = eigvals*[1,-1,1]
Igray3 = np.dot(eigvals3, np.dot(eigvecs, Izycc)).sum(0).reshape(Iycc.shape[:2])
eigvals4 = eigvals*[1,-1,-1]
Igray4 = np.dot(eigvals4, np.dot(eigvecs, Izycc)).sum(0).reshape(Iycc.shape[:2])

# Rescale Igray to [0, 255]. This is a fancy way to do this.
from scipy.interpolate import interp1d
Igray = np.floor((interp1d([Igray.min(), Igray.max()],
                            [0.0, 256.0 - 1e-4]))(Igray))
Igray2 = np.floor((interp1d([Igray2.min(), Igray2.max()],
                            [0.0, 256.0 - 1e-4]))(Igray2))
Igray3 = np.floor((interp1d([Igray3.min(), Igray3.max()],
                            [0.0, 256.0 - 1e-4]))(Igray3))
Igray4 = np.floor((interp1d([Igray4.min(), Igray4.max()],
                            [0.0, 256.0 - 1e-4]))(Igray4))

# Make sure we don't accidentally produce a photographic negative (flip image
# intensities). N.B.: `norm` is often expensive; in real life, try to see if
# there's a more efficient way to do this.
if norm(Iycc[:,:,0] - Igray) > norm(Iycc[:,:,0] - (255.0 - Igray)):
    Igray = 255 - Igray
if norm(Iycc[:,:,0] - Igray2) > norm(Iycc[:,:,0] - (255.0 - Igray2)):
    Igray2 = 255 - Igray2
if norm(Iycc[:,:,0] - Igray3) > norm(Iycc[:,:,0] - (255.0 - Igray3)):
    Igray3 = 255 - Igray3
if norm(Iycc[:,:,0] - Igray4) > norm(Iycc[:,:,0] - (255.0 - Igray4)):
    Igray4 = 255 - Igray4

# Display result
if True:
    import pylab
    pylab.ion()
fGray =    pylab.imshow(Igray, cmap='gray')
# Save result
cv2.imwrite('peppers-gray.png', Igray.astype(np.uint8))

fGray2 = pylab.imshow(Igray2, cmap='gray')
# Save result
cv2.imwrite('peppers-gray2.png', Igray2.astype(np.uint8))

fGray3 =pylab.imshow(Igray3, cmap='gray')
# Save result
cv2.imwrite('peppers-gray3.png', Igray3.astype(np.uint8))

fGray4 =pylab.imshow(Igray4, cmap='gray')
# Save result
cv2.imwrite('peppers-gray4.png', Igray4.astype(np.uint8))

****РЕДАКТИРОВАТЬ*****

После запроса Назлока о нестабильности направления собственных векторов (направление, в котором ориентированы любые собственные векторы, является произвольным, поэтому нет гарантии, что разные алгоритмы (или отдельные алгоритмы без воспроизводимого шага стандартизации для ориентации) дадут одинаковый результат. добавлено в двух дополнительных примерах, где я просто поменял знак собственных векторов (номер 2 и номера 2 и 3). Результаты снова разные, с переключением только PC2 дает гораздо более светлый тон, в то время как переключение 2 и 3 это аналогично (неудивительно, поскольку экспоненциальное масштабирование сводит к минимуму влияние PC3). Я оставлю этот последний для тех, кто потрудился запустить код.

Заключение

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