Added nii_nii2objMesh

nii_nii2objMesh.m

@@ -0,0 +1,267 @@
+function outnm = nii_nii2objMesh (fnm, thresh, clusterVox, isSmooth, reduce, outnm, floodfill)
+%convert NIfTI image to mesh, all arguments are optional
+% fnm : nifti image to meshify
+% thresh : air/surface threshold, e.g. 2=voxels darker than 2 are air
+%          n.b. for Altases, set thresh=0: one image will be created per region
+% clusterVox : only clusters of with more voxels than this value survive
+% isSmooth : if true, input is blurred, otherwise raw data
+% reduce : reduction ratio, e.g. 0.2 will decimate 80% of triangles
+% outnm : name for output image, e.g. "mesh.obj"
+% floodfill : if true, remove bubbles inside objects (e.g. ventricles embedded below cortex)
+%Examples
+% nii_nii2objMesh %use GUI
+% nii_nii2objMesh('motor.nii', 3, 30, false, 0.3);
+% nii_nii2objMesh('MNI152_T1_1mm_brain.nii.gz', 5500, 1000000, true, 0.1);
+% nii_nii2objMesh('HarvardOxford-cort-maxprob-thr0-1mm.nii.gz', 0, 1, true, 0.1); %convert smoothed atlas
+% nii_nii2objMesh('natbrainlab.nii.gz', 0, 1, false, 0.1); %convert unsmoothed atlas
+if ~exist('fnm','var')  %no files specified
+    [files,pth] = uigetfile({'*.gz;*.nii;*.hdr;';'*.*'},'Choose images to merge', 'MultiSelect', 'off');
+    fnm = strcat(pth,char(files));
+end
+%load image
+[hdr,img, fnm] = loadNiiSub(fnm);
+img(isnan(img)) = 0;
+%set preferences
+if ~exist('outnm','var') 
+    [pth nm] = spm_fileparts(fnm);
+    outnm = fullfile(pth, [nm '.obj']);
+end;
+if ~exist('thresh','var') thresh = (0.5 * max(img(:))-min(img(:))) + min(img(:)); end;
+if ~exist('clusterVox','var') clusterVox = 32; end;
+if ~exist('isSmooth','var') isSmooth = false; end;
+if  ~exist('floodfill','var') floodfill = false; end;
+fprintf('Image intensity range %g..%g\n', min(img(:)), max(img(:)));
+if ~exist('reduce','var') %provide user interface if prefences not specified
+    prompt = {'Enter threshold (2=only voxels brighter than 2, 0=image is atlas)','Enter minimum cluster size (voxels)','Smooth (blur) input(1=yes, 0=no)','Reduction ratio (0.5=remove 50% of vertices)'};
+    dlg_title = 'Values for thresholding';
+    num_lines = 1;
+    def = {num2str(thresh),num2str(clusterVox), num2str(isSmooth), '0.3'};
+    answer = inputdlg(prompt,dlg_title,num_lines,def);
+    thresh = str2double(answer{1});
+    clusterVox = str2double(answer{2});
+    isSmooth = str2double(answer{3});
+    reduce = str2double(answer{4});
+end;
+if thresh == 0
+    if clusterVox > 1, fprintf('nb: cluster size ignored for Atlases'); end;
+    if floodfill, fprintf('nb: floodfill ignored for Atlases'); end;
+    saveAtlasSub(hdr, img, isSmooth, reduce, outnm);
+    return;
+end
+if isSmooth
+    presmooth = img+0; %+0 forces new matrix
+    spm_smooth(presmooth,img,1,0); %smooth data
+end
+isInvert = thresh < 0;
+if isInvert
+   img = -img;
+   thresh = -thresh;
+end
+[img, nVoxPost] = clusterSub(hdr, img, thresh, clusterVox, isInvert);
+if nVoxPost < 1, fprintf('No voxels exceed threshold\n'); return; end;
+if floodfill
+    img = floodSub(img, thresh, 2); %optional: fill bubbles inside
+end;
+img2meshSub(hdr, img, outnm, reduce, thresh);
+%end nii_nii2obj()
+
+function saveAtlasSub(hdr, img, isSmooth, reduce, outnm)
+if ~(spm_type(hdr.dt,'intt')) %integer output
+    error('Threshold=0 suggests Atlas, but image is not saved as discrete integers.');
+end;
+mn = min(img(:))+1;
+mx = max(img(:));
+if (mn >= mx), error('No variability in input image'); end;
+if isSmooth, fprintf('n.b. each region will be smoothed'); end;
+[pth nm] = spm_fileparts(outnm);
+for i = mn : mx
+    outnm = fullfile(pth, [nm '_' num2str(i) '.obj']);
+    bw = (img == i)+0; %black and white image, +0 converts logical->double
+    if isSmooth
+        presmooth = bw+0; %+0 forces new matrix
+        spm_smooth(presmooth,bw,1,0); %smooth data
+    end
+    vx = sum(bw(:) >= 0.5);
+    if vx > 0
+        fprintf('Region %d has %d voxels\n', i, vx);
+        img2meshSub(hdr, bw, outnm, reduce, 0.5);
+    end;
+end
+%end saveAtlasSub()
+
+function img2meshSub(hdr, img, outnm, reduce, thresh)
+if (max(img(:))< thresh) || (min(img(:)) > thresh) return; end;
+FV = isosurface(img,thresh);
+FV = reducepatch(FV,reduce);
+FV.vertices = FV.vertices(:,[2,1,3]); %isosurface swaps X/Y
+vx = [ FV.vertices ones(size(FV.vertices,1),1)];
+vx = mtimes(hdr.mat,vx')';
+FV.vertices = vx(:,1:3);
+FV.faces = fliplr(FV.faces); %reverse winding
+writeObjSub(FV.vertices,FV.faces, outnm);
+%end img2meshSub()
+
+function [img, nVoxPost] = clusterSub(hdr, img, thresh, clusterVox, isInvert)
+bw=img;
+bw = (bw >= abs(thresh)) * 1.0;
+nVox = sum(bw(:));
+nVoxPost = nVox;
+if (clusterVox < 2) || (nVox < 1), return; end;
+[bw,nCluster] = spm_bwlabel(bw,18);
+nClusterPost = nCluster;
+minClusterVox = nCluster;
+maxClusterVox = 0;
+for i = 1:nCluster
+    vox = sum(bw(:)== i);
+    %fprintf('Cluster %d has %d voxels\n', i,vox);
+    minClusterVox = min(minClusterVox,vox);
+    maxClusterVox = max(maxClusterVox,vox);
+    if vox < clusterVox
+        nVoxPost = nVoxPost - vox;
+        nClusterPost = nClusterPost - 1;
+        bw(bw == i) = 0;        
+    end
+end
+%any voxels that do not survive threshold need to be set just below threshold
+bw = (bw < 1); %mask now inverted and logical: TRUEs are voxels that can be reduced
+maxsubthresh = thresh - (2 * eps);%max(img(img(:) < thresh));
+mask =  (img >= thresh) .* bw; %mask voxels that are over thresh undersized
+img(mask == 1) = maxsubthresh;
+if isInvert, thresh = - thresh; end;
+if minClusterVox ~= maxClusterVox
+    fprintf('Cluster sizes range from %d to %d voxels\n',minClusterVox, maxClusterVox);
+end
+outmm3=prod(abs(hdr.mat(1:3, 1:3)*[1;1;1]));
+fprintf('%d voxels (%d clusters) exceed %g, of which %d (%d) are in clusters larger than %d voxels (%gmm^3)\n', ...
+    nVox, nCluster, thresh, nVoxPost, nClusterPost, clusterVox, clusterVox * outmm3);
+%end clusterSub()
+
+function [hdr,img, fnm] = loadNiiSub(fnm)
+isGz = false;
+if ~exist('spm_fileparts','file'), fprintf('Please install SPM12 or later'); end;
+[pth,nam,ext] = spm_fileparts(fnm);
+if strcmpi(ext,'.gz') %.nii.gz
+    fnm = char(gunzip(fnm));
+    isGz = true;
+elseif strcmpi(ext,'.voi') %.voi ->
+    onam = char(gunzip(fnm));
+    fnm = fullfile(pth, [nam '.nii']);
+    movefile(onam,fnm);
+    isGz = true;
+end;
+hdr = spm_vol(fnm);
+img = spm_read_vols(hdr);
+if isGz, delete(fnm); end;
+%end loadNiiSub()
+function writeObjSub(vertex,face,filename)
+% --- save Face/Vertex data as WaveFront Object format file
+%inputs:
+%	vertex: vertices matrix where cols are xyz and each row a vertix
+%	face: face matrix where cols are xyz and each row is face
+%	fileName: the Wavefront Object file to create
+%notes
+% https://en.wikipedia.org/wiki/Wavefront_.obj_file
+[nF nFd] =size(face);
+[nV nVd] =size(vertex);
+if (nF <1) || (nV <3 || (nFd ~=3) || (nVd ~=3)), warning('Problem with writeObj'); return; end;  %#ok<WNTAG>
+fid = fopen(filename, 'wt');
+fprintf(fid, '# WaveFront Object format image created with MRIcroS\n');
+fprintf(fid, 'v %.12g %.12g %.12g\n', vertex');
+fprintf(fid, 'f %d %d %d\n', (face)');
+fclose(fid);
+%end writeObjSub()
+
+function img = floodSub(img, thresh, smoothVox)
+%similar to http://www.mathworks.com/matlabcentral/fileexchange/12184-floodfill3d
+imgBin = (img < thresh);
+imgBin = double(imgBin);      % In case a logical matrix comes in.
+imgBin(1,:,:) = NaN; %pad LEFT 
+imgBin(end,:) = NaN; %pad RIGHT
+imgBin(:,1,:) = NaN; %pad ANTERIOR 
+imgBin(:,end,:) = NaN; %pad POSTERIOR
+imgBin(:,:,1) = NaN; %pad INFERIOR 
+imgBin(:,:,end) = NaN; %pad SUPERIOR
+imgBin = floodFill3DSub(imgBin, [2,2,2]);
+%fill bubbles
+mx = max(img(:));
+img(isfinite(imgBin))= mx;
+%next: optional - blur to feather edges - only useful if marching cubes uses subvoxel smoothing
+maskIn = img + 0.0;
+if exist('smoothVox', 'var')
+    spm_smooth(maskIn, img, smoothVox); %feather the edges a lot: weak blur
+end
+%end floodSub()
+
+%Copyrights for Matlab Central Code
+% stlwriteSub Copyright (c) 2015, Sven Holcombe
+% floodFill3DSub Copyright (c) 2006,  F Dinath
+% 
+% All rights reserved.
+% 
+% Redistribution and use in source and binary forms, with or without
+% modification, are permitted provided that the following conditions are
+% met:
+% 
+%     * Redistributions of source code must retain the above copyright
+%       notice, this list of conditions and the following disclaimer.
+%     * Redistributions in binary form must reproduce the above copyright
+%       notice, this list of conditions and the following disclaimer in
+%       the documentation and/or other materials provided with the distribution
+% 
+% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+% POSSIBILITY OF SUCH DAMAGE.
+
+function [A] = floodFill3DSub(A, point)
+% http://www.mathworks.com/matlabcentral/fileexchange/12184-floodfill3d
+% [B] = FloodFill3D(A, slice);
+% This program flood fills a 6-connected 3D region. The input matrix MUST
+% be a binary image. The user will select a seed (point) in the matrix to
+% initiate the flood fill. You must specify the matrix slice in which you
+% wish to place the seed.
+% 
+% A = binary matrix
+% slice = a chosen slice in the matrix where you wish to place the seed.
+%
+% Enjoy,
+% F. Dinath
+if A(point(1), point(2), point(3));
+    A(point(1), point(2), point(3)) = NaN;
+    a{1} = sub2ind(size(A), point(1), point(2), point(3));
+    i = 1;
+    while 1
+        i = i+1;
+        a{i} = []; %#ok<AGROW>
+        [x, y, z] = ind2sub(size(A), a{i-1});
+        ob = nonzeros((A(sub2ind(size(A), x, y, z-1)) == 1).*sub2ind(size(A), x, y, z-1));
+        A(ob) = NaN;
+        a{i} = [a{i} ob']; %#ok<AGROW>
+        ob = nonzeros((A(sub2ind(size(A), x, y, z+1)) == 1).*sub2ind(size(A), x, y, z+1));
+        A(ob) = NaN;
+        a{i} = [a{i} ob']; %#ok<AGROW>
+        ob = nonzeros((A(sub2ind(size(A), x-1, y, z)) == 1).*sub2ind(size(A), x-1, y, z));
+        A(ob) = NaN;
+        a{i} = [a{i} ob']; %#ok<AGROW>
+        ob = nonzeros((A(sub2ind(size(A), x+1, y, z)) == 1).*sub2ind(size(A), x+1, y, z));
+        A(ob) = NaN;
+        a{i} = [a{i} ob']; %#ok<AGROW>
+        ob = nonzeros((A(sub2ind(size(A), x, y-1, z)) == 1).*sub2ind(size(A), x, y-1, z));
+        A(ob) = NaN;
+        a{i} = [a{i} ob']; %#ok<AGROW>
+        ob = nonzeros((A(sub2ind(size(A), x, y+1, z)) == 1).*sub2ind(size(A), x, y+1, z));
+        A(ob) = NaN;
+        a{i} = [a{i} ob']; %#ok<AGROW>
+        if isempty(a{i});
+            break;
+        end
+    end
+end
+%end floodFill3DSub()