From 32e93dc81e7a4b7a9a4f110a07aa38170a0a9be5 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Mon, 1 Aug 2022 15:41:09 +0100 Subject: [PATCH 01/18] Fix README instructions --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 5201026..c324b3e 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ To know more about the segmentation routine check out our wiki page: [epitools.g load('ProjIm.mat') %% crop image for testing (right click -> "Crop Image") -[crop, rect] = imcrop(ProjIm,[]); +[crop, rect] = imcrop(ProjIm); close(); %% compute segmentation and output segmentation feedback From fa80b042116f9c8bbaa661dfb54e83a8bcfe6c30 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Mon, 1 Aug 2022 15:41:20 +0100 Subject: [PATCH 02/18] Ignore direnv --- .gitignore | 1 + 1 file changed, 1 insertion(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7a6353d --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +.envrc From 4eaa1e949ed55554bddcae380564ee5103f80ced Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Mon, 1 Aug 2022 15:53:36 +0100 Subject: [PATCH 03/18] Add `python` `pre-commit` and apply --- .github/workflows/checks.yml | 43 +++ .pre-commit-config.yaml | 41 ++ ProjIm_example.m | 2 +- SegmentIm.m | 194 +++++----- calculateCellPositions.m | 44 +-- findcellsfromregiongrowing.cc | 540 +++++++++++++-------------- findcellsfromregiongrowing.mexmaci64 | Bin growcellsfromseeds3.cc | 424 ++++++++++----------- growcellsfromseeds3.mexmaci64 | Bin utilah.h | 17 +- 10 files changed, 692 insertions(+), 613 deletions(-) create mode 100644 .github/workflows/checks.yml create mode 100644 .pre-commit-config.yaml mode change 100755 => 100644 SegmentIm.m mode change 100755 => 100644 calculateCellPositions.m mode change 100755 => 100644 findcellsfromregiongrowing.cc mode change 100755 => 100644 findcellsfromregiongrowing.mexmaci64 mode change 100755 => 100644 growcellsfromseeds3.cc mode change 100755 => 100644 growcellsfromseeds3.mexmaci64 diff --git a/.github/workflows/checks.yml b/.github/workflows/checks.yml new file mode 100644 index 0000000..edd02bf --- /dev/null +++ b/.github/workflows/checks.yml @@ -0,0 +1,43 @@ +name: checks + +on: + push: + pull_request: + +jobs: + build: + runs-on: ubuntu-latest + steps: + - name: checkout source + uses: actions/checkout@v3 + + - name: set up python + uses: actions/setup-python@v3 + with: + python-version: "3.10" + + - name: set PY + run: echo "PY=$(python -VV | sha256sum | cut -d' ' -f1)" >> $GITHUB_ENV + + - name: cache stuff + uses: actions/cache@v3 + with: + path: | + ${{ env.pythonLocation }} + ~/.cache/pre-commit + key: | + pre-commit-${{ env.PY }}-${{ hashFiles('.pre-commit-config.yaml') }} + + - name: install dependencies + run: pip install pre-commit + + - name: Install pre-commit hooks + run: pre-commit install + + # This will run on all files in the repo not just those that have been + # committed. Once formatting has been applied once globally, from then on + # the files being changed by pre-commit should be just those that are + # being committed - provided that people are using the pre-commit hook to + # format their code. + - name: run pre-commit + run: pre-commit run --all-files --color always diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..dba1c99 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,41 @@ +repos: + - repo: https://github.com/pycqa/isort + rev: 5.10.1 + hooks: + - id: isort + types: [text] + types_or: [python, cython] + args: ["--profile", "black"] + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.1.0 + hooks: + - id: check-case-conflict + - id: check-docstring-first + - id: check-executables-have-shebangs + - id: check-merge-conflict + - id: check-toml + - id: end-of-file-fixer + - id: mixed-line-ending + args: [--fix=lf] + - id: requirements-txt-fixer + - id: trailing-whitespace + - repo: https://github.com/Lucas-C/pre-commit-hooks + rev: v1.3.0 + hooks: + - id: forbid-tabs + - repo: https://gitlab.com/pycqa/flake8 + rev: 3.9.2 + hooks: + - id: flake8 + - repo: https://github.com/psf/black + rev: 22.3.0 + hooks: + - id: black + - repo: https://github.com/pre-commit/mirrors-mypy + rev: v0.942 + hooks: + - id: mypy + - repo: https://github.com/pre-commit/mirrors-prettier + rev: v2.6.1 + hooks: + - id: prettier diff --git a/ProjIm_example.m b/ProjIm_example.m index 8937637..5c24cac 100644 --- a/ProjIm_example.m +++ b/ProjIm_example.m @@ -15,4 +15,4 @@ boundT=0.1; % mimimal seed/boundary ratio verbose=1; % segmentation feedback plot -[segmentation, seeds, labels] = SegmentIm(crop,s1,minA,minI,mergeT,s2,maxA,boundT,verbose); \ No newline at end of file +[segmentation, seeds, labels] = SegmentIm(crop,s1,minA,minI,mergeT,s2,maxA,boundT,verbose); diff --git a/SegmentIm.m b/SegmentIm.m old mode 100755 new mode 100644 index 485770f..1b87f41 --- a/SegmentIm.m +++ b/SegmentIm.m @@ -7,7 +7,7 @@ IBoundMax_pcnt,... show) % SegmentIm segments a single frame extracting the cell outlines -% IN: +% IN: % Im - [uint8] Image [increasing intesity for membrane] % sigma1 - [0+] size px of gaussian for smoothing image % mincellsize - [0+] size px of smallest cell expected @@ -18,7 +18,7 @@ % IBoundMax_pcnt - [0-1] minimum ratio for seed and membrane intensity % show - [0/1] show feedback for segmentation % -% OUT: +% OUT: % CellSeeds - [uint8] Rescaled Image to fit the range [0,252] % 253/254/255 are used for seed information % CellLables - [uint16] bitmap of cells colored with 16bit id @@ -40,28 +40,28 @@ CellLabels = zeros(ImSize,'uint16'); %todo: check casting: why using 16 bit for labels? %structuring element, SE, used for morphological operations -se = strel('disk',2); +se = strel('disk',2); %% Operations -% [0] Create starting seeds +% [0] Create starting seeds % Find the initial cell seeds (parameters: sigma1, threshold) DoInitialSeeding(); if show figure('Name','Intermediate steps'); end -if show +if show ax1 = subplot(4,6,[1 2 7 8]); %p = get(h, 'pos'); %p(3) = p(3) - 0.05; %set(h, 'pos', p); imshow(CellSeeds(:,:),[]); - title(['[Seeds] smooth1=' num2str(sigma1) '; min.int=' num2str(threshold)]); + title(['[Seeds] smooth1=' num2str(sigma1) '; min.int=' num2str(threshold)]); end - + % Remove initial cell regions which touch & whose boundary is insufficient % (parameters: params.MergeCriteria) MergeSeedsFromLabels() -if show +if show ax2 = subplot(4,6,[5 6 11 12]); imshow(CellSeeds(:,:),[]); title(['[Merging] min.ratio.thr=' num2str(MergeCriteria)]); @@ -69,11 +69,11 @@ % [1] Growing cells from seeds (parameter: sigma3) TODO: add paramters in Name description! GrowCellsInFrame() -if show - CreateColorBoundaries(); - ax3 = subplot(4,6,[13 14 19 20]); - imshow(ColIm,[]); - title(['[Boundaries] smooth2=' num2str(sigma3)]); +if show + CreateColorBoundaries(); + ax3 = subplot(4,6,[13 14 19 20]); + imshow(ColIm,[]); + title(['[Boundaries] smooth2=' num2str(sigma3)]); end % [2] Eliminate labels from seeds which have poor boundary intensity @@ -87,18 +87,18 @@ CreateColorBoundaries() if show ax4 = subplot(4,6,[17 18 23 24]); - imshow(ColIm,[]); + imshow(ColIm,[]); title(['[Final] boundary.thr=' num2str(IBoundMax_pcnt)]); - + linkaxes([ax1,ax2,ax3,ax4],'xy') - + % plot histogram in subplot(3,4,[2 3 4]) region_property = regionprops(CellLabels,'Area'); region_areas = cat(1,region_property.Area); subplot(4,6,[3 4]); hist(region_areas, 100); xlabel('Area of cells'); - + % add some general statistics subplot(4,6,[21 22]); description = {... @@ -111,7 +111,7 @@ load('zerogreen','mycmap') set(gcf, 'Colormap', mycmap) - + end @@ -120,7 +120,7 @@ function CreateColorBoundaries() % create nice pic with colors for cells - + % cellBoundaries = zeros(ImSize,'uint8'); % ColIm = zeros([ImSize(1) ImSize(2) 3],'double'); % fs=fspecial('laplacian',0.9); @@ -136,16 +136,16 @@ function CreateColorBoundaries() % ColIm(:,:,2) = .2*double(cellBoundaries(:,:)) + ColIm(:,:,2).*(1-double(cellBoundaries(:,:))); % ColIm(:,:,3) = .2*double(cellBoundaries(:,:)) + ColIm(:,:,3).*(1-double(cellBoundaries(:,:))); % ColIm = cast(ColIm*255, 'uint8'); %todo: typecasting -% +% %given that every cell has a different label - %we can compute the boundaries by computing + %we can compute the boundaries by computing %where the gradient changes cell_lables = double(CellLabels(:,:)); [gx,gy] = gradient(cell_lables); cell_outlines = (cell_lables > 0) & ((gx.^2+gy.^2)>0); - + cell_seeds = CellSeeds(:,:) > 253; - + ColIm = Im; ColIm(cell_outlines) = 0; ColIm(cell_seeds) = 255; @@ -153,50 +153,50 @@ function CreateColorBoundaries() end function DoInitialSeeding() - + % Create gaussian filter if sigma1 > 0.01 f1=fspecial( 'gaussian', [ImSize(1) ImSize(2)], sigma1); - + % Gaussian smoothing for the segmentation of individual cells SmoothedIm = real(fftshift(ifft2(fft2(Im(:,:)).*fft2(f1)))); else SmoothedIm = double(Im(:,:)); end %if show figure; imshow(SmoothedIm(:,:,1),[]); input('press to continue','s'); end - + SmoothedIm = SmoothedIm/max(max(SmoothedIm))*252.; - + % Use external c-code to find initial seeds InitialLabelling = findcellsfromregiongrowing(SmoothedIm , mincellsize, threshold); %if show figure; imshow(InitialLabelling(:,:),[]); input('press to continue','s'); end - + InitialLabelling(InitialLabelling==1) = 0; % set unallocated pixels to 0 - + % Generate CellLabels from InitalLabelling CellLabels(:,:) = uint16(InitialLabelling); - + % eliminate very large areas DelabelVeryLargeAreas(); % DelabelFlatBackground() - + % Use true centre of cells as labels centroids = round(calculateCellPositions(SmoothedIm,CellLabels(:,:), false)); centroids = centroids(~isnan(centroids(:,1)),:); for n=1:length(centroids); SmoothedIm(centroids(n,2),centroids(n,1))=255; end - - % CellSeeds contains the position of the true cell center. + + % CellSeeds contains the position of the true cell center. CellSeeds(:,:) = uint8(SmoothedIm); - + end % % Initial specification was encoding background pixels as zero values in cell images. % % DelabelFlatBackground() removes such background pixels from the cell label image, -% % i.e. it is applying a mask. -% function DelabelFlatBackground() +% % i.e. it is applying a mask. +% function DelabelFlatBackground() % L = CellLabels; % D = Im(:,:); % L(D==0) = 0; @@ -206,35 +206,35 @@ function DoInitialSeeding() function GrowCellsInFrame() bw=double(CellSeeds(:,:) > 252); % find labels - + if sigma3 > 0.01 f1=fspecial( 'gaussian', [ImSize(1) ImSize(2)], sigma3); SmoothedIm = real(fftshift(ifft2(fft2(Im(:,:)).*fft2(f1)))); else SmoothedIm = double(Im(:,:)); end - + ImWithSeeds = double(SmoothedIm).*(1-bw)+255*bw; % mark labels on image CellLabels = uint16(growcellsfromseeds3(ImWithSeeds,253)); - + end function UnlabelPoorSeedsInFrame() - + L = CellLabels; - + if sigma3 > 0.01 f1=fspecial( 'gaussian', [ImSize(1) ImSize(2)], sigma3); smoothedIm = real(fftshift(ifft2(fft2(Im(:,:)).*fft2(f1)))); else smoothedIm = double(Im(:,:)); end - - labelList = unique(L); %i.e. every cell is marked by one unique integer label + + labelList = unique(L); %i.e. every cell is marked by one unique integer label labelList = labelList(labelList~=0); IBounds = zeros(length(labelList),1); decisions = [ 0 0 0 0 0 ]; - + for c = 1:length(labelList) mask = L==labelList(c); [cpy cpx]=find(mask > 0); @@ -254,16 +254,16 @@ function UnlabelPoorSeedsInFrame() IEr = reducedIm(erodedMask>0); IBound = mean(boundaryIntensities); IBounds(c) = IBound; - + % cell seed information is retrieved as comparison F2 = CellSeeds; F2(~mask) = 0; [cpy cpx]=find(F2 > 252); ICentre = smoothedIm(cpy , cpx); - + IBoundMax = 255 * IBoundMax_pcnt; - - + + %Figure out which conditions make the label invalid %1. IBoundMax, gives the Lower bound to the mean intensity % 1.b condition upon that the cell seed has less than 20% intensity difference to the mean @@ -273,16 +273,16 @@ function UnlabelPoorSeedsInFrame() second_condition = (IBound < IBoundMax *25./30.); %3. If the minimum retrieved in the boundary mask is 0 (dangerous!) third_condition = (min(boundaryIntensities)==0); - %4. If the amount of low intensity signal (i.e. < 20) is more than 10% + %4. If the amount of low intensity signal (i.e. < 20) is more than 10% fourth_condition = (sum(H 0.1); if first_condition... || second_condition ... || third_condition... || fourth_condition - + %The label is cancelled (inverted mask multiplication.) CellLabels = CellLabels.*uint16(mask==0); - + % record the removal decisions if first_condition decisions(1) = decisions(1) + 1; @@ -295,23 +295,23 @@ function UnlabelPoorSeedsInFrame() else % should not happen decisions(5) = decisions(5) + 1; - end - + end + end - + end %The following debug figure shows the distribution of mean cell boundary intensity %if the threshold parameter IBoundMax is too high, valid cells might be delabeled if show subplot(4,6,[15 16]); - hist(IBounds/255,100); + hist(IBounds/255,100); xlabel(['mean cell boundary strength -[' num2str(decisions) ']']); % title(['[Cell boundary intensity] lower bound = ' num2str(IBoundMax_pcnt)]); end end function DelabelVeryLargeAreas() - + % remove cells which are bigger than LargeCellSizeThres L = CellLabels; dimInitL = length(L); @@ -320,7 +320,7 @@ function DelabelVeryLargeAreas() ls = unique(L); for i = 1:size(ls); l = ls(i); - if l == 0 + if l == 0 continue; end A = As(l); @@ -341,50 +341,50 @@ function MergeSeedsFromLabels() else smoothedIm = double(Im(:,:)); end - + labelList = unique(CellLabels); labelList = labelList(labelList~=0); c = 1; - + merge_intensity_distro = []; merge_decisions = 0; - + % loop over labels - while 1==1 + while 1==1 labelMask = CellLabels==labelList(c); label = labelList(c); - + [cpy cpx]=find(labelMask > 0); - + % find region of that label minx = min(cpx); maxx = max(cpx); miny = min(cpy); maxy = max(cpy); minx = max(minx-5,1); miny = max(miny-5,1); maxx = min(maxx+5,ImSize(2)); maxy = min(maxy+5,ImSize(1)); - + % reduce data to that region reducedLabelMask = labelMask(miny:maxy, minx:maxx); reducedIm = smoothedIm(miny:maxy, minx:maxx); reducedLabels = CellLabels(miny:maxy, minx:maxx); - + % now find boundaries ... dilatedMask = imdilate(reducedLabelMask, se); erodedMask = imerode(reducedLabelMask, se); borderMask = dilatedMask - erodedMask; borderIntensities = reducedIm(borderMask>0); centralIntensity = reducedIm(erodedMask>0); - + F2 = CellSeeds; F2(~labelMask) = 0; [cpy cpx]=find(F2 > 253); ICentre = smoothedIm(cpy , cpx); - + background_std = std(double(centralIntensity)); - + % get labels of surrounding cells (neighbours) neighbourLabels = unique(reducedLabels( dilatedMask > 0 )); neighbourLabels = neighbourLabels(neighbourLabels~=label); - + low_intensity_ratios = []; for i = 1:size(neighbourLabels) neighbLabel = neighbourLabels(i); @@ -392,65 +392,65 @@ function MergeSeedsFromLabels() neighbor_border(reducedLabels~=neighbLabel)=0; % slice of neighbour around cell cell_border = imdilate(neighbor_border,se); cell_border(reducedLabels~=label) = 0; % slice of cell closest to neighbour - + joint_border = ... (cell_border + neighbor_border) > 0; % combination of both creating boundary region border_intensities = reducedIm; border_intensities(~joint_border) = 0; % intensities at boundary - - % average number of points in boundary where intensity is + + % average number of points in boundary where intensity is % of low quality (dodgy) low_intensity_threshold = ICentre + (background_std/2.); low_intensity_pixels = ... border_intensities(joint_border) < low_intensity_threshold; - + low_intensity_ratio = ... sum(low_intensity_pixels)/size(border_intensities(joint_border),1); - + low_intensity_ratios = [low_intensity_ratios low_intensity_ratio]; end - - + + %Find out which is border with the lowest intensity ratio [worst_intensity_ratio,worst_neighbor_index] = max(low_intensity_ratios); neighbLabel = neighbourLabels(worst_neighbor_index); - - + + % if the label value is of poor quality, then recursively check % the merge criteria in order to add it as a potential label in - % the label set. - + % the label set. + merge_intensity_distro = [merge_intensity_distro; worst_intensity_ratio]; - + if ... worst_intensity_ratio > MergeCriteria && ... label~=0 && ... - neighbLabel~=0 - + neighbLabel~=0 + MergeLabels(label,neighbLabel); labelList = unique(CellLabels); labelList = labelList(labelList~=0); - c = c-1; % reanalyze the same cell for more + c = c-1; % reanalyze the same cell for more % possible mergings merge_decisions = merge_decisions + 1; - + end - + c = c+1; - + % Condition to break the while cycle -> as soon as all the % labels are processed, then exit if c > length(labelList); break; end end - + if show % plot the distro subplot(4,6,[9 10]) hist(merge_intensity_distro, 100); xlabel(['ratio of low intensity boundary px (merges=' num2str(merge_decisions) ')']); - % ylabel('percentage of cells'); + % ylabel('percentage of cells'); end - + end function MergeLabels(l1,l2) @@ -461,17 +461,17 @@ function MergeLabels(l1,l2) Il1 = Il; Il1(~m1) = 0; Il2 = Il; Il2(~m2) = 0; [cpy1 cpx1]=find( Il1 > 253); - [cpy2 cpx2]=find( Il2 > 253); - cpx = round((cpx1+cpx2)/2); + [cpy2 cpx2]=find( Il2 > 253); + cpx = round((cpx1+cpx2)/2); cpy = round((cpy1+cpy2)/2); - + CellSeeds(cpy1,cpx1) = 20; %background level - CellSeeds(cpy2,cpx2) = 20; + CellSeeds(cpy2,cpx2) = 20; if CellLabels(cpy,cpx)==l1 || CellLabels(cpy,cpx)==l2 CellSeeds(cpy,cpx) = 255; else % center is not actually under any of the previous labels ... - if sum(m1(:)) > sum(m2(:)) + if sum(m1(:)) > sum(m2(:)) CellSeeds(cpy1,cpx1) = 255; else CellSeeds(cpy2,cpx2) = 255; @@ -494,5 +494,5 @@ function NeutralisePtsNotUnderLabelInFrame() % end CellSeeds(:,:) = F; end - -end \ No newline at end of file + +end diff --git a/calculateCellPositions.m b/calculateCellPositions.m old mode 100755 new mode 100644 index 8265c0b..c3c4ab9 --- a/calculateCellPositions.m +++ b/calculateCellPositions.m @@ -1,5 +1,5 @@ function pos=calculateCellPositions(I1, labelledcells, type) -% This function calculates the min intensity position of each labelled cell +% This function calculates the min intensity position of each labelled cell % or the centroid position of each labelled region no_cells=max(labelledcells(:)); @@ -9,63 +9,59 @@ I2=255-I1; for n=1:no_cells, - + [sy , sx]=find(labelledcells==n); - + if (type == 1) % looking for the lowest intensity - + PlaceAtLowestInt(n) - + else % calculating the centroid from intensities - + PlaceAtCentroid(n) - + if ~isnan(pos(n,2)) && labelledcells(pos(n,2),pos(n,1)) ~= n % every so often the centroid is actually not in the label! PlaceAtLowestInt(n) end - + end - + end - - + + function PlaceAtLowestInt(n) val=255; for m=1:length(sy), - + if (I1(sy(m),sx(m)) < val) - + val=I1(sy(m),sx(m)); pos(n,2)=sy(m); pos(n,1)=sx(m); - + end - + end end - + function PlaceAtCentroid(n) sumX=0; sumY=0; sumI=0; - + for m=1:length(sy), - + sumX=sumX+sx(m)*I2(sy(m),sx(m)); sumY=sumY+sy(m)*I2(sy(m),sx(m)); sumI=sumI+I2(sy(m),sx(m)); end - + pos(n,2)=round(sumY/sumI); pos(n,1)=round(sumX/sumI); end - - - - -end +end diff --git a/findcellsfromregiongrowing.cc b/findcellsfromregiongrowing.cc old mode 100755 new mode 100644 index 7401dd1..1ffb3d4 --- a/findcellsfromregiongrowing.cc +++ b/findcellsfromregiongrowing.cc @@ -1,270 +1,270 @@ -/* - - 1. finding non allocated regions based on a intensity threshold - 2. finding new cells based on clustered unallocated regions - 3. grow cells with a maximum number of iterations and below the intensity threshold - 4. merge regions on the basis of cell size and cell intenity - - */ - - -#include "mex.h" -#include "utilah.h" - - -//pixel list -typedef struct -{ - int x; - int y; - int l; // label - -}pix; - - -// growing all existing cells (or rather, assigning newly available pixels to present regions) -void growingCells(matx* A, matx* B, int* cellsizes) //imageMatrix A, label B, cellsizes = how many pixels per cell -{ - int mrows=B->rows; - int ncols=B->cols; - int n,x,y; - int no_pixels=0; - - pix* pixlist = (pix*)malloc(ncols*mrows*sizeof(pix)); - - - for (y=0; y < A->rows; y++) - { - for (x=0; x < A->cols; x++) - { - if (B->rptr[x][y] == 1.0) - { - //Found an available ( try to assign it to a cell region - - //check neighborhood for possible assingnment - if(x < ncols-1 && B->rptr[x+1][y] > 1) - { - pixlist[no_pixels].x=x; - pixlist[no_pixels].y=y ; - pixlist[no_pixels].l=B->rptr[x+1][y]; - no_pixels++; - continue; - } - if(x > 0 && B->rptr[x-1][y] > 1) - { pixlist[no_pixels].x=x; pixlist[no_pixels].y=y ;pixlist[no_pixels].l=B->rptr[x-1][y]; no_pixels++; continue; } - if(y > 0 && B->rptr[x][y-1] > 1) - { pixlist[no_pixels].x=x; pixlist[no_pixels].y=y ;pixlist[no_pixels].l=B->rptr[x][y-1]; no_pixels++; continue; } - if(y < mrows-1 && B->rptr[x][y+1] > 1) - { pixlist[no_pixels].x=x; pixlist[no_pixels].y=y ;pixlist[no_pixels].l=B->rptr[x][y+1]; no_pixels++; continue; } - } - } - } - - - //Post Update of Label Matrix B (i.e. Prevent cell from growing indefinitely) - for (n=0; n < no_pixels; n++) - { - B->rptr[pixlist[n].x][pixlist[n].y]=pixlist[n].l; - cellsizes[ pixlist[n].l]++; - } - - free( pixlist); -} - - - -// search for cluster of unallocated but AVAILABLE pixels -int findingNewCells(matx* A, matx* B, int no_cells) -{ - - int x, y, x2, y2; - int no_pixels=0; - - //whole image with border - for (y2=3; y2 < A->rows-3; y2++) - for (x2=3; x2 < A->cols-3; x2++) - { - no_pixels=0; - - //going through rectangular region 5x5 - for (y=-2; y <= 2; y++) - for (x=-2; x <= 2; x++) - { - if (B->rptr[x2+x][y2+y] == 1.0) //i.e. only AVAILABLE pixels count! - no_pixels++; - - } - - int min_cluster_size = 10; - - // found a new cell when the area consists of more than 10 pixels - if (no_pixels > min_cluster_size) - { - no_cells++; // cell label == cell number - for (y=-2; y <= 2; y++) - for (x=-2; x <= 2; x++) - { - // mark the new cell on the label matrix - //if (B->rptr[x2+x][y2+y] == 1.0) TODO, should be corrected with here! - B->rptr[x2+x][y2+y] = no_cells; - } - } - } - - return no_cells; -} - - -// merge cell1 and cell2 which becomes cell2 -void mergeTwoCells( matx* B, int* cellsizes, int cell1, int cell2) -{ - - int x,y; - int counter=0; - int label1; - - for (y=0; y < B->rows; y++) - for (x=0; x < B->cols; x++) - { - label1=B->rptr[x][y]; - if (label1 == cell1) { B->rptr[x][y]=cell2; counter++; } - - } - - cellsizes[cell2]=cellsizes[cell2]+counter; - cellsizes[cell1]=0; -} - - -//merge touching regions if conditions are met -void mergeCells(matx* B, int* cellsizes, double mincellsize, bool belowthreshold) -{ - int x,y; - int label1,label2; - - for (y=1; y < B->rows-1; y++) - { - for (x=1; x < B->cols-1; x++) - { - // check whether two cells are touching - label1=B->rptr[x][y]; - label2=B->rptr[x+1][y]; - - if (label1 > 1 && label2 > 1 && label1 != label2) - { - if ((cellsizes[label1] < mincellsize || cellsizes[label2] < mincellsize) || belowthreshold) - { - // merging two cells - mergeTwoCells( B, cellsizes, label2, label1); - } - } - - label2=B->rptr[x][y+1]; - if (label1 > 1 && label2 > 1 && label1 != label2) - { - if ((cellsizes[label1] < mincellsize && cellsizes[label2] < mincellsize) || belowthreshold) - { - // merging two cells - mergeTwoCells( B, cellsizes, label2, label1); - } - } - } - } -} - - - - - -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) -{ - - matx *A, *B; - matx *pixel_list; - int mrows, ncols; - int n,x,y,count, no_pixels, no_cells, i; - double mincellsize,level, maxval, threshold; - double *img; - int* cellsizes; - - - // --- testing input and output values --- - if (nrhs != 3) { mexErrMsgTxt(" usage: B=findscellsfromregiongrowing(CellImage(A), min cell size, threshold)"); } - if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) - mexErrMsgTxt("Image input must be a double real matrix"); - if (!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1])) - mexErrMsgTxt("Image input must be a double real matrix"); - - // -- initialising arrays and stuff --- //matlab to c++ conversion - mrows = mxGetM(prhs[0]); - ncols = mxGetN(prhs[0]); - A=mx_from_vector( mxGetPr(prhs[0]), mrows, ncols); - - img=mxGetPr(prhs[0]); - - mincellsize = mxGetScalar( prhs[1]); - threshold = mxGetScalar( prhs[2]); - - plhs[0] = mxCreateDoubleMatrix( mrows,ncols, mxREAL); - B=mx_from_vector( mxGetPr(plhs[0]), mrows, ncols); - - // estimate number of pixels per image for the pixelist - // find the max value - cellsizes=(int*)malloc(mrows*ncols*sizeof(int)); - - no_pixels=0; - maxval=0; - - for (n=0; n < mrows*ncols; n++) - { - if (img[n] > maxval) maxval=img[n]; - cellsizes[n]=0; - } - - n=0; no_cells=1; - - // initialising label array B - for (y=0; y < A->rows; y++) - for (x=0; x < A->cols; x++) - { - B->rptr[x][y]=0.0; // I.e. the pixel is UNAVAILABLE - - } - - - // stepsize - int step_size = 1; - for (level=1; level <= maxval + 1; level=level+step_size) // grow regions until max intensity is reached - { - - for (y=0; y < A->rows; y++) - for (x=0; x < A->cols; x++) - { - // identify unallocated pixels at this level which have not been assigned to a cell region - // i.e. only certain pixels can be 'used' at every round (=> beneath current _level_) - if (A->rptr[x][y] <= level && A->rptr[x][y] > 0 && B->rptr[x][y] < 2) - B->rptr[x][y]=1.0; // i.e. pixel is now AVAILABLE for assignment! - - } - - // growing loop for ten iterations (i.e. the AVAILABLE pixel can be assigned to confining regions) - int region_growing_limit = 10; - for (i=0; i < region_growing_limit; i++) - growingCells(A,B, cellsizes); - - // finding regions (AVAILABLE pixel could cluster without ever touching => lookout for new regions) - // => This function creates the first seeding. - no_cells=findingNewCells(A, B, no_cells); - - //merge regions that touch that are below minimal cell size or if the intenisity threshold - //for small cells hasn't yet been reached. - mergeCells(B, cellsizes, mincellsize, (threshold > level)); - - } - - - free( cellsizes); - mx_free( A); mx_free(B); - -} +/* + + 1. finding non allocated regions based on a intensity threshold + 2. finding new cells based on clustered unallocated regions + 3. grow cells with a maximum number of iterations and below the intensity threshold + 4. merge regions on the basis of cell size and cell intenity + + */ + + +#include "mex.h" +#include "utilah.h" + + +//pixel list +typedef struct +{ + int x; + int y; + int l; // label + +}pix; + + +// growing all existing cells (or rather, assigning newly available pixels to present regions) +void growingCells(matx* A, matx* B, int* cellsizes) //imageMatrix A, label B, cellsizes = how many pixels per cell +{ + int mrows=B->rows; + int ncols=B->cols; + int n,x,y; + int no_pixels=0; + + pix* pixlist = (pix*)malloc(ncols*mrows*sizeof(pix)); + + + for (y=0; y < A->rows; y++) + { + for (x=0; x < A->cols; x++) + { + if (B->rptr[x][y] == 1.0) + { + //Found an available ( try to assign it to a cell region + + //check neighborhood for possible assingnment + if(x < ncols-1 && B->rptr[x+1][y] > 1) + { + pixlist[no_pixels].x=x; + pixlist[no_pixels].y=y ; + pixlist[no_pixels].l=B->rptr[x+1][y]; + no_pixels++; + continue; + } + if(x > 0 && B->rptr[x-1][y] > 1) + { pixlist[no_pixels].x=x; pixlist[no_pixels].y=y ;pixlist[no_pixels].l=B->rptr[x-1][y]; no_pixels++; continue; } + if(y > 0 && B->rptr[x][y-1] > 1) + { pixlist[no_pixels].x=x; pixlist[no_pixels].y=y ;pixlist[no_pixels].l=B->rptr[x][y-1]; no_pixels++; continue; } + if(y < mrows-1 && B->rptr[x][y+1] > 1) + { pixlist[no_pixels].x=x; pixlist[no_pixels].y=y ;pixlist[no_pixels].l=B->rptr[x][y+1]; no_pixels++; continue; } + } + } + } + + + //Post Update of Label Matrix B (i.e. Prevent cell from growing indefinitely) + for (n=0; n < no_pixels; n++) + { + B->rptr[pixlist[n].x][pixlist[n].y]=pixlist[n].l; + cellsizes[ pixlist[n].l]++; + } + + free( pixlist); +} + + + +// search for cluster of unallocated but AVAILABLE pixels +int findingNewCells(matx* A, matx* B, int no_cells) +{ + + int x, y, x2, y2; + int no_pixels=0; + + //whole image with border + for (y2=3; y2 < A->rows-3; y2++) + for (x2=3; x2 < A->cols-3; x2++) + { + no_pixels=0; + + //going through rectangular region 5x5 + for (y=-2; y <= 2; y++) + for (x=-2; x <= 2; x++) + { + if (B->rptr[x2+x][y2+y] == 1.0) //i.e. only AVAILABLE pixels count! + no_pixels++; + + } + + int min_cluster_size = 10; + + // found a new cell when the area consists of more than 10 pixels + if (no_pixels > min_cluster_size) + { + no_cells++; // cell label == cell number + for (y=-2; y <= 2; y++) + for (x=-2; x <= 2; x++) + { + // mark the new cell on the label matrix + //if (B->rptr[x2+x][y2+y] == 1.0) TODO, should be corrected with here! + B->rptr[x2+x][y2+y] = no_cells; + } + } + } + + return no_cells; +} + + +// merge cell1 and cell2 which becomes cell2 +void mergeTwoCells( matx* B, int* cellsizes, int cell1, int cell2) +{ + + int x,y; + int counter=0; + int label1; + + for (y=0; y < B->rows; y++) + for (x=0; x < B->cols; x++) + { + label1=B->rptr[x][y]; + if (label1 == cell1) { B->rptr[x][y]=cell2; counter++; } + + } + + cellsizes[cell2]=cellsizes[cell2]+counter; + cellsizes[cell1]=0; +} + + +//merge touching regions if conditions are met +void mergeCells(matx* B, int* cellsizes, double mincellsize, bool belowthreshold) +{ + int x,y; + int label1,label2; + + for (y=1; y < B->rows-1; y++) + { + for (x=1; x < B->cols-1; x++) + { + // check whether two cells are touching + label1=B->rptr[x][y]; + label2=B->rptr[x+1][y]; + + if (label1 > 1 && label2 > 1 && label1 != label2) + { + if ((cellsizes[label1] < mincellsize || cellsizes[label2] < mincellsize) || belowthreshold) + { + // merging two cells + mergeTwoCells( B, cellsizes, label2, label1); + } + } + + label2=B->rptr[x][y+1]; + if (label1 > 1 && label2 > 1 && label1 != label2) + { + if ((cellsizes[label1] < mincellsize && cellsizes[label2] < mincellsize) || belowthreshold) + { + // merging two cells + mergeTwoCells( B, cellsizes, label2, label1); + } + } + } + } +} + + + + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + + matx *A, *B; + matx *pixel_list; + int mrows, ncols; + int n,x,y,count, no_pixels, no_cells, i; + double mincellsize,level, maxval, threshold; + double *img; + int* cellsizes; + + + // --- testing input and output values --- + if (nrhs != 3) { mexErrMsgTxt(" usage: B=findscellsfromregiongrowing(CellImage(A), min cell size, threshold)"); } + if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) + mexErrMsgTxt("Image input must be a double real matrix"); + if (!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1])) + mexErrMsgTxt("Image input must be a double real matrix"); + + // -- initialising arrays and stuff --- //matlab to c++ conversion + mrows = mxGetM(prhs[0]); + ncols = mxGetN(prhs[0]); + A=mx_from_vector( mxGetPr(prhs[0]), mrows, ncols); + + img=mxGetPr(prhs[0]); + + mincellsize = mxGetScalar( prhs[1]); + threshold = mxGetScalar( prhs[2]); + + plhs[0] = mxCreateDoubleMatrix( mrows,ncols, mxREAL); + B=mx_from_vector( mxGetPr(plhs[0]), mrows, ncols); + + // estimate number of pixels per image for the pixelist + // find the max value + cellsizes=(int*)malloc(mrows*ncols*sizeof(int)); + + no_pixels=0; + maxval=0; + + for (n=0; n < mrows*ncols; n++) + { + if (img[n] > maxval) maxval=img[n]; + cellsizes[n]=0; + } + + n=0; no_cells=1; + + // initialising label array B + for (y=0; y < A->rows; y++) + for (x=0; x < A->cols; x++) + { + B->rptr[x][y]=0.0; // I.e. the pixel is UNAVAILABLE + + } + + + // stepsize + int step_size = 1; + for (level=1; level <= maxval + 1; level=level+step_size) // grow regions until max intensity is reached + { + + for (y=0; y < A->rows; y++) + for (x=0; x < A->cols; x++) + { + // identify unallocated pixels at this level which have not been assigned to a cell region + // i.e. only certain pixels can be 'used' at every round (=> beneath current _level_) + if (A->rptr[x][y] <= level && A->rptr[x][y] > 0 && B->rptr[x][y] < 2) + B->rptr[x][y]=1.0; // i.e. pixel is now AVAILABLE for assignment! + + } + + // growing loop for ten iterations (i.e. the AVAILABLE pixel can be assigned to confining regions) + int region_growing_limit = 10; + for (i=0; i < region_growing_limit; i++) + growingCells(A,B, cellsizes); + + // finding regions (AVAILABLE pixel could cluster without ever touching => lookout for new regions) + // => This function creates the first seeding. + no_cells=findingNewCells(A, B, no_cells); + + //merge regions that touch that are below minimal cell size or if the intenisity threshold + //for small cells hasn't yet been reached. + mergeCells(B, cellsizes, mincellsize, (threshold > level)); + + } + + + free( cellsizes); + mx_free( A); mx_free(B); + +} diff --git a/findcellsfromregiongrowing.mexmaci64 b/findcellsfromregiongrowing.mexmaci64 old mode 100755 new mode 100644 diff --git a/growcellsfromseeds3.cc b/growcellsfromseeds3.cc old mode 100755 new mode 100644 index 77497c3..49b3efe --- a/growcellsfromseeds3.cc +++ b/growcellsfromseeds3.cc @@ -1,212 +1,212 @@ -/* - - 1. Given the seeds with a certain region - 2. Find the boundary of the cells by fully "growing them" - - */ - -#include "mex.h" -#include "utilah.h" - -//pixel strucure -typedef struct -{ - int x; - int y; -}pix; - - - -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) -{ - //plhs => left hand side/output; prhs => right hand side/ input - matx *A, *B; - matx *pixel_list; - int mrows, ncols; - int n,x,y,count, no_pixels1, no_pixels2, i; - double level, celllevel, maxval; - double *img; - pix* pixlist; - pix* pixlist2; - pix* tmp; - double cellID; - - - // --- testing input and output values --- - if (nrhs != 2) { mexErrMsgTxt(" usage: B=growcellsfromseeds2(CellImage with seeds at 255, celllevel)"); } - if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) - mexErrMsgTxt("Image input must be a double real matrix"); - if (!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1])) - mexErrMsgTxt("Image input must be a double real matrix"); - - // -- initialising arrays and stuff --- - mrows = mxGetM(prhs[0]); - ncols = mxGetN(prhs[0]); - A=mx_from_vector( mxGetPr(prhs[0]), mrows, ncols); - - img=mxGetPr(prhs[0]); - - //threshold that defines seed location (intensity range of seeds) - celllevel = mxGetScalar( prhs[1]); - - - plhs[0] = mxCreateDoubleMatrix( mrows,ncols, mxREAL); - B=mx_from_vector( mxGetPr(plhs[0]), mrows, ncols); - - // estimate number of pixels per image for the pixellist - // find the max value - - no_pixels1=0; - maxval=0; - - for (n=0; n < mrows*ncols; n++) - { - if (img[n] > maxval) maxval=img[n]; - } - - // creating initial pixel list - - pixlist = (pix*)malloc(mrows*ncols*sizeof(pix)); - pixlist2 = (pix*)malloc(mrows*ncols*sizeof(pix)); - - n=0; - - // initial populating pixlist and label matrix B - // pixlist just contains the seed points - for (y=0; y < A->rows; y++) - for (x=0; x < A->cols; x++) - { - if (A->rptr[x][y] > celllevel) - { - pixlist[n].x=x; - pixlist[n].y=y; - B->rptr[x][y]=n+1; // assign new cell ID to output matrix - n++; no_pixels1++; - } - else - B->rptr[x][y]=0; - } - - //printf("maxval: %d, n: %d\n", (int)maxval, no_pixels1); - - for (level=1; level <= maxval; level=level+0.05) //maxval - { - no_pixels2=0; - - // growing cells - - //2 pixel lists - //alternative storing to avoid indefinitive growth - //i.e. the label matrix is updated only at the end to not influence the current run - - for (n=0; n < no_pixels1; n++) - { - x=pixlist[n].x; y=pixlist[n].y; - - //found new pixels and using - pixlist2[no_pixels2].x=x; - pixlist2[no_pixels2].y=y; - - no_pixels2++; - cellID=B->rptr[x][y]; - - // Are we touching another cell? If yes stop growing - - if(x < ncols-1 && B->rptr[x+1][y] > 0 && B->rptr[x+1][y] != cellID) { continue; } - if(x > 0 && B->rptr[x-1][y] > 0 && B->rptr[x-1][y] != cellID) { continue; } - if(y > 0 && B->rptr[x][y-1] > 0 && B->rptr[x][y-1] != cellID) { continue; } - if(y < mrows-1 && B->rptr[x][y+1] > 0 && B->rptr[x][y+1] != cellID) { continue; } - - - // trying to grow - are we having pixels : - // - connecting unallocated pixels to cells [ B->rptr[x+1][y] < 1 ] - // - below the current threshold level [ A->rptr[x+1][y] < level ] - - - if(x < ncols-1 && B->rptr[x+1][y] < 1 && A->rptr[x+1][y] < level) //grow by level again //todo avoid background explosion - { - B->rptr[x+1][y] = cellID; //growth! new label: B->rptr[x][y] - pixlist2[no_pixels2].x=x+1; - pixlist2[no_pixels2].y=y; - no_pixels2++; - } - - - if(x > 0 && B->rptr[x-1][y] < 1 && A->rptr[x-1][y] < level) - { - B->rptr[x-1][y]=cellID; - pixlist2[no_pixels2].x=x-1; - pixlist2[no_pixels2].y=y; - no_pixels2++; - } - - - if(y > 0 && B->rptr[x][y-1] < 1 && A->rptr[x][y-1] < level) - { - B->rptr[x][y-1]=cellID; - pixlist2[no_pixels2].x=x; - pixlist2[no_pixels2].y=y-1; - no_pixels2++; - } - - if(y < mrows-1 && B->rptr[x][y+1] < 1 && A->rptr[x][y+1] < level) - { - B->rptr[x][y+1]=cellID; - pixlist2[no_pixels2].x=x; - pixlist2[no_pixels2].y=y+1; - no_pixels2++; - } - - } - - //printf("level: %f, nopixels1 %d and nopixels2 %d\n", level, no_pixels1, no_pixels2); - - - // populating pixellist from old pixellist with allocated pixels only - no_pixels1 =0; - - for (n=0; n < no_pixels2; n++) - { - x=pixlist2[n].x; - y=pixlist2[n].y; - - //look for labels [B->rptr[x][y] > 0] - if (B->rptr[x][y] > 0) - { - pixlist[no_pixels1].x=x; - pixlist[no_pixels1].y=y; - no_pixels1++; - } - - } - } - - - // ------ post processing of completed segmentation to remove small gaps left over by the initial growing - for (n=0; n < no_pixels1; n++) - { - x=pixlist[n].x; y=pixlist[n].y; - - if(x < ncols-1 && B->rptr[x+1][y] < 1 && A->rptr[x+1][y] < level) //grow by level again - { - B->rptr[x+1][y] = B->rptr[x][y]; //growth! new label: B->rptr[x][y] - } - if(x > 0 && B->rptr[x-1][y] < 1 && A->rptr[x-1][y] < level) - { - B->rptr[x-1][y]=B->rptr[x][y]; - } - if(y > 0 && B->rptr[x][y-1] < 1 && A->rptr[x][y-1] < level) - { - B->rptr[x][y-1]=B->rptr[x][y]; - } - if(y < mrows-1 && B->rptr[x][y+1] < 1 && A->rptr[x][y+1] < level) - { - B->rptr[x][y+1]=B->rptr[x][y]; - } - - } - - - free( pixlist); free( pixlist2); - mx_free( A); mx_free(B); -} +/* + + 1. Given the seeds with a certain region + 2. Find the boundary of the cells by fully "growing them" + + */ + +#include "mex.h" +#include "utilah.h" + +//pixel strucure +typedef struct +{ + int x; + int y; +}pix; + + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + //plhs => left hand side/output; prhs => right hand side/ input + matx *A, *B; + matx *pixel_list; + int mrows, ncols; + int n,x,y,count, no_pixels1, no_pixels2, i; + double level, celllevel, maxval; + double *img; + pix* pixlist; + pix* pixlist2; + pix* tmp; + double cellID; + + + // --- testing input and output values --- + if (nrhs != 2) { mexErrMsgTxt(" usage: B=growcellsfromseeds2(CellImage with seeds at 255, celllevel)"); } + if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) + mexErrMsgTxt("Image input must be a double real matrix"); + if (!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1])) + mexErrMsgTxt("Image input must be a double real matrix"); + + // -- initialising arrays and stuff --- + mrows = mxGetM(prhs[0]); + ncols = mxGetN(prhs[0]); + A=mx_from_vector( mxGetPr(prhs[0]), mrows, ncols); + + img=mxGetPr(prhs[0]); + + //threshold that defines seed location (intensity range of seeds) + celllevel = mxGetScalar( prhs[1]); + + + plhs[0] = mxCreateDoubleMatrix( mrows,ncols, mxREAL); + B=mx_from_vector( mxGetPr(plhs[0]), mrows, ncols); + + // estimate number of pixels per image for the pixellist + // find the max value + + no_pixels1=0; + maxval=0; + + for (n=0; n < mrows*ncols; n++) + { + if (img[n] > maxval) maxval=img[n]; + } + + // creating initial pixel list + + pixlist = (pix*)malloc(mrows*ncols*sizeof(pix)); + pixlist2 = (pix*)malloc(mrows*ncols*sizeof(pix)); + + n=0; + + // initial populating pixlist and label matrix B + // pixlist just contains the seed points + for (y=0; y < A->rows; y++) + for (x=0; x < A->cols; x++) + { + if (A->rptr[x][y] > celllevel) + { + pixlist[n].x=x; + pixlist[n].y=y; + B->rptr[x][y]=n+1; // assign new cell ID to output matrix + n++; no_pixels1++; + } + else + B->rptr[x][y]=0; + } + + //printf("maxval: %d, n: %d\n", (int)maxval, no_pixels1); + + for (level=1; level <= maxval; level=level+0.05) //maxval + { + no_pixels2=0; + + // growing cells + + //2 pixel lists + //alternative storing to avoid indefinitive growth + //i.e. the label matrix is updated only at the end to not influence the current run + + for (n=0; n < no_pixels1; n++) + { + x=pixlist[n].x; y=pixlist[n].y; + + //found new pixels and using + pixlist2[no_pixels2].x=x; + pixlist2[no_pixels2].y=y; + + no_pixels2++; + cellID=B->rptr[x][y]; + + // Are we touching another cell? If yes stop growing + + if(x < ncols-1 && B->rptr[x+1][y] > 0 && B->rptr[x+1][y] != cellID) { continue; } + if(x > 0 && B->rptr[x-1][y] > 0 && B->rptr[x-1][y] != cellID) { continue; } + if(y > 0 && B->rptr[x][y-1] > 0 && B->rptr[x][y-1] != cellID) { continue; } + if(y < mrows-1 && B->rptr[x][y+1] > 0 && B->rptr[x][y+1] != cellID) { continue; } + + + // trying to grow - are we having pixels : + // - connecting unallocated pixels to cells [ B->rptr[x+1][y] < 1 ] + // - below the current threshold level [ A->rptr[x+1][y] < level ] + + + if(x < ncols-1 && B->rptr[x+1][y] < 1 && A->rptr[x+1][y] < level) //grow by level again //todo avoid background explosion + { + B->rptr[x+1][y] = cellID; //growth! new label: B->rptr[x][y] + pixlist2[no_pixels2].x=x+1; + pixlist2[no_pixels2].y=y; + no_pixels2++; + } + + + if(x > 0 && B->rptr[x-1][y] < 1 && A->rptr[x-1][y] < level) + { + B->rptr[x-1][y]=cellID; + pixlist2[no_pixels2].x=x-1; + pixlist2[no_pixels2].y=y; + no_pixels2++; + } + + + if(y > 0 && B->rptr[x][y-1] < 1 && A->rptr[x][y-1] < level) + { + B->rptr[x][y-1]=cellID; + pixlist2[no_pixels2].x=x; + pixlist2[no_pixels2].y=y-1; + no_pixels2++; + } + + if(y < mrows-1 && B->rptr[x][y+1] < 1 && A->rptr[x][y+1] < level) + { + B->rptr[x][y+1]=cellID; + pixlist2[no_pixels2].x=x; + pixlist2[no_pixels2].y=y+1; + no_pixels2++; + } + + } + + //printf("level: %f, nopixels1 %d and nopixels2 %d\n", level, no_pixels1, no_pixels2); + + + // populating pixellist from old pixellist with allocated pixels only + no_pixels1 =0; + + for (n=0; n < no_pixels2; n++) + { + x=pixlist2[n].x; + y=pixlist2[n].y; + + //look for labels [B->rptr[x][y] > 0] + if (B->rptr[x][y] > 0) + { + pixlist[no_pixels1].x=x; + pixlist[no_pixels1].y=y; + no_pixels1++; + } + + } + } + + + // ------ post processing of completed segmentation to remove small gaps left over by the initial growing + for (n=0; n < no_pixels1; n++) + { + x=pixlist[n].x; y=pixlist[n].y; + + if(x < ncols-1 && B->rptr[x+1][y] < 1 && A->rptr[x+1][y] < level) //grow by level again + { + B->rptr[x+1][y] = B->rptr[x][y]; //growth! new label: B->rptr[x][y] + } + if(x > 0 && B->rptr[x-1][y] < 1 && A->rptr[x-1][y] < level) + { + B->rptr[x-1][y]=B->rptr[x][y]; + } + if(y > 0 && B->rptr[x][y-1] < 1 && A->rptr[x][y-1] < level) + { + B->rptr[x][y-1]=B->rptr[x][y]; + } + if(y < mrows-1 && B->rptr[x][y+1] < 1 && A->rptr[x][y+1] < level) + { + B->rptr[x][y+1]=B->rptr[x][y]; + } + + } + + + free( pixlist); free( pixlist2); + mx_free( A); mx_free(B); +} diff --git a/growcellsfromseeds3.mexmaci64 b/growcellsfromseeds3.mexmaci64 old mode 100755 new mode 100644 diff --git a/utilah.h b/utilah.h index 6ff6361..7bed0f5 100644 --- a/utilah.h +++ b/utilah.h @@ -59,7 +59,7 @@ void mx_free( matx* m) { if (m->newmatrix == 1) free( m->rptr[0]); free(m->rptr); m->rptr=NULL; - free( m); + free( m); } matx* mx_clone(matx* a) @@ -260,9 +260,9 @@ double mx_find(matx* a, double value) for (y=0; y < a->rows; y++) { if (a->rptr[x][y]==value) sum=sum+1; - } - - return sum; + } + + return sum; } void mx_clear(matx* a) @@ -359,12 +359,12 @@ double mx_corr(matx* a, matx* b) { register int x,y; double meanA, meanB, sumAB=0.0, sumA2=0.0, sumB2=0.0; - + if (a->rows != b->rows || a->cols != b->cols) return 0.0; - + meanA = mx_mean(a); meanB= mx_mean(b); - + for(x=0; x < a->cols; x++) for (y=0; y < a->rows; y++) { @@ -372,7 +372,6 @@ double mx_corr(matx* a, matx* b) sumA2=sumA2+((a->rptr[x][y]-meanA)*(a->rptr[x][y]-meanA)); sumB2=sumB2+((b->rptr[x][y]-meanB)*(b->rptr[x][y]-meanB)); } - + return sumAB/(sqrt(sumA2*sumB2)); } - From 1fe113eb69fc915f8606333bd8a8f50658ed48d8 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Mon, 1 Aug 2022 17:04:02 +0100 Subject: [PATCH 04/18] Make start on porting --- example_segmentation.py | 44 +++++++++++++++++++++++++++++++++++++++++ pyproject.toml | 2 ++ requirements.txt | 2 ++ 3 files changed, 48 insertions(+) create mode 100644 example_segmentation.py create mode 100644 pyproject.toml create mode 100644 requirements.txt diff --git a/example_segmentation.py b/example_segmentation.py new file mode 100644 index 0000000..cdad7ec --- /dev/null +++ b/example_segmentation.py @@ -0,0 +1,44 @@ +from pathlib import Path + +from scipy import io as sio +from segment_image import segment_image + +_file_location = Path(__file__).resolve() +_matfile = _file_location.parent / "ProjIm.mat" + + +def main() -> None: + """ + implements the demo in the README in python + """ + # load example image (Projected Drosophila Wing Disc - Ecad:GFP) + proj_im = sio.loadmat(_matfile)["ProjIm"] + + # crop image for testing + crop = proj_im[300:500, 400:600] + + # compute segmentation and output segmentation feedback + first_smoothing_seeding = 0.5 + min_cell_area = 2 + min_membrane_intensity = 20 + min_border_intensity_ratio = 0.0 + second_smoothing_segmentation = 0.5 + max_cell_area = 3000 + min_seed_per_boundary = 0.1 + segmentation_feedback_plot = True + + segmentation, seeds, labels = segment_image( + crop, + first_smoothing_seeding, + min_cell_area, + min_membrane_intensity, + min_border_intensity_ratio, + second_smoothing_segmentation, + max_cell_area, + min_seed_per_boundary, + segmentation_feedback_plot, + ) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..0ba3bbd --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,2 @@ +[tool.mypy] +plugins = ["numpy.typing.mypy_plugin"] diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..f9f587b --- /dev/null +++ b/requirements.txt @@ -0,0 +1,2 @@ +numpy==1.23.1 +scipy==1.9.0 From fe073201e673bdfc882dee756aa27a80892bc27f Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Tue, 2 Aug 2022 14:14:19 +0100 Subject: [PATCH 05/18] Fix pre-commit --- pyproject.toml | 2 -- setup.cfg | 10 ++++++++++ 2 files changed, 10 insertions(+), 2 deletions(-) delete mode 100644 pyproject.toml create mode 100644 setup.cfg diff --git a/pyproject.toml b/pyproject.toml deleted file mode 100644 index 0ba3bbd..0000000 --- a/pyproject.toml +++ /dev/null @@ -1,2 +0,0 @@ -[tool.mypy] -plugins = ["numpy.typing.mypy_plugin"] diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..04d6e3c --- /dev/null +++ b/setup.cfg @@ -0,0 +1,10 @@ +import numpy + +[flake8] +ignore = E203,E741,W503,W605 +max-complexity = 18 +max-line-length = 88 +select = B,B9,BLK,C,E,F,I,S,T4,W + +[mypy] +plugins = numpy.typing.mypy_plugin From 9f40b2cdb3762f78ee3d5c6668dd29ef33c397ea Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Tue, 2 Aug 2022 14:27:23 +0100 Subject: [PATCH 06/18] Move mypy plugin for now --- example_segmentation.py | 31 +++++++++++-------------------- setup.cfg | 5 ----- 2 files changed, 11 insertions(+), 25 deletions(-) diff --git a/example_segmentation.py b/example_segmentation.py index cdad7ec..a4093f1 100644 --- a/example_segmentation.py +++ b/example_segmentation.py @@ -1,6 +1,7 @@ from pathlib import Path from scipy import io as sio + from segment_image import segment_image _file_location = Path(__file__).resolve() @@ -15,28 +16,18 @@ def main() -> None: proj_im = sio.loadmat(_matfile)["ProjIm"] # crop image for testing - crop = proj_im[300:500, 400:600] - - # compute segmentation and output segmentation feedback - first_smoothing_seeding = 0.5 - min_cell_area = 2 - min_membrane_intensity = 20 - min_border_intensity_ratio = 0.0 - second_smoothing_segmentation = 0.5 - max_cell_area = 3000 - min_seed_per_boundary = 0.1 - segmentation_feedback_plot = True + cropped_image = proj_im[300:500, 400:600] segmentation, seeds, labels = segment_image( - crop, - first_smoothing_seeding, - min_cell_area, - min_membrane_intensity, - min_border_intensity_ratio, - second_smoothing_segmentation, - max_cell_area, - min_seed_per_boundary, - segmentation_feedback_plot, + cropped_image, + sigma_1=0.5, + min_cell_size=2, + threshold=20, + merge_criteria=0.0, + sigma_3=0.5, + large_cell_size_thres=3000, + i_bound_max_pcnt=0.1, + show_feedback=True, ) diff --git a/setup.cfg b/setup.cfg index 04d6e3c..87c3241 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,10 +1,5 @@ -import numpy - [flake8] ignore = E203,E741,W503,W605 max-complexity = 18 max-line-length = 88 select = B,B9,BLK,C,E,F,I,S,T4,W - -[mypy] -plugins = numpy.typing.mypy_plugin From 7f8b0c7a050fac216c2da2a0e6e39bfa42db51ae Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Tue, 2 Aug 2022 16:59:41 +0100 Subject: [PATCH 07/18] Fix pre-commit --- .pre-commit-config.yaml | 1 + setup.cfg | 3 +++ 2 files changed, 4 insertions(+) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index dba1c99..a8c6cac 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -35,6 +35,7 @@ repos: rev: v0.942 hooks: - id: mypy + additional_dependencies: ["numpy"] - repo: https://github.com/pre-commit/mirrors-prettier rev: v2.6.1 hooks: diff --git a/setup.cfg b/setup.cfg index 87c3241..8129364 100644 --- a/setup.cfg +++ b/setup.cfg @@ -3,3 +3,6 @@ ignore = E203,E741,W503,W605 max-complexity = 18 max-line-length = 88 select = B,B9,BLK,C,E,F,I,S,T4,W + +[mypy] +plugins = numpy.typing.mypy_plugin From 217ace90d78910bf81f978e11349a87b4250f373 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Tue, 2 Aug 2022 17:58:08 +0100 Subject: [PATCH 08/18] Pseudocode --- segment_image.py | 78 ++++++++++++++++++++++++++++++++++++++++++++++++ utils.py | 75 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 153 insertions(+) create mode 100644 segment_image.py create mode 100644 utils.py diff --git a/segment_image.py b/segment_image.py new file mode 100644 index 0000000..9f590b7 --- /dev/null +++ b/segment_image.py @@ -0,0 +1,78 @@ +import numpy as np +import numpy.typing as npt + +from utils import ( + create_color_boundaries, + do_initial_seeding, + grow_cells_in_frame, + merge_seeds_from_labels, + neutralise_pts_not_under_label_in_frame, + unlabel_poor_seeds_in_frame, +) + + +def segment_image( + image: npt.NDArray[np.uint8], + *, + sigma_1: float, + min_cell_size: float, + threshold: int, + merge_criteria: float, + sigma_3: float, + large_cell_size_thres: float, + i_bound_max_pcnt: float, + show_feedback: bool = False +) -> tuple[npt.NDArray[np.uint8], npt.NDArray[np.uint16], npt.NDArray[np.uint8]]: + """ + Segments a single frame extracting the cell outlines. + + Args: + image: + increasing intesity for membrane + sigma_1: + size px of gaussian for smoothing image [0+] + min_cell_size: + size px of smallest cell expected [0+] + threshold: + minimum value for membrane signal [0-255] + merge_criteria: + minimum ratio of low intensity pxs upon which to merge cells [0-1] + sigma_3: + size px of gaussian for smoothing image [0+] + large_cell_size_thres: + size px of largest cell expected [0+] + i_bound_max_pcnt: + minimum ratio for seed and membrane intensity [0-1] + show_feedback: + show feedback for segmentation + + Returns: + segmented_image: + Im with seeds [255] and cell outlines [0] + cell_seeds: + Rescaled Image to fit the range [0,252] + 253/254/255 are used for seed information + cell_labels: + bitmap of cells colored with 16bit id + """ + # initialise + cell_seeds = np.zeros(image.shape, dtype=np.uint8) + # TODO: why using 16 bit for labels? + cell_labels = np.zeros(image.shape, dtype=np.uint16) + + # TODO: check image casting + image = image.astype(float) + image *= 252 / image.max() + image = image.astype(np.uint8) + + do_initial_seeding(sigma_1, min_cell_size, threshold) + + merge_seeds_from_labels(merge_criteria, sigma_3) + + grow_cells_in_frame(sigma_3) + + unlabel_poor_seeds_in_frame(threshold, sigma_3) + + neutralise_pts_not_under_label_in_frame(cell_labels, cell_seeds, value=253) + + create_color_boundaries() diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..756479b --- /dev/null +++ b/utils.py @@ -0,0 +1,75 @@ +import numpy as np +import numpy.typing as npt + + +def create_color_boundaries(image: npt.NDArray[np.uint8], cell_labels: npt.NDArray[np.uint8]) -> None: + """ + Generate final colored image (RGB) to represent the segmentation results + """ + # given that every cell has a different label we can compute + # the boundaries by computing where the gradient changes + cell_lables = cell_lables.astype(float) + [gx,gy] = gradient(cell_lables); + cell_outlines = (cell_lables > 0) & ((gx**2+gy**2)>0) + + cell_seeds = CellSeeds(:,:) > 253; + + ColIm = Im + ColIm(cell_outlines) = 0 + ColIm(cell_seeds) = 255 + + +def do_initial_seeding(sigma_1: float, min_cell_size: float, threshold: int) -> None: + """ + Find the initial cell seeds + """ + pass + + +def grow_cells_in_frame(sigma_3: float) -> None: + """ + Growing cells from seeds TODO: add paramters in Name description! + """ + pass + + +def unlabel_poor_seeds_in_frame(threshold: int, sigma_3: float) -> None: + """ + Eliminate labels from seeds which have poor boundary intensity + """ + pass + + +def delabel_very_large_areas() -> None: + """ """ + pass + + +def merge_seeds_from_labels(merge_criteria: float, sigma_3: float) -> None: + """ + Remove initial cell regions which touch & whose boundary is insufficient + """ + pass + + +def merge_labels() -> None: + """ """ + pass + + +def neutralise_pts_not_under_label_in_frame( + cell_labels: npt.NDArray[np.uint16], + cell_seeds: npt.NDArray[np.uint8], + *, + value: int = 253 +) -> npt.NDArray[np.uint8]: + """ + Seeds whose label has been eliminated are converted to neutral seeds + + the idea here is to set seeds not labelled to a value + ie invisible to retracking (and to growing, caution!) + """ + cell_seeds_copy = cell_seeds.copy() + cell_seeds_copy[cell_labels != 0] = 0 + cell_seeds[cell_seeds_copy > 252] = value + return cell_seeds From c2ef2fb4662ca5a04ca45ae07766928f13a0c525 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Wed, 3 Aug 2022 17:31:25 +0100 Subject: [PATCH 09/18] Some more porting --- example_segmentation.py | 2 +- requirements.txt | 1 + segment_image.py | 20 ++- utils.py | 387 +++++++++++++++++++++++++++++++++++++--- 4 files changed, 372 insertions(+), 38 deletions(-) diff --git a/example_segmentation.py b/example_segmentation.py index a4093f1..22b5857 100644 --- a/example_segmentation.py +++ b/example_segmentation.py @@ -26,7 +26,7 @@ def main() -> None: merge_criteria=0.0, sigma_3=0.5, large_cell_size_thres=3000, - i_bound_max_pcnt=0.1, + I_bound_max_pcnt=0.1, show_feedback=True, ) diff --git a/requirements.txt b/requirements.txt index f9f587b..81814bc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,3 @@ numpy==1.23.1 +scikit-image==0.19.3 scipy==1.9.0 diff --git a/segment_image.py b/segment_image.py index 9f590b7..97248ef 100644 --- a/segment_image.py +++ b/segment_image.py @@ -20,9 +20,9 @@ def segment_image( merge_criteria: float, sigma_3: float, large_cell_size_thres: float, - i_bound_max_pcnt: float, + I_bound_max_pcnt: float, show_feedback: bool = False -) -> tuple[npt.NDArray[np.uint8], npt.NDArray[np.uint16], npt.NDArray[np.uint8]]: +) -> tuple[npt.NDArray[np.uint8], npt.NDArray[np.uint8], npt.NDArray[np.uint16]]: """ Segments a single frame extracting the cell outlines. @@ -41,7 +41,7 @@ def segment_image( size px of gaussian for smoothing image [0+] large_cell_size_thres: size px of largest cell expected [0+] - i_bound_max_pcnt: + I_bound_max_pcnt: minimum ratio for seed and membrane intensity [0-1] show_feedback: show feedback for segmentation @@ -65,14 +65,16 @@ def segment_image( image *= 252 / image.max() image = image.astype(np.uint8) - do_initial_seeding(sigma_1, min_cell_size, threshold) + do_initial_seeding(image, sigma_1, min_cell_size, threshold, large_cell_size_thres) - merge_seeds_from_labels(merge_criteria, sigma_3) + merge_seeds_from_labels(image, cell_seeds, cell_labels, merge_criteria, sigma_3) - grow_cells_in_frame(sigma_3) + grow_cells_in_frame(image, cell_seeds, sigma_3) - unlabel_poor_seeds_in_frame(threshold, sigma_3) + unlabel_poor_seeds_in_frame( + image, cell_seeds, cell_labels, threshold, sigma_3, I_bound_max_pcnt + ) - neutralise_pts_not_under_label_in_frame(cell_labels, cell_seeds, value=253) + neutralise_pts_not_under_label_in_frame(cell_seeds, cell_labels) - create_color_boundaries() + create_color_boundaries(image, cell_seeds, cell_labels) diff --git a/utils.py b/utils.py index 756479b..a187f39 100644 --- a/utils.py +++ b/utils.py @@ -1,75 +1,406 @@ import numpy as np import numpy.typing as npt +from skimage.measure import regionprops -def create_color_boundaries(image: npt.NDArray[np.uint8], cell_labels: npt.NDArray[np.uint8]) -> None: +def create_color_boundaries( + image: npt.NDArray[np.uint8], + cell_seeds: npt.NDArray[np.uint8], + cell_labels: npt.NDArray[np.uint16], +) -> npt.NDArray[np.uint8]: """ Generate final colored image (RGB) to represent the segmentation results """ # given that every cell has a different label we can compute # the boundaries by computing where the gradient changes - cell_lables = cell_lables.astype(float) - [gx,gy] = gradient(cell_lables); - cell_outlines = (cell_lables > 0) & ((gx**2+gy**2)>0) + cell_labels = cell_labels.astype(float) + gx, gy = np.gradient(cell_labels) + cell_outlines = (cell_labels > 0) & ((gx**2 + gy**2) > 0) - cell_seeds = CellSeeds(:,:) > 253; + cell_seeds = cell_seeds > 253 - ColIm = Im - ColIm(cell_outlines) = 0 - ColIm(cell_seeds) = 255 + segmented_image = image.copy() + segmented_image[cell_outlines] = 0 + segmented_image[cell_seeds] = 255 + return segmented_image -def do_initial_seeding(sigma_1: float, min_cell_size: float, threshold: int) -> None: +def do_initial_seeding( + image: npt.NDArray[np.uint8], + sigma_1: float, + min_cell_size: float, + threshold: int, + large_cell_size_thres: float, +) -> None: """ Find the initial cell seeds """ - pass + # Create gaussian filter + if sigma_1 > 0.01: + f1 = _matlab_style_gauss2D(image.shape, sigma_1) + + # Gaussian smoothing for the segmentation of individual cells + smoothed_image = np.fft.fftshift( + np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1)) + ).real + else: + smoothed_image = image.astype(float) + + smoothed_image /= smoothed_image.max() * 252 + + # Use external c-code to find initial seeds + initial_labelling = findcellsfromregiongrowing( + smoothed_image, min_cell_size, threshold + ) + + # set unallocated pixels to 0 + initial_labelling[initial_labelling == 1] = 0 + + # Generate cell_labels from inital_labelling + cell_labels = initial_labelling.astype(np.uint16) + + # eliminate very large areas + delabel_very_large_areas(cell_labels, large_cell_size_thres) + # Use true centre of cells as labels + centroids = np.round(calculateCellPositions(smoothed_image, cell_labels, False)) + centroids = centroids[~np.isnan(centroids.T[0])] + for n in range(len(centroids)): + smoothed_image[centroids[n, 1], centroids[n, 0]] = 255 -def grow_cells_in_frame(sigma_3: float) -> None: + # cell_seeds contains the position of the true cell center. + cell_seeds = smoothed_image.astype(np.uint8) + + +def grow_cells_in_frame( + image: npt.NDArray[np.uint8], cell_seeds: npt.NDArray[np.uint8], sigma_3: float +) -> None: """ Growing cells from seeds TODO: add paramters in Name description! """ - pass + # find labels + bw = (cell_seeds > 252).astype(float) + + if sigma_3 > 0.01: + f1 = _matlab_style_gauss2D(image.shape, sigma_3) + smoothed_image = np.fft.fftshift( + np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1)) + ).real + else: + smoothed_image = image.astype(float) + + # mark labels on image + image_with_seeds = (smoothed_image).astype(float) * (1 - bw) + 255 * bw + cell_labels = growcellsfromseeds3(image_with_seeds, 253).astype(np.uint16) -def unlabel_poor_seeds_in_frame(threshold: int, sigma_3: float) -> None: +def unlabel_poor_seeds_in_frame( + image: npt.NDArray[np.uint8], + cell_seeds: npt.NDArray[np.uint8], + cell_labels: npt.NDArray[np.uint16], + threshold: int, + sigma_3: float, + I_bound_max_pcnt: float, +) -> None: """ Eliminate labels from seeds which have poor boundary intensity """ - pass + L = cell_labels + + if sigma_3 > 0.01: + f1 = _matlab_style_gauss2D(image.shape, sigma_3) + smoothed_image = np.fft.fftshift( + np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1)) + ).real + else: + smoothed_image = image.astype(float) + # i.e. every cell is marked by one unique integer label + label_list = np.unique(L) + label_list = label_list[label_list != 0] + IBounds = np.zeros(len(label_list)) + decisions = [0, 0, 0, 0, 0] -def delabel_very_large_areas() -> None: + for c in range(len(label_list)): + mask = L == label_list[c] + cpy, cpx = np.argwhere(mask > 0) + # find region of that label + minx = np.min(cpx) + maxx = np.max(cpx) + miny = np.min(cpy) + maxy = np.max(cpy) + minx = np.max(minx - 5, 1) + miny = np.max(miny - 5, 1) + maxx = np.min(maxx + 5, image.shape[1]) + maxy = np.min(maxy + 5, image.shape[0]) + # reduced to region of the boundary + reduced_mask = mask[miny:maxy, minx:maxx] + reduced_image = smoothed_image[miny:maxy, minx:maxx] + dilated_mask = imdilate(reduced_mask, se) + eroded_mask = imerode(reduced_mask, se) + boundary_mask = dilated_mask - eroded_mask + boundary_intensities = reduced_image[boundary_mask > 0] + H = reduced_image[boundary_mask > 0] + IEr = reduced_image[eroded_mask > 0] + I_bound = np.mean(boundary_intensities) + IBounds[c] = I_bound + + # cell seed information is retrieved as comparison + F2 = cell_seeds + F2[~mask] = 0 + cpy, cpx = np.argwhere(F2 > 252) + ICentre = smoothed_image[cpy, cpx] + + I_bound_max = 255 * I_bound_max_pcnt + + # Figure out which conditions make the label invalid + # 1. I_bound_max, gives the Lower bound to the mean intensity + # 1.b condition upon that the cell seed has less + # than 20% intensity difference to the mean + # => If the cell boundary is low and not very different + # from the seed, cancel the region + first_condition = I_bound < I_bound_max & I_bound / ICentre < 1.2 + # 2. W/o (1.b) the lower bound is reduced by ~17% (1 - 0.833) to be decisive + second_condition = I_bound < I_bound_max * 5 / 6 + # 3. If the minimum retrieved in the boundary mask is 0 (dangerous!) + third_condition = np.min(boundary_intensities) == 0 + # 4. If the amount of low intensity signal (i.e. < 20) is more than 10% + fourth_condition = np.sum(H < threshold) / len(H) > 0.1 + if first_condition | second_condition | third_condition | fourth_condition: + + # The label is cancelled (inverted mask multiplication.) + cell_labels *= (mask == 0).astype(np.uint16) + + # record the removal decisions + if first_condition: + decisions[1] = decisions[1] + 1 + elif second_condition: + decisions[2] = decisions[2] + 1 + elif third_condition: + decisions[3] = decisions[3] + 1 + elif fourth_condition: + decisions[4] = decisions[4] + 1 + else: + # should not happen + decisions[5] = decisions[5] + 1 + + +def delabel_very_large_areas( + cell_labels: npt.NDArray[np.uint16], large_cell_size_thres: float +) -> None: """ """ - pass + # remove cells which are bigger than large_cell_size_thres + L = cell_labels + dim_init_L = len(L) + A = regionprops(L, "area") + As = cat(1, A.Area) + ls = np.unique(L) + for i in range(len(ls)): + l = ls[i] + if l == 0: + continue + + A = As[l] + if A > large_cell_size_thres: + L[L == l] = 0 + dim_final_L = len(L) -def merge_seeds_from_labels(merge_criteria: float, sigma_3: float) -> None: + cell_labels = L + + +def merge_seeds_from_labels( + image: npt.NDArray[np.uint8], + cell_seeds: npt.NDArray[np.uint8], + cell_labels: npt.NDArray[np.uint16], + merge_criteria: float, + sigma_3: float, +) -> None: """ Remove initial cell regions which touch & whose boundary is insufficient """ - pass + # smoothing + if sigma_3 > 0.01: + f1 = _matlab_style_gauss2D(image.shape, sigma_3) + smoothed_image = np.fft.fftshift( + np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1)) + ).real + else: + smoothed_image = image.astype(float) + + label_list = np.unique(cell_labels) + label_list = label_list[label_list != 0] + c = 1 + + merge_intensity_distro = [] + merge_decisions = 0 + + # loop over labels + while True: + label_mask = cell_labels == label_list[c] + label = label_list[c] + + cpy, cpx = np.argwhere(label_mask > 0) + # find region of that label + minx = np.min(cpx) + maxx = np.max(cpx) + miny = np.min(cpy) + maxy = np.max(cpy) + minx = np.max(minx - 5, 1) + miny = np.max(miny - 5, 1) + maxx = np.min(maxx + 5, image.shape[1]) + maxy = np.min(maxy + 5, image.shape[0]) -def merge_labels() -> None: + # reduce data to that region + reduced_label_mask = label_mask[miny:maxy, minx:maxx] + reduced_image = smoothed_image[miny:maxy, minx:maxx] + reduced_labels = cell_labels[miny:maxy, minx:maxx] + + # now find boundaries + dilated_mask = imdilate(reduced_label_mask, se) + eroded_mask = imerode(reduced_label_mask, se) + border_mask = dilated_mask - eroded_mask + border_intensities = reduced_image[border_mask > 0] + central_intensity = reduced_image[eroded_mask > 0] + + F2 = cell_seeds + F2[~label_mask] = 0 + cpy, cpx = np.argwhere(F2 > 253) + ICentre = smoothed_image[cpy, cpx] + + background_std = np.std(central_intensity.astype(float)) + + # get labels of surrounding cells (neighbours) + neighbour_labels = np.unique(reduced_labels[dilated_mask > 0]) + neighbour_labels = neighbour_labels[neighbour_labels != label] + + low_intensity_ratios = [] + for i in range(len(neighbour_labels)): + neighb_label = neighbour_labels(i) + neighbor_border = dilated_mask + # slice of neighbour around cell + neighbor_border[reduced_labels != neighb_label] = 0 + cell_border = imdilate(neighbor_border, se) + # slice of cell closest to neighbour + cell_border[reduced_labels != label] = 0 + + # combination of both creating boundary region + joint_border = (cell_border + neighbor_border) > 0 + border_intensities = reduced_image + # intensities at boundary + border_intensities[~joint_border] = 0 + + # average number of points in boundary where intensity is + # of low quality (dodgy) + low_intensity_threshold = ICentre + (background_std / 2) + low_intensity_pixels = ( + border_intensities[joint_border] < low_intensity_threshold + ) + + low_intensity_ratio = np.sum(low_intensity_pixels) / np.shape( + border_intensities[joint_border] + ) + + low_intensity_ratios = [low_intensity_ratios, low_intensity_ratio] + + # Find out which is border with the lowest intensity ratio + worst_intensity_ratio, worst_neighbor_index = np.max(low_intensity_ratios) + neighb_label = neighbour_labels(worst_neighbor_index) + + # if the label value is of poor quality, then recursively check the merge + # criteria in order to add it as a potential label in the label set. + + merge_intensity_distro = [merge_intensity_distro, worst_intensity_ratio] + + if worst_intensity_ratio > merge_criteria & label != 0 & neighb_label != 0: + + merge_labels(cell_seeds, cell_labels, label, neighb_label) + label_list = np.unique(cell_labels) + label_list = label_list[label_list != 0] + c -= 1 + # reanalyze the same cell for more possible mergings + merge_decisions += 1 + + c += 1 + + # Condition to break the while cycle -> as soon as all the + # labels are processed, then exit + if c > len(label_list): + break + + +def merge_labels( + cell_seeds: npt.NDArray[np.uint8], cell_labels: npt.NDArray[np.uint16], l1, l2 +) -> None: """ """ - pass + Cl = cell_labels + Il = cell_seeds + m1 = Cl == l1 + m2 = Cl == l2 + Il1 = Il + Il1[~m1] = 0 + Il2 = Il + Il2[~m2] = 0 + cpy1, cpx1 = np.argwhere(Il1 > 253) + cpy2, cpx2 = np.argwhere(Il2 > 253) + cpx = round((cpx1 + cpx2) / 2) + cpy = round((cpy1 + cpy2) / 2) + + # background level + cell_seeds[cpy1, cpx1] = 20 + cell_seeds[cpy2, cpx2] = 20 + if cell_labels[cpy, cpx] == l1 | cell_labels[cpy, cpx] == l2: + cell_seeds[cpy, cpx] = 255 + elif np.sum(m1) > np.sum(m2): + cell_seeds[cpy1, cpx1] = 255 + else: + cell_seeds[cpy2, cpx2] = 255 + + Cl[m2] = l1 + cell_labels = Cl def neutralise_pts_not_under_label_in_frame( - cell_labels: npt.NDArray[np.uint16], cell_seeds: npt.NDArray[np.uint8], - *, - value: int = 253 -) -> npt.NDArray[np.uint8]: + cell_labels: npt.NDArray[np.uint16], +) -> None: """ Seeds whose label has been eliminated are converted to neutral seeds the idea here is to set seeds not labelled to a value ie invisible to retracking (and to growing, caution!) """ - cell_seeds_copy = cell_seeds.copy() - cell_seeds_copy[cell_labels != 0] = 0 - cell_seeds[cell_seeds_copy > 252] = value - return cell_seeds + L = cell_labels + F = cell_seeds + F2 = F + F2[L != 0] = 0 + F[F2 > 252] = 253 + cell_seeds = F + + +def growcellsfromseeds3(): + pass + + +def findcellsfromregiongrowing(): + pass + + +def calculateCellPositions(): + pass + + +def _matlab_style_gauss2D(shape: tuple[int, int], sigma: float): + """ + 2D gaussian mask - should give the same result as MATLAB's + fspecial('gaussian', [shape], [sigma]) + """ + m, n = [(ss - 1) / 2 for ss in shape] + y,x = np.ogrid[-m:m+1,-n:n+1] + h = np.exp(-(x * x + y * y) / (2 * sigma * sigma)) + h[h < np.finfo(h.dtype).eps*h.max()] = 0 + sumh = h.sum() + if sumh != 0: + h /= sumh + return h From d167e953df4b3b4c963f8e742f18a145d2dbf047 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Thu, 4 Aug 2022 15:09:58 +0100 Subject: [PATCH 10/18] Fix `delabel_very_large_areas` --- .gitignore | 4 ++++ utils.py | 29 ++++++++++++----------------- 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/.gitignore b/.gitignore index 7a6353d..3bbb967 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,5 @@ +__pycache__ +!ProjIm.mat +!zerogreen.mat .envrc +*.mat diff --git a/utils.py b/utils.py index a187f39..2df05ea 100644 --- a/utils.py +++ b/utils.py @@ -159,7 +159,7 @@ def unlabel_poor_seeds_in_frame( # than 20% intensity difference to the mean # => If the cell boundary is low and not very different # from the seed, cancel the region - first_condition = I_bound < I_bound_max & I_bound / ICentre < 1.2 + first_condition = (I_bound < I_bound_max) & (I_bound / ICentre < 1.2) # 2. W/o (1.b) the lower bound is reduced by ~17% (1 - 0.833) to be decisive second_condition = I_bound < I_bound_max * 5 / 6 # 3. If the minimum retrieved in the boundary mask is 0 (dangerous!) @@ -187,26 +187,21 @@ def unlabel_poor_seeds_in_frame( def delabel_very_large_areas( cell_labels: npt.NDArray[np.uint16], large_cell_size_thres: float -) -> None: - """ """ - # remove cells which are bigger than large_cell_size_thres - L = cell_labels - dim_init_L = len(L) - A = regionprops(L, "area") - As = cat(1, A.Area) - ls = np.unique(L) - for i in range(len(ls)): - l = ls[i] +) -> npt.NDArray[np.uint16]: + """ + remove cells which are bigger than large_cell_size_thres + """ + As = [r.area for r in regionprops(cell_labels, cache=False)] + ls = np.unique(cell_labels) + for l in range(len(ls)): if l == 0: continue A = As[l] if A > large_cell_size_thres: - L[L == l] = 0 - - dim_final_L = len(L) + cell_labels[cell_labels == l] = 0 - cell_labels = L + return cell_labels def merge_seeds_from_labels( @@ -313,7 +308,7 @@ def merge_seeds_from_labels( merge_intensity_distro = [merge_intensity_distro, worst_intensity_ratio] - if worst_intensity_ratio > merge_criteria & label != 0 & neighb_label != 0: + if (worst_intensity_ratio > merge_criteria) & label != 0 & neighb_label != 0: merge_labels(cell_seeds, cell_labels, label, neighb_label) label_list = np.unique(cell_labels) @@ -350,7 +345,7 @@ def merge_labels( # background level cell_seeds[cpy1, cpx1] = 20 cell_seeds[cpy2, cpx2] = 20 - if cell_labels[cpy, cpx] == l1 | cell_labels[cpy, cpx] == l2: + if (cell_labels[cpy, cpx] == l1) | (cell_labels[cpy, cpx] == l2): cell_seeds[cpy, cpx] = 255 elif np.sum(m1) > np.sum(m2): cell_seeds[cpy1, cpx1] = 255 From b4aea33f6ad6c7f1fc8af3526378ec53b5e46565 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Thu, 4 Aug 2022 15:15:06 +0100 Subject: [PATCH 11/18] Fix `neutralise_pts_not_under_label_in_frame` --- utils.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/utils.py b/utils.py index 2df05ea..95f88c9 100644 --- a/utils.py +++ b/utils.py @@ -359,19 +359,17 @@ def merge_labels( def neutralise_pts_not_under_label_in_frame( cell_seeds: npt.NDArray[np.uint8], cell_labels: npt.NDArray[np.uint16], -) -> None: +) -> npt.NDArray[np.uint8]: """ Seeds whose label has been eliminated are converted to neutral seeds the idea here is to set seeds not labelled to a value ie invisible to retracking (and to growing, caution!) """ - L = cell_labels - F = cell_seeds - F2 = F - F2[L != 0] = 0 - F[F2 > 252] = 253 - cell_seeds = F + cell_seeds_copy = cell_seeds.copy() + cell_seeds_copy[cell_labels != 0] = 0 + cell_seeds[cell_seeds_copy > 252] = 253 + return cell_seeds def growcellsfromseeds3(): @@ -392,9 +390,9 @@ def _matlab_style_gauss2D(shape: tuple[int, int], sigma: float): fspecial('gaussian', [shape], [sigma]) """ m, n = [(ss - 1) / 2 for ss in shape] - y,x = np.ogrid[-m:m+1,-n:n+1] + y, x = np.ogrid[-m : m + 1, -n : n + 1] h = np.exp(-(x * x + y * y) / (2 * sigma * sigma)) - h[h < np.finfo(h.dtype).eps*h.max()] = 0 + h[h < np.finfo(h.dtype).eps * h.max()] = 0 sumh = h.sum() if sumh != 0: h /= sumh From b3abf83fe66970cd74efe2d19cc9925388e4fc7f Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Thu, 4 Aug 2022 15:46:09 +0100 Subject: [PATCH 12/18] Fix `create_color_boundaries` --- utils.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/utils.py b/utils.py index 95f88c9..9eca1d5 100644 --- a/utils.py +++ b/utils.py @@ -14,15 +14,14 @@ def create_color_boundaries( # given that every cell has a different label we can compute # the boundaries by computing where the gradient changes cell_labels = cell_labels.astype(float) - gx, gy = np.gradient(cell_labels) - cell_outlines = (cell_labels > 0) & ((gx**2 + gy**2) > 0) + gy, gx = np.gradient(cell_labels) + cell_outlines = (cell_labels > 0) & ((gx ** 2 + gy ** 2) > 0) - cell_seeds = cell_seeds > 253 + cell_seeds_mask = cell_seeds > 253 - segmented_image = image.copy() - segmented_image[cell_outlines] = 0 - segmented_image[cell_seeds] = 255 - return segmented_image + image[cell_outlines] = 0 + image[cell_seeds_mask] = 255 + return image def do_initial_seeding( From d54d1f10671718d9569944577d351aa8736aa1cd Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Thu, 4 Aug 2022 18:18:46 +0100 Subject: [PATCH 13/18] Some more porting --- utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils.py b/utils.py index 9eca1d5..a0b9154 100644 --- a/utils.py +++ b/utils.py @@ -15,7 +15,7 @@ def create_color_boundaries( # the boundaries by computing where the gradient changes cell_labels = cell_labels.astype(float) gy, gx = np.gradient(cell_labels) - cell_outlines = (cell_labels > 0) & ((gx ** 2 + gy ** 2) > 0) + cell_outlines = (cell_labels > 0) & ((gx**2 + gy**2) > 0) cell_seeds_mask = cell_seeds > 253 @@ -383,7 +383,7 @@ def calculateCellPositions(): pass -def _matlab_style_gauss2D(shape: tuple[int, int], sigma: float): +def _matlab_style_gauss2D(shape: tuple[int, ...], sigma: float): """ 2D gaussian mask - should give the same result as MATLAB's fspecial('gaussian', [shape], [sigma]) From 056ecd934c78adae2003b91eb065fb76a335c652 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Fri, 5 Aug 2022 17:58:46 +0100 Subject: [PATCH 14/18] Some more porting --- .gitignore | 1 + segment_image.py | 10 +++- utils.py | 120 ++++++++++++++++++++++------------------------- 3 files changed, 65 insertions(+), 66 deletions(-) diff --git a/.gitignore b/.gitignore index 3bbb967..940f94a 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,5 @@ __pycache__ !ProjIm.mat !zerogreen.mat .envrc +*.asv *.mat diff --git a/segment_image.py b/segment_image.py index 97248ef..50d33db 100644 --- a/segment_image.py +++ b/segment_image.py @@ -1,5 +1,6 @@ import numpy as np import numpy.typing as npt +from skimage.morphology import disk from utils import ( create_color_boundaries, @@ -65,16 +66,21 @@ def segment_image( image *= 252 / image.max() image = image.astype(np.uint8) + # structuring element, SE, used for morphological operations + se = disk(2) + do_initial_seeding(image, sigma_1, min_cell_size, threshold, large_cell_size_thres) - merge_seeds_from_labels(image, cell_seeds, cell_labels, merge_criteria, sigma_3) + merge_seeds_from_labels(image, cell_seeds, cell_labels, se, merge_criteria, sigma_3) grow_cells_in_frame(image, cell_seeds, sigma_3) unlabel_poor_seeds_in_frame( - image, cell_seeds, cell_labels, threshold, sigma_3, I_bound_max_pcnt + image, cell_seeds, cell_labels, se, threshold, sigma_3, I_bound_max_pcnt ) neutralise_pts_not_under_label_in_frame(cell_seeds, cell_labels) create_color_boundaries(image, cell_seeds, cell_labels) + + return image, cell_seeds, cell_labels diff --git a/utils.py b/utils.py index a0b9154..7796067 100644 --- a/utils.py +++ b/utils.py @@ -1,6 +1,13 @@ import numpy as np import numpy.typing as npt from skimage.measure import regionprops +from skimage.morphology import binary_dilation, binary_erosion + +SIGMA_BOUND = 0.01 +# TODO: work out what these are for and give better names +TWO_FIVE_FIVE = 255 +TWO_FIVE_THREE = 253 +TWO_FIVE_TWO = 252 def create_color_boundaries( @@ -17,10 +24,10 @@ def create_color_boundaries( gy, gx = np.gradient(cell_labels) cell_outlines = (cell_labels > 0) & ((gx**2 + gy**2) > 0) - cell_seeds_mask = cell_seeds > 253 + cell_seeds_mask = cell_seeds > TWO_FIVE_THREE image[cell_outlines] = 0 - image[cell_seeds_mask] = 255 + image[cell_seeds_mask] = TWO_FIVE_FIVE return image @@ -30,12 +37,12 @@ def do_initial_seeding( min_cell_size: float, threshold: int, large_cell_size_thres: float, -) -> None: +) -> npt.NDArray[np.uint8]: """ Find the initial cell seeds """ # Create gaussian filter - if sigma_1 > 0.01: + if sigma_1 > SIGMA_BOUND: f1 = _matlab_style_gauss2D(image.shape, sigma_1) # Gaussian smoothing for the segmentation of individual cells @@ -45,7 +52,7 @@ def do_initial_seeding( else: smoothed_image = image.astype(float) - smoothed_image /= smoothed_image.max() * 252 + smoothed_image /= smoothed_image.max() * TWO_FIVE_TWO # Use external c-code to find initial seeds initial_labelling = findcellsfromregiongrowing( @@ -65,22 +72,22 @@ def do_initial_seeding( centroids = np.round(calculateCellPositions(smoothed_image, cell_labels, False)) centroids = centroids[~np.isnan(centroids.T[0])] for n in range(len(centroids)): - smoothed_image[centroids[n, 1], centroids[n, 0]] = 255 + smoothed_image[centroids[n, 1], centroids[n, 0]] = TWO_FIVE_FIVE # cell_seeds contains the position of the true cell center. - cell_seeds = smoothed_image.astype(np.uint8) + return smoothed_image.astype(np.uint8) def grow_cells_in_frame( image: npt.NDArray[np.uint8], cell_seeds: npt.NDArray[np.uint8], sigma_3: float -) -> None: +) -> npt.NDArray[np.uint16]: """ Growing cells from seeds TODO: add paramters in Name description! """ # find labels - bw = (cell_seeds > 252).astype(float) + bw = (cell_seeds > TWO_FIVE_TWO).astype(float) - if sigma_3 > 0.01: + if sigma_3 > SIGMA_BOUND: f1 = _matlab_style_gauss2D(image.shape, sigma_3) smoothed_image = np.fft.fftshift( np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1)) @@ -89,14 +96,15 @@ def grow_cells_in_frame( smoothed_image = image.astype(float) # mark labels on image - image_with_seeds = (smoothed_image).astype(float) * (1 - bw) + 255 * bw - cell_labels = growcellsfromseeds3(image_with_seeds, 253).astype(np.uint16) + image_with_seeds = (smoothed_image).astype(float) * (1 - bw) + TWO_FIVE_FIVE * bw + return growcellsfromseeds3(image_with_seeds, TWO_FIVE_THREE).astype(np.uint16) def unlabel_poor_seeds_in_frame( image: npt.NDArray[np.uint8], cell_seeds: npt.NDArray[np.uint8], cell_labels: npt.NDArray[np.uint16], + se: npt.NDArray[np.uint8], threshold: int, sigma_3: float, I_bound_max_pcnt: float, @@ -106,7 +114,7 @@ def unlabel_poor_seeds_in_frame( """ L = cell_labels - if sigma_3 > 0.01: + if sigma_3 > SIGMA_BOUND: f1 = _matlab_style_gauss2D(image.shape, sigma_3) smoothed_image = np.fft.fftshift( np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1)) @@ -122,35 +130,30 @@ def unlabel_poor_seeds_in_frame( for c in range(len(label_list)): mask = L == label_list[c] - cpy, cpx = np.argwhere(mask > 0) + cpy, cpx = np.argwhere(mask > 0).T # find region of that label - minx = np.min(cpx) - maxx = np.max(cpx) - miny = np.min(cpy) - maxy = np.max(cpy) - minx = np.max(minx - 5, 1) - miny = np.max(miny - 5, 1) - maxx = np.min(maxx + 5, image.shape[1]) - maxy = np.min(maxy + 5, image.shape[0]) + minx = max(cpx.min() - 5, 0) + miny = max(cpy.min() - 5, 0) + maxx = min(cpx.max() + 5, image.shape[1] - 1) + maxy = min(cpy.max() + 5, image.shape[0] - 1) # reduced to region of the boundary - reduced_mask = mask[miny:maxy, minx:maxx] - reduced_image = smoothed_image[miny:maxy, minx:maxx] - dilated_mask = imdilate(reduced_mask, se) - eroded_mask = imerode(reduced_mask, se) - boundary_mask = dilated_mask - eroded_mask + reduced_mask = mask[miny : maxy + 1, minx : maxx + 1] + reduced_image = smoothed_image[miny : maxy + 1, minx : maxx + 1] + dilated_mask = binary_dilation(reduced_mask, se) + eroded_mask = binary_erosion(reduced_mask, se) + boundary_mask = dilated_mask ^ eroded_mask boundary_intensities = reduced_image[boundary_mask > 0] H = reduced_image[boundary_mask > 0] - IEr = reduced_image[eroded_mask > 0] - I_bound = np.mean(boundary_intensities) + I_bound = (boundary_intensities).mean() IBounds[c] = I_bound # cell seed information is retrieved as comparison - F2 = cell_seeds + F2 = cell_seeds.copy() F2[~mask] = 0 - cpy, cpx = np.argwhere(F2 > 252) - ICentre = smoothed_image[cpy, cpx] + cpy, cpx = np.argwhere(F2 > TWO_FIVE_TWO).T + ICentre = smoothed_image[cpy, cpx][0] - I_bound_max = 255 * I_bound_max_pcnt + I_bound_max = TWO_FIVE_FIVE * I_bound_max_pcnt # Figure out which conditions make the label invalid # 1. I_bound_max, gives the Lower bound to the mean intensity @@ -172,16 +175,16 @@ def unlabel_poor_seeds_in_frame( # record the removal decisions if first_condition: - decisions[1] = decisions[1] + 1 + decisions[1] += 1 elif second_condition: - decisions[2] = decisions[2] + 1 + decisions[2] += 1 elif third_condition: - decisions[3] = decisions[3] + 1 + decisions[3] += 1 elif fourth_condition: - decisions[4] = decisions[4] + 1 + decisions[4] += 1 else: # should not happen - decisions[5] = decisions[5] + 1 + decisions[5] += 1 def delabel_very_large_areas( @@ -207,6 +210,7 @@ def merge_seeds_from_labels( image: npt.NDArray[np.uint8], cell_seeds: npt.NDArray[np.uint8], cell_labels: npt.NDArray[np.uint16], + se: npt.NDArray[np.uint8], merge_criteria: float, sigma_3: float, ) -> None: @@ -214,7 +218,7 @@ def merge_seeds_from_labels( Remove initial cell regions which touch & whose boundary is insufficient """ # smoothing - if sigma_3 > 0.01: + if sigma_3 > SIGMA_BOUND: f1 = _matlab_style_gauss2D(image.shape, sigma_3) smoothed_image = np.fft.fftshift( np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1)) @@ -226,7 +230,7 @@ def merge_seeds_from_labels( label_list = label_list[label_list != 0] c = 1 - merge_intensity_distro = [] + merge_intensity_distro: list = [] merge_decisions = 0 # loop over labels @@ -234,7 +238,7 @@ def merge_seeds_from_labels( label_mask = cell_labels == label_list[c] label = label_list[c] - cpy, cpx = np.argwhere(label_mask > 0) + cpy, cpx = np.argwhere(label_mask > 0).T # find region of that label minx = np.min(cpx) @@ -252,15 +256,15 @@ def merge_seeds_from_labels( reduced_labels = cell_labels[miny:maxy, minx:maxx] # now find boundaries - dilated_mask = imdilate(reduced_label_mask, se) - eroded_mask = imerode(reduced_label_mask, se) + dilated_mask = binary_dilation(reduced_label_mask, se) + eroded_mask = binary_erosion(reduced_label_mask, se) border_mask = dilated_mask - eroded_mask border_intensities = reduced_image[border_mask > 0] central_intensity = reduced_image[eroded_mask > 0] F2 = cell_seeds F2[~label_mask] = 0 - cpy, cpx = np.argwhere(F2 > 253) + cpy, cpx = np.argwhere(F2 > TWO_FIVE_THREE).T ICentre = smoothed_image[cpy, cpx] background_std = np.std(central_intensity.astype(float)) @@ -269,13 +273,13 @@ def merge_seeds_from_labels( neighbour_labels = np.unique(reduced_labels[dilated_mask > 0]) neighbour_labels = neighbour_labels[neighbour_labels != label] - low_intensity_ratios = [] + low_intensity_ratios: list = [] for i in range(len(neighbour_labels)): neighb_label = neighbour_labels(i) neighbor_border = dilated_mask # slice of neighbour around cell neighbor_border[reduced_labels != neighb_label] = 0 - cell_border = imdilate(neighbor_border, se) + cell_border = binary_dilation(neighbor_border, se) # slice of cell closest to neighbour cell_border[reduced_labels != label] = 0 @@ -336,8 +340,8 @@ def merge_labels( Il1[~m1] = 0 Il2 = Il Il2[~m2] = 0 - cpy1, cpx1 = np.argwhere(Il1 > 253) - cpy2, cpx2 = np.argwhere(Il2 > 253) + cpy1, cpx1 = np.argwhere(Il1 > TWO_FIVE_THREE).T + cpy2, cpx2 = np.argwhere(Il2 > TWO_FIVE_THREE).T cpx = round((cpx1 + cpx2) / 2) cpy = round((cpy1 + cpy2) / 2) @@ -345,11 +349,11 @@ def merge_labels( cell_seeds[cpy1, cpx1] = 20 cell_seeds[cpy2, cpx2] = 20 if (cell_labels[cpy, cpx] == l1) | (cell_labels[cpy, cpx] == l2): - cell_seeds[cpy, cpx] = 255 + cell_seeds[cpy, cpx] = TWO_FIVE_FIVE elif np.sum(m1) > np.sum(m2): - cell_seeds[cpy1, cpx1] = 255 + cell_seeds[cpy1, cpx1] = TWO_FIVE_FIVE else: - cell_seeds[cpy2, cpx2] = 255 + cell_seeds[cpy2, cpx2] = TWO_FIVE_FIVE Cl[m2] = l1 cell_labels = Cl @@ -367,22 +371,10 @@ def neutralise_pts_not_under_label_in_frame( """ cell_seeds_copy = cell_seeds.copy() cell_seeds_copy[cell_labels != 0] = 0 - cell_seeds[cell_seeds_copy > 252] = 253 + cell_seeds[cell_seeds_copy > TWO_FIVE_TWO] = TWO_FIVE_THREE return cell_seeds -def growcellsfromseeds3(): - pass - - -def findcellsfromregiongrowing(): - pass - - -def calculateCellPositions(): - pass - - def _matlab_style_gauss2D(shape: tuple[int, ...], sigma: float): """ 2D gaussian mask - should give the same result as MATLAB's From 2489455bd1abdc7113b3b071b53540bac9448020 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Fri, 5 Aug 2022 17:59:14 +0100 Subject: [PATCH 15/18] Demo of how to port --- example_segmentation.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/example_segmentation.py b/example_segmentation.py index 22b5857..5b00dd3 100644 --- a/example_segmentation.py +++ b/example_segmentation.py @@ -1,8 +1,10 @@ from pathlib import Path from scipy import io as sio +from skimage.morphology import disk from segment_image import segment_image +from utils import unlabel_poor_seeds_in_frame _file_location = Path(__file__).resolve() _matfile = _file_location.parent / "ProjIm.mat" @@ -32,4 +34,8 @@ def main() -> None: if __name__ == "__main__": - main() + se = disk(2) + Im = sio.loadmat(_file_location.parent / "Im.mat")["Im"] + CellSeeds = sio.loadmat(_file_location.parent / "CellSeeds.mat")["CellSeeds"] + CellLabels = sio.loadmat(_file_location.parent / "CellLabels.mat")["CellLabels"] + unlabel_poor_seeds_in_frame(Im, CellSeeds, CellLabels, se, 20, 0.5, 0.1), From fe623d1c66888d0b5fbd55453171ba1c11f166b3 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Fri, 12 Aug 2022 16:05:47 +0100 Subject: [PATCH 16/18] Start porting calculate cell positions --- calculate_cell_positions.py | 64 +++++++++++++++++++++++++++++++++++++ utils.py | 4 ++- 2 files changed, 67 insertions(+), 1 deletion(-) create mode 100644 calculate_cell_positions.py diff --git a/calculate_cell_positions.py b/calculate_cell_positions.py new file mode 100644 index 0000000..1c3b036 --- /dev/null +++ b/calculate_cell_positions.py @@ -0,0 +1,64 @@ +import numpy as np +import numpy.typing as npt + +TWO_FIVE_FIVE = 255 + + +def calculate_cell_positions( + image: npt.NDArray[np.float64], labelled_cells: npt.NDArray[np.uint16], type: bool +) -> None: + """ + This function calculates the min intensity position of each + labelled cell or the centroid position of each labelled region + """ + no_cells = labelled_cells.max() + + pos = np.zeros((no_cells, 2)) + + I2 = TWO_FIVE_FIVE - image + + for n in range(no_cells): + + sy, sx = np.argwhere(labelled_cells == n).T + + if type: + _place_at_lowest_int(pos, I2, sx, sy, n) + + else: + _place_at_centroid(pos, image, sx, sy, n) + + if ~np.isnan(pos[n, 1]) and labelled_cells[pos[n, 1], pos[n, 0]] != n: + # every so often the centroid is actually not in the label! + _place_at_lowest_int(pos, image, sx, sy, n) + + +def _place_at_centroid(pos, image, sx, sy, n) -> None: + """ + calculating the centroid from intensities + """ + sum_x = 0 + sum_y = 0 + sum_I = 0 + + for m in range(len(sy)): + sum_x += sx[m] * image[sy[m], sx[m]] + sum_y += sy[m] * image[sy[m], sx[m]] + sum_I += image[sy[m], sx[m]] + + pos[n, 1] = round(sum_y / sum_I) + pos[n, 0] = round(sum_x / sum_I) + + +def _place_at_lowest_int(pos, image, sx, sy, n) -> None: + """ + looking for the lowest intensity + """ + val = TWO_FIVE_FIVE + + for m in range(len(sy)): + + if image[sy[m], sx[m]] < val: + + val = image[sy[m], sx[m]] + pos[n, 2] = sy[m] + pos[n, 1] = sx[m] diff --git a/utils.py b/utils.py index 7796067..eef0ed5 100644 --- a/utils.py +++ b/utils.py @@ -3,6 +3,8 @@ from skimage.measure import regionprops from skimage.morphology import binary_dilation, binary_erosion +from calculate_cell_positions import calculate_cell_positions + SIGMA_BOUND = 0.01 # TODO: work out what these are for and give better names TWO_FIVE_FIVE = 255 @@ -69,7 +71,7 @@ def do_initial_seeding( delabel_very_large_areas(cell_labels, large_cell_size_thres) # Use true centre of cells as labels - centroids = np.round(calculateCellPositions(smoothed_image, cell_labels, False)) + centroids = np.round(calculate_cell_positions(smoothed_image, cell_labels, False)) centroids = centroids[~np.isnan(centroids.T[0])] for n in range(len(centroids)): smoothed_image[centroids[n, 1], centroids[n, 0]] = TWO_FIVE_FIVE From 8e78e56d2f9e381d9afa69683c23f65f9e617144 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Fri, 12 Aug 2022 17:26:15 +0100 Subject: [PATCH 17/18] Add return type --- calculate_cell_positions.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/calculate_cell_positions.py b/calculate_cell_positions.py index 1c3b036..e40eafa 100644 --- a/calculate_cell_positions.py +++ b/calculate_cell_positions.py @@ -31,6 +31,8 @@ def calculate_cell_positions( # every so often the centroid is actually not in the label! _place_at_lowest_int(pos, image, sx, sy, n) + return pos + def _place_at_centroid(pos, image, sx, sy, n) -> None: """ @@ -47,6 +49,7 @@ def _place_at_centroid(pos, image, sx, sy, n) -> None: pos[n, 1] = round(sum_y / sum_I) pos[n, 0] = round(sum_x / sum_I) + return pos def _place_at_lowest_int(pos, image, sx, sy, n) -> None: From 67e201fcb6e872726ccad56eaa30def800e294d2 Mon Sep 17 00:00:00 2001 From: paddyroddy Date: Mon, 15 Aug 2022 12:25:45 +0100 Subject: [PATCH 18/18] Rename `I2` --- calculate_cell_positions.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/calculate_cell_positions.py b/calculate_cell_positions.py index e40eafa..3ec14bf 100644 --- a/calculate_cell_positions.py +++ b/calculate_cell_positions.py @@ -15,14 +15,14 @@ def calculate_cell_positions( pos = np.zeros((no_cells, 2)) - I2 = TWO_FIVE_FIVE - image + subtracted_image = TWO_FIVE_FIVE - image for n in range(no_cells): sy, sx = np.argwhere(labelled_cells == n).T if type: - _place_at_lowest_int(pos, I2, sx, sy, n) + _place_at_lowest_int(pos, subtracted_image, sx, sy, n) else: _place_at_centroid(pos, image, sx, sy, n)