ソースコード

% ------------------------------------------------------

% 120224

% AlignStack作成

% ------------------------------------------------------

% test_FRETImageProcess → FRETImageProcess2

% Showを廃止、代わりにAdjustを作成

% Movieの画像補正の部分をAdjustにする

% SortNLSNESCAAXでNESかCAAXかを判定する基準のSDを引数で指定させる

% 核内intensityの微分のSDをNES判定材料に加える

% MovieFromTiff作成

% RatioImageFromTiff作成

% ------------------------------------------------------

% FRETImageProcess → test_FRETImageProcess

% parforに対応させる

% parfor内でtiffオブジェクトを読み込むとバグが起こるのでReadTiffFileでtiffオブジェクトではなく画像とタグとして読み込む

% それにあわせてGetImageとGetTimeofImageを変更、RegionForOutputのImageLengthとImageWidthの読み込み方を変更

% 画像表示（自動補正）

% レシオ画像作成

% ムービー

% コメント化された旧アルゴリズムを削除

% -------------------------------------------------------

classdef FRETImageProcess2

methods (Static)

% イメージング情報が記述されたndファイルを読み込んでimaging_info構造体に格納

function imaging_info = ReadNDFile(directory, filename)

% ファイルアドレスをimaging_data構造体に格納

imaging_info.directory = directory;

imaging_info.filename = filename;

% ndファイルを読み込んで文字列リストとして格納

