Author: Kevin-Mattheus-Moerman

HELP_patch2Im.m

@@ -171,7 +171,7 @@
 patch('Faces',F,'Vertices',V,'FaceColor','flat','CData',C,'EdgeColor','none','FaceAlpha',faceAlpha2);
 camlight('headlight'); lighting flat;
 axis equal; view(3); axis tight;  grid on;  set(gca,'FontSize',fontSize);
-colormap gjet; colorbar; 
+colormap gjet; colorbar; caxis([0 max(C(:))]);
 
 subplot(1,2,2);
 title('Patch data derived image data (3 slices)','FontSize',fontSize);
@@ -189,7 +189,7 @@
 Vm=Vm+imOrigin(ones(size(Vm,1),1),:);
 patch('Faces',Fm,'Vertices',Vm,'FaceColor','flat','CData',Cm,'EdgeColor','k','FaceAlpha',faceAlpha1);
 
-colormap gjet; colorbar;
+colormap gjet; colorbar; caxis([0 max(C(:))]);
 camlight('headlight'); lighting flat;
 axis equal; view(3); axis tight;  grid on;  set(gca,'FontSize',fontSize);
 drawnow;

html/helpsearch-v2/_0.cfe

@@ -165,8 +165,8 @@
 hf.UserData.ButtonHandles.Play=hPlay;
 hf.UserData.ButtonHandles.hCycle=hCycle;
 hf.UserData.ButtonHandles.hTextTime=hTextTime;
-hf.UserData.pauseTime=2/numel(animStruct.Time);
-hf.UserData.shiftMag=1; 
+hf.UserData.pauseTime=mean(diff(animStruct.Time));
+hf.UserData.shiftMag=ceil(numel(animStruct.Time)/20); 
 hf.UserData.icons.play=iconPlay;
 hf.UserData.icons.stop=iconStop;
 

html/helpsearch-v2/_0.cfs

@@ -83,9 +83,16 @@
 for q=1:1:numContours
 
     Vg=[Q(1,IND1(q):IND2(q))' Q(2,IND1(q):IND2(q))']; %Get coordinates
-    cLevel(q)=Q(1,IND1(q)-1); %Get level
+    cLevel(q)=Q(1,IND1(q)-1); %Get level    
+    
+    %Check for non-unique points
+    [~,indUni,~]=unique(pround(Vg,5),'rows');
+    logicKeep=false(size(Vg,1),1);
+    logicKeep(indUni)=1;
+    Vg=Vg(logicKeep,:);
 
     if size(Vg,1)>1
+        
         %Check if contour is closed loop or not
         d=sqrt(sum((Vg(1,:)-Vg(end,:)).^2)); %Distance between first and last point
         if d<(min(v)/10) %Smaller than 1/10th of a pixel dimension

html/helpsearch-v2/_0.si

@@ -1,31 +1,70 @@
-function [V_inner]=getInnerPoint(F_input,V_input,searchRadius,voxelSize,plotOn)
+function [varargout]=getInnerPoint(varargin)
 
 %%
 
-if isa(F_input,'cell')
-    [F_input,V_input,C_input]=joinElementSets(F_input,V_input);
+switch nargin
+    case 2
+        F=varargin{1};
+        V=varargin{2};
+        searchRadius=[];
+        voxelSize=[];
+        plotOn=0;
+    case 3        
+        F=varargin{1};
+        V=varargin{2};
+        searchRadius=varargin{3};
+        voxelSize=[];
+        plotOn=0;
+    case 4            
+        F=varargin{1};
+        V=varargin{2};
+        searchRadius=varargin{3};
+        voxelSize=varargin{4};
+        plotOn=0;
+    case 5
+        F=varargin{1};
+        V=varargin{2};
+        searchRadius=varargin{3};
+        voxelSize=varargin{4};
+        plotOn=varargin{5};
+    otherwise
+        error('Wrong number of input arguments');
+end
+
+if isempty(voxelSize)
+   D=patchEdgeLengths(F,V); 
+   voxelSize=mean(D)/2;
+end
+
+if isempty(searchRadius)
+   searchRadius=voxelSize*3;
+end
+
+%%
+
+if isa(F,'cell')
+    [F,V,C]=joinElementSets(F,V);
 else
