|
| 1 | +function [resMat,resMatNames,numImPts] = getTifBoundary(coords,img,object,imgName,distThresh,fibKey,endLength,fibProcMeth,distMini) |
| 2 | + |
| 3 | +% getTifBoundary.m - This function takes the coordinates from the boundary file, associates them with curvelets, and produces relative angle measures. |
| 4 | +% |
| 5 | +% Inputs: |
| 6 | +% coords - the locations of the endpoints of each line segment making up the boundary |
| 7 | +% img - the image being measured |
| 8 | +% object - a struct containing the center and angle of each measured curvelet, generated by the newCurv function |
| 9 | +% distThresh - number of pixels from boundary we should evaluate curvelets |
| 10 | +% boundaryImg - tif file with boundary outlines, must be a mask file |
| 11 | +% fibKey - list indicating the beginning of each new fiber in the object struct, allows for fiber level processing |
| 12 | +% distMini - minimum distrance to the boundary, to get rid of the fiber on or very close to the boundary |
| 13 | +% |
| 14 | +% Output: |
| 15 | +% measAngs - all relative angle measurements, not filtered by distance |
| 16 | +% measDist - all distances between the curvelets and the boundary points, not filtered |
| 17 | +% inCurvsFlag - curvelets that are considered |
| 18 | +% outCurvsFlag - curvelets that are not considered |
| 19 | +% measBndry = points on the boundary that are associated with each curvelet |
| 20 | +% inDist = distance between boundary and curvelet for each curvelet considered |
| 21 | +% numImPts = number of points in the image that are less than distThresh from boundary |
| 22 | +% insCt = number of curvelets inside an epithelial region |
| 23 | +% |
| 24 | +% |
| 25 | +% By Jeremy Bredfeldt, LOCI, Morgridge Institute for Research, 2013 |
| 26 | + |
| 27 | +%Note: a "curv" could be a curvelet or a fiber segment, depending on if CT or FIRE is used |
| 28 | + |
| 29 | +imHeight = size(img,1); |
| 30 | +imWidth = size(img,2); |
| 31 | +sz = [imHeight,imWidth]; |
| 32 | + |
| 33 | +% figure(600); |
| 34 | +% imshow(img); |
| 35 | + |
| 36 | +% figure(500); |
| 37 | +% hold on; |
| 38 | +% for k = 1:length(coords) |
| 39 | +% boundary = coords{k}; |
| 40 | +% plot(boundary(:,2), boundary(:,1), 'y', 'LineWidth', 2); |
| 41 | +% end |
| 42 | + |
| 43 | +%collect all fiber points |
| 44 | +allCenterPoints = vertcat(object.center); |
| 45 | +%collect all boundary points |
| 46 | +% coords = vertcat(coords{2:end,1}); |
| 47 | +coords = vertcat(coords{1:end,1}); |
| 48 | + |
| 49 | +%collect all region points |
| 50 | +linIdx = sub2ind(sz, allCenterPoints(:,1), allCenterPoints(:,2)); |
| 51 | + |
| 52 | +[idx_dist,dist] = knnsearch(coords,allCenterPoints); %closest point to a boundary |
| 53 | +reg_dist = img(linIdx); |
| 54 | + |
| 55 | + |
| 56 | +%YL: test the boundary association |
| 57 | +% figure(1002);clf;set(gcf,'pos',[200 300 imWidth imHeight ]); |
| 58 | +% plot(coords(:,2),coords(:,1),'k.'); axis ij |
| 59 | +% axis([1 imWidth 1 imHeight ]);hold on |
| 60 | + |
| 61 | +%[idx_reg,reg_dist] = knnsearch([reg_col,reg_row],allCenterPoints); %closest point to a filled in region |
| 62 | + |
| 63 | +%Make a list of points in the image (points scattered throughout the image) |
| 64 | +C = floor(imWidth/20); %use at least 20 per row in the image, this is done to speed this process up |
| 65 | +[I, J] = ind2sub(size(img),1:C:imHeight*imWidth); |
| 66 | +allImPoints = [I; J]'; |
| 67 | +%Get list of image points that are a certain distance from the boundary |
| 68 | +[~,dist_im] = knnsearch(coords(1:3:end,:),allImPoints); %returns nearest dist to each point in image |
| 69 | +%threshold distance |
| 70 | +if isempty(distMini) % keep all the fibers within the distance threshold |
| 71 | + inIm = dist_im <= distThresh; |
| 72 | +else % get rid of fibers on or very close to the boundary |
| 73 | + inIm = (dist_im <= distThresh & dist_im > distMini); |
| 74 | +end |
| 75 | +%count number of points |
| 76 | +inPts = allImPoints(inIm); |
| 77 | +numImPts = length(inPts)*C; |
| 78 | +% numImPts = 0; |
| 79 | + |
| 80 | + |
| 81 | +%process all curvs, at this point |
| 82 | +curvsLen = length(object); |
| 83 | +nbDist = nan(1,curvsLen); %nearest boundary distance |
| 84 | +nrDist = nan(1,curvsLen); %nearest region distance |
| 85 | +nbAng = nan(1,curvsLen); %nearest boundary relative angle |
| 86 | +epDist = nan(1,curvsLen); %distance to extension point intersection |
| 87 | +epAng = nan(1,curvsLen); %relative angle of extension point intersection |
| 88 | +measBndry = nan(curvsLen,2); |
| 89 | +inCurvsFlag = ~logical(1:curvsLen); |
| 90 | +outCurvsFlag = ~logical(1:curvsLen); |
| 91 | +% %ROI properties |
| 92 | +% bwROI.coords = coords; |
| 93 | +% bwROI.imageWidth = imWidth; |
| 94 | +% bwROI.imageHeight = imHeight; |
| 95 | +if isempty(distMini) |
| 96 | + for i = 1:curvsLen |
| 97 | + %-- inside region? |
| 98 | + nrDist(i) = reg_dist(i)==255|reg_dist(i)== 1; %YL: mask can be 1-0(matlab lab) or 255-0(ImageJ) |
| 99 | + %-- distance to nearest epithelial boundary |
| 100 | + nbDist(i) = dist(i); |
| 101 | + %-- relative angle at nearest boundary point |
| 102 | + if dist(i) <= distThresh |
| 103 | + [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 |
| 104 | + % %% |
| 105 | + % bwROI.index2object = idx_dist(i); |
| 106 | + % fiberobject.center = object(i).center; |
| 107 | + % fiberobject.angle = object(i).angle; |
| 108 | + % angleOption = 0; % caclulate all angles |
| 109 | + % [relativeAngles, ROImeasurements] = getRelativeangles(bwROI,fiberobject,angleOption); |
| 110 | + % fprintf(' original relative angle is %3.2f \n re-calculated relative angle is %3.2f \n ',nbAng(i),relativeAngles.angle2boundaryEdge); |
| 111 | + % fprintf(' relative angle to ROI orientation is %3.2f \n',relativeAngles.angle2boundaryCenter); |
| 112 | + % fprintf('relative angle to ROI-fiber centers line is %3.2f \n',relativeAngles.angle2centersLine); |
| 113 | + %% |
| 114 | + else |
| 115 | + nbAng(i) = nan; % if out of the region |
| 116 | + bPt = nan(1,2); % if out of the region |
| 117 | + end |
| 118 | + %-- extension point features |
| 119 | + [lineCurv orthoCurv] = getPointsOnLine(object(i),imWidth,imHeight,distThresh); |
| 120 | + [intLine, iLa, iLb] = intersect([lineCurv(:,2) lineCurv(:,1)],coords,'rows'); |
| 121 | + if (~isempty(intLine)) |
| 122 | + %get the closest distance from the curvelet center to the |
| 123 | + %intersection (get rid of the farther one(s)) |
| 124 | + [idxLineDist, lineDist] = knnsearch(intLine,object(i).center); |
| 125 | + boundaryPtIdx = iLb(idxLineDist); |
| 126 | + %%tentatively turn the extension feature off |
| 127 | + % %-- extension point distance |
| 128 | + % epDist(i) = lineDist; |
| 129 | + % %-- extension point angle |
| 130 | + % [epAng(i) bPt1] = GetRelAng([coords(:,2),coords(:,1)],boundaryPtIdx,object(i).angle,imHeight,imWidth,i); |
| 131 | + else |
| 132 | + epDist(i) = nan;% no intersection exists |
| 133 | + epAng(i) = nan; % no angle exists |
| 134 | + bPt1 = nan(1,2); % if no intersection set boundary to be [NaN NaN] |
| 135 | + end |
| 136 | + measBndry(i,:) = bPt; % nearest boundary |
| 137 | + end %of for loop |
| 138 | +else % get rid of fibers on or within the minimum distance to the boundary |
| 139 | + for i = 1:curvsLen |
| 140 | + %-- inside region? |
| 141 | + nrDist(i) = reg_dist(i)==255|reg_dist(i)== 1; %YL: mask can be 1-0(matlab lab) or 255-0(ImageJ) |
| 142 | + %-- distance to nearest epithelial boundary |
| 143 | + nbDist(i) = dist(i); |
| 144 | + %-- relative angle at nearest boundary point |
| 145 | + if (dist(i) <= distThresh & dist(i) > distMini) |
| 146 | + [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 |
| 147 | + else |
| 148 | + nbAng(i) = nan; % if out of the region |
| 149 | + bPt = nan(1,2); % if out of the region |
| 150 | + end |
| 151 | + %-- extension point features |
| 152 | + [lineCurv orthoCurv] = getPointsOnLine(object(i),imWidth,imHeight,distThresh); |
| 153 | + [intLine, iLa, iLb] = intersect([lineCurv(:,2) lineCurv(:,1)],coords,'rows'); |
| 154 | + if (~isempty(intLine)) |
| 155 | + %get the closest distance from the curvelet center to the |
| 156 | + %intersection (get rid of the farther one(s)) |
| 157 | + [idxLineDist, lineDist] = knnsearch(intLine,object(i).center); |
| 158 | + boundaryPtIdx = iLb(idxLineDist); |
| 159 | + %%tentatively turn the extension feature off |
| 160 | + % %-- extension point distance |
| 161 | + % epDist(i) = lineDist; |
| 162 | + % %-- extension point angle |
| 163 | + % [epAng(i) bPt1] = GetRelAng([coords(:,2),coords(:,1)],boundaryPtIdx,object(i).angle,imHeight,imWidth,i); |
| 164 | + else |
| 165 | + epDist(i) = nan;% no intersection exists |
| 166 | + epAng(i) = nan; % no angle exists |
| 167 | + bPt1 = nan(1,2); % if no intersection set boundary to be [NaN NaN] |
| 168 | + end |
| 169 | + measBndry(i,:) = bPt; % nearest boundary |
| 170 | + end %of for loopend |
| 171 | +end |
| 172 | + |
| 173 | + |
| 174 | + |
| 175 | +resMat = [nbDist', ... %nearest dist to a boundary |
| 176 | + nrDist', ... %flag, 0 for outside boundary, 1 for inside boundary |
| 177 | + nbAng', ... %nearest relative boundary angle |
| 178 | + epDist', ... %extension point distance |
| 179 | + epAng', ... %extension point relative boundary angle |
| 180 | + measBndry]; %list of boundary points associated with fibers |
| 181 | +resMatNames = { |
| 182 | + 'nearestBoundDist', ... |
| 183 | + 'nearestRegionDist', ... |
| 184 | + 'nearestBoundAng', ... |
| 185 | + 'extensionPointDist', ... |
| 186 | + 'extensionPointAng', ... |
| 187 | + 'bndryPtRow', ... |
| 188 | + 'bndryPtCol'}; |
| 189 | + |
| 190 | +end %of main function |
| 191 | + |
| 192 | +function [relAng, boundaryPt] = GetRelAng(coords,idx,fibAng,imHeight,imWidth,fnum) |
| 193 | +boundaryAngle = FindOutlineSlope([coords(:,2),coords(:,1)],idx); |
| 194 | +boundaryPt = coords(idx,:); |
| 195 | +if (boundaryPt(1) == 1 || boundaryPt(2) == 1 || boundaryPt(1) == imHeight || boundaryPt(2) == imWidth) |
| 196 | + %don't count fiber if boundary point is along edge of image |
| 197 | + tempAng = 0; |
| 198 | +else |
| 199 | + %--compute relative angle here-- |
| 200 | + %There is a 90 degree phase shift in fibAng and boundaryAngle due to image orientation issues in Matlab. |
| 201 | + % -therefore no need to invert (ie. 1-X) circ_r here. |
| 202 | + tempAng = circ_r([fibAng*2*pi/180; boundaryAngle*2*pi/180]); |
| 203 | + tempAng = 180*asin(tempAng)/pi; |
| 204 | + % %YL debug the NaN angle |
| 205 | + % if isnan(tempAng) |
| 206 | + % figure(1002),plot(coords(idx,1),coords(idx,2),'ro','MarkerSize',10) |
| 207 | + % text(coords(idx,1),coords(idx,2),sprintf('%d',fnum)); |
| 208 | + % disp(sprintf('fiber %d relative angle is Nan, fibAng = %f, boundaryAngle = %f, idx_dist = %d',fnum,fibAng,boundaryAngle,idx)) |
| 209 | + % % pause(3) |
| 210 | + % end |
| 211 | +end |
| 212 | +relAng = tempAng; |
| 213 | +end |
| 214 | + |
| 215 | +function [lineCurv orthoCurv] = getPointsOnLine(object,imWidth,imHeight,boxSz) |
| 216 | +center = object.center; |
| 217 | +angle = object.angle; |
| 218 | +slope = -tand(angle); |
| 219 | +%orthoSlope = -tand(angle + 90); %changed from tand(obj.ang) to -tand(obj.ang + 90) 10/12 JB |
| 220 | +intercept = center(1) - (slope)*center(2); |
| 221 | +%orthoIntercept = center(1) - (orthoSlope)*center(2); |
| 222 | + |
| 223 | +%[p1 p2] = getBoxInt(slope, intercept, imWidth, imHeight, center, boxSz); |
| 224 | +if isinf(slope) |
| 225 | + dist_y = 0; |
| 226 | + dist_x = boxSz; |
| 227 | +else |
| 228 | + dist_y = boxSz/sqrt(1+slope*slope); |
| 229 | + dist_x = dist_y*slope; |
| 230 | +end |
| 231 | +p1 = [center(2) - dist_y, center(1) - dist_x]; |
| 232 | +p2 = [center(2) + dist_y, center(1) + dist_x]; |
| 233 | +[lineCurv ~] = GetSegPixels(p1,p2); |
| 234 | + |
| 235 | +%Not using the orthogonal slope for anything now |
| 236 | +%[p1 p2] = getIntImgEdge(orthoSlope, orthoIntercept, imWidth, imHeight, center); |
| 237 | +%[orthoCurv, ~] = GetSegPixels(p1,p2); |
| 238 | +orthoCurv = []; |
| 239 | + |
| 240 | +end |
| 241 | + |
0 commit comments