Committing vectorization improvements

Deconvolution/calculateSmallOTF.m

@@ -50,9 +50,12 @@
 originPSF = ceil((sizePSF+ones(1,3))/2); % (for even sizes the origin occurs above the center, i.e. the origin of an image with size 4x4 occurs at (3,3) )
 PSF = PSF-min(PSF(:));
 PSF = PSF./max(PSF(:));
-smallPSF = PSF((originPSF(1)-ceil((newImageSize(1)-1)/2)):(originPSF(1)+floor((newImageSize(1)-1)/2)), ...
+smallPSF = PSF( ...
+    (originPSF(1)-ceil((newImageSize(1)-1)/2)):(originPSF(1)+floor((newImageSize(1)-1)/2)), ...
     (originPSF(2)-ceil((newImageSize(2)-1)/2)):(originPSF(2)+floor((newImageSize(2)-1)/2)), ...
-    (originPSF(3)-ceil((newImageSize(3)-1)/2)):(originPSF(3)+floor((newImageSize(3)-1)/2))); clear PSF;
+    (originPSF(3)-ceil((newImageSize(3)-1)/2)):(originPSF(3)+floor((newImageSize(3)-1)/2))); 
+
+clear PSF;
 
 % find the OTF
 smallOTF = fftshift(fftn(smallPSF)); clear smallPSF ;

Mesh/calculatePatchLength.m

@@ -22,12 +22,16 @@
 % 
 
 % calculate the minimum patchLength of the two regions
+% THESE FIRST TWO LINES ARE MOST TIME CONSUMING:
 firstFaces = faceIndex(watersheds == firstLabel);
 secondFaces = faceIndex(watersheds == secondLabel);
+
 firstSize = max([max(positions(firstFaces,1))-min(positions(firstFaces,1)), ... 
     max(positions(firstFaces,2))-min(positions(firstFaces,2)), max(positions(firstFaces,3))-min(positions(firstFaces,3))]);
+
 secondSize = max([max(positions(secondFaces,1))-min(positions(secondFaces,1)), ... 
     max(positions(secondFaces,2))-min(positions(secondFaces,2)), max(positions(secondFaces,3))-min(positions(secondFaces,3))]);
+
 patchLengthSmall = min([firstSize, secondSize, 0.2*meshLength]);
 patchLengthSmall = max([patchLengthSmall, 8]);
 patchLengthBig = max([firstSize, secondSize]);

Mesh/calculateTriangleMeasurePair.m

@@ -1,4 +1,4 @@
-function triangleMeasure = calculateTriangleMeasurePair(mesh, watersheds, watershedLabels, neighbors, closureSurfaceArea, firstRegionIndex, secondRegionIndex, patchLength, meshLength)
+function triangleMeasure = calculateTriangleMeasurePair(mesh, watersheds, watershedLabels, neighbors, closureSurfaceArea, firstRegionIndex, secondRegionIndex, patchLength, meshLength, labelIndex)
 
 % check to make sure that the patchLength isn't too large
 %
@@ -25,11 +25,13 @@
     return
 end
 
+%tic;
 % find the graph labels of the first two watersheds on the list
-labelIndex = 1:length(watershedLabels);
 gLabel1 = labelIndex(watershedLabels == firstRegionIndex);
 gLabel2 = labelIndex(watershedLabels == secondRegionIndex);  
+%toc;
 
+%tic;
 % make a list of watershed regions in which the two regions are merged 
 watershedsCombined = watersheds;
 mergeLabel = min([firstRegionIndex secondRegionIndex]);
@@ -39,9 +41,14 @@
 else
     watershedsCombined(watersheds == firstRegionIndex) = mergeLabel;
 end
+%toc;
 
+%tic;
 % find the closure surface area of the combined region
-[~, closureSurfaceAreaCombinedRegion, ~] = closeMesh(mergeLabel, mesh, watershedsCombined, neighbors);
+[~, closureSurfaceAreaCombinedRegion, ~, ~] = closeMesh(mergeLabel, mesh, watershedsCombined, neighbors);
+%toc;
 
+%tic;
 % calculate the value of the triangle measure (inspired by the law of cosines)
 triangleMeasure = (closureSurfaceArea(gLabel1)+closureSurfaceArea(gLabel2)-closureSurfaceAreaCombinedRegion)/(sqrt(closureSurfaceArea(gLabel1)*closureSurfaceArea(gLabel2)));
+%toc;

Mesh/measureCurvature.m

@@ -29,16 +29,21 @@
 
 % construct a graph of the faces
 try % this section is buggy because of irregularities in the mesh
+    tic;
     neighbors = findEdgesFaceGraph(mesh); % construct an edge list for the dual graph where the faces are nodes
