可視光像の生成

シミュレートした移動データから可視光像を生成

この例では、64 フレームでフレームサイズが 64 x 64 ピクセル (毎秒 10 フレーム) のムービーを生成する方法を示します。このムービーには、構造化背景 (これ自体も動いている) の中を動いているターゲットのシミュレーションが含まれています。ランダムな振動によるジッターの動きも (Simulink® モデル "aero_vibrati" で) 生成され、ジッターの動きはセンサー全体の動きに加えられます。最後に、ガウス光学点像分布関数を通じてイメージにブラーが加えられます。

メモ: ここでデルタを変更する場合、Simulink モデル "vibration" のパラメーター設定ダイアログボックスでも変更が必要です。

delt = 0.1;      % Sample time of the generated sequence
num_frames= 64;  % Number of frames to generate
framesize = 64;  % Square frame size in pixels

out = zeros(framesize,framesize,num_frames);    % Initialize movie storage as a 3D Array

ターゲットの生成とその動きの定義

第 1 段階として、ターゲットオブジェクトの形状と動きを定義します。選択した形状は大きな正符号であり、イメージは、各ピクセル位置でのイメージ強度を表す行列によって定義されます。このターゲットは、イメージの中央から右下へ移動するように定義されています。

target = [zeros(3,11)
          zeros(1,5) 6 zeros(1,5)
          zeros(1,5) 6 zeros(1,5)
          zeros(1,3) 6 6 6 6 6 zeros(1,3) % Target is a plus sign 5 by 5 pixels across
          zeros(1,5) 6 zeros(1,5)         %  with an intensity of 6 (S/N ratio is ~4).
          zeros(1,5) 6 zeros(1,5)         % The total target image is made on an 11x11 grid to
          zeros(3,11)];                   %  allow the image to be interpolated without error.

target_velx = 1;                 % target velocity in x direction in pixels per second
target_vely = 1;                 % target velocity in y direction in pixels per second
target_x_initially = framesize/2; % the target is initially in the center of the frame in x
target_y_initially = framesize/2; % and in y

figure(1);
colormap('gray');
image(target*32);
title('Target Image')

背景とターゲットの統合イメージの作成

正弦関数的に相関する背景を生成し、ドリフトの動きを与えます。次に、ターゲットを背景イメージに重ね合わせます。

backsize = framesize+36;  % Make the background bigger than the frame so when it
                          % drifts there are new pixels available to drift into.
xygrid = (1:backsize)/backsize;
B=2*sin(2*pi*xygrid).^2'*cos(2*pi*xygrid).^2;

psd = fft2(B);
psd = real(psd.*conj(psd));

background = B + 0.5*randn(backsize);    % Add a specular Gaussian white
					 % sequence to the structure with
					 % variance of 0.25  (sigma of 0.5).

xoff = 10;
yoff = 10;     % Sensor location is offset from the 0,0 of the background
driftx = 1;
drifty = 1;    % drift rate of the background in a and y directions pix/sec.

minout = min(min(min(background)));
maxout = max(max(max(background)));
colormap('gray');
image((background-minout)*64/(maxout-minout))
title('Background image with additive white specular noise')

トラッカーの回転振動のシミュレーション

トラッカーの回転振動は、モデル aero_vibrati を使用してシミュレートされます。トラッカーの振動のシミュレーションに必要なデータは、Simulink モデル "aero_vibrati" を実行することによって生成されます。

sim コマンドを使用して Simulink 振動モデルを実行します (デルタを 0.1 秒から変更した場合は、Simulink モデルも変更することにより、振動のサンプル時間がこのトラッカーイメージモデルのサンプル時間に一致するようにしなければなりません)。

結果として得られるランダムな回転を図 1 に示します。

omega = 2*pi*5;       % The structural frequencies are 5, 10 and 15 Hz in the model.
zeta  = 0.01;         % Damping ratio for all modes

open_system('aero_vibrati')
simout = sim('aero_vibrati','SrcWorkspace','current');

vibdat = simout.get('vibdat');          % The Simulink model "aero_vibrati"
                                        % generates the vibration data at
                                        % a sample time of 0.01 sec.
vibx = vibdat(1:10:1000);               % The output of simulation is
                                        % returned as the variable simout
                                        % The variable simout contains
viby = vibdat(1001:10:2000);            % the in array vibdat that contains
                                        % the vibration data

levarmx = 10;   % Rotational lever arm for vibration noise in x
levarmy = 10;   %  and in y.

subplot(211);
plot(0.01*(1:10:1000),vibx);grid;
title('Time history of the random Tracker rotations')
xlabel('Time');ylabel('x direction')

subplot(212);
plot(0.01*(1:10:1000),viby);grid;
xlabel('Time');ylabel('y direction')

