renaming directory

code/morphology/CalcInt.m

@@ -0,0 +1,47 @@
+function [results,Int]=CalcInt(MGstack,MGvolume,res,PunctaWanted,Pvolume_total,spMarker)
+
+
+MGPuncta=  MGstack + PunctaWanted;                        
+MGPuncta(MGPuncta == 1) = 0;                                      
+MGPuncta(MGPuncta == 2) = 1;     
+
+PunctaWanted_labeled  =  bwlabeln(PunctaWanted,8);  
+MGPuncta =  MGPuncta .* PunctaWanted_labeled; 
+[MGPuncta_size, ~]  = histcounts(MGPuncta,1 : max(MGPuncta(:)));
+MGPuncta_size(MGPuncta_size<=1)=0;
+
+[PunctaWanted_size, ~]  = histcounts(PunctaWanted_labeled ,1 : max(PunctaWanted_labeled(:)));
+PunctaWanted_size(PunctaWanted_size<=1)=0;
+
+
+
+%if parts extending Microglias
+
+if spMarker<1
+    MGPuncta_size(MGPuncta_size<=res.pixel_width_um/10)=0;
+    PunctaWanted_size(PunctaWanted_size<=res.pixel_width_um/10)=0;
+    
+    [Psize, ~]  =  histcounts(PunctaWanted_labeled, 1 : max(PunctaWanted_labeled(:)));
+    MGPuncta_size(end : length(Psize)) = 0;
+    size_ratio = MGPuncta_size ./ Psize;
+    MGPuncta_size = find(size_ratio == 1);
+    MGPuncta = ismember(MGPuncta,MGPuncta_size);
+end
+
+%% count number of inclusions and number of pixels in the intersection in the desired microglia
+
+results.number_inclusions  = nnz(MGPuncta_size);        
+results.number_inclusions_per_um = results.number_inclusions / MGvolume;
+results.volume_inclusions_um = sum(sum(sum(MGPuncta)))* (res.pixel_width_um*res.pixel_height_um*res.pixel_depth_um);
+results.volume_inclusions_per_um = results.volume_inclusions_um / MGvolume;
+results.microglia_volume=MGvolume;
+results.number_inclusions_total=nnz(PunctaWanted_size); 
+results.volume_inclusions_total=Pvolume_total;
+
+Int=MGPuncta;
+
+
+
+
+
+

code/morphology/MasterPruning.m