+    toc;
 catch
     disp('         Warning: The graph could not be constructed!')
     neighbors = [];
     return
 end
 
 % median filter the curvature in real space
+tic;
 medianFilteredCurvature = medianFilterKD(mesh, meanCurvatureUnsmoothed, medianFilterRadius);
+toc;
 
+tic;
 % check for lingering infinities and replace them
 if max(medianFilteredCurvature) > 1000  
     maxFiniteMeanCurvature = max(medianFilteredCurvature.*isfinite(medianFilteredCurvature));
@@ -52,6 +57,9 @@
 
 % replace any NaN's
 medianFilteredCurvature(~isfinite(medianFilteredCurvature)) = 0;
+toc;
 
 % diffuse curvature on the mesh geometry
+tic;
 meanCurvatureSmoothed = smoothDataOnMesh(mesh, neighbors, medianFilteredCurvature, smoothOnMeshIterations);
+toc;

Mesh/surfaceCurvatureFast.m

@@ -95,7 +95,6 @@
 
 % should probably vectorize/arrayfun this at some point...
 for i = 1:nTri
-
     % get the coordinates of this triangle's vertices
     X = S.vertices(S.faces(i,:),:);
 

Mesh/surfaceNormalsFast.m

@@ -72,7 +72,8 @@
 end
 
 %Number of faces
-nTri = size(S.faces,1);
+nTri = size(S.faces, 1);
+nVert = size(S.vertices, 1);
 
 if nargout > 1
     %Number of vertices
@@ -82,30 +83,54 @@
 % ------ Calculate the Face Normals -------- %
 
 faceN = zeros(nTri,3);