-    C_input=ones(size(F_input,1),1);
+    C=ones(size(F,1),1);
 end
 
 %% Get an inner point
 
 %Convert patch data to image
-[M,G,~]=patch2Im(F_input,V_input,C_input,voxelSize);
+[M,G,~]=patch2Im(F,V,C,voxelSize);
 L_test=(M==1); %Logic for interior voxels of surface mesh
-M=double(L_test); 
 imOrigin=G.origin; %Image origin can be used to allign image with surface
 
-
 %Construct search blurring kernel
-[x,y,z]=ndgrid(-searchRadius:1:searchRadius); %kernel coordinates
+[x,y,z]=ndgrid(-searchRadius:voxelSize:searchRadius); %kernel coordinates
 d=sqrt(x.^2+y.^2+z.^2); %distance metric
 k=(d<=searchRadius); %Define kernel based on radius
 k=k./sum(k(:)); %Normalize kernel
 
 %Convolve interior image with kernel
 ML=convn(double(L_test),k,'same');
-[~,indInternal]=max(ML(:)); %Kernel should yield max at "deep" (related to search radius) point
+ML(M~=1)=NaN; %Set other sites to NaN
+[~,indInternal]=nanmax(ML(:)); %Kernel should yield max at "deep" (related to search radius) point
 [I_in,J_in,K_in]=ind2sub(size(M),indInternal); %Convert to subscript coordinates
 
 %Derive point coordinates
@@ -35,47 +74,51 @@
 
 %% Plotting image, mesh and inner point
 if plotOn==1    
-    figColor='w'; figColorDef='white';
+    
     fontSize=20;
     markerSize1=50;%round(voxelSize*25);
-    faceAlpha1=0.7;
-    faceAlpha2=0.2;
+    faceAlpha1=1;
+    faceAlpha2=0.5;
 
-    hf=figuremax(figColor,figColorDef);
-    xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
+    cFigure;    
     hold on;
 
-    hp= patch('Faces',F_input,'Vertices',V_input,'FaceColor','flat','CData',C_input,'EdgeColor','none','FaceAlpha',faceAlpha2);
+    gpatch(F,V,0.5*ones(1,3),'none',faceAlpha2);
     plotV(V_inner,'k.','MarkerSize',markerSize1);
 
-    L_plot=false(size(M));
+    L_plot=false(size(ML));
     L_plot(:,:,K_in)=1;
-    L_plot=L_plot&M>0;
-    [Fm,Vm,Cm]=ind2patch(L_plot,double(M),'sk');
+    L_plot=L_plot & ~isnan(ML);
+    [Fm,Vm,Cm]=ind2patch(L_plot,double(ML),'sk');
     [Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3));
     Vm=Vm+imOrigin(ones(size(Vm,1),1),:);
-    patch('Faces',Fm,'Vertices',Vm,'FaceColor',0.5*ones(1,3),'EdgeColor','k','FaceAlpha',faceAlpha1);
 
-    L_plot=false(size(M));
+    gpatch(Fm,Vm,Cm,'k',faceAlpha1);    
+    
+    L_plot=false(size(ML));
     L_plot(I_in,:,:)=1;
-    L_plot=L_plot&M>0;
-    [Fm,Vm,Cm]=ind2patch(L_plot,M,'si');
+    L_plot=L_plot & ~isnan(ML);
+    [Fm,Vm,Cm]=ind2patch(L_plot,ML,'si');
     [Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3));
     Vm=Vm+imOrigin(ones(size(Vm,1),1),:);
-    patch('Faces',Fm,'Vertices',Vm,'FaceColor',0.5*ones(1,3),'EdgeColor','k','FaceAlpha',faceAlpha1);
+    gpatch(Fm,Vm,Cm,'k',faceAlpha1);    
 