背景、ターゲット、ジッターからのモーションエフェクトのシミュレーション

ムービーを構成することになるフレームは、多次元配列で作成、保存されるようになっています。各フレームの背景とターゲットは、ターゲットの動き、背景のドリフト、およびトラッカーの振動のため、位置が異なります。ムービーの最初のフレームを図 1 に示します。

clf; drawnow;

for t = 1:num_frames

  % Drift the Background at the rate driftx and drifty
  % (in pixels/second) and add in the vibration:
  xshift = driftx*delt*t+levarmx*vibx(t,1);
  yshift = drifty*delt*t+levarmy*viby(t,1);

  % Interpolate the 2D image using the MATLAB(R) function interp2:
  [xgrid, ygrid]   = meshgrid(1:backsize);
  [xindex, yindex] = meshgrid(xshift:1:xshift+backsize,yshift:1:yshift+backsize);
  outtemp = interp2(xgrid,ygrid,background,xindex,yindex);

  % Truncate the drifted image down from backsize to framesize:
  out(:,:,t) = outtemp(xoff:xoff+framesize-1,xoff:xoff+framesize-1);

  % Now let the target move also:
  tpixinx = floor(target_velx*delt*t);
  tpixiny = floor(target_vely*delt*t);  % Before interpolating extract the number of pixels moved
  txi = target_velx*delt*t - tpixinx;
  tyi = target_vely*delt*t - tpixiny;   % Interpolate on sub-pixels around the origin only
  [txgrid tygrid] = meshgrid(1:11);     % meshgrid here generates a matrix of grid elements
  [txi tyi] = meshgrid(txi+1:txi+11,tyi+1:tyi+11); % meshgrid generates 2 matrices with the x and y grids

  % Interpolate the intensity values first using interp2 -- a built in MATLAB command
  temp = interp2(txgrid,tygrid,target,txi,tyi);

  % Insert the target at the location determined by the initial offset, and the number of whole pixels moved
  tx = tpixinx + target_x_initially-1;
  ty = tpixiny + target_y_initially-1;
  out(tx:tx+6,ty:ty+6,t) = temp(9:-1:3,9:-1:3) + out(tx:tx+6,ty:ty+6,t);

end

minout = min(min(min(out)));
maxout = max(max(max(out)));
colormap('gray');
image((out(:,:,1)-minout) * 64/(maxout-minout));
title('First frame of combined target and background image.')

光学系を通じてイメージを渡す -- ガウス "開口関数" の使用

このコードセグメントでは、実測開口関数も簡単に使用できます。次の 5 行を "load measured_aperture" で置換するだけです。measured_aperture は ASCII に保存されている実測関数であり、ファイル measured_aperture.mat に保存されているデータは、行列 apfunction を含んでいる MATLAB .mat ファイルです (MATLAB .mat ファイルの読み書き方法を示す C および Fortran コードを読み込んで表示する方法を参照するには、MATLAB で「help load」と入力してください)。

(メモ: 点像分布関数がガウスの場合は、開口関数もガウスです。)

トラッカー光学系の効果をシミュレートするために、ムービーの各フレームには 2-D FFT (高速フーリエ変換) でブラーが加えられています。結果として得られるイメージの最初のフレームを図 1 に示します。

x = 1:framesize;
y = 1:framesize;
sigma      = 120;
apfunction = exp(-(x-framesize/2).^2/(2*sigma))' * exp(-(y-framesize/2).^2/(2*sigma));
apfunction = fftshift(apfunction);      % Rotate so it conforms with FFT convention

for j = 1:num_frames
  out(:,:,j) = real(ifft2(apfunction.*fft2(out(:,:,j))));
end

minout = min(min(min(out)));
maxout = max(max(max(out)));
colormap('gray');
image((out(:,:,1)-minout)*64/(maxout-minout));
title('First frame of blurred image.')

MATLAB ムービーの生成と再生

ムービーフレームを最小値から最大値まで 64 の強度値をもつようにスケーリングした後、イメージとしての結果を表示します。moviein と getframe の仕組みの説明は、MATLAB のヘルプを参照してください。

minout = min(min(min(out)));
maxout = max(max(max(out)));

M = moviein(num_frames);
for j = 1:num_frames
  image((out(:,:,j)-minout)*64/(maxout-minout))
  drawnow
  M(:,j) = getframe;
end

% colormap('gray')
% movie(M);

オプション: ムービーを .mat ファイルとして保存

生成されたトラッカームービーをオプションで mat ファイルとして保存できます。また、背景の PSD も、後でムービーと共に使用するために保存できます。

save trackerimage out
save psdback psd
save moviedat M

bdclose('aero_vibrati');