Consistent test data structure for processImage, newCurv

.github/workflows/ci.yml

@@ -1,14 +1,11 @@
 name: ci
 
 on:
-  # Allow manual runs from the Actions tab
   workflow_dispatch:
-  # Run on pushes to any branch (including feature branches)
   push:
-    branches: [ "**" ]
-  # And on PRs targeting main
+    branches: ["**"]
   pull_request:
-    branches: [ main ]
+    branches: [main]
 
 jobs:
   test:
@@ -19,19 +16,28 @@ jobs:
     steps:
       - uses: actions/checkout@v4
 
+      # Install system dependencies (like in your Docker container)
+      - name: Install system dependencies
+        run: sudo apt-get update && sudo apt-get install -y libgl1 libglib2.0-0 python3-venv
+
+      # Set up Python
       - uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
 
-      - name: Install dependencies
+      # Create virtual environment
+      - name: Create virtual environment
         run: |
-          python -m pip install --upgrade pip
-          # Install without curvelops to avoid FFTW/CurveLab build on CI
-          python -m pip install --break-system-packages -e .
-          python -m pip install --break-system-packages pytest ruff
+          python -m venv venv
+          source venv/bin/activate
+          python -m pip install --upgrade pip setuptools wheel
+          python -m pip install -e .
+          python -m pip install pytest ruff
 
+      # Run tests in the venv
       - name: Run tests
         env:
           QT_QPA_PLATFORM: "offscreen"
         run: |
-          pytest -q
+          source venv/bin/activate
+          pytest -rs

.gitignore

@@ -14,6 +14,7 @@ local.properties
 .classpath
 .settings/
 .loadpath
+*.idea
 
 # External tool builders
 .externalToolBuilders/
@@ -59,6 +60,7 @@ local.properties
 *.tli
 *.tlh
 *.tmp
+*.xlsx
 *.vspscc
 .builds
 *.dotCover