-    L_plot=false(size(M));
+    L_plot=false(size(ML));
     L_plot(:,J_in,:)=1;
-    L_plot=L_plot&M>0;
-    [Fm,Vm,Cm]=ind2patch(L_plot,M,'sj');
+    L_plot=L_plot & ~isnan(ML);
+    [Fm,Vm,Cm]=ind2patch(L_plot,ML,'sj');
     [Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3));
     Vm=Vm+imOrigin(ones(size(Vm,1),1),:);
-    patch('Faces',Fm,'Vertices',Vm,'FaceColor',0.5*ones(1,3),'EdgeColor','k','FaceAlpha',faceAlpha1);
+    gpatch(Fm,Vm,Cm,'k',faceAlpha1);    
 
-    colormap(gjet(numel(unique(C_input(:))))); colorbar; 
-    
-    axis equal; view(3); axis tight;  grid on;  set(gca,'FontSize',fontSize);
-    set(gca,'FontSize',fontSize);
+    colormap(gjet(250)); colorbar;     
+    axisGeom(gca,fontSize);
     drawnow;    
 end
 
+varargout{1}=V_inner;
+varargout{2}=M;
+varargout{3}=G;
+varargout{4}=ML;
+
+

html/helpsearch-v2/segments.gen

@@ -19,6 +19,7 @@
         voxelSize=[1 1 1];
         [I_in,J_in,K_in]=ind2sub(size(L),indInternal); %Convert to subscript coordinates
         [X_in,Y_in,Z_in]=im2cart(I_in,J_in,K_in,voxelSize);
+        V_in=[X_in Y_in Z_in];
 
         %Plot settings
         figColor='w'; figColorDef='white';
@@ -31,7 +32,7 @@
         hold on;
 
         %Found voxel location
-        plot3(X_in,Y_in,Z_in,'r.','MarkerSize',markerSize1);
+        plotV(V_in,'r.','MarkerSize',markerSize1);
 
         %Local slices
         L_plot=false(size(L));

html/helpsearch-v2/segments_1

@@ -514,7 +514,7 @@
 hf.UserData.sliceIndices=[sliceIndexI sliceIndexJ sliceIndexK];
 hf.UserData.contourSlice=sliceIndexK;
 hf.UserData.fontColor=fontColor;
-hf.UserData.faceAlpha=0.8;
+hf.UserData.faceAlpha=1;
 hf.UserData.lineColors=gjet(4);
 hf.UserData.logicThreshold=true(size(M));
 hf.UserData.ButtonHandles.Sample=hSample;

html/helpsearch-v3/_5.cfe

@@ -7,9 +7,10 @@
 indSeed=nan(n,1);
 indSeed(1:numStarts)=indStart; 
 
-for q=2:1:n
-    [marchCount,~]=patchMarchCount(F,V,indSeed(~isnan(indSeed)),optStruct);
+for q=numStarts:1:n
+    [marchCount,seedIndex]=patchMarchCount(F,V,indSeed(~isnan(indSeed)),optStruct);
     [~,indAdd]=max(marchCount);
-    indSeed(q)=indAdd;
+    if q<n
+        indSeed(q)=indAdd;
+    end
 end
-[marchCount,seedIndex]=patchMarchCount(F,V,indSeed(~isnan(indSeed)),optStruct);

html/helpsearch-v3/_5.cfs

