Add demo.

Author: mu40

.gitignore

@@ -0,0 +1,7 @@
+.dropbox
+.DS_Store
+Icon*
+out.*
+test.*
+test
+tmp

README.md

@@ -0,0 +1,67 @@
+% Reslice brain such that it's aligned with the anatomical planes, based on
+% the coordinates of its centre and those of the eyes (voxel coordinates).
+function out = alignbrain(mri, fetaltoras, outfile)
+
+if nargin() < 3
+    outfile = [];
+end
+
+if ischar(mri)
+    mri = MRIread(mri); % FreeSurfer.
+    mri.vol = permute(mri.vol, [2 1 3 4]);
+end
+voxtoras = mri.vox2ras1;
+dat = mri.vol(:,:,:,1);
+dat = dat ./ gaussblur(dat, 0.33*size(dat));
+
+% Resample in anatomical coordinates.
+fov = -70:1:70;
+[xf,yf,zf] = ndgrid(fov, fov, fov);
+fetaltovox = voxtoras \ fetaltoras;
+[xi,yi,zi] = coords(fetaltovox, xf, yf, zf);
+res = interpn(dat, xi, yi, zi, 'linear', 0);
+
+% Plot central slices.
+reset(clf());
+figure(gcf());
+colormap gray
+
+hs(1) = subplot(2,2,1);
+im = squeeze(res(:,round(end/2),:));
+im = flip(im', 1);
+imagesc(stretchcon(im));
+title('Coronal');
+
+hs(2) = subplot(2,2,2);
+im = squeeze(res(round(end/2),:,:));
+im = flip(flip(im', 1),2);
+imagesc(stretchcon(im));
+title('Sagittal');
+
+hs(3) = subplot(2,2,3);
+im = squeeze(res(:,:,round(end/2)));
+im = flip(im', 1);
+imagesc(stretchcon(im));
+title('Axial');
+
+hl = [];
+for i = 1:3
+    hl(end+1) = line(hs(i), xlim(), [1 1].*mean(ylim())); %#ok<*AGROW>
+    hl(end+1) = line(hs(i), [1 1].*mean(xlim()), ylim());
+end
+axis(hs, 'image', 'off');
+set(hl, 'color', 'y', 'linewidth', 1.5);
+drawnow();
+
+out = struct();
+out.vol = permute(res, [2 1 3 4]);
+voxsize = [xf(2,1,1)-xf(1,1,1) yf(1,2,1)-yf(1,1,1) zf(1,1,2)-zf(1,1,1)];
+out.volres = voxsize;
+scaling = diag([voxsize 1]);
+subone = [eye(4,3) [-1 -1 -1 1]'];
+voxtofetal = [eye(4,3) [xf(1) yf(1) zf(1) 1]'] * scaling * subone;
+out.vox2ras1 = fetaltoras * voxtofetal;
+out.vox2ras0 = vox2ras_1to0(out.vox2ras1);
+if ~isempty(outfile)
+    MRIwrite(out, outfile);
+end

alignbrain.m

@@ -0,0 +1,19 @@
+% Return prior knowledge of fetal anatomy in mm for specific GA in wk.
+% Snijders et al. (1994) and Robinson et al. (2008):
+% https://www.pedrad.org/Portals/5/Fetal%20MRI/fetal%20eyes.pdf
+function [ofd,bpd,odiam,odist] = anatomy(ga)
+
+if ga < 14 || ga > 39
+    error('ERROR: GA out of range 14-39 weeks.');
+end
+
+odiam = 0.47*ga - 0.7;
+bod = 1.47*ga + 1.76;
+odist = bod - odiam;
+xga = 14:39;
+bpd = [28 31 34 36 39 42 45 48 51 54 57 60 63 66 69 72 74 77 79 81 83 ...
+    85 86 87 88 89];
+ofd = [35 39 42 46 50 54 57 61 65 69 73 77 81 84 87 91 94 96 99 101 103 ...
+    105 106 107 107 107];
+bpd = interp1(xga, bpd, ga);
+ofd = interp1(xga, ofd, ga);

anatomy.m

@@ -0,0 +1,19 @@
+% Iteratively add voxels to each slice of a 3D mask that are included on
+% the closest non-zero slices on both sides.
+function volmask = closegaps(volmask)
+
+tmp = volmask;
+slcind = 1:size(volmask, 3);
+change = 1;
+while change
+    for i = slcind
+        masksize = squeeze(sum(sum(tmp)))';
+        j = find(slcind<i & masksize>0, 1, 'last');
+        k = find(slcind>i & masksize>0, 1, 'first');
+        if not(isempty(j)) && not(isempty(k))
+            tmp(:,:,i) = tmp(:,:,i) | ( tmp(:,:,j) & tmp(:,:,k) );
+        end
+    end
+    change = not(isequal(tmp, volmask));
+    volmask = tmp;
+end

closegaps.m

@@ -0,0 +1,7 @@
+% Multiply the 4-by-4 matrix M by the vector [x1, x2, x3, 1].
+function [y1,y2,y3,y4] = coords(M, x1, x2, x3)
+
+y1 = M(1,1)*x1 + M(1,2)*x2 + M(1,3)*x3 + M(1,4);
+y2 = M(2,1)*x1 + M(2,2)*x2 + M(2,3)*x3 + M(2,4);
+y3 = M(3,1)*x1 + M(3,2)*x2 + M(3,3)*x3 + M(3,4);
+y4 = M(4,1)*x1 + M(4,2)*x2 + M(4,3)*x3 + M(4,4);

coords.m

@@ -0,0 +1,16 @@
+% Remove protruding voxels (once) from each slice of a 3D mask if they are not
+% included on the closest slice on at least one side.
+function out = cutprot(volmask)
+
+numvox = squeeze(sum(sum(volmask)))';
+ind = find(numvox);
+out = volmask;
+for i = ind
+    if i == ind(1)
+        out(:,:,i) = volmask(:,:,i) & volmask(:,:,i+1);
+    elseif i == ind(end)
+        out(:,:,i) = volmask(:,:,i) & volmask(:,:,i-1);
+    else
+        out(:,:,i) = volmask(:,:,i) & (volmask(:,:,i-1) | volmask(:,:,i+1));
+    end
+end

cutprot.m

@@ -0,0 +1,39 @@
+% Demo visualizing fetal-brain localization and orientation estimation.
+
+% Plot stages: 1 brain localization, 2 masking, 3 eye detection (0 is off).
+doplot = struct();
+doplot.clust1 = 1;
+doplot.line1 = 1;
+doplot.mser1 = 1;
+doplot.mask1 = 1;
+doplot.mser2 = 1;
+doplot.line2 = 1;
+doplot.poly2 = 1;
+doplot.mask2 = 1;
+doplot.mser3 = 1;
+doplot.mask3 = 1;
+doplot.geom3 = 1;
+
+% Load data paths and gestational age.
+addpath freesurfer
+[data,ga] = loaddata();
+wait('This demo showcases the different stages of fetal-brain', ...
+    'localization and orientation estimation. You will be prompted for', ...
+    'input before each step is run.');
+
+% Import image.
+imnum = 1;
+mri = MRIread(data{imnum}); % FreeSurfer.
+mri.vol = permute(mri.vol, [2 1 3 4]);
+wait('Slices of input image', ['"' data{imnum} '"']);
+showmask(stretchcon(mri.vol));
+
+% Detect landmarks and reconstruct geometry.
+[bcenvox,ecenvox,bmask] = landmarks(mri, ga(imnum), doplot);
+fetaltoras = recongeo(mri.vox2ras1, bcenvox, ecenvox, bmask);
+
+% Resample in anatomical space and save for inspection. 
+outfile = 'out.nii';
+wait('Resampling in anatomical space');
+alignbrain(mri, fetaltoras, outfile);
+fprintf('Done\n');

data/ep2d_01.nii.gz

@@ -0,0 +1,88 @@
+% Fit an ellipse/ellipsoid to 2D/3D image blobs.
+function [bcen,semilen,rotmat,ellmask,fill,over] = fitellipse(mask, scalemask)
+
+if nargin() < 2
+    scalemask = 1;
+end
+
+iscellinput = iscell(mask);
+if ~iscellinput
+    mask = {mask};
+end
+sz = size(mask{1});
+numdim = numel(sz);
+if numel(sz) == 2
+    sz(3) = 1;
+end
+[x,y,z] = ndgrid(1:sz(1), 1:sz(2), 1:sz(3));
+
+nummask = numel(mask);
+bcen = zeros(nummask, numdim);
+semilen = zeros(nummask, numdim);
+rotmat = cell(nummask, 1);
+ellmask = cell(nummask, 1);
+fill = zeros(nummask, 1);
+over = zeros(nummask, 1);
+for i = 1:nummask
+    % Centre of mass and barycentric coordinates.
+    bw = mask{i};
+    ind = find(bw);
+    m000 = numel(ind);
+    cen = [sum(x(ind)) sum(y(ind)) sum(z(ind))] / m000;
+    xb = x(ind) - cen(1);
+    yb = y(ind) - cen(2);
+    zb = z(ind) - cen(3);
+    
+    % Second moments - weighted by intensity, which is 1.
+    m200 = sum(xb.^2);
+    m020 = sum(yb.^2);
+    m002 = sum(zb.^2);
+    m110 = sum(xb.*yb);
+    m011 = sum(yb.*zb);
+    m101 = sum(xb.*zb);
+    mi = [m200 m110 m101; m110 m020 m011; m101 m011 m002] / m000;
+    
+    % Semi axes lengths and rotation matrix from eigendecomposition.
+    [evec,eval] = eig(mi, 'vector');
+    [eval,order] = sort(eval, 'descend');
+    rot = evec(:,order);
+    semi = 2 * sqrt(eval);
+    
+    % Mask of ellipse.
+    halflen = semi(1);
+    lower = max(1, floor(cen-halflen));
+    upper = min(sz, ceil(cen+halflen));
+    indx = lower(1):upper(1);
+    indy = lower(2):upper(2);
+    indz = lower(3):upper(3);
+    xb = x(indx,indy,indz) - cen(1);
+    yb = y(indx,indy,indz) - cen(2);
+    zb = z(indx,indy,indz) - cen(3);
+    xr = rot(1,1)*xb + rot(2,1)*yb + rot(3,1)*zb;
+    yr = rot(1,2)*xb + rot(2,2)*yb + rot(3,2)*zb;
+    zr = rot(1,3)*xb + rot(2,3)*yb + rot(3,3)*zb;
+    ell = false(sz);
+    scale = semi .* scalemask;
+    scale(scale == 0) = 1; % Avoid 2D division by zero.
+    isell = (xr./scale(1)).^2 + (yr./scale(2)).^2 + (zr./scale(3)).^2 < 1;
+    ell(indx,indy,indz) = isell;
+    
+    if sz(3) == 1
+        cen = cen(1:2);
+        semi = semi(1:2);
+        rot = rot(1:2,1:2);
+    end
+    rot = sign(det(rot)) * rot; % Proper rotation matrix.
+    
+    bcen(i,:) = cen;
+    semilen(i,:) = semi;
+    rotmat{i} = rot;
+    ellmask{i} = ell;
+    fill(i) = sum(bw(ind) & ell(ind)) / nnz(isell);
+    over(i) = sum(~ell(ind) & bw(ind)) / m000;
+end
+
+if ~iscellinput
+    rotmat = rotmat{1};
+    ellmask = ellmask{1};
+end

data/ep2d_02.nii.gz

@@ -0,0 +1,16 @@
+% Fit line to ND point cloud: return centroid, direction and point distances.
+% http://www.ctralie.com/Teaching/COMPSCI290/Lectures/10_PCA/slides.pdf
+function [centroid, dir, dist] = fitline(points)
+
+centroid = mean(points, 1);
+points = points - centroid;
+
+covmat = points' * points;
+[evec,eval] = eig(covmat, 'vector');
+
+[~,order] = sort(abs(eval), 'descend');
+evec = evec(:,order);
+dir = evec(:,1)';
+
+p = points - centroid;
+dist = real(sqrt(sum(p.*p,2)-(sum(p.*dir,2)).^2));