二値化画像から主成分​分析するにはどうした​らよいですか？

二値化された画像の白部分に対して主成分分析をしたいです。

流れとしては、白部分1ピクセルごとの座標を取得して、それをもとに主成分分析を行い、第一主軸、第二主軸まで出すことを考えています。

また、第一主軸、第二主軸の長さも出したいと考えています。

座標取得のやり方から分からないので、ご教授頂きたいと存じます。

Answer 1

[Hiro Yoshino]

主成分分析 / KL 変換 はデータの集合に対して、その集合の情報量 (分散) が大きくなる部分空間 (基底) を順番に探していく手法ですので、画像一枚からではそれは求まりません。画像を沢山あつめて、その画像から基底となる画像を見つけるイメージです。

ここではやり方だけ示しておきます。

1. 画像を n 枚集めます。画像のピクセルを1行に並べて X とします
2. 行ベクトルは平均 0 になるようにしなくても pca でやってくれます。
3. 以下を実行:

[coeff,score,latent] = pca(X)

これで得られた、coeff が X の分散共分散行列の固有値で、PCF loadingとか因果負荷量と呼ばれる物です。これが先の基底に相当します。固有値が大きい順（影響が大きい順) に並んでいるはずです。score がこの部分空間中での元の画像の座標を示しているので、score と coeff の線形結合で元の画像に復元されます。

基本的に固有ベクトルは長さ1に正規化されるので、主軸の長さの議論は意味が有りません。主軸に画像を射影したときの長さが、coeff に対応しています。

pca の計算 (固有ベクトルの計算) に使っていない画像データも

img'*score(1,:) % 第一主軸の成分を抽出

こんな感じで、内積をとってあげれば 主軸の score (成分) を得ることが出来ます。

画像の注目するとこだけ抜き出す一般的な方法（画像処理の流れ）を以下に示しておきます：

I = imread('coins.png');
%BW=imbinarize(I); % そのままやると輝度によって綺麗に取れない
BW = bwareaopen(I>80,100); % モルフォロジー処理
imshowpair(I,BW,"montage");

% 余計な場所を削除
IoI = I;
IoI(~BW)=0;
imshow(IoI)

[R]

回答のほどありがとうございます。

説明が下手ですみません。

想定していることは、図に示すようなことです。

二値化画像の白色部分に注目し、1ピクセル1データ点として扱いたいです。

この場合のプログラムをご教授頂きたいです。

手間をとらせてしまって申し訳ございませんが、回答していただけると幸いです。

[Hernia Baby]

自分の勘違いだと思って消してしまいました…

scatter情報が画像になっていて、それを一度座標系に直してpcaかけたいんですね

ここでは二値化画像を座標データとして換算する方法を書いておきます。

findをご使用ください。

I = rgb2gray(imread('peppers.png'));
idx = I > 100;
imshow(idx)

[r, c] = find(idx);
scatter(c,max(r) - r,'+k')

[Hernia Baby]

ちなみに散布図の主成分分析についてはちょうど良さそうなのがありました。

Principal Component Analysis by Singular Value Decomposition

clear; clc; close all
%Generate Data
%--------------
u1=5*randn(1,10000);
u2=randn(1,10000);
x1=1/2^0.5*u1+1/2^0.5*u2;
x2=1/2^0.5*u1-1/2^0.5*u2;
X=[x1 ; x2];
%PCA by SVD
%--------------
X=X-mean(X,2);   %zero-mean data
[U,S,V]=svd(X);  %svd
[m, n]=size(X);
S2=S(1:m,1:m); Values=S2^2/n;  %Eigenvalues
Vectors=U                      %Eigenvectors

Vectors = 2×2

-0.7087 0.7056 -0.7056 -0.7087

Amplitudes=S*conj(V');         %Orthognal amplitudes

%Plot Output
%--------------
subplot(1,2,1);
plot(X(1,:),X(2,:),'.',Values(1,1)*[0 ; Vectors(1,1)],Values(1,1)*[0 ; Vectors(2,1)],'m',Values(1,1)*[0 ; Vectors(1,2)],Values(1,1)*[0 ; Vectors(2,2)],'g');
daspect([1 1 1]);
xlim([-20 20]); ylim([-20 20]);
title('Original Data');
subplot(1,2,2);
plot(Amplitudes(1,:),Amplitudes(2,:),'r.');
daspect([1 1 1]);
xlim([-20 20]); ylim([-20 20]);
title('Orthognal Amplitudes');

[Hiro Yoshino]

@R 失礼しました。てっきり画像の低ランク近似の話だと思ってしまいました。

単純に分散最大軸を見つけたいという事ですね。

@Hernia Baby さんのように、(x,y) 座標を取ってきて、行方向に並べて X に代入します。

mu = [2 3];
Sigma = [1 1.5; 1.5 3];
X = mvnrnd(mu,Sigma,50) % (x,y) 座標

X = 50×2

1.5184 1.6437 1.7754 3.2763 4.6480 6.4286 4.0800 5.6210 1.3291 1.9194 1.8776 1.5732 0.1347 0.8662 2.9640 4.0611 0.5880 1.1263 4.6938 6.4525

plot(X(:,1),X(:,2),'r.')

[coef, score] = pca(X)

coef = 2×2

0.4699 0.8827 0.8827 -0.4699

score = 50×2

-1.5330 0.1885 0.0289 -0.3518 4.1613 0.7027 3.1816 0.5809 -1.3786 -0.1081 -1.4264 0.5387 -2.8695 -0.6676 1.2802 0.3287 -2.4269 -0.3897 4.2039 0.7320

楽な方法は biplot 関数を使う事です。biplot

biplot(coef,'Scores',score);

[Hiro Yoshino]

@R すみません、タイミングの問題で @Hernia Baby さんと被ってしまいましたが、svd でも理論的には同じ固有ベクトルを獲得出来ます。

分散共分散行列を使って計算するよりもsvd の方が計算が速いので、恐らく pca でも行列サイズによっては内部的に svd が動いていると思います。若干用途は異なりますが、非常に似た考え方で、兄弟と思っても良いです。

2つは同じでは無いですが、主成分を抽出すると言う点では同じです。

[R]

お二方ともありがとうございます。

とても詳しく説明があり、分かりやすいです。

どちらも参考にし、試してみたいと思います。