latest version of matlab code

Author: eminizer

calcIndTr.m

@@ -0,0 +1,54 @@
+% This function takes the registered locations, ccl, and the original
+% locations, cl, and calculates the local transformations for these points.
+% This is done in a weighted manner by taking each point in turn and
+% finding its 10 nearest neighbors. Then, we draw rays from the point of
+% interest to each of its nearest neighbors and divide that ray into 10
+% segments.We then use all of these points to calculate the local affine
+% transformation. By interpolating points along the rays between points, we
+% create a higher density of points nearest our point of interest.
+%
+% Author: Joshua Doyle, 08/2019 - e.rational.2.718@gmail.com
+
+
+function T = calcIndTr(ccl,cl,c,conf,np,nn)
+
+T = zeros(length(np(conf>c)),6);
+[Idx,Dx] = knnsearch(ccl(conf>c,1:2),ccl(conf>c,1:2),'K',nn);
+% ind = (conf>c);
+dex = find(conf>c);
+Dxref = zeros(length(dex),size(Idx,2));
+Rreg = zeros(nn+1,size(Idx,2));
+Rref = zeros(nn+1,size(Idx,2));
+REG = [];%zeros(size(Rreg,1)*size(Rreg,2));
+REF = [];
+
+for i1 = 1:length(np(conf>c))
+    Dxref(i1,:) = ((cl(np(dex(Idx(i1,:))),1)-cl(np(dex(Idx(i1,1))),1)).^2 +...
+                   (cl(np(dex(Idx(i1,:))),2)-cl(np(dex(Idx(i1,1))),2)).^2).^0.5';
+    thetareg = atan2(ccl(dex(Idx(i1,:)),2)-ccl(dex(Idx(i1,1)),2),...
+                     ccl(dex(Idx(i1,:)),1)-ccl(dex(Idx(i1,1)),1));
+    thetaref = atan2(cl(np(dex(Idx(i1,:))),2)-cl(np(dex(Idx(i1,1))),2),...
+                     cl(np(dex(Idx(i1,:))),1)-cl(np(dex(Idx(i1,1))),1));
+    for i2 = 2:length(thetareg)
+    Rreg(:,i2) = 0:Dx(i1,i2)/nn:Dx(i1,i2);
+    Rref(:,i2) = 0:Dxref(i1,i2)/nn:Dxref(i1,i2);
+    end
+    regx = reshape(Rreg,[],1).*cos(reshape(repmat(thetareg',size(Rreg,1),1),[],1));
+    regy = reshape(Rreg,[],1).*sin(reshape(repmat(thetareg',size(Rreg,1),1),[],1));
+    refx = reshape(Rref,[],1).*cos(reshape(repmat(thetaref',size(Rref,1),1),[],1));
+    refy = reshape(Rref,[],1).*sin(reshape(repmat(thetaref',size(Rref,1),1),[],1));
+    
+    ref = [refx refy] + [cl(np(dex(Idx(i1,1))),1) cl(np(dex(Idx(i1,1))),2)];
+    reg = [regx regy] + [ccl(dex(Idx(i1,1)),1) ccl(dex(Idx(i1,1)),2)];
+  
+    [ab,cd,e,f] = newtform(ref,reg);
+    T(i1,:) = [ab' cd' [e;f]'];
+    REG = [REG;reg];
+    REF = [REF;ref];
+%     
+%         figure,scatter(ref(:,1),ref(:,2)),hold on, scatter(reg(:,1),reg(:,2),'x')
+%     hold on, scatter(cl(np(dex(Idx(i1,1))),1),cl(np(dex(Idx(i1,1))),2),'x')
+%     
+end
+
+end

checkVars.m

@@ -1,77 +1,78 @@
-function [regData,ihcrs,ifpaths,b,h,w,scale,ihcscale] = checkVars(ihc,ifdir)
+%% checkVars
+%% Created by: Joshua Doyle 
+%% Edited by: Benjamin Green
+%% -----------------------------------------------------------
+%% Description
+% part of the registration protocol for slides from different microscopes
+% reads in the images to memory and gets image parameters
+% -- Possible Improvement in readYCBCr2RGB for parloop reading in subset of
+% main IHC image to improve read time -- use imread 'Region Property'
+%% -----------------------------------------------------------
+function [regData,ihcrs,ifpaths,b,h,w,scale,ihcscale] ...
+    = checkVars(ihc,ifdir)
+%
 % pixels/u - Vectra scale (mft provided)
+%
 scale = 1.9981018;
+%
 % Number of bands in IF images
+%
 ncomponents = 8;
+%
 % pixels/u - Hamamatsu scale (mft provided)
+%
 ihcscale = 2.173913;
-
-if ~exist('regData', 'var')
-    %myDir = uigetdir('Load Componenet Tiff Directory'); %gets directory
-    ifpaths = dir(fullfile(ifdir,'*component_data.tif'));
-    
-    regData = struct('filename',{{}},'micronLoc',{{}},'regLoc',{{}});
-    for i1 = 1:length(ifpaths)
-        baseFileName = ifpaths(i1).name;
-%         B = strsplit(baseFileName,'_');
-%         if length(B) > 3
-%             if strcmp(B{4},'binary')
-%                 imstr.binaryfile = baseFileName;
-%             end
-%             if strcmp(B{4},'composite')
-%                 imstr.componentfile = baseFileName;
-%             end
-%             if strcmp(B{4},'phenotype')
-%                 imstr.phenofile = baseFileName;
-%             end
-%             
-%             if strcmp(B{8},'component')
-%                 imstr.compositefile = baseFileName;
-                regData.filename{i1} = baseFileName;
-                C = strsplit(baseFileName,{'[',',',']'});
-                % TODO
-                % add if wholeslide || if TMA
-                regData.micronLoc{i1} = [str2double(C{2}) str2double(C{3})];
-                %regData.micronLoc{i1} = [str2double(C{6}) str2double(C{7})];
-%             end
-%         end
-    end
-    b = cast(reshape(cell2mat(regData.micronLoc),2,[]),'uint32')';
+%
+% get if image names
+%
+ifpaths = dir(fullfile(ifdir, '*component_data.tif'));
+if isempty(ifpaths)
+    msg = ['Error: could not find IF images in - ',ifpaths];
+    error(msg)
 end
-
-% if ~exist('miData','var')
-%     load('miData');
-%     mirs = reshape(miData.mmi,[],1);
-%     midata2 = zeros(length(mirs),1);
-%     for i1 = 1:length(mirs)
-%         midata2(i1) = sum(mirs{i1}(:));
-%     end
-%     [ma,mind] = max(midata2);
-%     [I,J] = ind2sub([9,9],mind);
-% end
-
-if ~exist('ihcrs','var')
-    
-    %[ihcfile,ihcpath] = uigetfile('*/*.tif','Select a full slide image for registration');
-    %tobj = Tiff([ihcpath,ihcfile],'r');
-    tobj = Tiff(ihc,'r');
-    origh = getTag(tobj,'ImageLength');
-    origw = getTag(tobj,'ImageWidth');
-    ihcrs = imresize(readYCbCr2RGB(tobj),[origh*scale/ihcscale,origw*scale/ihcscale]);
-    
-    %ihcrs = imresize(readYCbCr2RGB(tobj),[origh*scale/miData.param.scale{I,J}(2),origw*scale/miData.param.scale{I,J}(1)]);
-    
-    %     tobjATP = Tiff('Z:\stitch\Hamamatsu Channel Images\M41-ATPase.tif','r');
-    %     tobjHematox = Tiff('Z:\stitch\Hamamatsu Channel Images\M41-Hematoxyln.tif','r');
-    %     origh = getTag(tobjATP,'ImageLength');
-    %     origw = getTag(tobjATP,'ImageWidth');
-    %     ihcrs(:,:,1) = imresize(read(tobjATP),[origh*scale/miData.param.scale{I,J}(2),origw*scale/miData.param.scale{I,J}(1)]);
-    %     ihcrs(:,:,2) = imresize(read(tobjHematox),[origh*scale/miData.param.scale{I,J}(2),origw*scale/miData.param.scale{I,J}(1)]);
+%
+regData = struct('filename', {{}}, 'micronLoc', {{}}, 'regLoc', {{}});
+for i1 = 1:length(ifpaths)
+    baseFileName = ifpaths(i1).name;
+    regData.filename{i1} = baseFileName;
+    C = strsplit(baseFileName,{'[',',',']'});
+    regData.micronLoc{i1} = [str2double(C{2}) str2double(C{3})];
 end
-
+%
+% variable for registered locations
+%
+b = cast(reshape(cell2mat(regData.micronLoc),2,[]),'uint32')';
+%
+% get IHC image and size params
+%
+if ~exist(ihc,'file')
+    msg = ['Error: could not find IHC image - ',ihc];
+    error(msg)
+end
+%
+fprintf(1, 'Reading %s\n', ihc);
+%
+tic
+imf = imfinfo(ihc);
+level = find(ismember({imf(:).ImageDescription}, 'x20 z0'));
+origh = imf(level).Height;
+origw = imf(level).Width;
+%
+img = imread(ihc,level);
+ihcrs = imresize(img,...
+    [origh*scale/ihcscale,origw*scale/ihcscale]);
+clear img
+fprintf('           ');
+toc
+%
+% get IF size params
+%
 baseFileName = ifpaths(1).name;
 fullFileName = fullfile(ifdir, baseFileName);
-fprintf(1, 'Now reading %s\n', fullFileName);
-tobj = Tiff(fullFileName);
-img = makeTifIm(tobj,ncomponents);
-[h,w,bands] = size(img);
+img = imfinfo(fullFileName);
+h = img.Height;
+w = img.Width;
+regData.h = h; %hpf height
+regData.w = w; %hpf width
+%
+end

cropandwritehpfs.m

@@ -0,0 +1,22 @@
+function cropandwritehpfs(rd,wsimf,wsimdir)
+
+if ~exist('wsimdir','dir')
+    mkdir(wsimdir);
+end
+S = strsplit(wsimdir,'\');
+suff = S{end-2};
+
+for i1 = 1:length(rd.filename)
+   
+   S = strsplit(rd.filename{i1}, '_');
+   fid = strjoin(S(1:end-2),'_');
+   fid = [fid,'_',suff];
+   
+   imr = wsimf(rd.regLoc(i1,2):rd.regLoc(i1,2)+rd.h-1,...
+               rd.regLoc(i1,1):rd.regLoc(i1,1)+rd.w-1,:);
+   
+   imwrite(imr,[wsimdir,'\',fid,'.tif']);
+   
+end
+save([wsimdir,'\rd'],'rd')
+end

doRegSet.m

@@ -0,0 +1,51 @@
+% topdir = '/Volumes/e$/Clinical_Specimen/';
+topdir = '\\bki04\e$\Clinical_Specimen\';
+A = dir(fullfile(topdir, 'M*'));
+specID = {A(:).name};%{'M1_1','M2_3'};%M9_1','M10_2','M11_1','M12_1'};
+% specID = {'M41_1'};
+for i1 = 1%round(length(specID)/2)+1:length(specID)
+    
+    tic;
+    he = ['\\bki04\e$\Clinical_Specimen\',specID{i1},'\HE\',specID{i1},'.ndpi'];
+    ihc = ['\\bki04\e$\Clinical_Specimen\',specID{i1},'\IHC\',specID{i1},'.ndpi'];
+    ifdir = ['\\bki04\e$\Clinical_Specimen\',specID{i1},'\inform_data\Component_Tiffs'];
+    if_hpfs = dir(fullfile(ifdir,'*component_data.tif'));
+    
+    if exist(ihc,'file')
+        S = strsplit(ihc,'\');
+        ihc_dir = ['\',strjoin(S(1:6),'\'),'\HPFs\'];
+        ihc_hpfs = dir(fullfile(ihc_dir, '*IHC.tif'));
+        if length(ihc_hpfs) == length(if_hpfs)
+            try
+                regFullSlide(ihc,ifdir);
+            catch ME
+                warning(['MATLAB returned the following error identifier - ',ME.identifier]);
+                warning(['MATLAB returned the following error message - ',ME.message]);
+                A(i1).error = 1;
+            end
+%         else
+%             fprintf(['Directory ',ihc_dir,' already exists with all HPFs. Delete directory and try again. \n'])
+            
+        end
+    end
+    toc;
+    tic;
+    if exist(he,'file')
+        S = strsplit(he,'\');
+        he_dir = ['\',strjoin(S(1:6),'\'),'\HPFs\'];
+        he_hpfs = dir(fullfile(he_dir, '*HE.tif'));
+        if length(he_hpfs) == length(if_hpfs)
+            try
+                regFullSlide(he,ifdir);
+            catch ME
+                warning(['MATLAB returned the following error identifier - ',ME.identifier]);
+                warning(['MATLAB returned the following error message - ',ME.message]);
+                A(i1).error = 1;
+            end
+%         else
+%             fprintf(['Directory ',he_dir,' already exists with all HPFs. Delete directory and try again. \n'])
+            
+        end
+    end
+    toc;
+end

emptyTask.m

@@ -0,0 +1,4 @@
+    function errStr = emptyTask(i1)
+        errStr = ['Job ',int2str(i1),' has no fields.'];
+        
+    end