Skip to content

Instantly share code, notes, and snippets.

@jcasado
Created October 25, 2018 15:30
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jcasado/333c8b4b0625f16a70e226bfdda4ce1b to your computer and use it in GitHub Desktop.
Save jcasado/333c8b4b0625f16a70e226bfdda4ce1b to your computer and use it in GitHub Desktop.
TMA Dearraying
%% 1. Input using Bio-formats
% Read image into 3D array and resize to be easy to handle
filePath = "tma_big.tif";
image = bfOpen3DVolume(filePath);
image= image{1}{1};
nuclei = image(:,:,1);
imagesub = imresize(nuclei,0.02);
%% 2. Morphological operations
% Edge detection
imageEdge = edge(imagesub,'Canny');
% Smooth - close gaps of 3px
mask = imclose(imageEdge,strel('disk',3));
% Smooth - fill blanks from closed shapes
mask = imfill(mask,'holes');
% Smooth - remove objects of less than 50 px
mask = bwareaopen(mask,50);
%% 3. Option A: Find circular objects
Rmin = 9;
Rmax = 14;
s = 0.95;
et = 0.05;
[centers, radii] = imfindcircles(mask,[Rmin Rmax],'ObjectPolarity','bright','Sensitivity',s,'EdgeThreshold',et);
% Plot bright circles in blue
imshowpair(imagesub,mask);
viscircles(centers, radii,'Color','b');
% Or plot centroids
hold on;
plot(centers(:,1), centers(:,2), 'r*');
% Create cropping areas as array of bounding boxes
coreSize = 28;
corners = centers-coreSize/2;
boundingBox = [corners coreSize+zeros(size(corners))];
numCores = length(centers);
%% 3. Option B: Use regionprops
stats=regionprops(mask,'Centroid','Area','BoundingBox');
centroids = cat(1, stats.Centroid);
% Plot centroids over image
imshowpair(imagesub,mask);
hold on;
plot(centroids(:,1), centroids(:,2), 'r*');
% Create cropping areas as array of bounding boxes
boundingBox = stats(iTMA).BoundingBox
numCores = numel(stats)
%% 4. Save cores as Tiff files
myFactor = 50;
outputFolder = "tma_data/cores/"
filePrefix = "TMA"
for iTMA = 1:2
TMA=[];
for iChan=1:size(image,3)
TMA(:,:,iChan) = imcrop(image(:,:,iChan),boundingBox(iTMA,:)*myFactor);
end
bfsave(TMA,[outputFolder filePrefix int2str(iTMA) '.tif'],'BigTiff',true);
end
@Yu-AnChen
Copy link

Yu-AnChen commented Oct 29, 2018

%% read input data
image = uiopen('D:\CAJ_AbVal10_TMA5_7\TMA6\TMA6.ome.tif',1);

%% resize
myFactor = 50;
imagesub = imresize(image,1/myFactor);

%% morphological operations
imageEdge = edge(imagesub,'Canny');
% Smooth - close gaps of 3px
mask = imclose(imageEdge,strel('disk',3));
% Smooth - fill blanks from closed shapes
mask = imfill(mask,'holes');
% Smooth - remove objects of less than 50 px
mask = bwareaopen(mask,50);

%% Detect circles
s = 0.95;
et = 0.05;
Rmin=21; Rmax=31;
[centers, radii] = imfindcircles(mask,[Rmin Rmax],'ObjectPolarity','bright','Sensitivity',s,'EdgeThreshold',et);
% Plot bright circles in blue
imshowpair(imagesub,mask);
viscircles(centers, radii,'Color','b');

%% Sort centers
sorted = sortrows(centers);
new = [];
minDist =10;
newx = 1;

new(1,:) = [sorted(1,:) newx];
for i = 2:size(sorted,1)
    if abs(sorted(i,1) - sorted(i-1,1)) > minDist
        newx = newx + 1;
    end
    new(i,:) = [sorted(i,:) newx];
end
sorted_centers = sortrows(new, [3 2]);

%% Save map
numCores = length(sorted_centers(:,1:2));
RGB = insertText(imagesub, sorted_centers(:,1:2) , 1:numCores,'AnchorPoint','Center');
imshow(RGB)

%% Calculate cropping area
coreSize = 62;
corners = sorted_centers(:,1:2)-coreSize/2;
boundingBox = [corners coreSize+zeros(size(corners))];



r = bfGetReader(char("TMA7.ome.tif"));

%% Crop boundingBox sizes if out of image
boundingBox = uint16(boundingBox*myFactor);
boundingBox(boundingBox == 0) = 1;
maxX = r.getSizeX();
maxY = r.getSizeY();
cornerX = boundingBox(:,1) + boundingBox(:,3);
cornerY = boundingBox(:,2) + boundingBox(:,4);
cornerX(cornerX > maxX) = maxX;
cornerY(cornerY > maxY) = maxY;
boundingBox(:,3) = cornerX - boundingBox(:,1);
boundingBox(:,4) = cornerY - boundingBox(:,2);

outputFolder="D:\CAJ_AbVal10_TMA5_7\TMA7\cores\";
filePrefix="tma7_core_";
tic
for iCore = 1:numCores
    core=[];
    display(iCore);
    for iChan=1:r.getImageCount()
        iBox = num2cell(boundingBox(iCore,:));
        core(:,:,iChan) =bfGetPlane(r, iChan, iBox{:});
    end
    bfsave(core,char(strcat(outputFolder, filePrefix,int2str(iCore),'.tif')),'BigTiff',true);
end
toc


%core(:,:,iChan) =bfGetPlane(r, iChan, uint16(boundingBox(iCore,1)*myFactor), uint16(boundingBox(iCore,2)*myFactor),uint16(boundingBox(iCore,3)*myFactor),uint16(boundingBox(iCore,4)*myFactor) );
        



@jcasado
Copy link
Author

jcasado commented Oct 31, 2018

%% Calculate cropping area
coreSize = 62;
corners = sorted_centers(:,1:2)-coreSize/2;
boundingBox = [corners coreSize+zeros(size(corners))];

numCores = length(sorted_centers(:,1:2));

r = bfGetReader(char("TMA5.ome.tif"));

%% Crop boundingBox sizes if out of image
boundingBox = uint16(boundingBox*myFactor);
maxX = r.getSizeX();
maxY = r.getSizeY();
cornerX = boundingBox(:,1) + boundingBox(:,3);
cornerY = boundingBox(:,2) + boundingBox(:,4);
cornerX(cornerX > maxX) = maxX;
cornerY(cornerY > maxY) = maxY;
boundingBox(:,3) = cornerX - boundingBox(:,1);
boundingBox(:,4) = cornerY - boundingBox(:,2);

outputFolder="D:\CAJ_AbVal10_TMA5_7\TMA5\cores\";
filePrefix="core_";
tic
for iCore = 48:numCores
    core=[];
    display(iCore)
    for iChan=1:r.getImageCount()
        iBox = num2cell(boundingBox(iCore,:));
        core(:,:,iChan) =bfGetPlane(r, iChan, iBox{:});
    end
    bfsave(core,char(strcat(outputFolder, filePrefix,int2str(iCore),'.tif')),'BigTiff',true);
end
toc

@Yu-AnChen
Copy link

Yu-AnChen commented Nov 13, 2018

I added one line in my comment because the minimum of boundingBox corners is 1 instead of 0.

boundingBox = uint16(boundingBox*myFactor);
boundingBox(boundingBox == 0) = 1;

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment