diff --git a/Deconvolution/calculateSmallOTF.m b/Deconvolution/calculateSmallOTF.m index cf4c6e9..84eb945 100644 --- a/Deconvolution/calculateSmallOTF.m +++ b/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 ; diff --git a/Mesh/calculatePatchLength.m b/Mesh/calculatePatchLength.m index 7c8a79d..215122a 100644 --- a/Mesh/calculatePatchLength.m +++ b/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]); \ No newline at end of file diff --git a/Mesh/calculateTriangleMeasurePair.m b/Mesh/calculateTriangleMeasurePair.m index 18bb917..cfaf9e0 100644 --- a/Mesh/calculateTriangleMeasurePair.m +++ b/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; \ No newline at end of file diff --git a/Mesh/measureCurvature.m b/Mesh/measureCurvature.m index 6ce5bec..0438b88 100644 --- a/Mesh/measureCurvature.m +++ b/Mesh/measureCurvature.m @@ -29,7 +29,9 @@ % 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 = []; @@ -37,8 +39,11 @@ 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; \ No newline at end of file diff --git a/Mesh/surfaceCurvatureFast.m b/Mesh/surfaceCurvatureFast.m index f7351f4..1e8a8b4 100644 --- a/Mesh/surfaceCurvatureFast.m +++ b/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,:),:); diff --git a/Mesh/surfaceNormalsFast.m b/Mesh/surfaceNormalsFast.m index 0861762..30593f3 100644 --- a/Mesh/surfaceNormalsFast.m +++ b/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); \ No newline at end of file diff --git a/SurfaceSegment/calculateMutualVisibilityPair.m b/SurfaceSegment/calculateMutualVisibilityPair.m index 5a4275c..6af2565 100644 --- a/SurfaceSegment/calculateMutualVisibilityPair.m +++ b/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 diff --git a/SurfaceSegment/closeMesh.m b/SurfaceSegment/closeMesh.m index b49f83f..22ae408 100644 --- a/SurfaceSegment/closeMesh.m +++ b/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) \ No newline at end of file diff --git a/SurfaceSegment/closeMeshcpp.cpp b/SurfaceSegment/closeMeshcpp.cpp new file mode 100644 index 0000000..d662fb3 --- /dev/null +++ b/SurfaceSegment/closeMeshcpp.cpp @@ -0,0 +1,515 @@ +#include "mex.h" +#include +#include +#include +#include +#include +#include + +extern "C"{ + #include "matrix.h" + #include "math.h" + #include +} + +using namespace std; + +#define MAX_ELEMENTS 100000 +#define SMALLER_MAX 1000 + +static vector *neighborsRegion_1 = NULL; +static vector *neighborsRegion_2 = NULL; +static vector *neighborsRegion_3 = NULL; + +static vector *faceIndex = NULL; +static vector *edgeList_1 = NULL; + +static vector *watersheds_equals_wLabel = NULL; + +static int * facesInRegion = NULL; + +static int * neighbors_ptr_int = NULL; + +static vector * facesAllNeighbors = NULL; + +static bool * boundaryEdgesMask = NULL; + +static vector * boundaryEdgeList_1 = NULL; +static vector * boundaryEdgeList_2 = NULL; + +static int * vertexPairs_1 = NULL; +static int * vertexPairs_2 = NULL; + +static int * faces_ptr_int = NULL; + +static vector * centered_v1 = NULL; +static vector * centered_v2 = NULL; +static vector * centered_v3 = NULL; + +static vector * centered_sum = NULL; + +static vector * cm_faces_1 = NULL; +static vector * cm_faces_2 = NULL; +static vector * cm_faces_3 = NULL; + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + + omp_set_num_threads(omp_get_max_threads()); + + //ALLOCATE MEMORY FIRST: + if(!neighborsRegion_1){ + cout << "allocating" << endl; + //length facesInRegion_len + neighborsRegion_1 = new vector(MAX_ELEMENTS); + neighborsRegion_2 = new vector(MAX_ELEMENTS); + neighborsRegion_3 = new vector(MAX_ELEMENTS); + + //length watersheds_size: + faceIndex = new vector(MAX_ELEMENTS); + + //length watersheds_size + watersheds_equals_wLabel = new vector(MAX_ELEMENTS); + + //length facesInRegion_len + facesInRegion = (int *) malloc(MAX_ELEMENTS * sizeof(int)); + + //length 3*neighbors_M + neighbors_ptr_int = (int*) malloc(3*MAX_ELEMENTS*sizeof(int)); + + //length 3*facesInRegion_len + edgeList_1 = new vector(3*MAX_ELEMENTS) ; + + //length 0 + facesAllNeighbors = new vector(); + + //length N + boundaryEdgesMask = (bool*) malloc(MAX_ELEMENTS * sizeof(bool)); + + //length BEM_greaterthan_zero_length + boundaryEdgeList_1 = new vector(SMALLER_MAX); + boundaryEdgeList_2 = new vector(SMALLER_MAX); + + //length BEM_greaterthan_zero_length + vertexPairs_1 = (int *) malloc(SMALLER_MAX * sizeof(int)); + vertexPairs_2 = (int *) malloc(SMALLER_MAX * sizeof(int)); + + //length 3*faces_M + faces_ptr_int = (int*) malloc(3*MAX_ELEMENTS*sizeof(int)); + + //length BEM_greaterthan_zero_length + centered_v1 = new vector(SMALLER_MAX); + centered_v2 = new vector(SMALLER_MAX); + centered_v3 = new vector(SMALLER_MAX); + + //length BEM_greaterthan_zero_length + centered_sum = new vector(SMALLER_MAX); + + //length facesInRegion_len + cm_faces_1 = new vector(MAX_ELEMENTS); + cm_faces_2 = new vector(MAX_ELEMENTS); + cm_faces_3 = new vector(MAX_ELEMENTS); + } + + // Inputs: + const double wLabel = (mxGetPr(prhs[0]))[0]; //Integer + const mxArray * mx_faces = prhs[1]; + const mxArray * mx_vertices = prhs[2]; + const mxArray * watersheds = prhs[3]; //Vector + const mxArray * neighbors = prhs[4]; + + /* Find the labels of the faces in the region: + * faceIndex = 1:length(watersheds); + */ + size_t watersheds_size = mxGetNumberOfElements(watersheds); + iota(begin(*faceIndex), begin(*faceIndex) + watersheds_size, 1); + + /* + * facesInRegion = faceIndex'.*(wLabel==watersheds); + * facesInRegion = facesInRegion(facesInRegion>0); + */ + double * watersheds_ptr = mxGetPr(watersheds); + + //replace_copy_if(std::execution::parallel_policy, watersheds_ptr, watersheds_ptr+watersheds_size,watersheds_equals_wLabel->begin(), [&](auto i){ return i == wLabel;}, 0); + + // Convert watersheds to a vector: + int facesInRegion_len = 0; + + for(int i = 0; i < watersheds_size; i++){ + watersheds_equals_wLabel->at(i) = watersheds_ptr[i] == wLabel; + if(watersheds_equals_wLabel->at(i)){ + ++facesInRegion_len; + } + } + + int k = 0; + for(int i = 0; i < watersheds_size; ++i) + if(watersheds_equals_wLabel->at(i)){ + facesInRegion[k] = faceIndex->at(i); + ++k; + } + + + /* + * 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)]; + */ + + double * neighbors_ptr = mxGetPr(neighbors); + size_t neighbors_M = mxGetM(neighbors); + + for (int i = 0; i < 3 * neighbors_M; ++i) { + neighbors_ptr_int[i] = (int)neighbors_ptr[i]; + } + + //vector *neighbors_1 = new vector(neighbors_ptr_int, neighbors_ptr_int + neighbors_M); + //vector *neighbors_2 = new vector(neighbors_ptr_int + neighbors_M, neighbors_ptr_int + 2*neighbors_M); + //vector *neighbors_3 = new vector(neighbors_ptr_int + 2*neighbors_M, neighbors_ptr_int + 3*neighbors_M); + + int * neighbors_1 = neighbors_ptr_int; + int * neighbors_2 = neighbors_ptr_int + neighbors_M; + int * neighbors_3 = neighbors_ptr_int + 2*neighbors_M; + + //vector *neighborsRegion_1 = new vector(facesInRegion_len); + //vector *neighborsRegion_2 = new vector(facesInRegion_len); + //vector *neighborsRegion_3 = new vector(facesInRegion_len); + + k = 0; + for(int i = 0; i < watersheds_size; ++i) + if(watersheds_equals_wLabel->at(i)){ + //length: facesInRegion_len + neighborsRegion_1->at(k) = neighbors_1[i]; + neighborsRegion_2->at(k) = neighbors_2[i]; + neighborsRegion_3->at(k) = neighbors_3[i]; + ++k; + } + + memcpy(edgeList_1->data() , facesInRegion, facesInRegion_len*sizeof(int)); + memcpy(edgeList_1->data() + facesInRegion_len , facesInRegion, facesInRegion_len*sizeof(int)); + memcpy(edgeList_1->data() + 2*facesInRegion_len , facesInRegion, facesInRegion_len*sizeof(int)); + + vector * edgeList_2 = new vector(*neighborsRegion_1); + edgeList_2->insert(edgeList_2->begin() + facesInRegion_len, neighborsRegion_2->begin(), neighborsRegion_2->begin()+facesInRegion_len); + edgeList_2->insert(edgeList_2->begin() + 2*facesInRegion_len, neighborsRegion_3->begin(), neighborsRegion_3->begin()+facesInRegion_len); + + /* + * Find a list of the neighboring faces + * facesAllNeighbors = setdiff(edgeList(:,2), facesInRegion); + */ + + + vector * facesInRegion_vect = new vector(facesInRegion, facesInRegion + facesInRegion_len); + vector * edgeList_2_sorted = new vector(*edgeList_2); + + // Need to sort input vectors and delete duplicates for this to work: + + sort(edgeList_2_sorted->begin(), edgeList_2_sorted->begin()+3*facesInRegion_len); + edgeList_2_sorted->erase( + unique( edgeList_2_sorted->begin(), edgeList_2_sorted->begin()+3*facesInRegion_len ), edgeList_2_sorted->begin()+3*facesInRegion_len + ); + + sort(facesInRegion_vect->begin(), facesInRegion_vect->begin()+facesInRegion_len); + facesInRegion_vect->erase( + unique( facesInRegion_vect->begin(), facesInRegion_vect->begin()+facesInRegion_len), facesInRegion_vect->begin()+facesInRegion_len + ); + + facesAllNeighbors->clear(); + // These "end's" are okay + set_difference(edgeList_2_sorted->begin(), edgeList_2_sorted->end(), + facesInRegion_vect->begin(), facesInRegion_vect->end(), + std::inserter(*facesAllNeighbors, facesAllNeighbors->begin())); + /* + * boundaryEdgesMask = ismembc(edgeList(:,2),facesAllNeighbors); + */ + + int * in1 = edgeList_2->data(); + int * in2 = facesAllNeighbors->data(); + + size_t N = 3*facesInRegion_len; + + for(int i = 0; i < N; ++i){ + boundaryEdgesMask[i] = find(facesAllNeighbors->begin(), facesAllNeighbors->end(), edgeList_2->at(i)) != facesAllNeighbors->end(); + } + + /* + * Find the edges that connect the regions to its neighors + * boundaryEdgeList = [edgeList(boundaryEdgesMask>0,1), edgeList(boundaryEdgesMask>0,2)]; + */ + + size_t BEM_greaterthan_zero_length = 0; + + for(int i = 0; i < N; ++i) + if(boundaryEdgesMask[i]) + ++BEM_greaterthan_zero_length; + + k = 0; + for(int i = 0; i < 3*facesInRegion_len; ++i){ + if(boundaryEdgesMask[i]){ + boundaryEdgeList_1->at(k) = edgeList_1->at(i); + boundaryEdgeList_2->at(k) = edgeList_2->at(i); + ++k; + } + } + + /* + * vertexPairs = zeros(size(boundaryEdgeList,1),2); + */ + + /* + % find the vertices that correspond to each boundary edge + for b = 1:size(boundaryEdgeList,1) + vertexPairs(b,:) = intersect(smoothedSurface.faces(boundaryEdgeList(b,1),:), smoothedSurface.faces(boundaryEdgeList(b,2),:), 'stable'); + + % maintain the order of the vertices in the pair so that the vertex directionality convention for normality will be preserved during closure + if vertexPairs(b,1) == smoothedSurface.faces(boundaryEdgeList(b,1),1) && vertexPairs(b,2) == smoothedSurface.faces(boundaryEdgeList(b,1),3) + vertexPairs(b,:) = fliplr(vertexPairs(b,:)); + end + end + */ + + double * faces_ptr = mxGetPr(mx_faces); + size_t faces_M = mxGetM(mx_faces); + + for (int i = 0; i < 3 * faces_M; ++i) { + faces_ptr_int[i] = (int)faces_ptr[i]; + } + + //vector * faces_1 = new vector(faces_ptr_int, faces_ptr_int + faces_M); + int * faces_1 = faces_ptr_int; + int * faces_2 = faces_ptr_int + faces_M; + int * faces_3 = faces_ptr_int + 2*faces_M; + + for(int b = 0; b < BEM_greaterthan_zero_length; ++b){ + unordered_set f1, f2, vertexPairs_set; + f1.insert(faces_1[boundaryEdgeList_1->at(b)-1]); + f2.insert(faces_1[boundaryEdgeList_2->at(b)-1]); + f1.insert(faces_2[boundaryEdgeList_1->at(b)-1]); + f2.insert(faces_2[boundaryEdgeList_2->at(b)-1]); + f1.insert(faces_3[boundaryEdgeList_1->at(b)-1]); + f2.insert(faces_3[boundaryEdgeList_2->at(b)-1]); + + copy_if(f1.begin(), f1.end(), inserter(vertexPairs_set, vertexPairs_set.begin()), [f2](const int element){return f2.count(element) > 0;} ); + + vertexPairs_1[b] = *(vertexPairs_set.begin()); + vertexPairs_2[b] = *next(vertexPairs_set.begin(), 1); + + if(vertexPairs_2[b] < vertexPairs_1[b]){ + int tmp = vertexPairs_1[b]; + vertexPairs_1[b] = vertexPairs_2[b]; + vertexPairs_2[b] = tmp; + } + + if ((vertexPairs_1[b] == faces_1[boundaryEdgeList_1->at(b)-1] && + vertexPairs_2[b] == faces_3[boundaryEdgeList_1->at(b)-1])){ + int tmp = vertexPairs_1[b]; + vertexPairs_1[b] = vertexPairs_2[b]; + vertexPairs_2[b] = tmp; + } + } + + /* + * Find the center of mass of the vertices + * nVertices = vertexPairs(:,1); + * 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))); + * + */ + + double * vertices_ptr = mxGetPr(mx_vertices); + size_t vertices_M = mxGetM(mx_vertices); + + double * vertices_1 = vertices_ptr; + double * vertices_2 = vertices_ptr + vertices_M; + double * vertices_3 = vertices_ptr + 2*vertices_M; + + double closeCenter_1 = 0, closeCenter_2 = 0, closeCenter_3 = 0; + + for(int b = 0; b < BEM_greaterthan_zero_length; ++b){ + centered_v1->at(b) = (double) vertices_1[vertexPairs_1[b]-1]; + centered_v2->at(b) = (double) vertices_2[vertexPairs_1[b]-1]; + centered_v3->at(b) = (double) vertices_3[vertexPairs_1[b]-1]; + + closeCenter_1 += (double) vertices_1[vertexPairs_1[b]-1]; + closeCenter_2 += (double) vertices_2[vertexPairs_1[b]-1]; + closeCenter_3 += (double) vertices_3[vertexPairs_1[b]-1]; + } + + closeCenter_1 /= (double) BEM_greaterthan_zero_length; + closeCenter_2 /= (double) BEM_greaterthan_zero_length; + closeCenter_3 /= (double) BEM_greaterthan_zero_length; + for (int i = 0; i < BEM_greaterthan_zero_length; ++i) { + centered_v1->at(i) = pow(centered_v1->at(i) - closeCenter_1, 2); + centered_v2->at(i) = pow(centered_v2->at(i) - closeCenter_2, 2); + centered_v3->at(i) = pow(centered_v3->at(i) - closeCenter_3, 2); + } + + // Subtract means: + //for(auto & element : *centered_v1) element -= closeCenter_1; + //for(auto & element : *centered_v2) element -= closeCenter_2; + //for(auto & element : *centered_v3) element -= closeCenter_3; + + // Square each element: + //for(auto & element : *centered_v1) element *= element; + //for(auto & element : *centered_v2) element *= element; + //for(auto & element : *centered_v3) element *= element; + + for (int i = 0; i < BEM_greaterthan_zero_length; ++i) { + centered_sum->at(i) = sqrt(centered_v1->at(i) + centered_v2->at(i) + centered_v3->at(i)); + } + + //for(auto & element : *centered_sum) element = sqrt(element); + + // Finally, get the mean: + double closeRadius_ = accumulate(centered_sum->begin(), centered_sum->begin()+BEM_greaterthan_zero_length, 0) / (double) BEM_greaterthan_zero_length; + + /* + * 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)]; + */ + + for(int i = 0; i < facesInRegion_len; ++i){ + cm_faces_1->at(i) = faces_1[facesInRegion[i]-1]; + cm_faces_2->at(i) = faces_2[facesInRegion[i]-1]; + cm_faces_3->at(i) = faces_3[facesInRegion[i]-1]; + } + + /* + * closedMesh.vertices = smoothedSurface.vertices; + * + * find a unique label for the new vertex: + * closeCenterLabel = size(closedMesh.vertices,1)+1; + * + * append the new vertex to the list of vertices + * closedMesh.vertices(closeCenterLabel,:) = closeCenter; + * + */ + + //HEREEEE + + + /* + * swap the order of the vertices in the pairs to maintain the directionality of the surface normal + * vertexPairs = fliplr(vertexPairs); + */ + + int * tmp = vertexPairs_1; + vertexPairs_1 = vertexPairs_2; + vertexPairs_2 = tmp; + + /* + * append the new faces to the list of faces + * newFaces = [vertexPairs, closeCenterLabel.*ones(size(vertexPairs,1),1)]; // Add a third column of all value closeCenterLabel + * closedMesh.faces = [closedMesh.faces; newFaces]; // Add the new faces to the bottom of this matrix: + * closureMesh.faces = newFaces; + * + * closureMesh.vertices = closedMesh.vertices; + * [closeCenter, closureSurfaceArea, closedMesh, closeRadius] + */ + + size_t fsize = facesInRegion_len + BEM_greaterthan_zero_length; + + mxArray * closeCenter = mxCreateDoubleMatrix(1, 3, mxREAL); + mxArray * closureSurfaceArea = mxCreateDoubleMatrix(1, 1, mxREAL); + mxArray * cmFaces = mxCreateDoubleMatrix(fsize, 3, mxREAL); + mxArray * cmVertices = mxCreateDoubleMatrix(vertices_M + 1, 3, mxREAL); + mxArray * closeRadius = mxCreateDoubleMatrix(1, 1, mxREAL); + + double * cmVertices_pr = mxGetPr(cmVertices); + memcpy(cmVertices_pr , vertices_1, vertices_M*sizeof(double)); + memcpy(cmVertices_pr + 1 + vertices_M , vertices_2, vertices_M*sizeof(double)); + memcpy(cmVertices_pr + 2 + 2*vertices_M, vertices_3, vertices_M*sizeof(double)); + + cmVertices_pr[vertices_M] = closeCenter_1; + cmVertices_pr[2*vertices_M + 1] = closeCenter_2; + cmVertices_pr[3*vertices_M + 2] = closeCenter_3; + + double * closeCenter_pr = mxGetPr(closeCenter); + + closeCenter_pr[0] = closeCenter_1; + closeCenter_pr[1] = closeCenter_2; + closeCenter_pr[2] = closeCenter_3; + + double * cmFaces_pr = mxGetPr(cmFaces); + + for(int i = 0; i < facesInRegion_len; ++i){ + cmFaces_pr[i] = (double)cm_faces_1->at(i); + cmFaces_pr[i + fsize] = (double)cm_faces_2->at(i); + cmFaces_pr[i + 2*fsize] = (double)cm_faces_3->at(i); + } + + for(int i = facesInRegion_len; i < fsize; ++i){ + cmFaces_pr[i] = vertexPairs_1[i - facesInRegion_len]; + cmFaces_pr[i + fsize] = vertexPairs_2[i - facesInRegion_len]; + cmFaces_pr[i + 2*fsize] = (double)vertices_M + 1; + + //cout << vertexPairs_1[i - facesInRegion_len] << endl; + } + + /* + * closureSurfaceArea = sum(measureAllFaceAreas(closureMesh)); + */ + + *(mxGetPr(closeRadius)) = closeRadius_; + + /* + unordered_set x, y, z; + x.insert(1126); + x.insert(1151); + x.insert(1127); + y.insert(1127); + y.insert(1151); + y.insert(1152); + + copy_if(x.begin(), x.end(), inserter(z, z.begin()), [y](const int element){return y.count(element) > 0;} ); + + for(auto & v : z) cout << v << endl; + */ + + plhs[0] = closeCenter; + plhs[1] = closureSurfaceArea; + plhs[2] = cmFaces; + plhs[3] = cmVertices; + plhs[4] = closeRadius; + + /* + delete faceIndex; + delete neighbors_1; + delete neighbors_2; + delete neighbors_3; + delete edgeList_2; + delete facesInRegion_vect; + delete facesAllNeighbors; + delete edgeList_2_sorted; + delete boundaryEdgeList_1; + delete boundaryEdgeList_2; + delete faces_1; + delete faces_2; + delete faces_3; + delete vertices_1; + delete vertices_2; + delete vertices_3; + delete centered_v1; + delete centered_v2; + delete centered_v3; + delete centered_sum; + delete cm_faces_1; + delete cm_faces_2; + delete cm_faces_3; + delete cm_vertices_1; + delete cm_vertices_2; + delete cm_vertices_3; + + free(facesInRegion); + free(boundaryEdgesMask); + free(vertexPairs_1); + free(vertexPairs_2); + free(neighbors_ptr_int); + free(watersheds_equals_wLabel); + */ + + return; +} \ No newline at end of file diff --git a/SurfaceSegment/closeMeshcpp.mexw64 b/SurfaceSegment/closeMeshcpp.mexw64 new file mode 100644 index 0000000..4b7d490 Binary files /dev/null and b/SurfaceSegment/closeMeshcpp.mexw64 differ diff --git a/SurfaceSegment/closeMeshcpp.mexw64.pdb b/SurfaceSegment/closeMeshcpp.mexw64.pdb new file mode 100644 index 0000000..34a5e0a Binary files /dev/null and b/SurfaceSegment/closeMeshcpp.mexw64.pdb differ diff --git a/SurfaceSegment/joinWatershedSpillDepth.m b/SurfaceSegment/joinWatershedSpillDepth.m index 93f64ce..39afc4a 100644 --- a/SurfaceSegment/joinWatershedSpillDepth.m +++ b/SurfaceSegment/joinWatershedSpillDepth.m @@ -77,7 +77,8 @@ ridgeHeights(mergeDestroyIndex) = 0; spillNeighbors(mergeDestroyIndex) = 0; spillNeighbors(spillNeighbors==mergeDestroy) = mergeLabel; - [spillDepths(mergeLabelIndex), spillNeighbors(mergeLabelIndex), ridgeHeights(mergeLabelIndex)] = measureDepthOneRegion(mergeLabelIndex, neighbors, watersheds, watershedLabels, watershedGraph, measure); + faceIndex = transpose(1:length(watersheds)); + [spillDepths(mergeLabelIndex), spillNeighbors(mergeLabelIndex), ridgeHeights(mergeLabelIndex)] = measureDepthOneRegion(mergeLabelIndex, neighbors, watersheds, watershedLabels, watershedGraph, measure, faceIndex); % check if there are more regions to be merged keepMerging = 0; diff --git a/SurfaceSegment/joinWatershedTriangle.m b/SurfaceSegment/joinWatershedTriangle.m index b845356..f3e37f8 100644 --- a/SurfaceSegment/joinWatershedTriangle.m +++ b/SurfaceSegment/joinWatershedTriangle.m @@ -79,6 +79,7 @@ end % find the closure surface area of the combined region + faceIndex = transpose(1:length(watershedsCombined)); [~, closureSurfaceAreaCombinedRegion, ~] = closeMesh(mergeLabel, mesh, watershedsCombined, neighbors); % calculate the value of the triangle measure (inspired by the law of cosines) diff --git a/SurfaceSegment/joinWatershedTriangleLOS.m b/SurfaceSegment/joinWatershedTriangleLOS.m index 87c23b3..66774cc 100644 --- a/SurfaceSegment/joinWatershedTriangleLOS.m +++ b/SurfaceSegment/joinWatershedTriangleLOS.m @@ -67,10 +67,27 @@ closureSurfaceArea = NaN(length(watershedLabels),1); for w = 1:length(watershedLabels) if watershedLabels(w) ~= 0 - [~, closureSurfaceArea(w), ~] = closeMesh(watershedLabels(w), mesh, watersheds, neighbors); + [~, closureSurfaceArea(w), ~, ~] = closeMesh(watershedLabels(w), mesh, watersheds, neighbors); end end +%save('closureSurfaceArea.mat', 'closureSurfaceArea'); +%load('closureSurfaceArea.mat'); + +%{ +closureSurfaceArea_ = NaN(length(watershedLabels),1); +faceIndex = transpose(1:length(watersheds)); +tic; +for w = 1:length(watershedLabels) + if watershedLabels(w) ~= 0 + closureSurfaceArea_(w) = closeMesh_(watershedLabels(w), mesh, watersheds, neighbors, faceIndex); + end +end +toc; +%} +%isequaln(closureSurfaceArea,closureSurfaceArea_) +%closureSurfaceArea_ = closureSurfaceArea; + % calculate the patch length of the edge pairs patchLengthSmall = zeros(size(edgesToCheck, 1), 1); patchLengthBig = zeros(size(edgesToCheck, 1), 1); @@ -78,12 +95,18 @@ for p = 1:size(edgesToCheck,1) [patchLengthSmall(p,1), patchLengthBig(p,1)] = calculatePatchLength(positions, watersheds, faceIndex, edgesToCheck(p,1), edgesToCheck(p,2), meshLength); end +%save('patchLengthSmall.mat', 'patchLengthSmall'); +%load('patchLengthSmall.mat'); % measure the closure surface area of all pairs of adjacent regions found triangleMeasure = NaN(size(edgesToCheck,1),1); +labelIndex = 1:length(watershedLabels); for p = 1:size(edgesToCheck,1) - triangleMeasure(p) = calculateTriangleMeasurePair(mesh, watersheds, watershedLabels, neighbors, closureSurfaceArea, edgesToCheck(p,1), edgesToCheck(p,2), patchLengthBig(p,1), meshLength); + triangleMeasure(p) = calculateTriangleMeasurePair(mesh, watersheds, watershedLabels, neighbors, closureSurfaceArea, edgesToCheck(p,1), edgesToCheck(p,2), patchLengthBig(p,1), meshLength, labelIndex); end +%save('triangleMeasure.mat', 'triangleMeasure'); +%load('triangleMeasure.mat'); + edgesToCheck(:,3) = NaN(1,length(triangleMeasure)); edgesToCheck(:,4) = triangleMeasure(:,1); edgesToCheck(:,5) = patchLengthSmall; diff --git a/SurfaceSegment/measureDepthOneRegion.m b/SurfaceSegment/measureDepthOneRegion.m index 1327136..029577e 100644 --- a/SurfaceSegment/measureDepthOneRegion.m +++ b/SurfaceSegment/measureDepthOneRegion.m @@ -1,5 +1,4 @@ -function [spillDepth, spillNeighbor, ridgeHeight] = measureDepthOneRegion(w, faceNeighbors, watersheds, watershedLabels, watershedGraph, measure) - +function [spillDepth, spillNeighbor, ridgeHeight] = measureDepthOneRegion(w, faceNeighbors, watersheds, watershedLabels, watershedGraph, measure, faceIndex) % measureDepthOneRegion - find the watershed depth and spill neighbor of a single region (here w is the label of the watershedRegion being measured) % % Copyright (C) 2019, Danuser Lab - UTSouthwestern @@ -44,37 +43,67 @@ return end +%tic; % remove the flat regions from the list of neighbors nLabels = nLabels(nLabels>0); % if only the flat region is a neighbor set the spillDepth to Inf and return if isempty(nLabels) - spillDepth = Inf; + spillDepth = Inf; spillNeighbor = 0; ridgeHeight = 0; return end +%toc; + +%tic; + +wLabel_equals_watersheds = wLabel==watersheds; % find the lowest value of the measure in the region -lowestValue = min(measure(wLabel==watersheds)); +lowestValue = min(measure(wLabel_equals_watersheds)); % find the labels of the faces in the region -faceIndex = 1:length(watersheds); -facesInRegion = faceIndex'.*(wLabel==watersheds); -facesInRegion = facesInRegion(facesInRegion>0); +% THIS METHOD WORKS BETTER THAN ORIGINAL: +facesInRegion = faceIndex(wLabel_equals_watersheds); +%facesInRegion_ = faceIndex.*(wLabel_equals_watersheds); +%facesInRegion_ = facesInRegion(facesInRegion>0); +%isequal(facesInRegion, facesInRegion_) +%tic; % make an edge list of edges associated with the region -neighborsRegion = faceNeighbors.*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 = faceNeighbors.*repmat(wLabel_equals_watersheds,1,3); +%toc; + +%tic; +% make an edge list of edges associated with the region +%neighborsRegion_ = bsxfun(@times, faceNeighbors, wLabel_equals_watersheds); +%toc; +%isequal(neighborsRegion, neighborsRegion_) + +%tic; +edgeList = [facesInRegion, neighborsRegion(neighborsRegion(:,1)>0,1);facesInRegion, neighborsRegion(neighborsRegion(:,2)>0,2);facesInRegion, neighborsRegion(neighborsRegion(:,3)>0,3)]; +%toc; + +%tic; +%edgeList = zeros(3 * size(facesInRegion,1), 2); +%edgeList(:,1) = repmat(facesInRegion, 3, 1); +%edgeList(:,2) = neighborsRegion(neighborsRegion > 0); +%toc; +%isequal(edgeList_, edgeList) % find the spillover depth for each neighbor depthsNeighbors = zeros(length(nLabels),1); heightNeighbors = zeros(length(nLabels),1); + for n=1:length(nLabels) - % find the labels of the faces in this neighboring region - facesInNeighbor = faceIndex'.*(nLabels(n)==watersheds); - facesInNeighbor = facesInNeighbor(facesInNeighbor>0); + + % THIS CALCULATION METHOD IS MORE EFFICIENT THAN ORIGINAL: + %facesInNeighbor_ = faceIndex.*(nLabels(n)==watersheds); + %facesInNeighbor_ = facesInNeighbor(facesInNeighbor>0); + facesInNeighbor = faceIndex(watersheds == nLabels(n)); + %isequal(facesInNeighbor, facesInNeighbor_) % find the edges that connect the watershed to the neighor boundaryEdgesMask = ismembc(edgeList(:,2),facesInNeighbor); @@ -85,7 +114,7 @@ for b=1:size(boundaryEdges,1) depthsEdges(b) = max(measure(boundaryEdges(b,1)), measure(boundaryEdges(b,2))); end - + % if something goes wrong, make the depth infinite if isempty(depthsEdges) depthsNeighbors(n,1) = Inf; @@ -94,10 +123,12 @@ depthsNeighbors(n,1) = min(depthsEdges)-lowestValue; heightNeighbors(n,1) = max(depthsEdges)-lowestValue; end - end + +%tic; % find the spill neighbor of the watershed region [spillDepth, spillNeighborIndex] = min(depthsNeighbors); spillNeighbor = nLabels(spillNeighborIndex(1)); -ridgeHeight = heightNeighbors(spillNeighborIndex,1); \ No newline at end of file +ridgeHeight = heightNeighbors(spillNeighborIndex,1); +%toc; \ No newline at end of file diff --git a/SurfaceSegment/measureDepthsAll.m b/SurfaceSegment/measureDepthsAll.m index 177045a..fe88425 100644 --- a/SurfaceSegment/measureDepthsAll.m +++ b/SurfaceSegment/measureDepthsAll.m @@ -28,10 +28,18 @@ spillNeighbors = zeros(numWatersheds, 1); ridgeHeights = zeros(numWatersheds, 1); +faceIndex = transpose(1:length(watersheds)); + % iterate through the regions +tic; for w = 1:numWatersheds - % measure the spill depth and spill neighbor for the region - [spillDepths(w,1), spillNeighbors(w,1), ridgeHeights(w,1)] = measureDepthOneRegion(w, faceNeighbors, watersheds, watershedLabels, watershedGraph, measure); - + [spillDepths(w,1), spillNeighbors(w,1), ridgeHeights(w,1)] = measureDepthOneRegion(w, faceNeighbors, watersheds, watershedLabels, watershedGraph, measure, faceIndex); end +toc; + +%tic; +%for w = 1:numWatersheds +% [spillDepths(w,1), spillNeighbors(w,1), ridgeHeights(w,1)] = measureDepthOneRegion(w, faceNeighbors, watersheds, watershedLabels, watershedGraph, measure); +%end +%toc; diff --git a/SurfaceSegment/mergeRegionPairTriangleLOS.m b/SurfaceSegment/mergeRegionPairTriangleLOS.m index 7e42963..191cd39 100644 --- a/SurfaceSegment/mergeRegionPairTriangleLOS.m +++ b/SurfaceSegment/mergeRegionPairTriangleLOS.m @@ -47,28 +47,34 @@ end % update the list of closure surface areas -[~, closureSurfaceAreaCombinedRegion, ~] = closeMesh(mergeLabel, mesh, watersheds, neighbors); +faceIndex = transpose(1:length(watersheds)); +[~, closureSurfaceAreaCombinedRegion, ~, ~] = closeMesh(mergeLabel, mesh, watersheds, neighbors); closureSurfaceArea(mergeLabelIndex) = closureSurfaceAreaCombinedRegion; closureSurfaceArea(mergeDestroyIndex) = NaN; % remove all instances of both watersheds in the list of edges to check toRemove = logical( (edgesToCheck(:,1)==mergeLabel) + (edgesToCheck(:,2)==mergeLabel) + ... (edgesToCheck(:,1)==mergeDestroy) + (edgesToCheck(:,2)==mergeDestroy) ); -edgesToCheckNew = []; -for p = 1:size(edgesToCheck,1) + +% SAVES APPROXIMATELY 88 SECONDS FOR ME: +edgesToCheck = edgesToCheck(find(toRemove==0), :); +%{ +edgesToCheckNew_ = []; +for p = 1:size(edgesToCheck_,1) if toRemove(p) == 0 - edgesToCheckNew = [edgesToCheckNew; edgesToCheck(p,:)]; + edgesToCheckNew_ = [edgesToCheckNew_; edgesToCheck(p,:)]; end end +isequaln(edgesToCheckNew, edgesToCheckNew_) edgesToCheck = edgesToCheckNew; +%} % calculate the triangle measure and mutual visibility for the pairs to be added and append the new pairs to the list nLabels = watershedGraph{mergeLabelIndex}; pairsToAdd = [mergeLabel.*ones(length(nLabels),1), nLabels']; -faceIndex = 1:length(watersheds); for p = 1:size(pairsToAdd, 1) [patchLengthSmall, patchLengthBig] = calculatePatchLength(positions, watersheds, faceIndex, pairsToAdd(p,1), pairsToAdd(p,2), meshLength); mutVis = calculateMutualVisibilityPair(mesh, positions, watersheds, pairsToAdd(p,1), pairsToAdd(p,2), patchLengthSmall, raysPerCompare, 1); - triMeas = calculateTriangleMeasurePair(mesh, watersheds, watershedLabels, neighbors, closureSurfaceArea, pairsToAdd(p,1), pairsToAdd(p,2), patchLengthBig, meshLength); + triMeas = calculateTriangleMeasurePair(mesh, watersheds, watershedLabels, neighbors, closureSurfaceArea, pairsToAdd(p,1), pairsToAdd(p,2), patchLengthBig, meshLength, labelIndex); edgesToCheck = [edgesToCheck; pairsToAdd(p,1:2), mutVis, triMeas, patchLengthSmall, patchLengthBig]; end \ No newline at end of file diff --git a/extern/TriangleRayIntersection.m b/extern/TriangleRayIntersection.m index 49d2d0f..5a477c2 100644 --- a/extern/TriangleRayIntersection.m +++ b/extern/TriangleRayIntersection.m @@ -73,25 +73,25 @@ % Jarek Tuszynski (jaroslaw.w.tuszynski@leidos.com) % % License: BSD license (http://en.wikipedia.org/wiki/BSD_licenses) -% -% Copyright (C) 2019, Danuser Lab - UTSouthwestern -% -% This file is part of Morphology3DPackage. -% -% Morphology3DPackage is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Morphology3DPackage is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Morphology3DPackage. If not, see . -% -% +% +% Copyright (C) 2019, Danuser Lab - UTSouthwestern +% +% This file is part of Morphology3DPackage. +% +% Morphology3DPackage is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Morphology3DPackage is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Morphology3DPackage. If not, see . +% +% %% Transpose inputs if needed @@ -108,7 +108,6 @@ if (size(vert0,1)==1 && N>1 && size(vert0,2)==3), vert0 = repmat(vert0, N, 1); end if (size(vert1,1)==1 && N>1 && size(vert1,2)==3), vert1 = repmat(vert1, N, 1); end if (size(vert2,1)==1 && N>1 && size(vert2,2)==3), vert2 = repmat(vert2, N, 1); end - %% Check if all the sizes match SameSize = (any(size(orig)==size(vert0)) && ... any(size(orig)==size(vert1)) && ... @@ -156,7 +155,6 @@ k = k+1; end end - %% Set up border parameter switch border case 'normal' @@ -172,7 +170,6 @@ %% initialize default output intersect = false(size(orig,1),1); % by default there are no intersections t = inf+zeros(size(orig,1),1); u=t; v=t; - %% Find faces parallel to the ray edge1 = vert1-vert0; % find vectors for two edges sharing vert0 edge2 = vert2-vert0; @@ -188,7 +185,6 @@ error('Triangle parameter must be either "one sided" or "two sided"'); end if all(~angleOK), return; end % if all parallel than no intersections - %% Different behavior depending on one or two sided triangles det(~angleOK) = nan; % change to avoid division by zero u = sum(tvec.*pvec,2)./det; % 1st barycentric coordinate @@ -216,7 +212,6 @@ % test if line/plane intersection is within the triangle ok = (ok & v>=-zero & u+v<=1.0+zero); end - %% Test where along the line the line/plane intersection occurs switch lineType case 'line' % infinite line @@ -236,4 +231,4 @@ xcoor(ok,:) = vert0(ok,:) ... + edge1(ok,:).*repmat(u(ok,1),1,3) ... + edge2(ok,:).*repmat(v(ok,1),1,3); -end +end \ No newline at end of file diff --git a/patchClassification/clickOnProtrusions.m b/patchClassification/clickOnProtrusions.m index d7b4089..26f63d0 100644 --- a/patchClassification/clickOnProtrusions.m +++ b/patchClassification/clickOnProtrusions.m @@ -248,4 +248,4 @@ function clickOnProtrusions(p) save(fullfile(savePath,'blebLocs.mat'), 'locations', 'frameIndex', 'pClick') end -end +end \ No newline at end of file diff --git a/patchClassification/trainProtrusionsClassifier.m b/patchClassification/trainProtrusionsClassifier.m index 0c49a4f..d598652 100644 --- a/patchClassification/trainProtrusionsClassifier.m +++ b/patchClassification/trainProtrusionsClassifier.m @@ -149,8 +149,10 @@ function trainProtrusionsClassifier(p) end % make a kd-tree to store the mesh faces - treeMesh = kdtree_build(positions); - + % The following no longer works for me: + %treeMesh = kdtree_build(positions); + treeMesh = KDTreeSearcher(positions); + % make a list of bleb indices that have been clicked on for the clickOnAll mode if strcmp(clickMode, 'clickOnAll') @@ -165,8 +167,9 @@ function trainProtrusionsClassifier(p) for b = 1:length(blebClickLocs) % find the closest face to the click - faceIndex = kdtree_k_nearest_neighbors(treeMesh, blebClickLocs(b,:), 1); - + %faceIndex = kdtree_k_nearest_neighbors(treeMesh, blebClickLocs(b,:), 1); + faceIndex = knnsearch(treeMesh, locations{f}.blebClickLocs(b,:), 'K', 1); + % find the bleb label of that watershed clickedOnProtrusionsFrame(b) = segment(faceIndex); @@ -181,8 +184,9 @@ function trainProtrusionsClassifier(p) for b = 1:length(locations{f}.notBlebs) % find the closest face to the click - faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.notBlebs(b,:), 1); - + %faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.notBlebs(b,:), 1); + faceIndex = knnsearch(treeMesh, locations{f}.notBlebs(b,:), 'K', 1); + % find the bleb label of that watershed clickedOnNotProtrusionsFrame(b) = segment(faceIndex); @@ -201,7 +205,9 @@ function trainProtrusionsClassifier(p) for b = 1:size(locations{f}.blebs,1) % find the closest face to the click - faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.blebs(b,:), 1); + % Matlab crashes everytime I run: + %faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.blebs(b,:), 1); + faceIndex = knnsearch(treeMesh, locations{f}.blebs(b,:), 'K', 1); % find the bleb label of that watershed clickedOnProtrusionsFrame(b) = segment(faceIndex); @@ -214,7 +220,8 @@ function trainProtrusionsClassifier(p) for b = 1:size(locations{f}.notBlebs,1) % find the closest face to the click - faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.notBlebs(b,:), 1); + %faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.notBlebs(b,:), 1); + faceIndex = knnsearch(treeMesh, locations{f}.notBlebs(b,:), 'K', 1); % find the bleb label of that watershed clickedOnNotProtrusionsFrame(b) = segment(faceIndex); @@ -230,8 +237,8 @@ function trainProtrusionsClassifier(p) for b = 1:size(locations{f}.protrusions{n},1) % find the closest face to the click - faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.protrusions{n}(b,:), 1); - + %faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.protrusions{n}(b,:), 1); + faceIndex = knnsearch(treeMesh, locations{f}.protrusions{n}(b,:), 'K', 1); % find the patch label of that watershed clickedOnProtrusionsFrame{n}(b) = segment(faceIndex); end @@ -244,8 +251,9 @@ function trainProtrusionsClassifier(p) for b = 1:size(locations{f}.notProtrusions,1) % find the closest face to the click - faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.notProtrusions(b,:), 1); - + %faceIndex = kdtree_k_nearest_neighbors(treeMesh, locations{f}.notProtrusions(b,:), 1); + faceIndex = knnsearch(treeMesh, locations{f}.notProtrusions(b,:), 'K', 1); + % find the patch label of that watershed clickedOnNotProtrusionsFrame(b) = segment(faceIndex); end @@ -347,4 +355,4 @@ function trainProtrusionsClassifier(p) save(fullfile(p.saveDirectory,[p.saveNameCells '.mat']), 'clickedOnNotProtrusions'); else save(fullfile(p.saveDirectory,[p.saveNameCells '.mat']), 'clickedOnProtrusions', 'clickedOnNotProtrusions'); -end +end \ No newline at end of file diff --git a/patchDescribe/calculateShapeDiameterFunctionRegion.m b/patchDescribe/calculateShapeDiameterFunctionRegion.m index 9b6ab43..1da21f6 100644 --- a/patchDescribe/calculateShapeDiameterFunctionRegion.m +++ b/patchDescribe/calculateShapeDiameterFunctionRegion.m @@ -62,7 +62,14 @@ % check every face in the mesh to see if the ray and face intersect using a one-sided ray-triangle intersection algorithm startMatrix = repmat([startPosition(1) startPosition(2) startPosition(3)], size(mesh.faces,1), 1); rayMatrix = repmat([ray(1) ray(2) ray(3)], size(mesh.faces,1), 1); - [intersect, dist] = TriangleRayIntersection(startMatrix, rayMatrix, mesh.vertices(mesh.faces(:,2),:), mesh.vertices(mesh.faces(:,1),:), mesh.vertices(mesh.faces(:,3),:), 'planeType', 'one sided'); + + [intersect, dist] = TriangleRayIntersection(startMatrix,... + rayMatrix, ... + mesh.vertices(mesh.faces(:,2),:), ... + mesh.vertices(mesh.faces(:,1),:), ... + mesh.vertices(mesh.faces(:,3),:), ... + 'planeType', ... + 'one sided'); intersectMask = dist.*(intersect & (dist > 2)); intersectMask(intersectMask == 0) = Inf; [minDist, faceIntersectIndex] = min(intersectMask); @@ -89,13 +96,11 @@ if numLoops > 100*raysPerCompare break end - end % % calculate the self visibility % selfVisibility = nanmean(selfVisArray); % if isempty(selfVisibility), selfVisibility = nan; end - % calculate the shape diameter function, remove values more than one standard deviation from the median %selfVisArray(isnan(sdfArray)) = []; sdfArray(isnan(sdfArray)) = []; @@ -108,7 +113,6 @@ sdf = mean(sdfArray); %selfVisibilityCentral = mean(selfVisArray); %if isempty(selfVisibilityCentral), selfVisibilityCentral = nan; end - % % debug code % climits = [0 Inf]; % cmap = jet(3); diff --git a/patchDescribe/patchDescriptionForMergeMeshMD.m b/patchDescribe/patchDescriptionForMergeMeshMD.m index 7eb7f6f..b8a1c5b 100644 --- a/patchDescribe/patchDescriptionForMergeMeshMD.m +++ b/patchDescribe/patchDescriptionForMergeMeshMD.m @@ -92,7 +92,7 @@ function patchDescriptionForMergeMeshMD(processOrMovieData, varargin) neighborsPath = fullfile(inFilePaths{1,c},'neighbors_%i_%i.mat'); surfaceSegPath = fullfile(inFilePaths{2,c},'surfaceSegment_%i_%i.mat'); - parfor t = 1:MD.nFrames_ % parfor + for t = 1:MD.nFrames_ % parfor % display progress disp([' image ' num2str(t) ' (channel ' num2str(c) ')']) diff --git a/patchDescribe/patchDescriptionMeshMD.m b/patchDescribe/patchDescriptionMeshMD.m index b61eafc..d75599e 100644 --- a/patchDescribe/patchDescriptionMeshMD.m +++ b/patchDescribe/patchDescriptionMeshMD.m @@ -119,7 +119,7 @@ function patchDescriptionMeshMD(processOrMovieData, varargin) % surfaceSegPath = fullfile(inFilePaths{2,c},'surfaceSegment_%i_%i.mat'); % end - parfor t = 1:MD.nFrames_ % parfor + for t = 1:MD.nFrames_ % parfor % display progress disp([' image ' num2str(t) ' (channel ' num2str(c) ')']) diff --git a/runKrasBlebSegment.m b/runKrasBlebSegment.m index 24d38e3..4aad99c 100644 --- a/runKrasBlebSegment.m +++ b/runKrasBlebSegment.m @@ -26,15 +26,15 @@ function runKrasBlebSegment() %% Set directories -imageDirectory = '/home2/mdrisc/Desktop/motif3DExampleDataFinal/testData/krasMV3'; % directory the image is in. The images to be analyzed should be the only thing in the directory. -motifModelDirectory = '/home2/mdrisc/Desktop/motif3DExampleDataFinal/svmModels'; % directory of SVM motif models -psfDirectory = '/home2/mdrisc/Desktop/motif3DExampleDataFinal/PSFs'; % directory of microscope PSFs -saveDirectory = '/home2/mdrisc/Desktop/Motif3DExampleDataFinal/Analysis/krasScript'; % directory for the analysis output +imageDirectory = 'C:\Users\bsterling\Desktop\ZebrafishScripted\Image'; % directory the image is in. The images to be analyzed should be the only thing in the directory. +motifModelDirectory = 'C:\Users\bsterling\Desktop\ZebrafishScripted\svmModels'; % directory of SVM motif models +%psfDirectory = '/home2/mdrisc/Desktop/motif3DExampleDataFinal/PSFs'; % directory of microscope PSFs +saveDirectory = 'C:\Users\bsterling\Desktop\ZebrafishScripted\Analysis\krasScript'; % directory for the analysis output %% Set movie parameters -pixelSizeXY = 160; % in nm -pixelSizeZ = 200; % in nm +pixelSizeXY = 1000; % in nm +pixelSizeZ = 500; % in nm timeInterval = 60; % in seconds @@ -57,19 +57,25 @@ function runKrasBlebSegment() %% Override Default Parameters -p.deconvolution.pathPSF = fullfile(psfDirectory,'meSpimPSF.mat'); +%p.deconvolution.pathPSF = fullfile(psfDirectory,'meSpimPSF.mat'); +% Turn off deconvolution: +p.control.deconvolution = 0; +p.mesh.useUndeconvolved = 1; p.mesh.imageGamma = 0.7; +%{ p.motifDetect.svmPath = {fullfile(motifModelDirectory,'SVMcertainBleb1'), ... % input more than one motif model to have the models vote fullfile(motifModelDirectory, 'SVMcertainBleb2'), ... fullfile(motifModelDirectory, 'SVMcertainBleb3')}; -%p.motifDetect.svmPath = fullfile(motifModelDirectory,'SVMtestKras'); -p.motifDetect.removePatchesSVMpath = fullfile(motifModelDirectory, 'SVMcertainRetractionFiber'); +%} +p.motifDetect.svmPath = motifModelDirectory;%fullfile(motifModelDirectory,'SVMZebrafish'); +p.motifDetect.removePatchesSVMpath = fullfile(motifModelDirectory, 'SVMZebrafishcertainRetractionFiber'); p.intensityBlebCompare.analyzeOnlyFirst = 1; p.intensityBlebCompare.calculateVonMises = 1; %% Analyze kras cells -imageList = 1:3; +%imageList = 1:3; +imageList = 1; parfor c = 1:length(imageList) % can be made a parfor loop if sufficient RAM is available. disp(['--------- Analysing KRAS Cell ' num2str(imageList(c))]) @@ -107,6 +113,7 @@ function runKrasBlebSegment() %% Plot analyses (or use runPlotIntensityBlebCompare) % (Note that more than three images would normally be analyzed to draw robust conclusions.) +%{ disp('--------- Plotting Analyses') p.savePath = fullfile(saveDirectory, 'analysisFigures'); p.mainDirectory = saveDirectory; @@ -115,3 +122,4 @@ function runKrasBlebSegment() p.cellsList{3} = ['Cell3']; p.analyzeDiffusion = 1; p.analyzeVonMises = 1; p.analyzeDistance = 1; plotIntensityBlebCompare(p); % see this function for more plots, many plots are commented out +%} \ No newline at end of file diff --git a/runKrasTrain.m b/runKrasTrain.m index a9ff80f..c2b8af3 100644 --- a/runKrasTrain.m +++ b/runKrasTrain.m @@ -35,15 +35,15 @@ function runKrasTrain() %% Set Directories -analysisDirectory = '/Users/meghan/Desktop/Morphology3DPackage/Analysis/krasScript'; -motifModelDirectory = '/Users/meghan/Desktop/Morphology3DPackage/svmModels'; +analysisDirectory = 'C:\Users\bsterling\Desktop\ZebrafishScripted\Analysis\'; +motifModelDirectory = 'C:\Users\bsterling\Desktop\ZebrafishScripted\Analysis\Cell1\Morphology\svmModels\'; %% Generate training data % set Kras cells to click on p.cellsList{1} = 'Cell1'; -p.cellsList{2} = 'Cell2'; -p.cellsList{3} = 'Cell3'; +%p.cellsList{2} = 'Cell2'; +%p.cellsList{3} = 'Cell3'; % set parameters for generating training data (see the documentation of clickOnBlebs.m for more information) p.mainDirectory = analysisDirectory; @@ -57,11 +57,13 @@ function runKrasTrain() % generate training data clickOnProtrusions(p) +%load('C:\Users\bsterling\Desktop\ZebrafishScripted\Analysis\Cell1\TrainingData\TestClicker\blebLocs.mat'); %% Train and validate an SVM motif classifier -p.MDsPathList{1} = p.cellsList{1}; p.clicksPathList{1} = fullfile(p.MDsPathList{1}, 'TrainingData', p.nameOfClicker); -p.MDsPathList{2} = p.cellsList{2}; p.clicksPathList{2} = fullfile(p.MDsPathList{2}, 'TrainingData', p.nameOfClicker); -p.MDsPathList{3} = p.cellsList{3}; p.clicksPathList{3} = fullfile(p.MDsPathList{3}, 'TrainingData', p.nameOfClicker); +p.MDsPathList{1} = p.cellsList{1}; +p.clicksPathList{1} = fullfile(p.MDsPathList{1}, 'TrainingData', p.nameOfClicker); +%p.MDsPathList{2} = p.cellsList{2}; p.clicksPathList{2} = fullfile(p.MDsPathList{2}, 'TrainingData', p.nameOfClicker); +%p.MDsPathList{3} = p.cellsList{3}; p.clicksPathList{3} = fullfile(p.MDsPathList{3}, 'TrainingData', p.nameOfClicker); p.saveNameModel = 'SVMtestKras'; p.saveNameCells = 'testKrasCells'; p.mainDirectory = analysisDirectory; @@ -72,4 +74,4 @@ function runKrasTrain() trainProtrusionsClassifier(p); % validate the classifier -validateBlebClassifier(p); \ No newline at end of file +%validateBlebClassifier(p); \ No newline at end of file diff --git a/utilities/medianFilterKD.m b/utilities/medianFilterKD.m index 94d1aef..7f128b8 100644 --- a/utilities/medianFilterKD.m +++ b/utilities/medianFilterKD.m @@ -25,6 +25,7 @@ % get the face center positions nFaces = size(surface.faces,1); faceCenters = zeros(nFaces,3); + for f = 1:nFaces faceCenters(f,:) = mean(surface.vertices(surface.faces(f,:),:),1); end