bwdistsc.m 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483
  1. function D=bwdistsc(bw,aspect)
  2. % D=BWDISTSC(BW,ASPECT)
  3. % BWDISTSC computes Euclidean distance transform of the binary 3D image
  4. % BW. For each pixel in BW, the distance transform assignes a number
  5. % that is the distance between that pixel and the nearest nonzero pixel
  6. % of BW. BW may be a single 2D image, 3D array or a cell array of
  7. % 2D slices. ASPECT is 3-component vector defining aspect ratio in
  8. % the dataset BW. If ASPECT is not specified, isotropic aspect
  9. % ratio [1 1 1] is assumed.
  10. %
  11. % BWDISTSC uses fast optimized scanning algorithm and cell-arrays to
  12. % represent internal data, so that it is less demanding to physical
  13. % memory. In many cases BWDISTSC is actually faster than MATLAB's
  14. % optimized kd-tree algorithm used for Euclidean distance
  15. % transform in 3D.
  16. %
  17. % BWDISTSC tries to use MATLAB bwdist for 2D scans if possible, which
  18. % is significantly faster. Otherwise BWDISTSC uses internal algorithm
  19. % to perform 2D scans.
  20. %
  21. % Yuriy Mishchenko JFRC HHMI Chklovskii Lab JUL 2007
  22. % This code is free for use or modifications, just please give credit
  23. % where appropriate. And if you modify code or fix bugs, please drop
  24. % me a message at gmyuriy@hotmail.com.
  25. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  26. % Scan algorithms below use following Lema: %
  27. % LEMA: let F(X,z) be lower envelope of a family of parabola: %
  28. % F(X,z)=min_{i} [G_i(X)+(z-k_i)^2]; %
  29. % and let H_k(X,z)=A(X)+(z-k)^2 be a parabola. %
  30. % Then for H_k(X,z)==F(X,z) at each X there exist at most %
  31. % two solutions k1<k2 such that H_k12(X,z)=F(X,z), and %
  32. % H_k(X,z)<F(X,z) is restricted to at most k1<k2. %
  33. % Here X is any-dimensional coordinate. %
  34. % %
  35. % Thus, simply scan away from any z such that H_k(X,z)<F(X,z) %
  36. % in either direction as long as H_k(X,z)<F(X,z) and update %
  37. % F(X,z). Note that need to properly choose starting point; %
  38. % starting point is any z such that H_k(X,z)<F(X,z); z==k is %
  39. % usually, but not always the starting point!!! %
  40. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  41. % parse inputs
  42. if(nargin<2 || isempty(aspect)) aspect=[1 1 1]; end
  43. % determine geometry of data
  44. if(iscell(bw)) shape=[size(bw{1}),length(bw)]; else shape=size(bw); end
  45. % fix to handle 1D & 2D images
  46. if(length(shape)<3) shape(length(shape)+1:3)=1; end
  47. if(length(aspect)<3) aspect(length(aspect)+1:3)=1; end
  48. % allocate space
  49. D=cell(1,shape(3)); for k=1:shape(3) D{k}=zeros(shape(1:2)); end
  50. %%%%%%%%%%%%% scan along XY %%%%%%%%%%%%%%%%
  51. for k=1:shape(3)
  52. if(iscell(bw)) bwXY=bw{k}; else bwXY=bw(:,:,k); end
  53. % initialize arrays
  54. DXY=zeros(shape(1:2));
  55. D1=zeros(shape(1:2));
  56. DK=zeros(shape(1:2));
  57. % if can, use 2D bwdist from image processing toolbox
  58. if(exist('bwdist') && aspect(1)==aspect(2))
  59. D1=aspect(1)^2*bwdist(bwXY).^2;
  60. else % if not, use full XY-scan
  61. %%%%%%%%%%%%%%% X-SCAN %%%%%%%%%%%%%%%
  62. % reference nearest bwXY "on"-pixel in x direction downward:
  63. % scan bottow-up, copy x-reference from previous row unless
  64. % there is bwXY "on"-pixel in that point in current row
  65. xlower=repmat(Inf,shape(1:2));
  66. xlower(1,find(bwXY(1,:)))=1; % fill in first row
  67. for i=2:shape(1)
  68. xlower(i,:)=xlower(i-1,:); % copy previous row
  69. xlower(i,find(bwXY(i,:)))=i;% unless there is pixel
  70. end
  71. % reference nearest bwXY "on"-pixel in x direction upward:
  72. xupper=repmat(Inf,shape(1:2));
  73. xupper(end,find(bwXY(end,:)))=shape(1);
  74. for i=shape(1)-1:-1:1
  75. xupper(i,:)=xupper(i+1,:);
  76. xupper(i,find(bwXY(i,:)))=i;
  77. end
  78. % find points for which distance needs to be updated
  79. idx=find(~bwXY); [x,y]=ind2sub(shape(1:2),idx);
  80. % set distance as the shortest to upward or to downward
  81. DXY(idx)=aspect(1)^2*min((x-xlower(idx)).^2,(x-xupper(idx)).^2);
  82. %%%%%%%%%%%%%%% Y-SCAN %%%%%%%%%%%%%%%
  83. % this will be the envelop
  84. % envelop is initialized at Inf to ensure single scan direction,
  85. % otherwise may end up in infinite loop when trying to find
  86. % starting point
  87. D1=repmat(Inf,shape(1:2));
  88. % these will be the references to parabolas defining the envelop
  89. DK=repmat(Inf,shape(1:2));
  90. % starting points
  91. i0=zeros(shape(1),1);
  92. % convenience x-coords array
  93. x=(1:shape(1))';
  94. for i=1:shape(2)
  95. % need to select starting point for each X:
  96. % * starting point should be below current envelop
  97. % * i0==i is not necessarily a starting point
  98. % * there is at most one starting point
  99. % * there may be no starting point
  100. % i0 is the starting points for each X: i0(X) is the first
  101. % y-index such that parabola from line i is below the envelop
  102. % first guess is the current y-line
  103. i0(:)=i;
  104. % some auxiliary datasets
  105. d0=DXY(:,i);
  106. % L0 indicates for which X starting point had been fixed
  107. L0=isinf(d0) | (d0==0);
  108. while(~isempty(find(~L0,1)))
  109. % reference starting points in DXY
  110. idx=sub2ind(shape(1:2),x(~L0),i0(~L0));
  111. % reduce out trivial points (DXY==0)
  112. L=(DXY(idx)==0);
  113. L0(~L0)=L;
  114. idx=idx(~L);
  115. if(isempty(idx)) continue; end
  116. % these are current best parabolas for starting points
  117. ik=DK(idx);
  118. % these are new values from parabola from line #i
  119. dtmp=d0(~L0)+aspect(2)^2*(i0(~L0)-i).^2;
  120. % these starting points are OK - below the envelop
  121. L=D1(idx)>dtmp; D1(idx(L))=dtmp(L);
  122. % points which are still above the envelop but ik==i0,
  123. % will not get any better, so fix them as well
  124. L=L | (ik==i0(~L0));
  125. % all other points are not OK, need new starting point:
  126. % starting point should be at least below parabola
  127. % beating us at current choice of i0
  128. % solve quadratic equation to find where this happens
  129. ik=(ik-i);
  130. di=(D1(idx(~L))-dtmp(~L))./ik(~L)/2/aspect(2)^2;
  131. % should select next highest index to the equality
  132. di=fix(di)+sign(di);
  133. % the new starting points
  134. idx=find(~L0);
  135. i0(idx(~L))=i0(idx(~L))+di;
  136. % update L0 to indicate which points we've fixed
  137. L0(~L0)=L; L0(idx(~L))=(di==0);
  138. % points that went out of boundaries can't get better;
  139. % fix them as well
  140. idx=idx(~L);
  141. idx=idx((i0(idx)<1) | (i0(idx)>shape(2)));
  142. i0(idx)=i;
  143. L0(idx)=1;
  144. end
  145. % reduce out trivial points DXY(idx)<DXY(:,i)
  146. idx=sub2ind(shape(1:2),x,i0);
  147. L=(DXY(idx)>0) | (i0==i);
  148. idx=idx(L);
  149. % these will keep track along which X should
  150. % keep updating distances
  151. map_lower=L;
  152. map_upper=L;
  153. idx_lower=idx;
  154. idx_upper=idx;
  155. % set trivial pixels D==0 in line #i:
  156. % this has to be done b/s we manually discarded them from L0
  157. D1(d0==0,i)=0;
  158. % scan from starting points for each X,i0 in increments of 1
  159. di=0; % distance from current y-line
  160. eols=2; % end-of-line-scan flag
  161. totlen=prod(shape(1:2));
  162. while(eols)
  163. eols=2;
  164. di=di+1;
  165. % select X which can be updated for di<0;
  166. % i.e. X which had been below envelop all way till now
  167. if(~isempty(idx_lower))
  168. % shift y by -1
  169. idx_lower=idx_lower-shape(1);
  170. % prevent index dropping below 1st
  171. L=(idx_lower>=1);
  172. map_lower(map_lower)=L;
  173. idx_lower=idx_lower(L);
  174. if(~isempty(idx_lower))
  175. dtmp=d0(map_lower)+...
  176. aspect(2)^2*(i0(map_lower)-di-i).^2;
  177. % these pixels are to be updated with i0-di
  178. L=D1(idx_lower)>dtmp & DXY(idx_lower)>0;
  179. map_lower(map_lower)=L;
  180. idx_lower=idx_lower(L);
  181. D1(idx_lower)=dtmp(L);
  182. DK(idx_lower)=i;
  183. end
  184. else % if this is empty, get ready to quit
  185. eols=eols-1;
  186. end
  187. % select X which can be updated for di>0;
  188. % i.e. X which had been below envelop all way till now
  189. if(~isempty(idx_upper))
  190. % shift y by +1
  191. idx_upper=idx_upper+shape(1);
  192. % prevent index from going over array limits
  193. L=(idx_upper<=totlen);
  194. map_upper(map_upper)=L;
  195. idx_upper=idx_upper(L);
  196. if(~isempty(idx_upper))
  197. dtmp=d0(map_upper)+...
  198. aspect(2)^2*(i0(map_upper)+di-i).^2;
  199. % check which pixels are to be updated with i0+di
  200. L=D1(idx_upper)>dtmp & DXY(idx_upper)>0;
  201. map_upper(map_upper)=L;
  202. idx_upper=idx_upper(L);
  203. D1(idx_upper)=dtmp(L);
  204. DK(idx_upper)=i;
  205. end
  206. else % if this is empty, get ready to quit
  207. eols=eols-1;
  208. end
  209. end
  210. end
  211. end
  212. D{k}=D1;
  213. end
  214. %%%%%%%%%%%%% scan along Z %%%%%%%%%%%%%%%%
  215. % this will be the envelop:
  216. % envelop has to be initialized at Inf to ensure single direction of scan,
  217. % otherwise may end up in infinite loop when trying to find starting point
  218. D1=cell(size(D));
  219. for k=1:shape(3) D1{k}=repmat(Inf,shape(1:2)); end
  220. % these will be the references to parabolas defining the envelop
  221. DK=cell(size(D));
  222. for k=1:shape(3) DK{k}=repmat(Inf,shape(1:2)); end
  223. % start building the envelope
  224. for k=1:shape(3)
  225. % need to select starting point for each X:
  226. % * starting point should be below current envelop
  227. % * k0==k is not necessarily a starting point
  228. % * there may be no starting point
  229. % k0 is the starting points for each XY: k0(XY) is the first
  230. % z-index such that parabola from line k is below the envelop
  231. % initial starting point guess is current slice
  232. k0=repmat(k,shape(1:2));
  233. % L0 indicates which starting points had been fixed
  234. L0=isinf(D{k}) | (D{k}==0);
  235. idxtot=find(~L0);
  236. while(~isempty(idxtot))
  237. % because of using cells need to explicitly scan in Z
  238. % to avoid repetitious searches in k0, parse first
  239. ss=getregions(k0(idxtot));
  240. sslen=length(ss);
  241. for kk=1:sslen
  242. % these are starting points @kk which had not been set
  243. idx=idxtot(ss(kk).PixelIdxList);
  244. % reduce out trivial points (D==0)
  245. if(kk~=k)
  246. L=(D{kk}(idx)==0);
  247. L0(idx)=L;
  248. idx=idx(~L);
  249. end
  250. if(isempty(idx)) continue; end
  251. % these are currently best parabolas for slice kk
  252. ik=DK{kk}(idx);
  253. % these are new values for slice kk from parabola from k
  254. dtmp=D{k}(idx)+aspect(3)^2*(kk-k)^2;
  255. % these points are OK - below current envelop
  256. L=D1{kk}(idx)>dtmp; D1{kk}(idx(L))=dtmp(L);
  257. % these points are not OK, but since ik==k0
  258. % can't get any better
  259. L=L | (ik==kk);
  260. % all other points are not OK, need new starting point:
  261. % starting point should be at least below parabola
  262. % beating us at current choice of k0
  263. % solve quadratic equation to find where this happens
  264. ik=(ik-k);
  265. dk=(D1{kk}(idx(~L))-dtmp(~L))./ik(~L)/2/aspect(3)^2;
  266. dk=fix(dk)+sign(dk);
  267. k0(idx(~L))=k0(idx(~L))+dk;
  268. % update starting points that had been set
  269. L0(idx)=L;
  270. L0(idx(~L))=(dk==0);
  271. % points that went out of boundaries can't get better
  272. idx=idx(~L);
  273. idx=idx((k0(idx)<1) | (k0(idx)>shape(3)));
  274. L0(idx)=1;
  275. k0(idx)=k;
  276. end
  277. idxtot=find(~L0);
  278. end
  279. % map_lower/map_upper keeps track of which pixels can be yet updated
  280. % with new distance, i.e. all such XY that had been below envelop for
  281. % all dk up to now for dk<0/dk>0 respectively
  282. map_lower=true(shape(1:2));
  283. map_upper=true(shape(1:2));
  284. % parse different values in k0 to avoid repetitious searching below
  285. ss=getregions(k0);
  286. sslen=length(ss);
  287. % reduce out trivially faulty starting points
  288. for kk=1:sslen
  289. if(kk==k) continue; end
  290. idx=ss(kk).PixelIdxList;
  291. L=D{kk}(idx)>D{k}(idx);
  292. map_lower(idx)=L;
  293. map_upper(idx)=L;
  294. end
  295. % these are maintained to keep fast track of whether maps are empty
  296. idx_lower=find(map_lower);
  297. idx_upper=find(map_upper);
  298. % set trivial pixels D==0 in slice k:
  299. % this has to be done b/s we manually discarded them from L0
  300. D1{k}(D{k}==0)=0;
  301. % scan away from starting points in increments of 1
  302. dk=0; % distance from current xy-slice
  303. eols=2; % end-of-scan flag
  304. while(eols)
  305. eols=2;
  306. dk=dk+1;
  307. if(~isempty(idx_lower))
  308. % prevent index from going over the boundaries
  309. L=(k0(map_lower)-dk>=1);
  310. map_lower(map_lower)=L;
  311. % need to explicitly scan in Z because of using cell-arrays
  312. for kk=1:sslen-dk
  313. % get all XY such that k0-dk==kk
  314. idx=ss(kk+dk).PixelIdxList;
  315. L=map_lower(idx);
  316. idx=idx(L);
  317. if(~isempty(idx))
  318. dtmp=D{k}(idx)+aspect(3)^2*(kk-k)^2;
  319. % these pixels are to be updated with k0-dk
  320. L=D1{kk}(idx)>dtmp & D{kk}(idx)>0;
  321. map_lower(idx)=L;
  322. D1{kk}(idx(L))=dtmp(L);
  323. % ridiculously, but this is faster than
  324. % direct assignment
  325. dtmp=idx(L);
  326. dtmp(:)=k;
  327. DK{kk}(idx(L))=k;
  328. end
  329. end
  330. idx_lower=idx_lower(map_lower(idx_lower));
  331. else
  332. eols=eols-1;
  333. end
  334. if(~isempty(idx_upper))
  335. % prevent index from going over the boundaries
  336. L=(k0(map_upper)+dk<=shape(3));
  337. map_upper(map_upper)=L;
  338. % need to explicitly scan in Z because of using cell-arrays
  339. for kk=dk+1:min(shape(3),sslen+dk)
  340. % get all XY such that k0+dk==kk
  341. idx=ss(kk-dk).PixelIdxList;
  342. L=map_upper(idx);
  343. idx=idx(L);
  344. if(~isempty(idx))
  345. dtmp=D{k}(idx)+aspect(3)^2*(kk-k)^2;
  346. % these pixels are to be updated with k0+dk
  347. L=D1{kk}(idx)>dtmp & D{kk}(idx)>0;
  348. map_upper(idx)=L;
  349. D1{kk}(idx(L))=dtmp(L);
  350. dtmp=idx(L);
  351. dtmp(:)=k;
  352. DK{kk}(idx(L))=dtmp;
  353. end
  354. end
  355. idx_upper=idx_upper(map_upper(idx_upper));
  356. else
  357. eols=eols-1;
  358. end
  359. end
  360. end
  361. % the answer
  362. if(iscell(bw))
  363. D=cell(size(bw));
  364. for k=1:shape(3) D{k}=sqrt(D1{k}); end
  365. else
  366. D=zeros(shape);
  367. for k=1:shape(3) D(:,:,k)=sqrt(D1{k}); end
  368. end
  369. function s=getregions(map)
  370. % this function is replacer for regionprops(map,'PixelIdxList');
  371. % it produces the list of different values along with list of
  372. % indexes of pixels in map with these values; 's' is struct-array
  373. % such that s(i).PixelIdxList contains list of pixels in map
  374. % with value i.
  375. % enable using regionprops if available, faster on 7.3
  376. fregionprops=1;
  377. % version control for using regionprops
  378. v=version; v=str2num(v(1:3));
  379. fregionprops=fregionprops & v>=7.3;
  380. % in later matlab regionprops is actually faster than this code
  381. if(exist('regionprops') & fregionprops)
  382. s=regionprops(map,'PixelIdxList');
  383. return
  384. end
  385. idx=(1:prod(size(map)))';
  386. dtmp=double(map(:));
  387. [dtmp,ind]=sort(dtmp);
  388. idx=idx(ind);
  389. ind=[0;find([diff(dtmp(:));1])];
  390. s=[];
  391. for i=2:length(ind)
  392. if(dtmp(ind(i)))==0 continue; end
  393. s(dtmp(ind(i))).PixelIdxList=idx(ind(i-1)+1:ind(i));
  394. end