demoROI_NP.m

%% 初期化
clear variables
close all;

%% 定数,ディレクトリ
NPWIDTH = 10;

rootDir = 'C:\Users\rysk\Documents\MATLAB\oskd';
roiRoot = [rootDir filesep 'RoiSet'];
meanImgName = [rootDir filesep 'averaged.tif'];

%% read mean Image
meanImg = imread(meanImgName);


%% read ROI
SIZEIMG = size(meanImg);
files = dir([roiRoot filesep '*.roi']);%ROIが入っているフォルダの中のroiと名前のついたファイルをリストアップ（*:ワイルドカード）
for c = 1:length(files)                                                 %いくつROIが入っているか
	%1つ1つのROIを順番に取得
	rois{c} = ReadImageJROI([roiRoot filesep files(c).name]); 
	
	% 中心座標の取得
	meanPos(c,1) = mean(rois{c}.mnCoordinates(:,1));
	meanPos(c,2) = mean(rois{c}.mnCoordinates(:,2));	
end


%% sort by x pos
[~, I] = sort(meanPos(:,1));
rois = rois(I);
meanPos(:,1) = meanPos(I,1);
meanPos(:,2) = meanPos(I,2);


%% imfill & show
maskROIsum= meanImg*0;
for c = 1:length(files)
	maskROItmp = meanImg*0;		
	for r = 1:length(rois{c}.mnCoordinates)
		maskROItmp(rois{c}.mnCoordinates(r, 2), rois{c}.mnCoordinates(r, 1)) = 1;
	end
	maskROItmp   = imfill(maskROItmp);
	maskROIsum   = maskROIsum + maskROItmp;
	rois{c}.Mask = maskROItmp;
end
maskMerge = zeros(SIZEIMG(1), SIZEIMG(2), 3);
maskMerge(:,:,1) = maskROIsum/max(maskROIsum(:));
maskMerge(:,:,2) = meanImg/max(meanImg(:));
maskMerge(:,:,3) = 0;


%% make npROIs
maskNPROIsum = maskROIsum*0;
for c = 1:length(files)
	[npmn, rois{c}.npReg] = ROI_neuropil(rois{c}.Mask, NPWIDTH);
	rois{c}.npmn = npmn{1};
	rois{c}.npReg(maskROIsum > 0) = 0;
	maskNPROIsum = maskNPROIsum + rois{c}.npReg;
end


%% Show Animation
figure; 
subplot(1,2,1);
imagesc(maskMerge); axis image;
for c = 1:length(files)
	subplot(1,2,2);
	imagesc(rois{c}.npReg); axis image; colormap gray;
	pause(1/30)
end
maskMerge(:,:,3) = maskNPROIsum/max(maskNPROIsum(:));
imagesc(maskMerge); axis image;


%% show Summary
subplot(2,2,1)
imagesc(meanImg); axis image;
subplot(2,2,2)
imagesc(maskMerge); axis image;
subplot(2,2,3)
maskOnly = maskMerge;
maskOnly(:,:,2) =maskOnly(:,:,3);
imagesc(maskOnly); axis image
subplot(2,2,4)
imagesc(meanImg); hold on;
axis image;
for c = 1:length(files)
%  	plot(rois{c}.mnCoordinates(:,1), rois{c}.mnCoordinates(:,2), '.r')		
% 	plot(rois{c}.npmn(:,1), rois{c}.npmn(:,2), 'g.')
	text(meanPos(c,1), meanPos(c,2), num2str(c), 'Color', [1 1 1])
end
scatter(meanPos(:,1), meanPos(:,2), '.r');

%% Show ROIs
% figure;
% imagesc(meanImg); hold on;
% axis image;
% colormap gray;
% for c = 1:length(files)
% 	plot(rois{c}.mnCoordinates(:,1), rois{c}.mnCoordinates(:,2), '.r')		
% 	plot(rois{c}.npmn(:,1), rois{c}.npmn(:,2), 'g.')
% 	text(meanPos(c,1), meanPos(c,2), num2str(c), 'Color', [1 1 1])
% end
% scatter(meanPos(:,1), meanPos(:,2), '.r');


%% Show imfill result
% figure; 
% subplot(1,2,1)
% imagesc(meanImg); axis image; colormap gray;


%%
function [NP, NPregion] = ROI_neuropil(ROImask, Psize)
% [NP, NPregion] = ROI_neuropil(ROI, Psize)
%
% making Neuropil-ROI from ROI by bmp2roi
%
% Hakumoto, Takeuchi 2016

L = ROImask;
labels = unique(L);

cellRegion = (L ~= labels(1));

se = strel('disk', Psize);
bw = imdilate(L, se);

NPregion = bw - cellRegion;
NPregion(cellRegion==1) = 0;

for i = 1:length(labels)
	tmpImg = (L == i);
	tmpImgDil = imdilate(tmpImg, se);
	tmpImgDil = tmpImgDil - cellRegion;
	tmpImgDil(tmpImgDil < 1) = 0;
	[NP{i}(:,2) NP{i}(:,1)] = ind2sub(size(L), find(tmpImgDil));
end

end