@@ -0,0 +1,146 @@
+function [varargout]=patchMarchDistMapIterative(varargin)
+
+
+%%
+
+switch nargin
+    case 2
+        F=varargin{1};
+        V=varargin{2};
+        indStart=1;
+        W=ones(size(V,1),1);
+        optStruct.waitBarOn=0;
+        optStruct.numIterationsMax=200;
+        optStruct.toleranceLevel=1e-3;
+    case 3
+        F=varargin{1};
+        V=varargin{2};
+        indStart=varargin{3};
+        W=ones(size(V,1),1);
+        optStruct.waitBarOn=0;
+        optStruct.numIterationsMax=200;
+        optStruct.toleranceLevel=1e-3;
+    case 4
+        F=varargin{1};
+        V=varargin{2};
+        indStart=varargin{3};
+        W=varargin{4};
+        optStruct.waitBarOn=0;
+        optStruct.numIterationsMax=200;
+        optStruct.toleranceLevel=1e-3;
+    case 5
+        F=varargin{1};
+        V=varargin{2};
+        indStart=varargin{3};
+        W=varargin{4};
+        optStruct=varargin{5};
+end
+
+if isempty(W)
+    W=ones(size(V,1),1);
+end
+
+if isempty(indStart)
+    indStart=1;
+end
+
+if isfield(optStruct,'numIterationsMax')
+    numIterationsMax=optStruct.numIterationsMax;
+else
+    numIterationsMax=200;
+end
+
+if isfield(optStruct,'toleranceLevel')
+    toleranceLevel=optStruct.toleranceLevel;
+else
+    toleranceLevel=1e-3;
+end
+
+%%
+
+i = F(:);
+j = [F(:,2); F(:,3); F(:,1)];
+k = [F(:,3); F(:,1); F(:,2)];
+numPoints = size(V,1);
+x  = V(i,:);
+x1 = V(j,:) - x;
+x2 = V(k,:) - x;
+
+Inv1 = @(M,d)[M(:,4)./d -M(:,3)./d -M(:,2)./d M(:,1)./d];
+Inv  = @(M)Inv1(M, M(:,1).*M(:,4) - M(:,3).*M(:,2));
+Mult = @(M,u)[M(:,1).*u(:,1) + M(:,3).*u(:,2)  M(:,2).*u(:,1) + M(:,4).*u(:,2)];
+
+U = getoptions(optStruct, 'U', []);
+if isempty(U)
+    U = zeros(numPoints,1);
+end
+
+% inner product matrix
+C = [dot(x1,x1,2) dot(x1,x2,2) dot(x2,x1,2) dot(x2,x2,2)];
+
+S = Inv(C);
+
+% a = <S*1,1>
+a = sum(S,2);
+
+w = W(i);
+
+% edge length
+L1 = sqrt(dot(x1,x1,2));
+L1 = L1(:).*w(:);
+L2 = sqrt(dot(x2,x2,2));
+L2 = L2(:).*w(:);
+
+E=nan(numPoints,1);
+
+ii=1:1:numel(i);
+
+for q=1:numIterationsMax
+    uj = U(j);
+    uk = U(k);
+    u = [uj uk];
+    b = sum([S(:,1)+S(:,3)  S(:,2)+S(:,4)].*u,2);
+    c = sum(Mult(S,u).*u,2)-w.^2;
+    delta = max( b.^2 - a.*c, 0);
+    d = (b + sqrt(delta) )./a; % solution
+
+    alpha = Mult(S,u-repmat(d, 1,2)); % direction of the update
+
+    logic_d =alpha(:,1)>0 | alpha(:,2)>0;
+
+    % update along edges
+    d1 = L1 + uj(:);
+    d2 = L2 + uk(:);
+    d = d(:);
+    d(logic_d) = min(d1(logic_d), d2(logic_d));
+
+    %Create U1 using sparse array to get minima
+    maxd=max(d(:));
+    df=abs(maxd-d);
+    U1 = sparse(i,ii,df,numPoints,max(ii),numel(d));
+    U1 =maxd-full(max(U1,[],2));
+    U1(U1==0) = Inf;
+    U1(indStart) = 0; % boundary condition
+
+    % enforce monotony
+    if min(U1-U)<-1e-5
+      warning('Monotony problem');
+    end
+
+    % Store error
+    if q==1
+        e1=norm(U-U1,'fro');
+    end
+
+    E(q) = norm(U-U1,'fro')/e1;
+    if E(q)<toleranceLevel
+        break;
+    end
+
+    % Update
+    U = U1;
+end
+
+%Output
+varargout{1}=U;
+varargout{2}=E;