+% vert_big = permute(repmat(S.vertices', 1, 1, nTri), [3,1,2]);
+% S.faces: nTri x 3
+% S.vertices: nVert x 3
+% 3 x nVert x nTri
+% vert_big: nTri x 3 x nVert
 
-for j = 1:nTri
-
-   %Get vertices of this face
-   X = S.vertices(S.faces(j,:),:);
-
-   %Face normal is the cross of two sides.
-   faceN(j,:) = -crossProduct(X(1,:)-X(2,:),X(2,:)-X(3,:));
-
-
+X = zeros(nTri, 3, 3);
+
+for i = 1:nTri
+    for j = 1:3
+        for k = 1:3
+            X(i,j,k) = S.vertices(S.faces(i,j), k);
+        end
+    end
 end
 
+faceN = squeeze(-cross(X(:,1,:)-X(:,2,:), X(:,2,:)-X(:,3,:),3));
+clear X;
+
 % ----- If requested, calculate the vertex normals ----- %
 
+
 if nargout > 1
 
-    vertN = zeros(nVert,3);
-
+    vertN = zeros(nVert, 3);
+    vertN_ = zeros(nVert, 3);
+    faces = sort(S.faces, 2);
+    remaining_indices = 1:nTri;
+    %tic;
     for j = 1:nVert
-
         %Average the normals of faces adjacent to this vertex        
-        vertN(j,:) = mean(faceN(any(S.faces == j,2),:));  % this is the slow line      
-
+        vertN(j,:) = mean( faceN(remaining_indices(any(faces == j,2)),:) );  % this is the slow line  
+        % Remove rows that are irrelevant:
+        if(mod(j,100) == 0)
+            remove_ind = find(all(faces <= j, 2));
+            remaining_indices(remove_ind) = [];
+            faces(remove_ind,:) = [];
+        end
     end
+    %toc;
+
+    % My method takes 45% as long as this:
+    %tic;
+    %for j = 1:nVert
+    %    vertN_(j,:) = mean(faceN(any(S.faces == j,2),:));  % this is the slow line      
+    %end
+    %toc;
+    %isequal(vertN, vertN_)
     if normalize
         %Normalize the vertex normals
         vertN = vertN ./ repmat(sqrt(dot(vertN,vertN,2)),1,3);
@@ -119,8 +144,11 @@
 
 
 % a faster cross product
+%{
 function z = crossProduct(x,y)
 z = x;
 z(:,1) = x(:,2).*y(:,3) - x(:,3).*y(:,2);
 z(:,2) = x(:,3).*y(:,1) - x(:,1).*y(:,3);
 z(:,3) = x(:,1).*y(:,2) - x(:,2).*y(:,1);
+%}
+%z = x;z(:,1) = x(:,2).*y(:,3) - x(:,3).*y(:,2);z(:,2) = x(:,3).*y(:,1) - x(:,1).*y(:,3);z(:,3) = x(:,1).*y(:,2) - x(:,2).*y(:,1);

SurfaceSegment/calculateMutualVisibilityPair.m

@@ -21,17 +21,35 @@
 % 
 % 
 
+%edgesToCheck = edgesToCheck(find(toRemove==0), :);
+%{
+edgesToCheckNew_ = [];
+for p = 1:size(edgesToCheck_,1)
+    if toRemove(p) == 0
+        edgesToCheckNew_ = [edgesToCheckNew_; edgesToCheck(p,:)];
+    end
+end
+isequaln(edgesToCheckNew, edgesToCheckNew_)
+edgesToCheck = edgesToCheckNew;
+%}
 
 % find the faces in each region
 faceIndex = 1:size(mesh.faces,1);
-firstFaces = []; secondFaces = [];
+firstFaces = faceIndex(regionLabels == firstRegionIndices);
+secondFaces = faceIndex(regionLabels == secondRegionIndices);
+%{
+firstFaces_ = []; secondFaces_ = [];
 for f = 1:length(firstRegionIndices)
-    firstFaces = [firstFaces, faceIndex(regionLabels == firstRegionIndices(f))];
+    firstFaces_ = [firstFaces_, faceIndex(regionLabels == firstRegionIndices(f))];
 end
 for f = 1:length(secondRegionIndices)
-    secondFaces = [secondFaces, faceIndex(regionLabels == secondRegionIndices(f))];
+    secondFaces_ = [secondFaces_, faceIndex(regionLabels == secondRegionIndices(f))];
 end
 
+isequal(firstFaces,firstFaces_)
+isequal(secondFaces,secondFaces_)
+%}
+
 % just go through the whole list, and if you run out return 0
 mutualVisibility = 0;
 
@@ -43,6 +61,7 @@
 else
     localSize = Inf;
 end
+
 mutualVisibilityArray = makeMutVisArray(maxPairLimitMult, localSize, mesh, firstFaces, secondFaces, patchLength, raysPerCompare);
 
 % if the array is partially empty, try again with a larger list of tests
@@ -73,7 +92,6 @@
     disp('  found too few rays to compare') % this is more problematic 
 end
 
-
 function mutualVisibilityList = makeMutVisArray(maxPairLimitMult, localSize, mesh, firstFaces, secondFaces, patchLength, raysPerCompare)
 
 % if there are many, many possible rays then subsample the pairs
@@ -97,8 +115,8 @@
 % iterate through the ray intersections
 mutualVisibilityList = NaN(1,raysPerCompare);
 rayCount = 0;
+
 for r = 1:size(pairsList,1)
-
     % find the position of the first face
     verticesFace = mesh.faces(pairsList(r,1),:);
     firstPosition = (mesh.vertices(verticesFace(1),:) + mesh.vertices(verticesFace(2),:) + mesh.vertices(verticesFace(3),:))/3;
@@ -110,7 +128,7 @@
     % find the direction and length of a ray from the first to the second face
     ray = secondPosition - firstPosition;
     rayLength = sqrt(sum(ray.^2));
-    
+
     % the local line-of-sight condition
     if (rayLength > localSize*patchLength)
         continue
@@ -120,6 +138,7 @@
     if (rayLength < 2)  
         continue
     end
+
     rayCount = rayCount + 1;
     ray = ray./rayLength;
 
@@ -129,49 +148,14 @@
     % check every face in the mesh to see if the ray and face intersect using a one-sided ray-triangle intersection algorithm
     firstMatrix = repmat([firstPosition(1) firstPosition(2) firstPosition(3)], size(mesh.faces,1), 1);
     rayMatrix = repmat([ray(1) ray(2) ray(3)], size(mesh.faces,1), 1);
+
     %[intersect, dist] = TriangleRayIntersectionFastOneSided(firstMatrix, rayMatrix, mesh.vertices(mesh.faces(:,2),:), mesh.vertices(mesh.faces(:,1),:), mesh.vertices(mesh.faces(:,3),:));
     [intersect, dist] = TriangleRayIntersection(firstMatrix, rayMatrix, mesh.vertices(mesh.faces(:,2),:), mesh.vertices(mesh.faces(:,1),:), mesh.vertices(mesh.faces(:,3),:), 'planeType', 'one sided');
+
     intersectMask = intersect & (dist > 0) & (dist <= rayLength);
     mutualVisibilityList(1,rayCount) = 1-max(intersectMask);
-
 %     % iteratively merge regions until all adjacent regions have mutual visibility below the cutoff
-     if rayCount == raysPerCompare
-%        mutualVisibility = mean(mutualVisibilityList);
-%           
-%         %% debug code
-%         %if rand(1) > 0.95 
-%         if (mutualVisibility > 0.4) && (mutualVisibility <= 0.8)
-%             %% debug code (plot the two patches with the mutual visibility displayed at the top
-%             climits = [0 Inf];
-%             cmap = jet(4);
-% 
-%             meshColor = ones(length(mesh.faces),1);
-%             meshColor(firstFaces) = 2;
-%             meshColor(secondFaces) = 4;
-% 
-%             % plot the mesh
-%             figure
-%             meshHandle = patch(mesh,'FaceColor','flat','EdgeColor','none','FaceAlpha',1);
-% 
-%             % color the mesh
-%             meshHandle.FaceVertexCData = meshColor;
-%             colormap(cmap);
-%             caxis(climits);
-% 
-%             % properly set the axis
-%             %axis([130 330 0 400 0 200]);
-%             daspect([1 1 1]); axis off;
-% 
-%             % light the scene
-%             %light_handle = camlight('headlight');
-%             camlookat(meshHandle);
-%             camlight(0,0); camlight(120,-60); camlight(240,60);
-%             lighting phong;
-% 
-%             title(num2str(mutualVisibility))
-%             1;
-%        end
-%        
+    if rayCount == raysPerCompare
         return
     end
 

SurfaceSegment/closeMesh.m

@@ -23,13 +23,16 @@
 
 
 % find the labels of the faces in the region
-faceIndex = 1:length(watersheds);
-facesInRegion = faceIndex'.*(wLabel==watersheds);
-facesInRegion = facesInRegion(facesInRegion>0);
+
+wLabel_equals_watersheds = wLabel==watersheds;
+faceIndex = (1:length(watersheds))';
+facesInRegion = faceIndex(wLabel_equals_watersheds);
 
 % make an edge list of edges associated with the region
-neighborsRegion = neighbors.*repmat(wLabel==watersheds,1,3);
-edgeList = [facesInRegion, neighborsRegion(neighborsRegion(:,1)>0,1); facesInRegion, neighborsRegion(neighborsRegion(:,2)>0,2); facesInRegion, neighborsRegion(neighborsRegion(:,3)>0,3)];
+neighborsRegion = neighbors.*repmat(wLabel_equals_watersheds,1,3);
+edgeList = [facesInRegion, neighborsRegion(neighborsRegion(:,1)>0,1); ...
+            facesInRegion, neighborsRegion(neighborsRegion(:,2)>0,2); 
+            facesInRegion, neighborsRegion(neighborsRegion(:,3)>0,3)];
 
 % find a list of the neighboring faces
 facesAllNeighbors = setdiff(edgeList(:,2), facesInRegion);
@@ -54,7 +57,7 @@
 closeCenter = mean([smoothedSurface.vertices(nVertices,1), smoothedSurface.vertices(nVertices,2), smoothedSurface.vertices(nVertices,3)], 1);
 
 % find the average distance from each edge vertex to the closeCenter
-closeRadius = mean(sqrt(sum((smoothedSurface.vertices(nVertices,:) - repmat(closeCenter, size(smoothedSurface.vertices(nVertices,:),1), 1)).^2, 2)));
+closeRadius = mean(vecnorm(smoothedSurface.vertices(nVertices,:) - repmat(closeCenter, size(smoothedSurface.vertices(nVertices,:),1), 1), 2, 2));
 
 % make a fv (faces-vertices) structure for the region  (note that the vertices are not relabeled and so the structure is large)
 closedMesh.faces = [smoothedSurface.faces(facesInRegion,1), smoothedSurface.faces(facesInRegion,2), smoothedSurface.faces(facesInRegion,3)];
@@ -78,12 +81,17 @@
 closureMesh.faces = newFaces;
 closureSurfaceArea = sum(measureAllFaceAreas(closureMesh));
 
+%{
+[closeCenter, closureSurfaceArea, closedMesh.faces, closedMesh.vertices, closeRadius] = closeMeshcpp(wLabel, smoothedSurface.faces, smoothedSurface.vertices, watersheds, neighbors);
+closureSurfaceArea = sum(measureAllFaceAreas(closedMesh));
+%}
+
+%closureSurfaceArea = sum(measureAllFaceAreas(closureMesh));
+
 % % if the closeCenter is not defined, set it to be the middle of the object
 % if isnan(closeCenter(1))
 %    nVertices = smoothedSurface.faces(facesInRegion(:,1));
 %    closeCenter = mean([smoothedSurface.vertices(nVertices,1), smoothedSurface.vertices(nVertices,2), smoothedSurface.vertices(nVertices,3)], 1); 
 % end
-
-% % Debug code
 % figure
-% patch(closedMesh)
+% patch(closedMesh)