@@ -174,3 +176,9 @@ lastPATH_CTF.mat
 *.tif
 
 pycurvelets/opt/*.*
+
+# test tmp folder
+tests/test_images/CA_Out/
+
+# referenced MATLAB files
+*.m

pyproject.toml

@@ -9,16 +9,17 @@ description = "Python translation of CurveAlign; curvelet-based quantification a
 authors = [{ name = "Eliceiri Lab" }]
 requires-python = ">=3.11"
 dependencies = [
-  "pandas",
-  "numpy",
-  "matplotlib",
-  "scipy",
-  "opencv-python",
-  "scikit-image",
-  "openpyxl",
+  "matplotlib>=3.10.1",
   "napari",
+  "numpy>=2.2.6",
+  "opencv-python-headless>=4.10.0.84",
+  "openpyxl>=3.1.5",
+  "pandas>=2.2.3",
   "pyqt6",
   "qtpy",
+  "scikit-image>=0.25.0",
+  "scikit-learn>=1.7.1",
+  "scipy>=1.15.2",
 ]
 
 [dependency-groups]

src/curvealign_matlab/getTifBoundary.m

@@ -0,0 +1,241 @@
+function [resMat,resMatNames,numImPts] = getTifBoundary(coords,img,object,imgName,distThresh,fibKey,endLength,fibProcMeth,distMini)
+
+% getTifBoundary.m - This function takes the coordinates from the boundary file, associates them with curvelets, and produces relative angle measures.
+%
+% Inputs:
+%   coords - the locations of the endpoints of each line segment making up the boundary
+%   img - the image being measured
+%   object - a struct containing the center and angle of each measured curvelet, generated by the newCurv function
+%   distThresh - number of pixels from boundary we should evaluate curvelets
+%   boundaryImg - tif file with boundary outlines, must be a mask file
+%   fibKey - list indicating the beginning of each new fiber in the object struct, allows for fiber level processing
+%   distMini - minimum distrance to the boundary, to get rid of the fiber on or very close to the boundary
+%
+% Output:
+%   measAngs - all relative angle measurements, not filtered by distance
+%   measDist - all distances between the curvelets and the boundary points, not filtered
+%   inCurvsFlag - curvelets that are considered
+%   outCurvsFlag - curvelets that are not considered
+%   measBndry = points on the boundary that are associated with each curvelet
+%   inDist = distance between boundary and curvelet for each curvelet considered
+%   numImPts = number of points in the image that are less than distThresh from boundary
+%   insCt = number of curvelets inside an epithelial region
+%
+%
+% By Jeremy Bredfeldt, LOCI, Morgridge Institute for Research, 2013
+
+%Note: a "curv" could be a curvelet or a fiber segment, depending on if CT or FIRE is used
+
+imHeight = size(img,1);
+imWidth = size(img,2);
+sz = [imHeight,imWidth];
+
+% figure(600);
+% imshow(img);
+
+% figure(500);
+% hold on;
+% for k = 1:length(coords)
+%    boundary = coords{k};
+%    plot(boundary(:,2), boundary(:,1), 'y', 'LineWidth', 2);
+% end
+
+%collect all fiber points
+allCenterPoints = vertcat(object.center);
+%collect all boundary points
+% coords = vertcat(coords{2:end,1});
+coords = vertcat(coords{1:end,1});
+
+%collect all region points
+linIdx = sub2ind(sz, allCenterPoints(:,1), allCenterPoints(:,2));
+
+[idx_dist,dist] = knnsearch(coords,allCenterPoints); %closest point to a boundary
+reg_dist = img(linIdx);
+
+
+%YL: test the boundary association
+% figure(1002);clf;set(gcf,'pos',[200 300 imWidth imHeight ]);
+% plot(coords(:,2),coords(:,1),'k.'); axis ij
+% axis([1 imWidth 1 imHeight ]);hold on
+
+%[idx_reg,reg_dist] = knnsearch([reg_col,reg_row],allCenterPoints); %closest point to a filled in region
+
+%Make a list of points in the image (points scattered throughout the image)
+C = floor(imWidth/20); %use at least 20 per row in the image, this is done to speed this process up
+[I, J] = ind2sub(size(img),1:C:imHeight*imWidth);
+allImPoints = [I; J]';
+%Get list of image points that are a certain distance from the boundary
+[~,dist_im] = knnsearch(coords(1:3:end,:),allImPoints); %returns nearest dist to each point in image
+%threshold distance
+if isempty(distMini)  % keep all the fibers within the distance threshold
+  inIm = dist_im <= distThresh;
+else                  %  get rid of fibers on or very close to the boundary
+  inIm = (dist_im <= distThresh & dist_im > distMini);
+end
+%count number of points
+inPts = allImPoints(inIm);
+numImPts = length(inPts)*C;
+% numImPts = 0;
+
+
+%process all curvs, at this point
+curvsLen = length(object);
+nbDist = nan(1,curvsLen); %nearest boundary distance
+nrDist = nan(1,curvsLen); %nearest region distance
+nbAng = nan(1,curvsLen); %nearest boundary relative angle
+epDist = nan(1,curvsLen); %distance to extension point intersection
+epAng = nan(1,curvsLen); %relative angle of extension point intersection
+measBndry = nan(curvsLen,2);
+inCurvsFlag = ~logical(1:curvsLen);
+outCurvsFlag = ~logical(1:curvsLen);
+% %ROI properties
+% bwROI.coords = coords;
+% bwROI.imageWidth = imWidth;
+% bwROI.imageHeight = imHeight;
+if isempty(distMini)
+  for i = 1:curvsLen
+    %-- inside region?
+    nrDist(i) = reg_dist(i)==255|reg_dist(i)== 1; %YL: mask can be 1-0(matlab lab) or 255-0(ImageJ)
+    %-- distance to nearest epithelial boundary
+    nbDist(i) = dist(i);
+    %-- relative angle at nearest boundary point
+    if dist(i) <= distThresh
+      [nbAng(i), bPt] = GetRelAng([coords(:,2),coords(:,1)],idx_dist(i),object(i).angle,imHeight,imWidth,i);    % add i as an input argument for debug
+      % %%
+      % bwROI.index2object  = idx_dist(i);
+      % fiberobject.center = object(i).center;
+      % fiberobject.angle = object(i).angle;
+      % angleOption = 0; % caclulate all angles
+      % [relativeAngles, ROImeasurements] = getRelativeangles(bwROI,fiberobject,angleOption);
+      % fprintf(' original relative angle is %3.2f \n  re-calculated relative angle is %3.2f \n ',nbAng(i),relativeAngles.angle2boundaryEdge);
+      % fprintf(' relative angle to ROI orientation is %3.2f \n',relativeAngles.angle2boundaryCenter);
+      % fprintf('relative angle to ROI-fiber centers line is %3.2f \n',relativeAngles.angle2centersLine);
+      %%
+    else
+      nbAng(i) = nan; % if out of the region
+      bPt = nan(1,2);  % if out of the region
+    end
+    %-- extension point features
+    [lineCurv orthoCurv] = getPointsOnLine(object(i),imWidth,imHeight,distThresh);
+    [intLine, iLa, iLb] = intersect([lineCurv(:,2) lineCurv(:,1)],coords,'rows');
+    if (~isempty(intLine))
+      %get the closest distance from the curvelet center to the
+      %intersection (get rid of the farther one(s))
+      [idxLineDist, lineDist] = knnsearch(intLine,object(i).center);
+      boundaryPtIdx = iLb(idxLineDist);
+      %%tentatively turn the extension feature off
+      %         %-- extension point distance
+      %         epDist(i) = lineDist;
+      %         %-- extension point angle
+      %         [epAng(i) bPt1] = GetRelAng([coords(:,2),coords(:,1)],boundaryPtIdx,object(i).angle,imHeight,imWidth,i);
+    else
+      epDist(i) = nan;% no intersection exists
+      epAng(i) = nan; % no angle exists
+      bPt1 = nan(1,2);      % if no intersection set boundary to be [NaN NaN]
+    end
+    measBndry(i,:) = bPt;  % nearest boundary
+  end %of for loop
+else  % get rid of fibers on or within the minimum distance to the boundary
+  for i = 1:curvsLen
+    %-- inside region?
+    nrDist(i) = reg_dist(i)==255|reg_dist(i)== 1; %YL: mask can be 1-0(matlab lab) or 255-0(ImageJ)
+    %-- distance to nearest epithelial boundary
+    nbDist(i) = dist(i);
+    %-- relative angle at nearest boundary point
+    if (dist(i) <= distThresh & dist(i) > distMini)
+      [nbAng(i), bPt] = GetRelAng([coords(:,2),coords(:,1)],idx_dist(i),object(i).angle,imHeight,imWidth,i);    % add i as an input argument for debug
+    else
+      nbAng(i) = nan; % if out of the region
+      bPt = nan(1,2);  % if out of the region
+    end
+    %-- extension point features
+    [lineCurv orthoCurv] = getPointsOnLine(object(i),imWidth,imHeight,distThresh);
+    [intLine, iLa, iLb] = intersect([lineCurv(:,2) lineCurv(:,1)],coords,'rows');
+    if (~isempty(intLine))
+      %get the closest distance from the curvelet center to the
+      %intersection (get rid of the farther one(s))
+      [idxLineDist, lineDist] = knnsearch(intLine,object(i).center);
+      boundaryPtIdx = iLb(idxLineDist);
+      %%tentatively turn the extension feature off
+      %         %-- extension point distance
+      %         epDist(i) = lineDist;
+      %         %-- extension point angle
+      %         [epAng(i) bPt1] = GetRelAng([coords(:,2),coords(:,1)],boundaryPtIdx,object(i).angle,imHeight,imWidth,i);
+    else
+      epDist(i) = nan;% no intersection exists
+      epAng(i) = nan; % no angle exists
+      bPt1 = nan(1,2);      % if no intersection set boundary to be [NaN NaN]
+    end
+    measBndry(i,:) = bPt;  % nearest boundary
+  end %of for loopend
+end
+
+
+
+resMat = [nbDist', ... %nearest dist to a boundary
+  nrDist', ... %flag, 0 for outside boundary, 1 for inside boundary
+  nbAng',  ... %nearest relative boundary angle
+  epDist', ... %extension point distance
+  epAng',  ... %extension point relative boundary angle
+  measBndry];  %list of boundary points associated with fibers
+resMatNames = {
+  'nearestBoundDist', ...
+  'nearestRegionDist', ...
+  'nearestBoundAng', ...
+  'extensionPointDist', ...
+  'extensionPointAng', ...
+  'bndryPtRow', ...
+  'bndryPtCol'};
+
+end %of main function
+
+function [relAng, boundaryPt] = GetRelAng(coords,idx,fibAng,imHeight,imWidth,fnum)
+boundaryAngle = FindOutlineSlope([coords(:,2),coords(:,1)],idx);
+boundaryPt = coords(idx,:);
+if (boundaryPt(1) == 1 || boundaryPt(2) == 1 || boundaryPt(1) == imHeight || boundaryPt(2) == imWidth)
+  %don't count fiber if boundary point is along edge of image
+  tempAng = 0;
+else
+  %--compute relative angle here--
+  %There is a 90 degree phase shift in fibAng and boundaryAngle due to image orientation issues in Matlab.
+  % -therefore no need to invert (ie. 1-X) circ_r here.
+  tempAng = circ_r([fibAng*2*pi/180; boundaryAngle*2*pi/180]);
+  tempAng = 180*asin(tempAng)/pi;
+  %         %YL debug the NaN angle
+  %         if isnan(tempAng)
+  %             figure(1002),plot(coords(idx,1),coords(idx,2),'ro','MarkerSize',10)
+  %             text(coords(idx,1),coords(idx,2),sprintf('%d',fnum));
+  %             disp(sprintf('fiber %d relative angle is Nan, fibAng = %f, boundaryAngle = %f, idx_dist = %d',fnum,fibAng,boundaryAngle,idx))
+  % %            pause(3)
+  %         end
+end
+relAng = tempAng;
+end
+
+function [lineCurv orthoCurv] = getPointsOnLine(object,imWidth,imHeight,boxSz)
+center = object.center;
+angle = object.angle;
+slope = -tand(angle);
+%orthoSlope = -tand(angle + 90); %changed from tand(obj.ang) to -tand(obj.ang + 90) 10/12 JB
+intercept = center(1) - (slope)*center(2);
+%orthoIntercept = center(1) - (orthoSlope)*center(2);
+
+%[p1 p2] = getBoxInt(slope, intercept, imWidth, imHeight, center, boxSz);
+if isinf(slope)
+  dist_y = 0;
+  dist_x = boxSz;
+else
+  dist_y = boxSz/sqrt(1+slope*slope);
+  dist_x = dist_y*slope;
+end
+p1 = [center(2) - dist_y, center(1) - dist_x];
+p2 = [center(2) + dist_y, center(1) + dist_x];
+[lineCurv ~] = GetSegPixels(p1,p2);
+
+%Not using the orthogonal slope for anything now
+%[p1 p2] = getIntImgEdge(orthoSlope, orthoIntercept, imWidth, imHeight, center);
+%[orthoCurv, ~] = GetSegPixels(p1,p2);
+orthoCurv = [];
+
+end
+

src/pycurvelets/__init__.py

@@ -1,10 +1,11 @@
-# Optional import: allow the package to be imported without native curvelet deps.
 try:
     from .new_curv import new_curv  # type: ignore
+
     HAS_CURVELETS = True
 except Exception:
     HAS_CURVELETS = False
 
-__all__ = ["new_curv", "HAS_CURVELETS"]
-
-
+__all__ = [
+    "HAS_CURVELETS",
+    "new_curv",
+]