%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%    画像分離シミュレーションのためのサンプルプログラム    %%%%%%%%%
%%%%%%%%                                         ver 3.0          %%%%%%%%%
%%%%%%%%                   i)元画像を表示するようにした           %%%%%%%%%
%%%%%%%%                  ii)ノイズも加えた                       %%%%%%%%%
%%%%%%%%                 iii)アルゴリズムを選択するようにした     %%%%%%%%%
%%%%%%%%                  iv)再現性確保のため乱数のたねを指定     %%%%%%%%%
%%%%%%%%                   v)画像表示の途中に停止(pause)を導入    %%%%%%%%%
%%%%%%%%                                        樺島祥介          %%%%%%%%%
%%%%%%%%                            modified by 池田思朗          %%%%%%%%%
%%%%%%%% samp1,samp2 (n1 × T1 データ) から混合画像を作って分離   %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 画像を読み込む
load samp1;
load samp2;

% 画像の割愛
samp1 = samp1(1:40,:);
samp2 = samp2(1:40,:);

% 画素スケールの設定
scal = 64;

% 行列のサイズ [n1, T1]=[n2, T2] でないとだめなことに注意.
[n1,T1]	= size(samp1);
[n2,T2]	= size(samp2);

% 画像のサイズが同じであるかチェック modified by kaba
if (n1~=n2)
  fprintf('Data sizes of the two images must be the same!\n');
  clear samp1 samp2;
  quit;
elseif (T1~=T2),
  fprintf('Data sizes of the two images must be the same!\n');
  clear samp1 samp2;
  quit;
end;

% ノイズのたねを入れる
rand('seed',65875);
% ノイズを作る modified by kaba
noise = scal * rand(n1,T1);

disp( ' Original images:');

% 元画像の表示 modified by kaba
clf;

%figure(1);
subplot(4,3,1);
colormap(gray);
image(samp1);
axis equal; axis off;
title( 'Original images:');

%figure(2);
subplot(4,3,2);
colormap(gray);
image(samp2);
axis equal; axis off;

%figure(3);
subplot(4,3,3);
colormap(gray);
image(noise);
axis equal; axis off;

% 画像ヒストグラムの表示 added by Doya

subplot(4,3,4);
hist( samp1(:), 20);
title( 'Pixel value histograms:');

subplot(4,3,5);
hist( samp2(:), 20);

subplot(4,3,6);
hist( noise(:), 20);

% 何かタイプするまで待つ
% disp( ' Hit a key to see mixed images:');
% pause;


% 混合比 modified by kaba and ysaka
f11=.4;f12=.3;f13=1.0-(f11+f12);
f21=.3;f22=.4;f23=1.0-(f21+f22);
f31=.3;f32=.3;f33=1.0-(f31+f32);

% 混合行列 added by Doya
A = [ f11,f12,f13; f21,f22,f23; f31,f32,f33];
disp( 'mixing matrix A:');
disp( A);

% 画像の混合 modified by ysaka

mix1=f11*samp1+f12*samp2+f13*noise;
mix2=f21*samp1+f22*samp2+f23*noise;
mix3=f31*samp1+f32*samp2+f33*noise;
NT=n1*T1;

% 混合画像の表示 modified by shiro 

%figure(4);
subplot(4,3,7);
colormap(gray);
image(mix1);
axis equal; axis off;
title( 'Mixed images:');

%figure(5);
subplot(4,3,8);
colormap(gray);
image(mix2);
axis equal; axis off;

%figure(6);
subplot(4,3,9);
colormap(gray);
image(mix3);
axis equal; axis off;

% disp( ' Hit a key to see separated images:');
% pause;
method = input( ' Select an algorithm...1:Bell&Sejnowski, 2:Cardoso, 3:Olshausen&Field ');

% 学習データへの加工 3 × n*T 行列へ modified by kaba and ysaka
X=[reshape(mix1,1,NT);reshape(mix2,1,NT);reshape(mix3,1,NT)];

% 分離行列の学習: どのアルゴリズムを用いるか選択．modified by kaba
switch( method)
case 1
	B=bell(X); % Bell & Sejnowsky のアルゴリズム  (チューニング未)
case 2
	B=jadeR(X); % Cardoso のICA アルゴリズム JADE (チューニング済)
case 3
	B=nlem(X); % Olshausen & Field の非線形EMアルゴリズム (チューニング未)
otherwise
	error( [' ...invalid choice' method]);
end

disp( 'separating matrix B:')
disp( B);

disp( 'B*A');
disp( B*A);

% 分離
Y=B*X;

% 分離画像
C1=Y(1,:);
C2=Y(2,:);
C3=Y(3,:);


% 分離画像１表示 modified by shiro
%figure(7);
subplot(4,3,10);
colormap(gray);
image(reshape(scal*(C1-min(C1)*ones(1,NT))/(max(C1)-min(C1)),n1,T1));
axis equal; axis off;
title( 'Separated images:');

% 分離画像２表示 modified by shiro
%figure(8);
subplot(4,3,11);
colormap(gray);
image(reshape(scal*(C2-min(C2)*ones(1,NT))/(max(C2)-min(C2)),n1,T1));
axis equal; axis off;

% 分離画像３表示 modified by shiro
%figure(9);
subplot(4,3,12);
colormap(gray);
image(reshape(scal*(C3-min(C3)*ones(1,NT))/(max(C3)-min(C3)),n1,T1));
axis equal; axis off;

% データ消去
clear samp1;
clear samp2;
clear noise;