@@ -0,0 +1,136 @@
+%Jastyn Poepplau 15.10.2021
+close all
+clear all
+
+%% EXPERIMENT SET
+%% Load basic information for groups, the information is stored in an excel file. An example template is uploaded in the repository (MG_pruning).
+Nexperiment=[];
+Path='Q:\Personal\Anne\FromJastyn\Matlab\MG_Pruning.xlsx';%change path->where you store the excel file
+
+experiments=get_experiment_list(Path,Nexperiment);
+
+load_main='Q:\Personal\Anne\FromJastyn\PruningStacks\';
+results_main='Q:\Personal\Anne\FromJastyn\Results\';
+
+groups={''}; %Load groups you want to analyze, ExpType in your excel file
+
+Num1=0;
+for f1=1:size(experiments,1)
+    if ~isempty(experiments(f1).Basic.General.P_age)%checks for entrys
+        if any(strcmp(experiments(f1).Basic.General.ExpType,groups))%compares ExpType
+            Num1=Num1+1;
+            experiments1(Num1,1)=experiments(f1);
+            groupsExp{Num1}=experiments(f1).Basic.General.ExpType;
+        end
+    end
+end
+experiments=experiments1;
+clear experiments1
+
+%% Custom Modis
+
+%Script for analyzing overlaying punkta of one or two seperate channels with main stack of complete cell
+%Stacks are preprocessed in imageJ. Ask for fiji scripts or use record function to find your parameters
+%Matlab expects a raw stack.tif for the main stack and thresholded stack.tif for the punkta 
+%Script can analyze up to 3 cells from the same stack. Specify tihis in the max column in your excel file
+
+loadMainStack='MG'; %name of your stack for which you want to have the whole extracted cell, i.e. of microglia cell, but could also work with neurons
+load1='H1'; %name of your first stack for overlay
+load2='GAD'; %name of your second stack for overlay
+
+NumChans=3; %number of channels you want to analyze, works for 2 or 3
+spMarker=0; %Do you have a second specific marker for your main stack? i.e. iba1 +CD68 for microglia
+
+
+recalc=true; %Do you want to calculate all cells again? ->true; Or only the new ones? ->false
+
+%% Batch analysis
+for f1=1:size(experiments,1)
+    Animal=num2str(experiments(f1).Basic.General.animalID);
+    Cell=num2str(experiments(f1).MicrogliaAnalysis.General.Cell);
+    FileCell=num2str(experiments(f1).MicrogliaAnalysis.General.FileCell);
+    max=experiments(f1).MicrogliaAnalysis.General.max;% Set to 2 if you have two cells you want to analyze in the same stack
+    range=str2num(experiments(f1).MicrogliaAnalysis.General.range);
+    interval=experiments(f1).MicrogliaAnalysis.General.interval;
+    savedir=strcat(results_main,Animal,'\',Cell,'\');
+    
+    if (~exist(strcat(savedir,'results.mat'),'file') || recalc)
+        fprintf(['---------\nanalyzing animal ' Animal  ' cell number ' Cell '\n']);
+        
+        %Preprocess microglia stack and extract cell for analysis
+        [Mainstack,MainVolume,res]=PreProMainStack(load_main,Animal,FileCell,max,range,interval,loadMainStack);
+        
+        %Preprocess puncta
+        if NumChans==3
+            [PunctaWanted,Pvolume,Punkta1Stack,Punkta2Stack]=PreProPunkta_2chan(load_main,Animal,FileCell,range,res,Mainstack,load1,load2);
+        elseif NumChans==2
+            [PunctaWanted,Pvolume,Punkta1Stack]=PreProPunkta_2chan(load_main,Animal,FileCell,range,res,Mainstack,load1);
+        end
+        %quantify intersections
+        [results,Int]=CalcInt(Mainstack,MainVolume,res,PunctaWanted,Pvolume,spMarker);
+
+        if ~exist(savedir,'dir')
+            mkdir(savedir);
+        end
+        save(strcat(savedir,'results.mat'),'results');
+        clear results;
+        %% write the images -> for image overlay generation as in the paper
+        Mainstack = double(Mainstack);
+        path1=(strcat(savedir,'Chan1_processed\'));
+        mkdir(path1)
+        path2=(strcat(savedir,'Chan2_processed\'));
+        mkdir(path2)
+        if NumChans==3
+            path3=(strcat(savedir,'Chan3_processed\'));
+            mkdir(path3)
+        end
+        for r=1:size(Mainstack,3)
+            cd (path1)
+            i = Mainstack(:,:,r);
+            saveImages = ['Chan1-', num2str(r), '.tif'];
+            imwrite(i,saveImages);
+            clear i saveImages
+
+            cd (path2)
+            i = Punkta1Stack(:,:,r);
+            saveImages = ['Chan2-', num2str(r), '.tif'];
+            imwrite(i,saveImages);
+            clear i saveImages
+            
+            if NumChans==3
+                cd (path3)
+                i = Punkta2Stack(:,:,r);
+                saveImages = ['Chan3-', num2str(r), '.tif'];
+                imwrite(i,saveImages);
+                clear i saveImages
+            end
+            
+        end
+    else
+        fprintf(['Cell number ' Cell 'of animal ' Animal ' already underwent pruning analysis\n']);
+    end
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

code/morphology/PreProMainStack.m

@@ -0,0 +1,49 @@
+function [MGwanted,MGvolume,res]=PreProMainStack(load_main,Animal,Cell,m,range,int,loadMainStack)
+
+info  =  imfinfo(strcat(load_main,Animal,'\',Cell,'\',loadMainStack,'.tif'));   
+
+MGstack = zeros(info(1).Width,info(1).Height, range(2) - range (1)+1);   
+
+count=0;
+for r=range(1):range(2)
+    count=count+1;
+    mg  =  imread(strcat(load_main,Animal,'\',Cell,'\',loadMainStack,'.tif'), r, 'Info', info);
+    background  =  imopen(mg, strel('disk',50));
+    mg  =  mg - background;
+    MGstack(:,:,count)  =  mg;
+    clear mg
+end
+
+[~, BWmg]  =  hysteresis3d(MGstack, 0.1, 0.5, 26); % thresholds the microglia image with lower threshold 0.1 and upper threshold 0.5 and segments it with connectivity of 26, only returns the hysteresis images (~negates first output)
+clear MGstack
+[MGstack]  =  ordfilt3D(BWmg,14); % filters the BW microglia stack for 26 neighbours
+clear BWmg
+
+[MGstack_label, MGnum]  =  bwlabeln (MGstack);% labels parts of one distinct microglia with one number
+size_mg = sum(MGstack_label(:)== 1:MGnum); % sums up pixels with the same number (so pixels of each microglial part or microglia)
+
+[MGvolume, idx] = max(size_mg);% returns the volume and the index of the highest pixel volume
+if m==2
+    size_mg(1,idx)=0;
+    clear idx MGvolume
+    [MGvolume, idx] = max(size_mg);
+elseif m==3
+    size_mg(1,idx)=0;
+    clear idx MGvolume
+    [MGvolume, idx] = max(size_mg);
+    size_mg(1,idx)=0;
+    clear idx MGvolume
+    [MGvolume, idx] = max(size_mg);
+end
+
+MGwanted = (MGstack_label == idx);% makes stack from the original microglia stack but with only the desired microglia displayed as logical one (so only it is shown)
+
+res.pixel_width_um = info(1).Width./info(1).XResolution;% enter your pixel width in um (get from info in .tif)
+res.pixel_height_um = info(1).Height./info(1).YResolution; % enter your pixel heigth in um
+res.pixel_depth_um = int; % enter the step size (slice thickness) in um
+
+MGvolume = MGvolume* (res.pixel_width_um*res.pixel_height_um*res.pixel_depth_um);% calculates the volume in um of the microglia from the volume in pixels
+                     
+                     
+                     
+                     

code/morphology/PreProPunkta_1chan.m

@@ -0,0 +1,25 @@
+function [PunctaWanted,Pvolume,PunctaS1]=PreProPunkta_1chan(load_main,Animal,Cell,range,res,MGstack,load1)
+
+CDinfo  =  imfinfo(strcat(load_main,Animal,'\',Cell,'\',load1,'.tif'));   
+
+CDstack = zeros(CDinfo(1).Width,CDinfo(1).Height, range(2) - range (1)+1);   
+
+count=0;
+for r=range(1):range(2)
+    count=count+1;
+    CD  =  imread(strcat(load_main,Animal,'\',Cell,'\',load1,'.tif'), r, 'Info', CDinfo); % returns the value of the k^th image
+    CD  =  logical(CD); % converts it to logical zero (if info is zero) or one (if info is nonzero)
+    CDstack(:, :, count)  =  CD; % fills the matrix with the logical zero or one
+    clear CD
+end
+
+PunctaWanted =CDstack;                        
+
+Pvolume = sum(sum(sum(PunctaWanted)))* (res.pixel_width_um*res.pixel_height_um*res.pixel_depth_um);
+
+PunctaS1=  PSDstack + MGstack;                        
+PunctaS1(PunctaS1 == 1) = 0;                                      
+PunctaS1(PunctaS1 == 2) = 1;     
+
+
+

code/morphology/PreProPunkta_2chan.m

@@ -0,0 +1,45 @@
+function [PunctaWanted,Pvolume,PunctaS1,PunctaS2]=PreProPunkta_2chan(load_main,Animal,Cell,range,res,MGstack,load1,load2)
+
+CDinfo  =  imfinfo(strcat(load_main,Animal,'\',Cell,'\',load1,'.tif'));   
+PSDinfo  =  imfinfo(strcat(load_main,Animal,'\',Cell,'\',load2,'.tif'));   
+
+CDstack = zeros(CDinfo(1).Width,CDinfo(1).Height, range(2) - range (1)+1);   
+PSDstack = zeros(PSDinfo(1).Width,PSDinfo(1).Height, range(2) - range (1)+1);   
+
+count=0;
+for r=range(1):range(2)
+    count=count+1;
+    CD  =  imread(strcat(load_main,Animal,'\',Cell,'\',load1,'.tif'), r, 'Info', CDinfo); % returns the value of the k^th image
+    CD  =  logical(CD); % converts it to logical zero (if info is zero) or one (if info is nonzero)
+    CDstack(:, :, count)  =  CD; % fills the matrix with the logical zero or one
+    clear CD
+    PSD  =  imread(strcat(load_main,Animal,'\',Cell,'\',load2,'.tif'), r, 'Info', PSDinfo); % returns the value of the k^th image
+    PSD  =  logical(PSD); % converts it to logical zero (if info is zero) or one (if info is nonzero)
+    PSDstack(:, :, count)  =  PSD; % fills the matrix with the logical zero or one
+    clear PSD
+end
+
+PunctaWanted =  PSDstack + CDstack;                        
+PunctaWanted(PunctaWanted == 1) = 0;                                      
+PunctaWanted(PunctaWanted == 2) = 1;     
+
+Pvolume = sum(sum(sum(PunctaWanted)))* (res.pixel_width_um*res.pixel_height_um*res.pixel_depth_um);
+
+PunctaS1=  PSDstack + MGstack;                        
+PunctaS1(PunctaS1 == 1) = 0;                                      
+PunctaS1(PunctaS1 == 2) = 1;     
+
+
+PunctaS2=  CDstack + MGstack;                        
+PunctaS2(PunctaS2 == 1) = 0;                                      
+PunctaS2(PunctaS2 == 2) = 1;   
+
+
+
+
+
+
+
+
+
+

code/morphology/get_experiment_list.m

@@ -0,0 +1,42 @@
+function experiments=get_experiment_list(Path,Nexperiment)
+% function creates the selected experiment plus correlated parameters from
+% the specified excel file
+% input variable optional (vector), if not defined: all experiments are selected,
+% if defined: selected experiments will be used
+
+% Path='Q:\Personal\Sebastian\ProjectOptoChronicPFC\ExperimentPlan.xlsx';
+ExcelSheet='InfoandDevMil';
+xlRange='A1:FZ500';
+[~,~,InfoandDevMil] = xlsread(Path,ExcelSheet,xlRange); % Import recording summary from excel sheet
+for f1=1:size(InfoandDevMil,1)
+    for f2=1:size(InfoandDevMil,2)
+        if strcmp(InfoandDevMil{f1,f2},'NaN')
+            InfoandDevMil{f1,f2}=NaN;
+        end
+    end
+end
+
+idxR_header1=3;
+idxR_header2=4;
+idxR_header3=5;
+
+
+%% assign InfoandDevMil to experiments structure
+idxC_header1=[find(cellfun(@ischar,InfoandDevMil(idxR_header1,:))) size(InfoandDevMil,2)+1];
+IdxR=0;
+idxC_Nexperiment=find(strcmp(InfoandDevMil(idxR_header3,:),'Nexperiment'));
+idxC_ExpType=find(strcmp(InfoandDevMil(idxR_header3,:),'ExpType'));
+for r1=6:size(InfoandDevMil,1)
+    if ~isnan(InfoandDevMil{r1,idxC_Nexperiment}) && ischar(InfoandDevMil{r1,idxC_ExpType}) && or(isempty(Nexperiment),ismember(InfoandDevMil{r1,idxC_Nexperiment},Nexperiment))
+        IdxR=IdxR+1;
+        for h1=1:length(idxC_header1)-1
+            idxC_header2=[find(cellfun(@ischar,InfoandDevMil(idxR_header2,idxC_header1(h1):idxC_header1(h1+1)-1)))+idxC_header1(h1)-1 idxC_header1(h1+1)];
+            for h2=1:length(idxC_header2)-1
+                idxC_header3=find(cellfun(@ischar,InfoandDevMil(idxR_header3,idxC_header2(h2):idxC_header2(h2+1)-1)))+idxC_header2(h2)-1;
+                for h3=1:length(idxC_header3)
+                        experiments(IdxR,1).(InfoandDevMil{idxR_header1,idxC_header1(h1)}).(InfoandDevMil{idxR_header2,idxC_header2(h2)}).(InfoandDevMil{idxR_header3,idxC_header3(h3)})=InfoandDevMil{r1,idxC_header3(h3)};
+                end
+            end
+        end
+    end
+end

code/morphology/hysteresis3d.m

@@ -0,0 +1,74 @@
+function [tri,hys]=hysteresis3d(img,t1,t2,conn)
+% function [tri,hys]=HYSTERESIS3D(img,t1,t2,conn)
+%
+% Hysteresis3d is a simple function that performs trinarisation and
+% hysteresis for 2D and 3D images. Hysteresis3d was inspired by Peter
+% Kovesi's 2D hysteresis function
+% (http://www.csse.uwa.edu.au/~pk/research/matlabfns/). This 3D function
+% takes advantage of the 3D connectivities of imfill instead of the 2D
+% connectivities of bwselect.
+%
+% Usage:        [tri,hys]=HYSTERESIS3D(img,t1,t2,conn)
+%
+% Arguments:    img - image for hysteresis (assumed to be non-negative)
+%               t1 - lower threshold value (fraction b/w 0-1, e.g.: 0.1)
+%               t2 - upper threshold value (fraction b/w 0-1, e.g.: 0.9)
+%                   (t1/t2 can be entered in any order, larger one will be 
+%                   set as the upper threshold)
+%               conn - number of connectivities (4 or 8 for 2D)
+%                                               (6, 18, or 26 for 3D)       
+% Returns:
+%               tri - the trinarisation image (values are 0, 1, or 2)
+%               hys - the hysteresis image (logical mask image)
+% 
+% Examples:     [tri,hys]=HYSTERESIS3D(img,0.25,0.8,26)
+%
+% 2012/07/10: written by Luke Xie 
+% 2013/12/09: defaults added 
+%
+% To see an example of hysteresis used to segment a kidney region, please 
+% refer to supplement in QSM of Kidney, NMR Biomed, 2013 Dec;26(12):1853-63 
+% (http://onlinelibrary.wiley.com/doi/10.1002/nbm.3039/abstract).
+% Supplemental material is also available on our CIVMspace: 
+% http://www.civm.duhs.duke.edu/lx201204/
+
+%% arguments
+if nargin<3
+    disp('function needs at least 3 inputs')
+    return;
+elseif nargin==3
+    disp('inputs=3')
+    if numel(size(img))==2;
+        disp('img=2D')
+        disp('conn set at 4 connectivies (number of neighbors)')
+        conn=4;
+    end
+    if numel(size(img))==3; 
+        disp('img=3D')
+        disp('conn set at 6 connectivies (number of neighbors)')
+        conn=6;
+    end
+end
+
+%% scale t1 & t2 based on image intensity range
+if t1>t2    % swap values if t1>t2 
+	tmp=t1;
+	t1=t2; 
+	t2=tmp;
+end
+minv=min(img(:));                % min image intensity value
+maxv=max(img(:));                % max image intensity value
+t1v=t1*(maxv-minv)+minv;
+t2v=t2*(maxv-minv)+minv;
+
+%% trinarisation
+tri=zeros(size(img));
+tri(img>=t1v)=1;
+tri(img>=t2v)=2;
+
+%% hysteresis
+abovet1=img>t1v;                                     % points above lower threshold
+seed_indices=sub2ind(size(abovet1),find(img>t2v));   % indices of points above upper threshold
+hys=imfill(~abovet1,seed_indices,conn);              % obtain all connected regions in abovet1 that include points with values above t2
+hys=hys & abovet1;
+

code/morphology/imshow3D.m

@@ -0,0 +1,291 @@
+function  imshow3D( Img, disprange )
+%IMSHOW3D displays 3D grayscale images in slice by slice fashion with mouse
+%based slice browsing and window and level adjustment control.
+%
+% Usage:
+% imshow3D ( Image )
+% imshow3D ( Image , [] )
+% imshow3D ( Image , [LOW HIGH] )
+%   
+%    Image:      3D image MxNxK (K slices of MxN images) 
+%    [LOW HIGH]: display range that controls the display intensity range of
+%                a grayscale image (default: the widest available range)
+%
+% Use the scroll bar or mouse scroll wheel to switch between slices. To
+% adjust window and level values keep the mouse right button pressed and
+% drag the mouse up and down (for level adjustment) or right and left (for
+% window adjustment). 
+% 
+% "Auto W/L" button adjust the window and level automatically 
+%
+% While "Fine Tune" check box is checked the window/level adjustment gets
+% 16 times less sensitive to mouse movement, to make it easier to control
+% display intensity rang.
+%
+% Note: The sensitivity of mouse based window and level adjustment is set
+% based on the user defined display intensity range; the wider the range
+% the more sensitivity to mouse drag.
+% 
+% 
+%   Example
+%   --------
+%       % Display an image (MRI example)
+%       load mri 
+%       Image = squeeze(D); 
+%       figure, 
+%       imshow3D(Image) 
+%
+%       % Display the image, adjust the display range
+%       figure,
+%       imshow3D(Image,[20 100]);
+%
+%   See also IMSHOW.
+
+%
+% - Maysam Shahedi (mshahedi@gmail.com)
+% - Released: 1.0.0   Date: 2013/04/15
+% - Revision: 1.1.0   Date: 2013/04/19
+% 
+
+sno = size(Img,3);  % number of slices
+S = round(sno/2);
+
+global InitialCoord;
+
+MinV = 0;
+MaxV = max(Img(:));
+LevV = (double( MaxV) + double(MinV)) / 2;
+Win = double(MaxV) - double(MinV);
+WLAdjCoe = (Win + 1)/1024;
+FineTuneC = [1 1/16];    % Regular/Fine-tune mode coefficients
+
+if isa(Img,'uint8')
+    MaxV = uint8(Inf);
+    MinV = uint8(-Inf);
+    LevV = (double( MaxV) + double(MinV)) / 2;
+    Win = double(MaxV) - double(MinV);
+    WLAdjCoe = (Win + 1)/1024;
+elseif isa(Img,'uint16')
+    MaxV = uint16(Inf);
+    MinV = uint16(-Inf);
+    LevV = (double( MaxV) + double(MinV)) / 2;
+    Win = double(MaxV) - double(MinV);
+    WLAdjCoe = (Win + 1)/1024;
+elseif isa(Img,'uint32')
+    MaxV = uint32(Inf);
+    MinV = uint32(-Inf);
+    LevV = (double( MaxV) + double(MinV)) / 2;
+    Win = double(MaxV) - double(MinV);
+    WLAdjCoe = (Win + 1)/1024;
+elseif isa(Img,'uint64')
+    MaxV = uint64(Inf);
+    MinV = uint64(-Inf);
+    LevV = (double( MaxV) + double(MinV)) / 2;
+    Win = double(MaxV) - double(MinV);
+    WLAdjCoe = (Win + 1)/1024;
+elseif isa(Img,'int8')
+    MaxV = int8(Inf);
+    MinV = int8(-Inf);
+    LevV = (double( MaxV) + double(MinV)) / 2;
+    Win = double(MaxV) - double(MinV);
+    WLAdjCoe = (Win + 1)/1024;
+elseif isa(Img,'int16')
+    MaxV = int16(Inf);
+    MinV = int16(-Inf);
+    LevV = (double( MaxV) + double(MinV)) / 2;
+    Win = double(MaxV) - double(MinV);
+    WLAdjCoe = (Win + 1)/1024;
+elseif isa(Img,'int32')
+    MaxV = int32(Inf);
+    MinV = int32(-Inf);
+    LevV = (double( MaxV) + double(MinV)) / 2;
+    Win = double(MaxV) - double(MinV);
+    WLAdjCoe = (Win + 1)/1024;
+elseif isa(Img,'int64')
+    MaxV = int64(Inf);
+    MinV = int64(-Inf);
+    LevV = (double( MaxV) + double(MinV)) / 2;
+    Win = double(MaxV) - double(MinV);
+    WLAdjCoe = (Win + 1)/1024;
+elseif isa(Img,'logical')
+    MaxV = 0;
+    MinV = 1;
+    LevV =0.5;
+    Win = 1;
+    WLAdjCoe = 0.1;
+end    
+
+SFntSz = 9;
+LFntSz = 10;
+WFntSz = 10;
+LVFntSz = 9;
+WVFntSz = 9;
+BtnSz = 10;
+ChBxSz = 10;
+
+if (nargin < 2)
+    [Rmin Rmax] = WL2R(Win, LevV);
+elseif numel(disprange) == 0
+    [Rmin Rmax] = WL2R(Win, LevV);
+else
+    LevV = (double(disprange(2)) + double(disprange(1))) / 2;
+    Win = double(disprange(2)) - double(disprange(1));
+    WLAdjCoe = (Win + 1)/1024;
+    [Rmin Rmax] = WL2R(Win, LevV);
+end
+
+axes('position',[0,0.2,1,0.8]), imshow(Img(:,:,S), [Rmin Rmax])
+
+FigPos = get(gcf,'Position');
+S_Pos = [50 45 uint16(FigPos(3)-100)+1 20];
+Stxt_Pos = [50 65 uint16(FigPos(3)-100)+1 15];
+Wtxt_Pos = [50 20 60 20];
+Wval_Pos = [110 20 60 20];
+Ltxt_Pos = [175 20 45 20];
+Lval_Pos = [220 20 60 20];
+BtnStPnt = uint16(FigPos(3)-250)+1;
+if BtnStPnt < 300
+    BtnStPnt = 300;
+end
+Btn_Pos = [BtnStPnt 20 100 20];
+ChBx_Pos = [BtnStPnt+110 20 100 20];
+
+if sno > 1
+    shand = uicontrol('Style', 'slider','Min',1,'Max',sno,'Value',S,'SliderStep',[1/(sno-1) 10/(sno-1)],'Position', S_Pos,'Callback', {@SliceSlider, Img});
+    stxthand = uicontrol('Style', 'text','Position', Stxt_Pos,'String',sprintf('Slice# %d / %d',S, sno), 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', SFntSz);
+else
+    stxthand = uicontrol('Style', 'text','Position', Stxt_Pos,'String','2D image', 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', SFntSz);
+end    
+ltxthand = uicontrol('Style', 'text','Position', Ltxt_Pos,'String','Level: ', 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', LFntSz);
+wtxthand = uicontrol('Style', 'text','Position', Wtxt_Pos,'String','Window: ', 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', WFntSz);
+lvalhand = uicontrol('Style', 'edit','Position', Lval_Pos,'String',sprintf('%6.0f',LevV), 'BackgroundColor', [1 1 1], 'FontSize', LVFntSz,'Callback', @WinLevChanged);
+wvalhand = uicontrol('Style', 'edit','Position', Wval_Pos,'String',sprintf('%6.0f',Win), 'BackgroundColor', [1 1 1], 'FontSize', WVFntSz,'Callback', @WinLevChanged);
+Btnhand = uicontrol('Style', 'pushbutton','Position', Btn_Pos,'String','Auto W/L', 'FontSize', BtnSz, 'Callback' , @AutoAdjust);
+ChBxhand = uicontrol('Style', 'checkbox','Position', ChBx_Pos,'String','Fine Tune', 'BackgroundColor', [0.8 0.8 0.8], 'FontSize', ChBxSz);
+
+set (gcf, 'WindowScrollWheelFcn', @mouseScroll);
+set (gcf, 'ButtonDownFcn', @mouseClick);
+set(get(gca,'Children'),'ButtonDownFcn', @mouseClick);
+set(gcf,'WindowButtonUpFcn', @mouseRelease)
+set(gcf,'ResizeFcn', @figureResized)
+
+
+% -=< Figure resize callback function >=-
+    function figureResized(object, eventdata)
+        FigPos = get(gcf,'Position');
+        S_Pos = [50 45 uint16(FigPos(3)-100)+1 20];
+        Stxt_Pos = [50 65 uint16(FigPos(3)-100)+1 15];
+        BtnStPnt = uint16(FigPos(3)-250)+1;
+        if BtnStPnt < 300
+            BtnStPnt = 300;
+        end
+        Btn_Pos = [BtnStPnt 20 100 20];
+        ChBx_Pos = [BtnStPnt+110 20 100 20];
+        if sno > 1
+            set(shand,'Position', S_Pos);
+        end
+        set(stxthand,'Position', Stxt_Pos);
+        set(ltxthand,'Position', Ltxt_Pos);
+        set(wtxthand,'Position', Wtxt_Pos);
+        set(lvalhand,'Position', Lval_Pos);
+        set(wvalhand,'Position', Wval_Pos);
+        set(Btnhand,'Position', Btn_Pos);
+        set(ChBxhand,'Position', ChBx_Pos);
+    end
+
+% -=< Slice slider callback function >=-
+    function SliceSlider (hObj,event, Img)
+        S = round(get(hObj,'Value'));
+        set(get(gca,'children'),'cdata',Img(:,:,S))
+        caxis([Rmin Rmax])
+        if sno > 1
+            set(stxthand, 'String', sprintf('Slice# %d / %d',S, sno));
+        else
+            set(stxthand, 'String', '2D image');
+        end
+    end
+
+% -=< Mouse scroll wheel callback function >=-
+    function mouseScroll (object, eventdata)
+        UPDN = eventdata.VerticalScrollCount;
+        S = S - UPDN;
+        if (S < 1)
+            S = 1;
+        elseif (S > sno)
+            S = sno;
+        end
+        if sno > 1
+            set(shand,'Value',S);
+            set(stxthand, 'String', sprintf('Slice# %d / %d',S, sno));
+        else
+            set(stxthand, 'String', '2D image');
+        end
+        set(get(gca,'children'),'cdata',Img(:,:,S))
+    end
+
+% -=< Mouse button released callback function >=-
+    function mouseRelease (object,eventdata)
+        set(gcf, 'WindowButtonMotionFcn', '')
+    end
+
+% -=< Mouse click callback function >=-
+    function mouseClick (object, eventdata)
+        MouseStat = get(gcbf, 'SelectionType');
+        if (MouseStat(1) == 'a')        %   RIGHT CLICK
+            InitialCoord = get(0,'PointerLocation');
+            set(gcf, 'WindowButtonMotionFcn', @WinLevAdj);
+        end
+    end
+
+% -=< Window and level mouse adjustment >=-
+    function WinLevAdj(varargin)
+        PosDiff = get(0,'PointerLocation') - InitialCoord;
+
+        Win = Win + PosDiff(1) * WLAdjCoe * FineTuneC(get(ChBxhand,'Value')+1);
+        LevV = LevV - PosDiff(2) * WLAdjCoe * FineTuneC(get(ChBxhand,'Value')+1);
+        if (Win < 1)
+            Win = 1;
+        end
+
+        [Rmin, Rmax] = WL2R(Win,LevV);
+        caxis([Rmin, Rmax])
+        set(lvalhand, 'String', sprintf('%6.0f',LevV));
+        set(wvalhand, 'String', sprintf('%6.0f',Win));
+        InitialCoord = get(0,'PointerLocation');
+    end
+
+% -=< Window and level text adjustment >=-
+    function WinLevChanged(varargin)
+
+        LevV = str2double(get(lvalhand, 'string'));
+        Win = str2double(get(wvalhand, 'string'));
+        if (Win < 1)
+            Win = 1;
+        end
+
+        [Rmin, Rmax] = WL2R(Win,LevV);
+        caxis([Rmin, Rmax])
+    end
+
+% -=< Window and level to range conversion >=-
+    function [Rmn Rmx] = WL2R(W,L)
+        Rmn = L - (W/2);
+        Rmx = L + (W/2);
+        if (Rmn >= Rmx)
+            Rmx = Rmn + 1;
+        end
+    end
+
+% -=< Window and level auto adjustment callback function >=-
+    function AutoAdjust(object,eventdata)
+        Win = double(max(Img(:))-min(Img(:)));
+        Win (Win < 1) = 1;
+        LevV = double(min(Img(:)) + (Win/2));
+        [Rmin, Rmax] = WL2R(Win,LevV);
+        caxis([Rmin, Rmax])
+        set(lvalhand, 'String', sprintf('%6.0f',LevV));
+        set(wvalhand, 'String', sprintf('%6.0f',Win));
+    end
+
+end
+% -=< Maysam Shahedi (mshahedi@gmail.com), April 19, 2013>=-

code/morphology/ordfilt3D.m

@@ -0,0 +1,100 @@
+function [Vr] = ordfilt3D(V0,ord,padoption);
+% ordfilt3D:    Perform 3-D order-statistic filtering on 26 neighbors
+%
+%   [Vr] = ordfilt3D(V0,ord,padoption)
+%          use 26 neighbors
+%       ord = 14 <=> median filtering
+%       ord = 1 <=> min
+%       ord = [1 27] <=> [min max]
+%       padoption: same as in padarray
+%
+% Olivier Salvado, Case Western Reserve University, 16Aug04
+
+
+
+if ~exist('padoption','var')
+    padoption = 'replicate';
+end
+
+
+%%
+% special care for uint8
+if isa(V0,'uint8')
+    V = uint8(padarray(V0,[1 1 1],padoption));
+    S = size(V);
+    Vn = uint8(zeros(S(1),S(2),S(3),26));  % all the neighbor
+else
+    V = single(padarray(V0,[1 1 1],padoption));
+    S = size(V);
+    Vn = single(zeros(S(1),S(2),S(3),26));  % all the neighbor
+end
+
+%%
+% build the neighboord
+Vn(:,:,:,1) = V;
+i = 1:S(1); ip1 = [i(2:end) i(end)]; im1 = [i(1) i(1:end-1)];
+j = 1:S(2); jp1 = [j(2:end) j(end)]; jm1 = [j(1) j(1:end-1)];
+k = 1:S(3); kp1 = [k(2:end) k(end)]; km1 = [k(1) k(1:end-1)];
+
+%%
+% left
+Vn(:,:,:,2)     = V(im1    ,jm1    ,km1);
+Vn(:,:,:,3)     = V(im1    ,j      ,km1);
+Vn(:,:,:,4)     = V(im1    ,jp1    ,km1);
+
+Vn(:,:,:,5)     = V(im1    ,jm1    ,k);
+Vn(:,:,:,6)     = V(im1    ,j      ,k);
+Vn(:,:,:,7)     = V(im1    ,jp1    ,k);
+
+Vn(:,:,:,8)     = V(im1    ,jm1    ,kp1);
+Vn(:,:,:,9)     = V(im1    ,j      ,kp1);
+Vn(:,:,:,10)    = V(im1    ,jp1    ,kp1);
+
+%%
+% right
+Vn(:,:,:,11)    = V(ip1    ,jm1    ,km1);
+Vn(:,:,:,12)    = V(ip1    ,j      ,km1);
+Vn(:,:,:,13)    = V(ip1    ,jp1    ,km1);
+
+Vn(:,:,:,14)    = V(ip1    ,jm1    ,k);
+Vn(:,:,:,15)    = V(ip1    ,j      ,k);
+Vn(:,:,:,16)    = V(ip1    ,jp1    ,k);
+
+Vn(:,:,:,17)    = V(ip1    ,jm1    ,kp1);
+Vn(:,:,:,18)    = V(ip1    ,j      ,kp1);
+Vn(:,:,:,19)    = V(ip1    ,jp1    ,kp1);
+
+%%
+% top
+Vn(:,:,:,20)    = V(i       ,jm1    ,kp1);
+Vn(:,:,:,21)    = V(i       ,j      ,kp1);
+Vn(:,:,:,22)    = V(i       ,jp1    ,kp1);
+
+%%
+% bottom
+Vn(:,:,:,23)    = V(i       ,jm1    ,km1);
+Vn(:,:,:,24)    = V(i       ,j      ,km1);
+Vn(:,:,:,25)    = V(i       ,jp1    ,km1);
+
+%%
+% front
+Vn(:,:,:,26)    = V(i       ,jp1    ,k);
+
+%%
+% back
+Vn(:,:,:,27)    = V(i       ,jm1    ,k);
+
+%%
+% perform the processing
+Vn = sort(Vn,4);
+Vr = Vn(:,:,:,ord);
+
+
+%%
+% remove padding on the 3 first dimensions
+Vr = Vr(2:end-1,2:end-1,2:end-1,:);
+
+
+
+
+