fid = fopen([directory, '\', filename, '.nd']);

cellarray = textscan(fid, '%s', 'delimiter', ','); % セル配列が返る

nd = cellarray{1}; % 行列に格納

fclose(fid);

% 読み込んだ文字列リストからimaging_info構造体に格納

for i = 1 : length(nd)

if strcmp(nd(i), '"StartTime1"') % start day and time (秒単位 seconds)

date = nd{i + 1}; % nd(i + 1)とするとセル配列が入る

imaging_info.start_time = datenum(date,'yyyymmdd HH:MM:SS') * 86400; % 文字列から秒単位の数値に変換

elseif strcmp(nd(i), '"DoTimelapse"') % Timelapse on or off

imaging_info.do_timelapse = nd{i + 1};

elseif strcmp(nd(i), '"NTimePoints"') % number of timepoints

imaging_info.timepoints = str2double(nd(i + 1));

elseif strcmp(nd(i), '"NStagePositions"') % number of positions

imaging_info.positions = str2double(nd(i + 1));

elseif strcmp(nd(i), '"NWavelengths"') % number of wavelengths

imaging_info.wavelengths = str2double(nd(i + 1));

elseif strcmp(nd{i}(2:length(nd{i}) -2) , 'WaveName') % 「" "」と数字を除く("WaveName1" → WaveName)

num = str2double(nd{i}(length(nd{i}) -1)); % nd(i)とすると文字列ではなくセル配列が返る

imaging_info.wavename{num} = nd{i + 1}(2:length(nd{i + 1}) -1); % wavenameはセル配列

end

end

end % ReadNDFile

% Tiffファイルの読み込み

function [Image_array Time_array] = ReadTiffFile(imaging_info) % 引数はReadNDFile関数の返り値の構造体

Image_array = cell(imaging_info.wavelengths, imaging_info.positions, imaging_info.timepoints);

Time_array = zeros(imaging_info.wavelengths, imaging_info.positions, imaging_info.timepoints);

warning off all; % 警告メッセージを非表示にする

for i = 1 : imaging_info.wavelengths

for j = 1 : imaging_info.positions

for k = 1 : imaging_info.timepoints

% Tiffファイルのパスを設定

if strcmp(imaging_info.do_timelapse, 'FALSE')

path = [imaging_info.directory, '\', imaging_info.filename, '_w', int2str(i), imaging_info.wavename{i}, '_s', int2str(j), '.TIF'];

else

path = [imaging_info.directory, '\', imaging_info.filename, '_w', int2str(i), imaging_info.wavename{i}, '_s', int2str(j), '_t', int2str(k), '.TIF'];

end

% Tiffオブジェクト（読み込みモード）として読み込み

tiff_temp = Tiff(path, 'r');

% readで画像を読み込み

Image_array{i,j,k} = read(tiff_temp);

% 撮影時刻を読み込み

DateTime = tiff_temp.getTag('DateTime');

% 文字列から秒単位の数値に変換

Time_array(i,j,k) = datenum(DateTime, 'yyyy:mm:dd HH:MM:SS') * 86400;

end

end

end

warning on all; % 警告メッセージ非表示の解除

end % ReadTiffFile

% 読み込んだTiffオブジェクトから指定した画像を取得

function Image = GetImage(Imaging_info, Tiff_array, timepoint, position, wavename)

% wavenameが何番目のwavelengthかを調べる

for i = 1:Imaging_info.wavelengths

if strcmp(Imaging_info.wavename{i} , wavename)

waveID = i;

break

end

if i == Imaging_info.wavelengths

fprintf('指定したwavenameは存在しません \n')

waveID = 0;

end

end

Image = Tiff_array{waveID, position, timepoint};

end % GetImage

% 読み込んだTiffオブジェクトから指定した画像の撮影時刻を取得

function TimeOfImage = GetTimeOfImage(Imaging_info, Time_array, timepoint, position, wavename)

% wavenameが何番目のwavelengthかを調べる

for i = 1:Imaging_info.wavelengths

if strcmp(Imaging_info.wavename{i} , wavename)

waveID = i;

break

end

if i == Imaging_info.wavelengths

fprintf('指定したwavenameは存在しません \n')

waveID = 0;

end

end

% 撮影時刻を取得

TimeOfImage = Time_array(waveID, position, timepoint);

end % GetTimeOfImage

% 画像をintensityが最大255,最小0になるように自動補正

function AdjustedImage = Adjust(Image)

% doubleに型変換

Image = double(Image);

% ratio画像（カラー）の場合、各色の画像を並べる

if size(Image,3) == 3

ref_image = [Image(:,:,1), Image(:,:,2), Image(:,:,3)];

else

ref_image = Image;

end

width = size(ref_image,2);

height = size(ref_image,1);

% intensity histogram作成

mi = max(max(ref_image));

dist = zeros(mi+1,2);

dist(:,1) = 0:mi;

for y = 1:size(ref_image,1)

for x = 1:size(ref_image,2)

dist(ref_image(y,x)+1,2) = dist(ref_image(y,x)+1,2) + 1;

end

end

% intensity windowをmin_intとmax_intをintensity histogramの上下0.3%になるように設定

min_adjust = 0;

n = 0;

while n < 0.003 * width * height

min_adjust = min_adjust + 1;

n = sum(dist(1:min_adjust+1, 2));

end

rows = size(dist,1);

max_adjust = rows - 1;

m = 0;

while m < 0.003 * width * height

max_adjust = max_adjust - 1;

m = sum(dist(max_adjust+1:rows, 2));

end

% intensityを最大255,最小0になるように補正

AdjustedImage = uint8((Image - min_adjust) / (max_adjust - min_adjust) * 255);

end % Adjust

% アニメーション表示、ムービーフレーム作成（フレームはmovie2avi関数でAVIファイルに変換することができる）

function Frames = Movie(ImageCellArray, interval)

% モノクロまたはカラーの画像のセル配列ImageCellArrayをアニメーション表示

figure;

for i = 1 : length(ImageCellArray)

% １枚の画像を取得、doubleに型変換

Image = double(ImageCellArray{i});

% 画像を自動補正

adjusted_image = FRETImageProcess2.Adjust(Image);

% 描画

imshow(adjusted_image);

% 画像更新

drawnow;

% ムービーフレームとして格納（フレームにすると上下が反転するためflipdimで修正）

Frames(i) = im2frame(flipdim(adjusted_image, 1), colormap);

% インターバル

pause(interval);

end

% フレームから再生

% movie(Frames);

% aviファイル作成

% movie2avi(Frames, filename, 'fps', 10);

end % Movie

% アニメーション表示、ムービーフレーム作成（フレームはmovie2avi関数でAVIファイルに変換することができる）

% ReadTiffFileで読み込んだ{wavelength, position, timepoint}の配列を引数とする

function Frames = MovieFromTiff(Imaging_info, TiffArray, Position, Wavename, interval)

figure;

% wavenameが何番目のwavelengthかを調べる

for i = 1:Imaging_info.wavelengths

if strcmp(Imaging_info.wavename{i} , Wavename)

waveID = i;

break

end

if i == Imaging_info.wavelengths

fprintf('指定したwavenameは存在しません \n')

waveID = 0;

end

end

for timepoint = 1 : size(TiffArray, 3)

% １枚の画像を取得、doubleに型変換

Image = double(TiffArray{waveID, Position, timepoint});

% 画像を自動補正

adjusted_image = FRETImageProcess2.Adjust(Image);

% 描画

imshow(adjusted_image);

% 画像更新

drawnow;

% ムービーフレームとして格納（フレームにすると上下が反転するためflipdimで修正）

Frames(timepoint) = im2frame(flipdim(adjusted_image, 1), colormap);

% インターバル

pause(interval);

end

% フレームから再生

% movie(Frames);

% aviファイル作成

% movie2avi(Frames, filename, 'fps', 10);

end % MovieWPT

% ratio画像作成(timepoint数分)、ムービー表示

function ResultImage_array = RatioImage(FRETImages, CFPImages, min_ratio, max_ratio)

% 引数のFRETImages,CFPImagesは(height,width)の２次元配列または(height,width,timepoints)の３次元配列

% 戻り値はカラー画像行列のセル配列

% 画像サイズとtimepoint数を取得

height = size(FRETImages,1);

width = size(FRETImages,2);

timepoints = size(FRETImages, 3);

% 戻り値宣言

ResultImage_array = cell(1,timepoints);

% 各timepointにおけるRatio画像作成

for t = 1:timepoints

% 各timepointにおける画像を取得

image_FRET = FRETImages(:,:,t);

image_CFP = CFPImages(:,:,t);

% ratioを計算

ratio = double(image_FRET) ./ double(image_CFP);

% 0で割ってNaNやInfになったものを0にする

ratio(isnan(ratio)) = 0;

ratio(isinf(ratio)) = 0;

% 最小0、最大510になるように補正

ratio = (ratio - min_ratio) / (max_ratio - min_ratio) * 510;

color_ratio = zeros(height, width, 3);

% 1～255 : 青低下、緑上昇 256～510 : 緑低下、赤上昇

color_ratio(:,:,1) = uint8(ratio - 255);

color_ratio(:,:,2) = uint8(-abs(ratio - 255) + 255);

color_ratio(:,:,3) = uint8(255 - ratio);

% intensityによって補正

ref_image = zeros(height, width, 3);

ref_image(:,:,1) = double(image_FRET);

ref_image(:,:,2) = double(image_FRET);

ref_image(:,:,3) = double(image_FRET);

max_int = max(max(max(ref_image)));

min_int = min(min(min(ref_image)));

adjusted_color_ratio = double(color_ratio) .* ((ref_image - min_int) / (max_int - min_int));

ResultImage_array{t} = uint8(adjusted_color_ratio);

end

% ムービー表示

FRETImageProcess2.Movie(ResultImage_array, 0.02);

end % RatioImage

% backgroundを引いてratio画像作成、ムービー表示

% ReadTiffFileで読み込んだ{wavelength, position, timepoint}の配列を引数とする

function ResultImage_array = RatioImageFromTiff(Imaging_info, TiffArray, position, min_ratio, max_ratio)

% FRETとCFPが何番目のwavelengthかを調べる

waveID_FRET = 0;

waveID_CFP = 0;

for i = 1:Imaging_info.wavelengths

if strcmp(Imaging_info.wavename{i} , 'FRET')

waveID_FRET = i;

end

if strcmp(Imaging_info.wavename{i} , 'CFP')

waveID_CFP = i;

end

if i == Imaging_info.wavelengths

if waveID_FRET == 0

fprintf('FRETのwavenameが存在しません \n')

return

end

if waveID_CFP == 0

fprintf('CFPのwavenameが存在しません \n')

return

end

end

end

width = size(TiffArray{1,1,1}, 2);

height = size(TiffArray{1,1,1}, 1);

% FRETとCFPの３次元画像行列(height,width,timepoint)を作成

FRETImages = uint16(zeros(height, width, Imaging_info.timepoints));

CFPImages = uint16(zeros(height, width, Imaging_info.timepoints));

for t = 1:Imaging_info.timepoints

FRETImages(:,:,t) = TiffArray{waveID_FRET, position, t};

CFPImages(:,:,t) = TiffArray{waveID_CFP, position, t};

end

% backgroundを引く

FRET_background = FRETImageProcess2.FigureBackground(FRETImages(:,:,1));

CFP_background = FRETImageProcess2.FigureBackground(CFPImages(:,:,1));

for t = 1:Imaging_info.timepoints

FRETImages(:,:,t) = FRETImages(:,:,t) - FRET_background;

CFPImages(:,:,t) = CFPImages(:,:,t) - CFP_background;

end

% RatioImageを呼び出してレシオ画像作成

ResultImage_array = FRETImageProcess2.RatioImage(FRETImages, CFPImages, min_ratio, max_ratio);

end % RatioImageWPT

% Intensityヒストグラム作成

function Distribution = IntensityHistogram(Image) % Image : Tiff_array{i,j,k} or read(Tiff_array{i,j,k})

if strcmp(class(Image), 'Tiff')

img = read(Image);

else

img = round(Image);

end

m = max(max(img));

dist = zeros(m+1,2);

dist(:,1) = 0:m;

for i = 1:size(img,1)

for j = 1:size(img,2)

dist(img(i,j)+1,2) = dist(img(i,j)+1,2) + 1;

end

end

dist(:,1) = 0:m;

figure, bar(dist(:,2)) % ヒストグラム描画(新規figure)

Distribution = dist;

end % IntensityHistogram

% Background intensityを求める

function Background_intensity = FigureBackground(Image)

% 8x8ピクセルの領域を4ピクセルずつ動かしていき、それぞれの領域における平均輝度のうち最小のものをbackgroundとする

if strcmp(class(Image), 'Tiff')

img = read(Image);

else

img = round(Image);

end

radius = (size(img,1) + size(img,2)) / 4;

mean_array = zeros(length(1:4:size(img,1)-8) * length(1:4:size(img,2)-8) , 3);

count = 1;

for y = 1:4:size(img,1)-8

for x = 1:4:size(img,2)-8

if (((x+3 - radius)^2 + (y+3 - radius)^2) < radius^2) % 四隅は暗いので除く

region = img(y:y+7, x:x+7);

mean_array(count,1) = x;

mean_array(count,2) = y;

mean_array(count,3) = round(mean(mean(region)));

count = count + 1;

end

end

end

% mean_arrayの余分な行を削除

mean_array = mean_array(1:count-1, 1:3);

% mean_arrayの最小値をbackgroundとする

Background_intensity = min(mean_array(:,3));

end

% 画像の位置が時間がたつにつれてずれていくのを修正

% ImageStackは{wavelength,position,timepoint}のセル配列

function AlignedStack = AlignStack(ImageStack)

width = size(ImageStack{1,1,1}, 2);

height = size(ImageStack{1,1,1}, 1);

num_wavelengths = size(ImageStack, 1);

num_positions = size(ImageStack, 2);

num_timepoints = size(ImageStack, 3);

ImageStack2 = ImageStack;

for p = 1 : num_positions

% 各background計算

background = zeros(num_wavelengths, num_positions);

for w = 1:num_wavelengths

background(w,p) = FRETImageProcess2.FigureBackground(ImageStack{w,p,1});

end

% 各timepointについてずれを修正

adjust_x = 0;

adjust_y = 0;

t = 2;

while t <= num_timepoints

subtract = zeros(5,5);

% timepoint=1の画像とtimepoint=tの画像を比較

% timepoint = 1, t の画像を縦横に-2～2ピクセルずらしながら差分を計算

for y = -2:2

for x = -2:2

% adjust_x,y はひとつ前のtについての補正値

cx = adjust_x + x;

cy = adjust_y + y;

subtract(y+3,x+3) = sum(sum(ImageStack{1,p,1}(max(1,1+cy):min(height,height+cy), max(1,1+cx):min(width,width+cx))...

- ImageStack{1,p,t}(max(1,1-cy):min(height,height-cy), max(1,1-cx):min(width,width-cx))));

end

end

% 差分が最小になるようなずらし幅x,yを取得

min_of_subtract = subtract * 0;

min_of_subtract(subtract == min(min(subtract))) = 1;

for x = 1:5

if sum(min_of_subtract(:,x) == 1) > 0

break;

end

end

for y = 1:5

if sum(min_of_subtract(y,:) == 1) > 0

break;

end

end

% 1～5 を -2～2 に補正

x = x - 3;

y = y - 3;

% 合計のadjust_x,y( = position1と合わせるためのadjustx,y)

adjust_x = adjust_x + x;

adjust_y = adjust_y + y;

% x,yが-2または2であればさらにずらして調べる

if (x == -2) || (x == 2) || (y == -2) || (y == 2)

continue;

end

% timepoint = tの画像をずらす

for w = 1:num_wavelengths

ImageStack2{w,p,t} = ImageStack2{w,p,t} * 0 + background(w,p);

ImageStack2{w,p,t}(max(1,1+adjust_y):min(height,height+adjust_y), max(1,1+adjust_x):min(width,width+adjust_x))...

= ImageStack{w,p,t}(max(1,1-adjust_y):min(height,height-adjust_y), max(1,1-adjust_x):min(width,width-adjust_x));

end

t = t + 1;

end

end

AlignedStack = ImageStack2;

end % AlignStack

% watershedアルゴリズムによる核染色画像のsegmentation

function SegmentationMatrix = Segmentation(NucleusImage) % Image : Tiff_array{i,j,k} or read(Tiff_array{i,j,k})

if strcmp(class(NucleusImage), 'Tiff')

img = read(NucleusImage);

else

img = round(NucleusImage);

end

% 前処理として細かくopening-by-reconstruction

se = strel('disk', 2);

le = imerode(img,se);

img = imreconstruct(le,img);

% intensityを平らにする

for i = 1 : 3

rm = imregionalmax(img);

rm = img .* uint16(rm);

rmlow = rm * 0.90; % 局所的最大値の少し下を取得

robr = imreconstruct(rmlow,img);

fgmr = imregionalmax(robr);

cc = bwconncomp(fgmr); % 領域分割

for j = 1 : cc.NumObjects

img(cc.PixelIdxList{1,j}) = mean(img(cc.PixelIdxList{1,j})); % 各領域内の値を平均化

end

end

% figure,imshow(img * 16);

fgm = imregionalmax(img);

% figure,imshow(fgm);

% watershedで領域分割

L_dist = -bwdist(~fgm);

L_dist(~fgm) = -Inf;

L_seg = watershed(L_dist);

% figure,imshow(label2rgb(L_seg,'jet'));

SegmentationMatrix = L_seg;

end

% segmentationした核の位置情報をもとに、核内と核周囲の平均intensityを計算

function intensity_array = CheckIntensityAroundNucleus(FRETImage, SegmentationMatrix)

if strcmp(class(FRETImage), 'Tiff')

img = read(FRETImage);

else

img = uint16(round(FRETImage));

end

se = strel('disk', 1); % imdilate(膨張処理)用

count = 3; % 何ピクセル外側まで調べるか

dilated_nucleusregion = cell(1,count+1);

external_nucleusregion = cell(1,count);

external_image = cell(1,count);

external_intensity = cell(1,count);

result_array = zeros(max(max(SegmentationMatrix)), count+1);

for i = 2:max(max(SegmentationMatrix))

% 核の領域を取得

nucleusregion = uint16(SegmentationMatrix == i);

nucleusregion = imdilate(nucleusregion, strel('disk',1)); % 境界部が含まれないので拡張

% 核内全体の平均intensityを求める

nucleus_image = img .* nucleusregion;

nucleus_intensity = sum(sum(nucleus_image)) / sum(sum(nucleusregion));

result_array(i, 2) = nucleus_intensity;

% 核の縁（0ピクセル外側）の平均intensityを求める

eroded_nucleusregion = imerode(nucleusregion, se);

border_nucleusregion = uint16(eroded_nucleusregion ~= nucleusregion);

border_image = img .* border_nucleusregion;

border_intensity = sum(sum(border_image)) / sum(sum(border_nucleusregion));

result_array(i, 3) = border_intensity;

% 核の外側の平均intensityを求める

dilated_nucleusregion{1} = nucleusregion;

for j = 1:count

dilated_nucleusregion{j+1} = imdilate(dilated_nucleusregion{j},se); % 1ピクセル分膨張処理

external_nucleusregion{j} = uint16(dilated_nucleusregion{j+1}~=dilated_nucleusregion{j}); % 外側の領域だけ取り出し

external_image{j} = img .* external_nucleusregion{j};

external_intensity{j} = sum(sum(external_image{j})) / sum(sum(external_nucleusregion{j}));

result_array(i, j+3) = external_intensity{j};

end

end

intensity_array = result_array;

end

% 領域拡張法 : seedbinaryimage（核）の２倍の直径程度に膨張処理した領域内での領域拡張

% 拡張条件はcriteiaの値（文字列）により異なる。 'NLS', 'NES', 'CAAX'のいずれか

% SeedBinaryImageは１細胞の核の領域のみであることを前提とする

% メモ：imresizeで画像サイズ変更

function GrownBinaryImage = RegionGrowing(image, SeedBinaryImage)

if strcmp(class(image), 'Tiff')

img = read(image);

else

img = uint16(round(image));

end

width = size(img, 2);

height = size(img, 1);

GrownBinaryImage = SeedBinaryImage; % 拡張したイメージ

FinishedBinaryImage = logical(SeedBinaryImage * 0); % region growingを終えたseed point(次のseed pointにはならない)

NextSeedBinaryImage = (GrownBinaryImage ~= FinishedBinaryImage);

% 核(SeedBinaryImage)をひとまわり大きくした領域を求める。GrownBinaryImageはこの領域内におさめる

numberofpoints = sum(sum(SeedBinaryImage));

strel_radius = round(sqrt(numberofpoints)); % 面積がnumberofpointsの円の半径

dilated_region = imdilate(SeedBinaryImage, strel('disk',strel_radius));

% criteria = 'CAAX'のとき用にSeedBinaryImage領域内の平均intensityを計算

%mean_seed_intensity = sum(img(SeedBinaryImage)) / sum(sum(SeedBinaryImage));

% 領域拡張ループ

while sum(sum(NextSeedBinaryImage)) > 0

% seed pointの座標を取得

seed_list = zeros(sum(sum(NextSeedBinaryImage)), 2); % 行数 = seed pointの数、列 = x,y

count = 0;

for y = 1:height

for x = 1:width

if NextSeedBinaryImage(y,x) == 1

count = count + 1;

seed_list(count, 1) = x;

seed_list(count, 2) = y;

end

end

end

% 各seed pointについて領域拡張

for i = 1:count

seed_point_x = seed_list(i,1);

seed_point_y = seed_list(i,2);

seed_intensity = img(seed_point_y, seed_point_x);

% seed pointの周囲８ピクセルを調べる

for y = 1:3

for x =1:3

grow_point_x = seed_point_x + x - 2;

grow_point_y = seed_point_y + y - 2;

if (grow_point_x < 1) || (grow_point_x > width) || (grow_point_y < 1) || (grow_point_y > height)

continue

end

% 領域拡張条件

if (img(grow_point_y, grow_point_x) > seed_intensity * 0.9)

GrownBinaryImage(grow_point_y, grow_point_x) = 1;

end

end

end

% GrownBinaryImageはdilated_region内におさめる

GrownBinaryImage = logical(GrownBinaryImage .* dilated_region);

FinishedBinaryImage(seed_point_y, seed_point_x) = 1;

end

NextSeedBinaryImage = (GrownBinaryImage ~= FinishedBinaryImage);

end

end

% 領域拡張法によるNLS,NES,CAAXの判別

function SortArray = SortNLSNESCAAX(FRETimage, SegmentationMatrix, SD_NES, SD_CAAX)

if strcmp(class(FRETimage), 'Tiff')

img = read(FRETimage);

else

img = uint16(round(FRETimage));

end

cells = zeros(max(max(SegmentationMatrix)) , 1); % 各細胞がNLS,NES,CAAXのどれであるかを格納する行列

nucleus_size_array = zeros(max(max(SegmentationMatrix)) , 1); % 各細胞の核のサイズを格納

nucleus_intensity_array = zeros(max(max(SegmentationMatrix)) , 1); % 各細胞の核の平均intensityを格納

% sd_array = zeros(max(max(SegmentationMatrix)) , 1); j = 1;

intensity_array = FRETImageProcess2.CheckIntensityAroundNucleus(FRETimage, SegmentationMatrix);

% prewittフィルタでエッジ（微分）抽出

horiz = filter2([-1,0,1;-1,0,1;-1,0,1],img);

vert = filter2([-1,-1,-1;0,0,0;1,1,1],img);

differential_img = (horiz .* horiz + vert .* vert) .^ 0.5;

for i = 2:length(cells) % i = 1は背景

nucleus_region = SegmentationMatrix == i;

nucleus_region = imdilate(nucleus_region, strel('disk',1)); % 境界部が含まれないので拡張

mean_nucleus_intensity = mean(img(nucleus_region));

nucleus_size_array(i) = sum(sum(nucleus_region));

nucleus_intensity_array(i) = mean_nucleus_intensity;

% NLSかどうかの判定（周辺数ピクセルを調べる）

r0 = intensity_array(i,3) / intensity_array(i,2); % 核の0ピクセル外側 / 核内全体

r1 = intensity_array(i,4) / intensity_array(i,3); % 核の1ピクセル外側 / 核の0ピクセル外側

r2 = intensity_array(i,5) / intensity_array(i,4); % 核の2ピクセル外側 / 核の1ピクセル外側

r3 = intensity_array(i,6) / intensity_array(i,5); % 核の3ピクセル外側 / 核の2ピクセル外側

min_ratio = min([r0,r1,r2,r3]);

if min_ratio < 0.85

cells(i) = 2; % NLSと断定

continue;

end

% (核の1ピクセル外側 / 核内全体) < 0.9 ならNES,CAAXではないと判断

if r0 * r1 < 0.90

cells(i) = 5; % 判定不可

continue;

end

% 細胞がNESあるいはCAAXだと仮定して細胞質領域(cytoplasm_region)を求め、判定

grown_region = FRETImageProcess2.RegionGrowing(img, nucleus_region);

cytoplasm_region = grown_region ~= nucleus_region;

% 細胞質の面積 < 核内の面積 * 1.5 ならNES,CAAXではないと判断

if sum(sum(cytoplasm_region)) < sum(sum(nucleus_region)) * 1.5

cells(i) = 5; % 判定不可

continue;

end

mean_cytoplasm_intensity = mean(img(cytoplasm_region));

% 細胞内のintensityのSDをとる

sd_of_intensity = std(double(img(grown_region)));

% 核内intensityのエッジ（微分）のSDをとる

sd_of_differential_nuc = std(double(differential_img(nucleus_region)));

% sd_array(j) = sd_of_intensity; j = j + 1;

% 以前はSD_NES = 120, SD_CAAX = 60で固定

% 細胞内のSDと核内の微分SDで判断

if (sd_of_intensity > SD_NES) && (mean_cytoplasm_intensity > mean_nucleus_intensity)

if sd_of_differential_nuc > SD_NES * 2

cells(i) = 3; % NESと断定

continue;

end

elseif sd_of_intensity < SD_CAAX && (mean_cytoplasm_intensity >= mean_nucleus_intensity)

cells(i) = 4; % CAAXと断定

continue;

end

% ここまでで断定できなければ判定不可

cells(i) = 5; % 判定不可

end

% sd_array = sd_array(1:j-1);

% 核の大きさやintensityによってノイズをふるいにかける

mean_size = mean(nucleus_size_array(2:length(nucleus_size_array))); % 行列の１番目は背景分(0)なので除く

NLS_mean_intensity = mean(nucleus_intensity_array(cells == 2));

NES_mean_intensity = mean(nucleus_intensity_array(cells == 3));

CAAX_mean_intensity = mean(nucleus_intensity_array(cells == 4));

for i = 2:length(cells)

if (nucleus_size_array(i) < mean_size / 1.7) || (nucleus_size_array(i) > mean_size * 1.7)

cells(i) = 5; % サイズが平均より大きく離れていたら判定不可

end

if (cells(i) == 2) && ((nucleus_intensity_array(i) < NLS_mean_intensity / 1.7) || (nucleus_intensity_array(i) > NLS_mean_intensity * 1.7))

cells(i) = 5; % NLSのうち、intensityがNLSの平均より大きく離れていたら判定不可

elseif (cells(i) == 3) && ((nucleus_intensity_array(i) < NES_mean_intensity / 1.7) || (nucleus_intensity_array(i) > NES_mean_intensity * 1.7))

cells(i) = 5; % NESのうち、intensityがNESの平均より大きく離れていたら判定不可

elseif (cells(i) == 4) && ((nucleus_intensity_array(i) < CAAX_mean_intensity / 1.7) || (nucleus_intensity_array(i) > CAAX_mean_intensity * 1.7))

cells(i) = 5; % CAAXのうち、intensityがCAAXの平均より大きく離れていたら判定不可

end

end

% SortMatrix = SegmentationMatrix;

% for i = 2:length(cells)

% SortMatrix(SortMatrix==i) = cells(i); % NLS,NES,CAAX,判定不可で色分け

% end

% figure, imshow(label2rgb(SortMatrix,'jet','w','shuffle'));

SortArray = cells;

end

% 各細胞のIntensityを出力する領域を取得（ノイズを除いた領域）

function RegionForOutput = RegionForOutput(Imaging_Info, Tiff_Array, Position, WaveName, SegmentationMatrix, SortArray)

% 指定したPosition, WaveNameの画像を全timepoint分読み込み

height = size(Tiff_Array{1,1,1}, 1);

width = size(Tiff_Array{1,1,1}, 2);

imgarray = zeros(height, width, Imaging_Info.timepoints, 'uint16');

for i = 1:Imaging_Info.timepoints

imgarray(:,:,i) = FRETImageProcess2.GetImage(Imaging_Info, Tiff_Array, i, Position, WaveName);

end

% maximum imageとminimum imageを求めてそのsubtractを計算（noiseが浮き彫りになる）

maximum = max(imgarray, [], 3);

minimum = min(imgarray, [], 3);

subtract = maximum - minimum;

% subtractのintensity分布を計算

m = max(max(subtract));

distribution = zeros(m+1,2);

distribution(:,1) = 0:m;

for i = 1:size(subtract,1)

for j = 1:size(subtract,2)

distribution(subtract(i,j)+1,2) = distribution(subtract(i,j)+1,2) + 1;

end

end

distribution(:,1) = 0:m;

% subtract intensityを100ごとに区切り、それぞれの範囲内のピクセル数を比較、比が1.5以下になればそれ以上のsubtract intensityはノイズとする

for i = 1:100:(m-199)

sum1 = sum(distribution(i:i+99, 2));

sum2 = sum(distribution(i+100:i+199, 2));

if sum1 / sum2 < 1.5

break;

end

end

noise_removed_region = subtract <= i+99;

% minimum imageとenabled_regionからRegionForOutputを求める

% イメージング中に細胞が動いたり分裂したりする部分を除く

resultregion = SegmentationMatrix * 0 + 1; % SegmentationMatrixと同じサイズの1の行列(1は背景)

for i = 2:length(SortArray) % SortArray(1)は背景

% 核の領域を取得（enabled_regionに限る）

nucleus_region = SegmentationMatrix == i;

nucleus_region = imdilate(nucleus_region, strel('disk',1)); % 境界部が含まれないので拡張

nucleus_region = nucleus_region & noise_removed_region; % ノイズ領域を除く

if SortArray(i) == 2 % NLS

% 核内のintensityの平均値を求める

mean_int = mean(minimum(nucleus_region));

% 平均値より1.2倍以上低いピクセルは除く（細胞が動いた時などの影響）

outputregion = (minimum .* uint16(nucleus_region)) > mean_int / 1.2;

resultregion(outputregion) = i;

elseif SortArray(i) == 3 % NES

% 核の１～２ピクセル周辺のintensityの平均値を求める

periphery_region = imdilate(nucleus_region, strel('disk',2)) ~= nucleus_region;

mean_int = mean(minimum(periphery_region));

% 平均値より1.2倍以上低いピクセルは除く（細胞が動いた時などの影響）

outputregion = (minimum .* uint16(periphery_region)) > mean_int / 1.2;

resultregion(outputregion) = i;

elseif SortArray(i) == 4 % CAAX

% 核内のintensityの平均値を求める

mean_int = mean(minimum(nucleus_region));

% 平均値より1.2倍以上低いピクセルは除く（細胞が動いた時などの影響）

outputregion = (minimum .* uint16(nucleus_region)) > mean_int / 1.2;

resultregion(outputregion) = i;

end

end

RegionForOutput = resultregion;

end % RegionForOutput

% RegionForOutputに基づいてSortArrayを修正する

function ModifiedSortArray = ModifySortArray(Imaging_Info, Tiff_Array, Position, WaveName, SegmentationMatrix, SortArray, RegionForOutput)

% 指定したPosition, WaveNameの画像を全timepoint分読み込み

height = size(Tiff_Array{1,1,1}, 1);

width = size(Tiff_Array{1,1,1}, 2);

imgarray = zeros(height, width, Imaging_Info.timepoints, 'uint16');

for i = 1:Imaging_Info.timepoints

imgarray(:,:,i) = FRETImageProcess2.GetImage(Imaging_Info, Tiff_Array, i, Position, WaveName);

end

% minimum imageを求める

minimum_image = min(imgarray, [], 3);

new_sortarray = SortArray;

for i = 2:length(SortArray) % SortArray(1)は背景

% １細胞についての出力領域が5ピクセル以下なら除外

if sum(sum(RegionForOutput == i)) <= 5

new_sortarray(i) = 5; % 判定不可

end

% % NLSをNESと誤判定したものを修正する（NLSの隣に明るい死細胞があるとが判定されうる）

% % minimum imageにおいて核内のintensityが出力領域のintensityより1.5倍以上明るければ修正

% nucleus_region = (SegmentationMatrix == i);

% output_region = (RegionForOutput == i);

% if mean(minimum_image(nucleus_region)) > mean(minimum_image(output_region)) * 1.5

% new_sortarray(i) = 5;

% end

end

ModifiedSortArray = new_sortarray;

% SortMatrix = SegmentationMatrix;

% for i = 2:length(new_sortarray)

% SortMatrix(SortMatrix==i) = new_sortarray(i); % NLS,NES,CAAX,判定不可で色分け

% end

% figure, imshow(label2rgb(SortMatrix,'jet','w','shuffle'));

end

% NLSの細胞のintesityを出力((NLSの細胞数)行 x 1列)

function NLSIntensities = OutputNLSIntensity(FRETimage, RegionForOutput, SortArray)

if strcmp(class(FRETimage), 'Tiff')

img = read(FRETimage);

else

img = uint16(round(FRETimage));

end

numofNLS = sum(SortArray == 2); % 2:NLS , 3:NES , 4:CAAX

output = zeros(numofNLS, 1);

j = 0;

for i = 2:length(SortArray)

if SortArray(i) == 2 % NLS

% 核内の平均intensityを出力

j = j + 1;

output(j) = mean(img(RegionForOutput == i));

end

end

NLSIntensities = output;

end

% NESの細胞のintesityを出力((NESの細胞数)行 x 1列)

function NESIntensities = OutputNESIntensity(FRETimage, RegionForOutput, SortArray)

if strcmp(class(FRETimage), 'Tiff')

img = read(FRETimage);

else

img = uint16(round(FRETimage));

end

numofNES = sum(SortArray == 3); % 2:NLS , 3:NES , 4:CAAX

output = zeros(numofNES, 1);

j = 0;

for i = 2:length(SortArray)

if SortArray(i) == 3 % NES

% 核の1～2ピクセル周辺の平均intensityを出力

j = j + 1;

output(j) = mean(img(RegionForOutput == i));

end

end

NESIntensities = output;

end

% CAAXの細胞のintesityを出力((CAAXの細胞数)行 x 1列)

function CAAXIntensities = OutputCAAXIntensity(FRETimage, RegionForOutput, SortArray)

if strcmp(class(FRETimage), 'Tiff')

img = read(FRETimage);

else

img = uint16(round(FRETimage));

end

numofCAAX = sum(SortArray == 4); % 2:NLS , 3:NES , 4:CAAX

output = zeros(numofCAAX, 1);

j = 0;

for i = 2:length(SortArray)

if SortArray(i) == 4 % CAAX

% 核内の平均intensityを出力

j = j + 1;

output(j) = mean(img(RegionForOutput == i));

end

end

CAAXIntensities = output;

end

end % method

end % classdef