From 25efae8269e5647b17dc69329f412f914d26751c Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 13 Jul 2012 00:53:51 -0400 Subject: [PATCH 01/21] Add morphological reconstruction with test and example. --- .../plot_peak_detection_comparison.py | 218 ++++++++++++++++++ skimage/data/noisy_circles.jpg | Bin 0 -> 26206 bytes skimage/morphology/__init__.py | 1 + skimage/morphology/_greyreconstruct.pyx | 85 +++++++ skimage/morphology/greyreconstruct.py | 150 ++++++++++++ skimage/morphology/setup.py | 3 + .../morphology/tests/test_reconstruction.py | 73 ++++++ 7 files changed, 530 insertions(+) create mode 100644 doc/examples/applications/plot_peak_detection_comparison.py create mode 100644 skimage/data/noisy_circles.jpg create mode 100644 skimage/morphology/_greyreconstruct.pyx create mode 100644 skimage/morphology/greyreconstruct.py create mode 100644 skimage/morphology/tests/test_reconstruction.py diff --git a/doc/examples/applications/plot_peak_detection_comparison.py b/doc/examples/applications/plot_peak_detection_comparison.py new file mode 100644 index 00000000..2cf7f933 --- /dev/null +++ b/doc/examples/applications/plot_peak_detection_comparison.py @@ -0,0 +1,218 @@ +""" +============== +Peak detection +============== + +Peak detection (a.k.a. spot detection or particle detection) is a common +image analysis step. For example, it's used to detect tracer particles in a +flow for particle image velocimetry (i.e. PIV) and to identify features in +the Hough transform. + +To simplify plotting code, let's define a simple function that creates a new +figure on each call, and removes tick labels. +""" + +import matplotlib.pyplot as plt + +plt.rcParams['axes.titlesize'] = 10 +plt.rcParams['font.size'] = 10 +def imshow(image, **kwargs): + plt.figure(figsize=(2.5, 2.5)) + plt.imshow(image, **kwargs) + plt.axis('off') + +""" +To explore different peak detection techniques, we use an image of circles with +added noise: +""" + +from skimage import data +img = data.load('noisy_circles.jpg') +imshow(img) + +""" + +.. image:: PLOT2RST.current_figure + +This image is noisy and has uneven background illumination. The peaks in the +image, while readily identified by eye, can be tricky to find algorithmically. +The first thing we need to do is remove the high-frequency noise; this can +be done with a simple Gaussian filter. +""" + +import scipy.ndimage as ndimg +img_smooth = ndimg.gaussian_filter(img, 3) + +imshow(img_smooth) + +""" + +.. image:: PLOT2RST.current_figure + +Thresholding +============ + +One way to extract the background is to threshold the image. +""" + +thresh_value = 100 +background = img_smooth.copy() +background[img_smooth > thresh_value] = 0 +peaks = img_smooth - background + +""" +Here, all pixels values below the threshold value are subtracted from the +image. The resulting background image and the extracted peaks are shown below. +""" + +imshow(background, vmin=0, vmax=255) +plt.title("background image (thresholding)") + +""" +.. image:: PLOT2RST.current_figure +""" + +imshow(peaks, vmin=0, vmax=255) +plt.title("peaks (thresholding)") + +""" +.. image:: PLOT2RST.current_figure + +Because of uneven illumination, peaks on the right bleed into each other. +Increasing the threshold will fix this problem, but it will also cause some +peaks on the left to go undetected. + +Morphological reconstruction +============================ +""" + +import numpy as np +img_r = np.int32(img_smooth) + +import skimage.morphology as morph +h = 20 +rec = morph.reconstruction(img_r-h, img_r) + +imshow(img_r, vmin=0, vmax=255) +plt.title("original (smoothed) image") + +""" +.. image:: PLOT2RST.current_figure +""" + +imshow(rec, vmin=0, vmax=255) +plt.title("background image (reconstruction)") + +""" +.. image:: PLOT2RST.current_figure + +This reconstructed image looks pretty much like the original, except that the +peaks in the image are truncated. The reconstructed image can then be +subtracted from the original image to reveal the peaks of the image. +""" + +imshow(img_r-rec) +plt.title("h-dome of image") + +""" +.. image:: PLOT2RST.current_figure + +The result is known as the h-dome transformation [2]_, which extracts peaks of +height `h` from the original image. To better understand what's going on, +let's take a 1D slice along the middle of the image (cutting through peaks in +the image). +""" + +img_slice = img_r[99:100, :] +rec_slice = morph.reconstruction(img_slice-h, img_slice) + +""" +Plotting the reconstructed image (slice) next to the original image and the +seed image shed light on the reconstruction process +""" +plt.figure(figsize=(4, 3)) +plt.plot(img_slice[0], 'k', label='original image') +plt.plot(img_slice[0]-h, '0.5', label='seed image') +plt.plot(rec_slice[0], 'r', label='reconstructed') +plt.title("image slice") +plt.xlabel('x') +plt.ylabel('intensity') +plt.legend() + +""" +.. image:: PLOT2RST.current_figure + +Here, you see that morphological reconstruction dilates the seed image (i.e. +the `h`-shifted image) until it intersects the mask (original image). Note that +the peaks in the original image have very different intensity values (e.g. the +peak at x=200 and x=100 differ by about 80). Subtracting the reconstructed +image from the original image gives peaks of roughly equal intensity. Thus, the +h-dome transformation is quiet effective at removing uneven, dark backgrounds +from bright features. The inverse operation---the h-basin +transformation---should be used when removing bright backgrounds from dark +features. + + +White tophat +============ +""" +selem = morph.disk(10) +img_t = np.uint8(img_smooth) +opening = morph.greyscale_open(img_t, selem) +top_hat = img_t - opening + +imshow(opening, vmin=0, vmax=255) +plt.title("Greyscale opening of image") + +""" +.. image:: PLOT2RST.current_figure +""" + + +imshow(top_hat) +plt.title("Tophat with disk of r = 10") + +""" +.. image:: PLOT2RST.current_figure +""" + +selem = morph.disk(5) +top_hat = morph.greyscale_white_top_hat(img_t, selem) + +imshow(top_hat) +plt.title("Tophat with disk of r = 5") + +""" +.. image:: PLOT2RST.current_figure +""" + +selem = morph.square(20) +opening = morph.greyscale_open(img_t, selem) +# scikit's top hat filter uses uint8 and doesn't check for over(under)flow. +mask = opening > img_t +opening[mask] = img_t[mask] +top_hat = img_t - opening + +imshow(opening, vmin=0, vmax=255) +plt.title("Greyscale opening of image") + +""" +.. image:: PLOT2RST.current_figure +""" + +imshow(top_hat) +plt.title("Tophat with square of w = 10") + +plt.show() + +""" +.. image:: PLOT2RST.current_figure + + +References +========== + +.. [1] Crocker and Grier, Journal of Colloid and Interface Science (1996) +.. [2] Vincent, L., IEEE Transactions on Image Processing (1993) + +""" diff --git a/skimage/data/noisy_circles.jpg b/skimage/data/noisy_circles.jpg new file mode 100644 index 0000000000000000000000000000000000000000..a0f49d6767ca4ca3cfbef6c54c36039223d243ad GIT binary patch literal 26206 zcmeIa2UJr_*FSthBIOF9T@^x$B*3+R^b$&lAP{LPC@7*-A#_M+B81{axJp+6m7*X; znhHo$1BwDFy(!WJq?&+$L|XDc!E*2OywAVA-}=`2*7`hI=bR~f_UxJ2)ABpBnQxgN zgclA45zYaCsVT4*002&a8v+1O&;X?uv z;rsE%Ie9xcBagUxle~$p-abe@3l*fPy(7VkP7TRy-koJGfc4_pMeIXOAGIJvpF zxcRsLww-?mFE=;u4k1Co9fE>F{M=jD_t!Vizds>7+qd(;dHCRPzTe?+`0uO-{P%B3 z`2SBBFzW$+P9Ph|hC(C(Hhu_{AHp003;7$f5rPG!+y>dYK!I5o6o9b7wy|??a{aa) z^6!xlfDQU1k`I8uAW$|K+cpk%&TUX`B`}g73fnE9unl+AUQp6ELXrLVq}vY;qR=~z zH8?0K`$cx@aY!XwEPVLvNL1N*=V6qzeu|~jp=UWH)m>f3$D9|-FFx<~A73(%xtjZD z#kIzsiRHVMFM2{bwfffD#i3T{)3Ncx>6LBp}m1tqk@2fxT?{E|q*b1EfkNc$ayoQe zmE=73Q-uGy1Tue%z#IU$p<8+K16W`=Js*+RttVo8T&1M+5zto|xANeMOQ++hzIE$| z57%;T=9uu`R%*{$y=QBVa%$+AtE#_{s+5{`dST=uAMC=Hjh?F`<)t=?{Ta3mQ+m`u z@7X(!{+B&2MCILttpqx^Vy3E7#acf{4T>)okJ}Vvp!J)3rcvwoJyeDA!4EVfeeqN2 zl@ev-%Rha;`13tOCZ6LJP53f+EyKm#s7pj
  • F%n8k9it%ULC>PP!Y|GKDK z=-BzU=3Ug)BAIIQd$JKV4uQrd9rsVD6s#b~1}=xW{d zK8fSL-cOctBV8g*$P=ktW;j>4llYRtvR){q_r%TKGddAb1CV}7rIdE5m9;E^o0P1K zI=MMN!~}BhE3OUQUhJnX!nh>*Ha^C>Sk%tv23!&&=_u8uhVHPvZaw)rHVzuo1K1gp zuO?@7_AADswW8i8De8V@%yu-+X$Yuclm+mN9y!9Fd7g;t#i!Psk}9`IlKdDx32gJbV&__sH0(SiCtd zwnsFy)Ca9khBE;}x4t_s-*ze*^}q159Gyi(q^fF~W_G!nc{VhOEu(81_#yrG^ih6B^qI?N0@D71i>0gfwqK>t@ z9X|Z2PSx19R+6hnylIFER(#RPw9yX;`T8DcVK+VL-HD|}y3&Yw!2ISPdGZkP(J4Q4 z7XLcs-AZm=;%fowXog*A$H>~AO&n4Xjn}*q_3~XOr8dE9pLL+HFAu*$E}g!&^wFbB zT;XW=xc=3*U$k!}<(g$AtQKzCwhKxk`82Xx-!ETNeX{;qFZ8~hkw9Qb<<)?wVs}hs zCt%L=p20g4B4FWFRgnPTfH$*YNwhpMxCm+d%;1LF-U*vJ99*m|AHRmzpb|kzQuT$3D2D*3``RbLvINaLl_8`#$$MLtp1WaR-q;@TBqsPj^kwf07_wgHQ za&mg*a}tl)jre*!2lCbUXI$u%k7K;Uh5fSOCWg^b<{K@ogIA;_va`G|jV0go?BKCm z_U=b{pi37cC`zCXzT4;zxt>vT{I|`hcq2r_FgaLC zT@2q4hX%%!Q$zX4YOZT?E7zhw-OCGPYwixWixYc4qw;>oaJhZRFrE*R3K>Uv_0djB zp0;las56VeSK&VG6ErrXXj zn$FWmP4D8GV7OFmeCZ~I-96nnrzOUAIKZ|`Xs^IquN!7r3QQo`TKjPHu|;K?&HJa6 z)YZx-fv1YP2ATpi&4a$m5er-%%qQz2<8;&evj|K;LrUddSPAL@Y^LBdec#)*kW2_7 zwr7qb(*LAs3n8pS_cQ5hMmHmHvZ%U7Uf-Rc=CdGSL_XQFXRbB?b(SW}=dYdXnvfCG z(amSuu7m(P&Yg0f^F@ey_V_;BdOl+#_F8XT+|&aFe#sCODdo(*`{D2|8iY?~x}js8 zlAQ5pU3iycq$c^~NNv%3-_UG>>`0@IsIP?&Uxs&`GD@kkalP^gkwd+Nz3uLWc4$;LVv`)8}c%X|^&Mx?>mmcc&&! zRY%bjid(vLQolSZy)-A%Ow8171yVJW9q+IA?+9H{SeTZR7AcFqw3c;9O1_2odSF!7 zbwLpQHTJl$=d~vfyK8kN14E+k7W*91@T|O>hf}J_HD4i5*YFw=q;lWrr@=o`6WYq& zF#)gKycG#cf15zDQ)75>R02V#ua1tTR=N#W+PlIU+5E*W`N!MZ__INc2K2)DNaxDe z8)j=#g}ghZ)x`A8>ZYb^D5g?sv{~MYHv&${V%a3Rw@5$@@OHca*1o$@_`{2wvv)Tt z_>b2n>cT|=iI9BvQ$yXW&2d2DTdi0mKJ5!xZ8SO6Nj>C?J6=|#LMf-(t^D}p9TDrk zkh$?sm~73!%0{vJIqVbctg z`Uyiw-aEQ5Zzj;<;Z|j1d_%_I<-CqgB4O=tj@!UAfsP4i&YCSZe-k!;bZ3>iTG+I^ z^0t>GtR9h~xxQylr+Z5`siJ{r9(&`%AH|hTu0GqE(&xuF&r$&WmXx9itf`nrFQAW42E+~B6}CUqpE=ALueeuPQFopXSW zgz@__hF!uny;H`zX%p$=P4vtuje}PVYwBG&9Jm|iMZ?@F*Bh@kRMQ?TKzzj7lBco8 zt{WrSU))O5Ltf5~Y94(|5WuP}1u%grclSKfeQLyjXB;nKb4aKbbKX-j8;Fz13}FHU zDXkq}ReL(oR8g;<>CpFN4aU%`YZZ@$i1v+HA8|4xq9LKfy#rqYRN?~Od51#crLqj^ zR;X+Cj<0Wwcz33zCFQ*Zn--nR<9WNYJ?e&2dkPbnbtg|1?HJdlDX^Jjj^T7@Se+f` zqXri=q<24mk6Eu^m@}>{_nUSsh}nq@VS+gsMFJL=S2X8h?0TkcbA_hpfmN%Ib}ty? z>1>V{+d z5eJ)97V(!zqaNQoD1>n0748T}EnO^uXO2!q^Py+7QY30tuY}EiMCbMCpa#)wxLCHJ z@z;YR%YU`M^E@}c`h2NgX?1`J9Q7u>!Hk^kghq)~Z?+E|YJXl@`}zum9Y1Q8HGntA zJew4X?9Po>#@>jtQHONsE5nU8EM*CBpA8xf8nya9V(pJ>G+1T1xG3+^Z7J2lf)*KT zx`2PU7I|@>W4=SA6BKG266U?RdbHmDK^KD5#le1f)<+>t;e!bBI2j8$y1_`!ll}Xk;+;Gxh(J=PWM9mC>PF*Ev3iY*ZffzIJuTokoq`xe zY7UOyj$J!5;<_>aDL4P#Oie^wI@%q%h*j2P{M8p|vF#!)24RxSwRHNB)ze3~fDRZW zANF$oW0>FW*#ok&GDq9Gl%X@|eB*ml+iji1OQ z2raXeT*lel33-=B? zC#8&vQs1(3o2RAiq9X>4Lkm8BOo*;kb@k6$qLecM&)oddSF@x2+}e-mm=sKn-|D-t zlDnln8Ab}@kCZoiEkB;23y?g#qTYP2piD%2$;=2)hr2%TbBxm%c_SZl@6sQ)7h+sz zVtdw^!2G=tb4vFk^DZ<$G@;46=E$QVY9b_E>b_I^lo1o4Fvc>Z>y7nsffy)Fe2UzK zH}~(xl%KSQh3n%4H@(kik}W+;Y+i`~g3`JJX+_Fa?)7LPjIdb6+TLrX#7E7EPEq67TmDvV(STkHl53Sf zOQGV;uFDSJV&BL=bAbswtkt=h+2qk#OUS$HcP+y?1*1>CCBhf5-j{D$;T@}q>mjr( ziI3!hs?^)F*4*#qixlujy4w5HR0?A7lF*tIT3j0{plSZh?Y{eZ4J*xp+RF>4&5t+y zz6+0h6C>6ZyJynIT4qo`L(d3ETeEqZaJT$D)sW{BMnx`jskm)q;x55!)fbwK7tGd@ zI-56MZqSlML!I52bE#tjC(wEFi>DEIY2I@4HTdpBIm>@$2MMBi3||jd6g2E?4}^#V zL1bqHE$ROH&6=DfOehnmPQF#`YDeg7G>Akkxzu%>p4z2S8j&5TNvv|y{z`L?(vyUb zsViY;8>K5wOnLdBsy0^wXe;Xhxp|&>r$p*p(bZf6`jy?_u415f*V9Mp4?1}RatJ|9V`cR5oO7O(5N$W1q!(j3~TxccWJNs>#%tnNvhfHMDLfv$U>awvA zfV`y2seov&0;kFLd#CK5+C$2AzYSM`pri z$IqAAWR7d`O1LUQKCj{TeDo`n4f2&1X^?3;LbD!$ zksYra^E)*qje0Z_Ha-&5GAykCv24ybv(DAcC*Wq3gLdNK`48`)gE)ks!1Tz-l7Pjv zzP+x&u03skc4y1Hmw2+AcE2SaBACm|`C(+>RC^BS^B!pp9}Y@qV_~RMEB-|nN7tH7 zDz3Y>gtU|9!?ASfqQrq&vTAZ>4pNav$HM zCk&&kW0;}7(HqgD7C0^08%m5AOW zc)Gaj@|pN5YVOz4$)Bn6nt4ZWUm-94k#2jm-PKi*H*DZ5Z|nKP1Q(S*r;d%P$5N@d zyP2a-?sbc-KSkn0J|SpRXH@uf5;uDHo-6>;N3*y-Z;Z6z&Du@x=0?>cqCO(j0xsWv zS#rheqM`JKuiz%XrPOg3stZ-qeKJ$KG>aN?Vp!bktAycEol~L0;3~a;21D(bnuC!` zldA{cyZ8h==Ipz00UPX}5nATo=596;sDmiAlHq=9FrPE@)J`)p`Bq3&dz}xPq!A$Y zn6OOOU_*#@&G^V7n*Q4P3U#9x`aoPF~Y@XxfYwvBnW-8=4AXx}aW1RpSn!@LQj z&;L)66QR!kMwSOMkR<9KqLQU_E3{hbMc-wQr+`O>_KzxwX`;pG&V-}eDjG%CC z@yYE;b6#DdDzWj;=}!};c0asn@PW<6I<5a<)uGUwIPf36K(0Mfrt@!NhNNYx7=9H# zcT4${M^_mElorb(*EM-wGPLw*Ux_Ft)+gc4uzB$2@G%Q!2NjbE^&l=!BzCuN-+lll z-!=qmnArAiV zLkPO0;g?;uSzp4|#f=D)x}PAc$fI9~yYq}yY3V%tFvHclV-_h1U3%{aG=Mp}JCz=G zzLs{e!6uv2qrEZg6S?H_-G=G*0-g~g_ri%d+h;X_rq}i1FSW+7MXN1)mgMesuSDvC z8|m~7+Q`j~yYaIk^MlVI;?&Sa0(EGaqIrJ#7859{I%T4afoYm+ob+DH30s&qfhr;m zQPwLAzr#VZGZi(LPRgJ8z;{|bvho4Oc7@L@(!uedSv)akMs`5w-f%xSwV*=}aEg!Q zY+S30?Z~nTsJU;5Uc{~D=3lPM9jHI;FV9XhQsCH_&n=`#bnULQi-4&^%2up$oJcM; z`3L8tZ|x6xS~6$^xpWNQk>O}Naw9qO)D3Z|>}F9;EpVx!HeJ+XzkZ`bjGq?iXr8c% z=zBpv@lc%`Dd3=K4(}hUl&|Hv;xnjSt`$VyBc*cQ8=pL+>6=ICbL$T@b{8`8$m~$f z2tH+<=6N?<3_oQ`p3WIEo4me&g6pctzwioZHk7GWmF#H~JO5 zcHzNRb6%Ys4-b=>@0e8NvA+=$B{_LO{UF{H#zI>FKvooT3V@188UB2O=nNacCdy_2 z@E;)rIeU28lf3C2>94z~bet=m!37BKPqTU;?)aVg-F54>4u z2LNpM@dE)GbAOVLKMC}*AR<7|!_UIm%ir^doAo|scuqFkXN~spmuT`oR{I!NKpR=a>u8^ruJ9&r9!bj~D50kBNf^@n?*g3n}2I z=a8p|{!b5B%)eP5baZzycW`qA^Z!o>hV^ zA&@pIHY`7K*k>p`hoZocnG z8sx`&t7WCP;^5z?Sx;>E_crI*3S;TJ_!gz`zxts|TYj#sC-{JDX$purvaa8LGywqF zl3;l&Q~*!`-w;3+e1g~@%e@s#-*Srr0P{ayP!QSNl7MrI?z;xE>4H5DyuZh>kWmn; zWIeWM5iA$f9DIDsiRA@YNbT0m19$=600AHZ!N6I-5p?;1PY=Ku4Ed=(tOP%=?<%$B z|H0*FB+Kznf*%oocNW0Q--Ff5SRgLDgSWqz6Y<}4)sdtC#%&dV)xaL!UM@e_uoU5+ z{(~fM&+mSgD8Mln*B|jLKReNb;OI;|;bF2>n19J<+g2P%#OkLUTOPdrKX`Mv_<8&L z{7VqKw;#cU;PqWgxer;exSMZ9Ks{};Od3*G;P?*Bsff1&%o(EVTN{x5X@7rOro z-T#H||3deFq5Hqk{a@(*FLeJGy8jE^|Ap@VLihjM=>8T$u^R*)0N^-ygO~*fjDx3Q zoIqI60fYgOz!AU|bd$iFxaINzDFpvF6i7f1u=oygZsqU|d*s+cJit7dA9ylcNhBZ4 zfdgJdIeXTrOgTqy&jUgBJ_i)!nM^)f|ItSm5Qmn zsgItsE5YczpR>hzGfT(w?v5Hxl2{$2RuCr0)5p`9WRDE;^zb5Lg0v;Km}5XUEBJsU za*KrIt}VHBS`uk(dIYKG?dObCk&~BolvkBUsw&7SC@HI{sqIIC2SVi!$SWRDP?S|r z!zd|Yv>B zdy(veWW9)}9~|_ZiH?4Rt?veqERG;1?oZN|1SQ=nf|KJAxj30!Ke#HMw=m>r&;b$-FXCEjF z;%p>>mkUcQ`JZBeIhdOM=cGUK@bvsiO(Y!(0L%ZsQV${v{X1~Xndt5B=jePW0IUd< zB+|(dbI#ij{N6-clHh6Y;(Wl($Js>^`Lom*6K5BD$6)*sOFeypi!%|VQ1@(LOXxPPYJLKjWFot(8L6KYSD5O&`&d0}t;0$&jOFw_-t?W!qF?cT`$==J+ z8LzJ`36@8WKybpyE2x7~swvB>Daxy=D=6Ufa0UmJ4OCSRDJv@}C@cO*tMBdT&r-l2 zX`TKr(tg9t?LGcmJivYjHZ;-S!43R~gH`8|;E->R8!3J5vJF9;=nf^m8M^}3<7gj4CklZRSYqR|Uf|Lq5Zp-pn=u&xG$y0>|JaxeXbWaL zC-~86GX%m9Y*~R3W;OzMz&+>x+zQP92MoP{hE_Ns;I0@nvQkvWO-PbSo7%LYo(v8( zErAQdQ^uRRP6~?=CK4_PqD@@GQjYuf7t8s`MJaZ0<9_Rtu_k868CNE3dl1*NJpw)| zj7n5lJG*>s(rKL};21j^HY>Jz{~bGLzaEMSD6Z|3Q7YbfiH{bBFsvJRpopP9R_O>% zp|Wv&NwMgo5noQ}W)byQ8E$W%&GzkLuj%57P9&wBk}3RDkZ8w)FcEWP6gp3RDVOIC z*M3k;KIw3>kLo_EMdi>{m$*qg;(g}Hix;ISZMjSL>Mr+tN!`3s>%_$b&bkdY%$5cL zwc>_iOkld?)Z)No@)RyMZnxFzSsqU>WHU59VJ)kT?rr_f_GGK1C>OgehtctOy(K$^ zeLGr2!4sVtCtuD6q_UZ(qVW>+`3czz7oLqlI!qGLP4jBgCF~ldmuM+s3}f)X3Gv~G zHF~BP`gV>9q=z|7DU}DB2#gnW<3Bn&O|+!bAuSzQvy53mFUX#H z+E^ucIxJBGG|Mr$Pw`W_j5RCv7tLiUAHmb_cDn)HNnB z8j~;=HpB&wlCdRMl>42(nbp*piNZst!BetUs;NSXMoNgrFybSI+rB^(L9~FOFLts> zXhr9gBr=xcK6bUg`;>Lfd_g2Pu0wqqtG4kmyRF>L9jD0E@g|%^N80w>n&P9PjgS|x zi-&LDud{NU!Nf+AH8>gYZS^T|Ayr<=6nZjy4c)FNuL%&Lb zM@r)Pa9Lkga{80|*o-=2edeVtpAN!sNmw4NTHxA~--dI~aYYp4`x&eX7@2!_Ys@GO zJtJ1JIX~g4`)~|D)WU~i#81etCTjYbLzj{cB}xs{PteoW@U8KHSkuu(hI!>~n&z9& zg-W*>W)tfs6@p5I?AUuT^_?&wS#}4!mfFzjtZ8MdCOkZTAi&T$C-!>A1x0d;nl&A>g8AF2hjJ2K80t|Uc38X+yCgnRqU z9qgdXhDytis2daKJqC0~-+UFa5fX;;@;NoSzc-~EA~U-NW+#nfr2U5^uuE35$!7FS z5hT~=&5?}{$3x?jQBljb)P%V|&X;*bjdj0ecunoTl#(wOXV}TGAgdf|C*Ckr42t!^ z;vpX=xH{mF7-}WlV?Re^CyiaXdx&bm`(Oa$VCVOHRWd#xH#0b4?QDz1PPTVxNd=fd z>~e~?_N3$m&Ocg1?%?B)O`FCP?Oi<`f&=_HcK0u#s=wwTTKxFTfzhWGPC^2?*e37ojw@KGL^xClcxuCz9R$D62N36#bR6 zLmY4DUVk47i^@&l4XhEq&MEDHgzy=z9JPifbjG3+b&Rcm$2HmdY%1d2>#v#JBeNAX zO+(oqA0iE>G>97eiDz)&StxUf42ECpeJUnE9`}{u{wVwW?Z&kJL=A1aie69;d$Tms zgm-;R#Xq$c=uZSsioW1;;4$lPGJC=VGzWV6#wys1#$*~StAG1~WCKley;Vj;A^(tn zYqBC69o4k>^fGxSr$YF8tjr8}kg<-+I4titc1j~{&+{|f*ONIhpC#s&f0ufG1|qLK zUp9G1WkK3bKlV^bd{Ab06RM6}T5Mr)Yt*KC4klPA0ooNPvN{=O?<#h2mb|oM`yUB+ zOjHsztgd_QBukTn)iwf zvv(HPcV~pfk9^#>&awwSxltHDeLOtyewXoqwIPx`QVoMHf|@Nkm+Q~>hG)4$r=U%0bJ>2gYw zfKaI$v#(7ch5fGmjLl{8(wHgGSU;s#7CmuJS$(c0j{PEekDJ{X&1cYw*rceF+G#5X z-*-xlJz~^+^T66&X74;vTVB$s%GO+>1vK(h<5l7ai=MSh#`iCgXF{w{J*9*P^*-Eu zB+$0g@^y<+_gmYqq>necCM5U_Dn~TO88ND8r78J|x-@)ryDt+6U9;(X8u^ebwueVW zqHpu#pS@IfShzuo%n2qSGb3lSw_jtCYVdd(8yGM$n;n5+%qYT=!2VM7J3lf5!Y|ge z`I)3G_mU0frvjeTTf5uyVZu^~E8tS3`ZkN%lh?%9qtWW_61}G%+1uZzBbKzECT#A! z*GfEPXE0422m99v*9Uh*L!T^9=u>Tm_O5V7@!!V4UuTyUoah6b2z<&hI=y2ev?ZTkGXopq+#LYkoCQunpXrMFg_Zm(lU8wpLLz_&FW_8nG4vZ9i7Yd&9GX? zV*19=NbgE%l!RAj`&DuosC<3XYo)$+VE$Z)!11NjJw3$CC0TSpNATvZ{oHrXfJ2Lz z?W;E8MUyufj+iBvQ5Zj6>wdl5sM~xUc7Fo8iri^hOYeuMN563MPy{&Gv|rKWoQs3eFCQ8bz|)VC zl@t$+6257Ez>KomRYE?l3Y9i#;Ls)ayWRGr$?;LyFJ`b)o^0M_uW4wENqQTHW88aS zdvSfxq%moH<64M>LMweiKNb;QrdC+JZkC)Ku44-1FDpxs-07>*u7a!3E!CXIb|MZ?nUUtqg4NRysgWh$> zBUOOP3$*25a}JNmG>5HVjhB+|p}Bu$YwWaYiz;Bh^-JaBL(SBG3IM^uPwyPnb94^l5*=hKaDi&DG-p#y37t zjtX!dB|Z<2=S8nd$zn2)uXD;TGJ(#4^ma&YJ)cuDsUKX}q*rn^`N}J2=uY^h=1u|9 zCW4}KdreEZt6LK4QdRH|X^qpI&?xLH#;cUv@ZcE1@m0#9!l>|$^znIy)oM#qwceG? z(=dS``iZ_8{Sy#avsQQz!0_9*mj{^v^Ucj#v*~$n%{MG_zfWxWs8APPloyhnv_4Wt zS&j*hDj?13gR71}Nj97uX6CaJ*ghLC{9V8nzmc-^pKO)kO=7#*5L0Hp(2} zUFH=b2`@>Tz0#>B{4fbBdD7GYY46LCObSHI)R-K#o51rcFagUXIMjgJPni;?%S^;M z1cOKXL&*1XUvqWl@kWeF_DzYML3+map&9hmf=^|>9lHET!zEExj|ue&0#Wkc)3@>{X5*aHYZCpi|lp^^2Bx!ij!#GU#w(P zalA-QQO*faF+(klWJwjV5`hfASnH7Xs*?-Ujgt2PvRWo=iyb0`qdYFxp2T({5T|u@ z&KV`l?K#^m^s_V*rN_F0)lC8(#-!DX~z!E(cq-L+gDyg3ZYT_T|}A^lQa|s=B@d(nd7wz25Tf`!^UQSh(I}-py4nt8RqNOG+)ae`arbP-IfU z8roF_15SUFFQB8g^GNoID_Fi{sZ!_5=GH9x3+4W)*~eC{%Jk5#fhF=ofhGfw*QkdnFh z8}Gr9e5WMwJREu>gCie9jz2z8F>AwyjzEOb!DPLqm-GG_%BF?IE=cq;jIUEqSoX1X zY`^0auCe>rZ2Ep~U&IEyHpCXxm>TiB#50R3sFNJeO({3QQDT1;;E&Kv69u<7dn_SA z{4#i0FwfJG^|9$%bpj~{DNF}D*3{YMUB?eqtfBezBaN&1x62gS9u-hl&>oLiw%V60 zbXq@H=o9wdNNy3W1D7PAkr`U>woPgj3ZFC>;SKDvwK5{Y?{J0Vv4)R6_8j+j&XlRv z?o+dDDNaTN=`aDo&9YYOM|@{FtlVG5z8#&GsUGzqw{G+X3YQ#k#RNvY-8`_Aitl=@ zHz3>+2az!!vQYsn5AFuAXKIzwz2EOZQx!?N3wipFtx*o~LPt zsFj^sDP{Jn!ym0_Q&kSipHC=Ixs{3>kC>IqFbdtr&P7+&Hn8ajPJu(%jSAtO1}HHX zK6yCjW=~)u!k_&@^~s%YZKawHs>poUxNAFxmjtA7*yYedU`j~Q#Uhul%-*>+O^L(Q zunA<)d8%Y~p*44!)NQ>npR;+*aGz3d3Tjj{fx5qH7sb$(D^|;W*-`~a4yFZlniX*d zOfIUR;|xLLX=|VQpS^iIBXblHhgzy4A9s&09eYP3fz9H4f}+yLk4RxRizJ7v-K!3s zJVz2NjM?ZV+%)2EK{VnbKqD);r3xQQv@e=H76YdV@_RRis$LpSKqZ_bt8SCrS(=c4 z?hPDT2M-zob*@*fZ{O*DKKPx`oP$<7uF`9!e3MO-gXXC8P;BdQH4eejY@c93D&v*k zJ0zT&Ph)&_&+~J@C-tv}p!J;2n37;SbQ;oKN8&~n)r^2lj{_DrMKQiK_EeC+Yu0+irl<2qn4t(<3fLW!S7glE1Rl za=UGHxHNYo7o(h#9P?_-PzEPEU&e6Ag((|66&%@^Z#}9AjZi3i&B$j&#zt3*yVT?U zxZ;X11l%rAQkDhPsmQ7$$_o;pnV8&Y_>xvomN0fTIj!VHqB5VqRz{%LF_e33hd4xg zQe>GDIFKV7$DWaMw3rF(n+dVBQ;$RYq9V|L)a15IMS|AXRGy9uGw%3wz;8iY=1Y4F zZ_vhE|wVEOUh-|!}%+k~32Nv1ub%Fga zZ0ynOna_CnVBpw+QaYONb!4MTdA3!Z+t)C9crwLG<4j&UaG0Kt<@LuPz5@^`dV_JK0?G?K>5?fj(~bJBs;^<`%b=K0WuQ# zk^(O3H8(vDozuRUs;s}Dbu0i?K{APj5%$6h-^41hnKKLndVn+NM#Ximqnn_@rG>XodbD)D3s>8m}%0E>}S&z*LxW2 zcvI#GZ*2FUviNYv2j1i_8!z079Tstj2&(p%wX^T6?ls=0uC#C@s&_9@W_yz^Ioc&r z1v%=QDmC4GX#FDGyz~9-tmKJ@?zh zcZ_3A?G4-+yPrA0;sZk{Z6|js%U`EPG-UC%9sGi_$oJDUT&p?!Q1pUHGjxjKM~RJ| zsuqySUCNJH0lTP8b%u~evWQ$pNQ%Ubn`*h0af7$4n*s8`sdza@h{vG|mLVPG9#7xn-nO*4%uu z4wkg~7@8YK8+t_(In@u4lnS+)7^K?xG*@>7hk|6V7C673D5Amh;Nfw|SsA*OX0Nsy zn-C;D+o3^>*%|~8e89hapgTX_fSscr_0XSzkZ*uvSdeC4ikIXT0n=P-2r&);c zB`tB#H1b4g8Rv<;=q6%a3=^=uS9auK_A>~u=+d-A(eSf9)O$%N2H`#kPKCy@We0kU zFioH#>`4>^nS_*Z6_>~Wr`9@n+z}DH0B}9hzNqNVtM_NS5|z>E*p`8N{vvY%(zc35 zT1)Fi>92eaL;x8u&_WV)yM8m>hyf3npD5ANg@$7QgG;t9bSr$gJ0u&>$1Uu$Pq6n` zsyoxgF-c{Zug)>MX4mp`?c~5Op_&9BNGhT7sBT0zeZhElkr`K zwBo_i+s+TgirrAQJN=)IBV)(3C`-9A^2!Z-T{35Qr$TPzTkr1hN8&}@H}2kzy|G8Y zyqf=IIg}VkmU8d!hgt)qvMc@<6a0Q!K`K~u)7a!^%DQiyScuT8)+W-}k zj%fOFO-5j9+KG>(6)-qITlgScC=Q>-KTdnUeviswpQ!5_RX9$k%wErwcE9;DE;fk` z7Gk^U?Kbr6)d7@G3O|SWr!r>u+qqpP4GB7qYm~`SL%FdjF(9=!y0`7+?td6`}!!{e)mYWUYI{4QMMF#AmU=a#~n+HB^ zqF=UcaAjd9vQyJin%9n*-8ND;iQ+Hwl+uhv@(M?dGwQAe?wlzQ(1IvzT2gD9IgXhn zI~cBz=chC4CV@zUl-mMil>v%>OT5^~DOR_6VB$`+ zK@RSiULp z3dL#jWoxw4Pwmjju#{R<(njiDxxAOc=$TQ}=IIDeeH?tSWIxaK6xa-4LU#8%SYDATz#W2*Q`wjz@0PKIZ8n0|ZoyrE`HZSl zV4u0*`n`kqIVFuGm8a+y(P@fkLv~aI&1@h7S%u@DAggVRxXI)~FZLS6Kim zLQPtdcwEwr2e{~lRi95Fnm8c!*c&8pg8OUVzEA#I{O#m6Z(0$li^?^T>Z+cMy_#^` zIa?rTQnS70WWn=uPJ)U~*IB;Na00K3$wk|wYC4mTM+a2K1yHs$v=su@OS7(Z#VQurez1Br(| zy(oLZ>(*3+V^BD+13CCK+0QO{`E<=5>QB|FWToVdx1R1I=yG5I%pI z^lrvqJ4jSEBxL;U+%mcC1R>ASOg}uBVGM4V0q4wa6(&vnV5I3}4^+iJjhuLgPrWJkrYj>o(+j4?WhZe0-yy zo7-n6Y{XAoAQuHd?^o&S4w6JePE1M|?wB^{3Ta$mpZ9yRu63PFUp)rt4w_sm>pmHr zv8xAGm6E*#!ge2vFEfe;Z2^;5U}J~H7qfS9wJ&c7|7IC_lHwEX25XrF*8}dn*-!67 zZ$buI_^>6*p`}&{1~5Q>ZiG<;Zg}jEyq-1>7O&S#MoWAup4vQJQJwr8V8bIMh7MDT z)O=2&jkIH|K0ht3iFU+6ao84{YFoM58=ATi!#b&zt()_HJOV3?CvcWkoJeX48`v^-vd@qMIrl+TIx Q%i7-Ox?=fb6z04C0pd5rzW@LL literal 0 HcmV?d00001 diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py index d639be1c..24982269 100644 --- a/skimage/morphology/__init__.py +++ b/skimage/morphology/__init__.py @@ -4,3 +4,4 @@ from .ccomp import label from .watershed import watershed, is_local_maximum from .skeletonize import skeletonize, medial_axis from .convex_hull import convex_hull_image +from .greyreconstruct import reconstruction diff --git a/skimage/morphology/_greyreconstruct.pyx b/skimage/morphology/_greyreconstruct.pyx new file mode 100644 index 00000000..18197812 --- /dev/null +++ b/skimage/morphology/_greyreconstruct.pyx @@ -0,0 +1,85 @@ +""" +`reconstruction_loop` originally part of CellProfiler, code licensed under both GPL and BSD licenses. + +Website: http://www.cellprofiler.org +Copyright (c) 2003-2009 Massachusetts Institute of Technology +Copyright (c) 2009-2011 Broad Institute +All rights reserved. +Original author: Lee Kamentsky + +""" + +from __future__ import division +import numpy as np + +cimport numpy as np +cimport cython + + +@cython.boundscheck(False) +def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, + negative_indices = False, + mode = 'c'] avalues, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices = False, + mode = 'c'] aprev, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices = False, + mode = 'c'] anext, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices = False, + mode = 'c'] astrides, + np.int32_t current, + int image_stride): + """The inner loop for reconstruction""" + cdef: + np.int32_t neighbor + np.uint32_t neighbor_value + np.uint32_t current_value + np.uint32_t mask_value + np.int32_t link + int i + np.int32_t nprev + np.int32_t nnext + int nstrides = astrides.shape[0] + np.uint32_t *values = (avalues.data) + np.int32_t *prev = (aprev.data) + np.int32_t *next = (anext.data) + np.int32_t *strides = (astrides.data) + + while current != -1: + if current < image_stride: + current_value = values[current] + if current_value == 0: + break + for i in range(nstrides): + neighbor = current + strides[i] + neighbor_value = values[neighbor] + # Only do neighbors less than the current value + if neighbor_value < current_value: + mask_value = values[neighbor + image_stride] + # Only do neighbors less than the mask value + if neighbor_value < mask_value: + # Raise the neighbor to the mask value if + # the mask is less than current + if mask_value < current_value: + link = neighbor + image_stride + values[neighbor] = mask_value + else: + link = current + values[neighbor] = current_value + # unlink the neighbor + nprev = prev[neighbor] + nnext = next[neighbor] + next[nprev] = nnext + if nnext != -1: + prev[nnext] = nprev + # link the neighbor after the link + nnext = next[link] + next[neighbor] = nnext + prev[neighbor] = link + if nnext >= 0: + prev[nnext] = neighbor + next[link] = neighbor + current = next[current] + diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py new file mode 100644 index 00000000..7f0def38 --- /dev/null +++ b/skimage/morphology/greyreconstruct.py @@ -0,0 +1,150 @@ +""" +`reconstruction` originally part of CellProfiler, code licensed under both GPL and BSD licenses. + +Website: http://www.cellprofiler.org +Copyright (c) 2003-2009 Massachusetts Institute of Technology +Copyright (c) 2009-2011 Broad Institute +All rights reserved. +Original author: Lee Kamentsky + +""" +import numpy as np + +from skimage.filter.rank_order import rank_order + + +def reconstruction(image, mask, selem=None, offset=None): + """Perform a morphological reconstruction of the image. + + Reconstruction requires a "seed" image and a "mask" image. The seed image + gets dilated until it is constrained by the mask. The "seed" and "mask" + images will be the minimum and maximum possible values of the reconstructed + image. + + Parameters + ---------- + image : ndarray + The seed image. + + mask : ndarray + The maximum allowed value at each point. + + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + + Returns + ------- + reconstructed : ndarray + The result of morphological reconstruction. + + Notes + ----- + The algorithm is taken from: + Robinson, "Efficient morphological reconstruction: a downhill filter", + Pattern Recognition Letters 25 (2004) 1759-1767. + + Applications for greyscale reconstruction are discussed in: + Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: + Applications and Efficient Algorithms", IEEE Transactions on Image + Processing (1993) + + Examples + -------- + Uses for greyscale reconstruction are described in Vincent (1993). For + example, let's try to extract the features of an image by subtracting a + background image created by reconstruction. + + First, create an image where the "bumps" are the features that + we want to extract: + + >>> import numpy as np + >>> from scikits.image.morphology.grey import grey_reconstruction + >>> y, x = np.mgrid[:20:0.5, :20:0.5] + >>> bumps = np.sin(x) + np.sin(y) + + To create the background image, set the mask image to the original image, + and the seed image to the original image with an intensity offset, `h`. + + >>> h = 0.3 + >>> seed = bumps - h + >>> rec = grey_reconstruction(seed, bumps) + + The resulting reconstructed image looks exactly like the original image, + but with the peaks of the bumps cut off. Subtracting this reconstructed + image from the original image leaves just the peaks of the bumps + + >>> hdome = bumps - rec + + This operation is known as the h-dome of the image, which leaves features + of height `h` in the subtracted image. The h-dome transform, and its + inverse h-basin, are analogous to the white top-hat and black top-hat + transforms, but don't require a structuring element. + + """ + assert tuple(image.shape) == tuple(mask.shape) + assert np.all(image <= mask) + try: + from ._morphrec import reconstruction_loop + except ImportError: + raise ImportError("_morphrec extension not available.") + + if selem is None: + selem = np.ones([3]*image.ndim, bool) + else: + selem = selem.copy() + + if offset == None: + assert all([d % 2 == 1 for d in selem.shape]),\ + "Footprint dimensions must all be odd" + offset = np.array([d/2 for d in selem.shape]) + # Cross out the center of the selem + selem[[slice(d,d+1) for d in offset]] = False + # + # Construct an array that's padded on the edges so we can ignore boundaries + # The array is a dstack of the image and the mask; this lets us interleave + # image and mask pixels when sorting which makes list manipulations easier + # + padding = (np.array(selem.shape)/2).astype(int) + dims = np.zeros(image.ndim+1,int) + dims[1:] = np.array(image.shape)+2*padding + dims[0] = 2 + inside_slices = [slice(p,-p) for p in padding] + values = np.ones(dims)*np.min(image) + values[[0]+inside_slices] = image + values[[1]+inside_slices] = mask + # + # Create a list of strides across the array to get the neighbors + # within a flattened array + # + value_stride = np.array(values.strides[1:]) / values.dtype.itemsize + image_stride = values.strides[0] / values.dtype.itemsize + selem_mgrid = np.mgrid[[slice(-o,d - o) + for d,o in zip(selem.shape,offset)]] + selem_offsets = selem_mgrid[:,selem].transpose() + strides = np.array([np.sum(value_stride * selem_offset) + for selem_offset in selem_offsets], np.int32) + values = values.flatten() + value_sort = np.lexsort([-values]).astype(np.int32) + # + # Make a linked list of pixels sorted by value. -1 is the list terminator. + # + prev = -np.ones(len(values), np.int32) + next = -np.ones(len(values), np.int32) + prev[value_sort[1:]] = value_sort[:-1] + next[value_sort[:-1]] = value_sort[1:] + # + # Create a rank-order value array so that the Cython inner-loop + # can operate on a uniform data type + # + values, value_map = rank_order(values) + current = value_sort[0] + + reconstruction_loop(values, prev, next, strides, current, image_stride) + # + # Reshape the values array to the shape of the padded image + # and return the unpadded portion of that result + # + values = value_map[values[:image_stride]] + values.shape = np.array(image.shape)+2*padding + return values[inside_slices] + diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py index fcf33f7f..fda13679 100644 --- a/skimage/morphology/setup.py +++ b/skimage/morphology/setup.py @@ -18,6 +18,7 @@ def configuration(parent_package='', top_path=None): cython(['_skeletonize.pyx'], working_path=base_path) cython(['_pnpoly.pyx'], working_path=base_path) cython(['_convex_hull.pyx'], working_path=base_path) + cython(['_greyreconstruct.pyx'], working_path=base_path) config.add_extension('ccomp', sources=['ccomp.c'], include_dirs=[get_numpy_include_dirs()]) @@ -31,6 +32,8 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_convex_hull', sources=['_convex_hull.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_greyreconstruct', sources=['_greyreconstruct.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/morphology/tests/test_reconstruction.py b/skimage/morphology/tests/test_reconstruction.py new file mode 100644 index 00000000..61b72033 --- /dev/null +++ b/skimage/morphology/tests/test_reconstruction.py @@ -0,0 +1,73 @@ +""" +These tests are originally part of CellProfiler, code licensed under both GPL and BSD licenses. + +Website: http://www.cellprofiler.org +Copyright (c) 2003-2009 Massachusetts Institute of Technology +Copyright (c) 2009-2011 Broad Institute +All rights reserved. +Original author: Lee Kamentsky +""" +import numpy as np + +from skimage.morphology.greyreconstruct import reconstruction + + +def test_zeros(): + """Test reconstruction with image and mask of zeros""" + assert np.all(reconstruction(np.zeros((5, 7)), np.zeros((5, 7))) == 0) + + +def test_image_equals_mask(): + """Test reconstruction where the image and mask are the same""" + assert np.all(reconstruction(np.ones((7, 5)), np.ones((7, 5))) == 1) + + +def test_image_less_than_mask(): + """Test reconstruction where the image is uniform and less than mask""" + image = np.ones((5, 5)) + mask = np.ones((5, 5)) * 2 + assert np.all(reconstruction(image,mask) == 1) + + +def test_one_image_peak(): + """Test reconstruction with one peak pixel""" + image = np.ones((5, 5)) + image[2, 2] = 2 + mask = np.ones((5, 5)) * 3 + assert np.all(reconstruction(image,mask) == 2) + + +def test_two_image_peaks(): + """Test reconstruction with two peak pixels isolated by the mask""" + image = np.array([[1, 1, 1, 1, 1, 1, 1, 1], + [1, 2, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 1, 3, 1], + [1, 1, 1, 1, 1, 1, 1, 1]]) + + mask = np.array([[4, 4, 4, 1, 1, 1, 1, 1], + [4, 4, 4, 1, 1, 1, 1, 1], + [4, 4, 4, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 4, 4, 4], + [1, 1, 1, 1, 1, 4, 4, 4], + [1, 1, 1, 1, 1, 4, 4, 4]]) + + expected = np.array([[2, 2, 2, 1, 1, 1, 1, 1], + [2, 2, 2, 1, 1, 1, 1, 1], + [2, 2, 2, 1, 1, 1, 1, 1], + [1, 1, 1, 1, 1, 3, 3, 3], + [1, 1, 1, 1, 1, 3, 3, 3], + [1, 1, 1, 1, 1, 3, 3, 3]]) + assert np.all(reconstruction(image,mask) == expected) + + +def test_zero_image_one_mask(): + """Test reconstruction with an image of all zeros and a mask that's not""" + result = reconstruction(np.zeros((10, 10)), np.ones((10, 10))) + assert np.all(result == 0) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite() From 70153abb834d79fe2ac49ed3bf2f257b2b6adbbf Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 9 Aug 2012 00:15:42 -0400 Subject: [PATCH 02/21] BUG: fix import of Cython extension. --- skimage/morphology/greyreconstruct.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index 7f0def38..bffdbe0d 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -19,16 +19,14 @@ def reconstruction(image, mask, selem=None, offset=None): Reconstruction requires a "seed" image and a "mask" image. The seed image gets dilated until it is constrained by the mask. The "seed" and "mask" images will be the minimum and maximum possible values of the reconstructed - image. + image, respectively. Parameters ---------- image : ndarray The seed image. - mask : ndarray The maximum allowed value at each point. - selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. @@ -84,9 +82,9 @@ def reconstruction(image, mask, selem=None, offset=None): assert tuple(image.shape) == tuple(mask.shape) assert np.all(image <= mask) try: - from ._morphrec import reconstruction_loop + from ._greyreconstruct import reconstruction_loop except ImportError: - raise ImportError("_morphrec extension not available.") + raise ImportError("_greyreconstruct extension not available.") if selem is None: selem = np.ones([3]*image.ndim, bool) From e4dd658daf138dec5b05527dc6076ee0dbe03f21 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 9 Aug 2012 00:21:59 -0400 Subject: [PATCH 03/21] STY: PEP8 and other clean up. --- skimage/morphology/greyreconstruct.py | 43 ++++++++----------- .../morphology/tests/test_reconstruction.py | 6 +-- 2 files changed, 22 insertions(+), 27 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index bffdbe0d..8556e444 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -87,62 +87,57 @@ def reconstruction(image, mask, selem=None, offset=None): raise ImportError("_greyreconstruct extension not available.") if selem is None: - selem = np.ones([3]*image.ndim, bool) + selem = np.ones([3]*image.ndim, dtype=bool) else: selem = selem.copy() if offset == None: - assert all([d % 2 == 1 for d in selem.shape]),\ - "Footprint dimensions must all be odd" - offset = np.array([d/2 for d in selem.shape]) + if not all([d % 2 == 1 for d in selem.shape]): + ValueError("Footprint dimensions must all be odd") + offset = np.array([d / 2 for d in selem.shape]) # Cross out the center of the selem - selem[[slice(d,d+1) for d in offset]] = False - # + selem[[slice(d, d + 1) for d in offset]] = False + # Construct an array that's padded on the edges so we can ignore boundaries # The array is a dstack of the image and the mask; this lets us interleave # image and mask pixels when sorting which makes list manipulations easier - # - padding = (np.array(selem.shape)/2).astype(int) - dims = np.zeros(image.ndim+1,int) + padding = (np.array(selem.shape) / 2).astype(int) + dims = np.zeros(image.ndim + 1, dtype=int) dims[1:] = np.array(image.shape)+2*padding dims[0] = 2 inside_slices = [slice(p,-p) for p in padding] - values = np.ones(dims)*np.min(image) - values[[0]+inside_slices] = image - values[[1]+inside_slices] = mask - # + values = np.ones(dims) * np.min(image) + values[[0] + inside_slices] = image + values[[1] + inside_slices] = mask + # Create a list of strides across the array to get the neighbors # within a flattened array - # value_stride = np.array(values.strides[1:]) / values.dtype.itemsize image_stride = values.strides[0] / values.dtype.itemsize - selem_mgrid = np.mgrid[[slice(-o,d - o) - for d,o in zip(selem.shape,offset)]] - selem_offsets = selem_mgrid[:,selem].transpose() + selem_mgrid = np.mgrid[[slice(-o, d - o) + for d, o in zip(selem.shape, offset)]] + selem_offsets = selem_mgrid[:, selem].transpose() strides = np.array([np.sum(value_stride * selem_offset) for selem_offset in selem_offsets], np.int32) values = values.flatten() value_sort = np.lexsort([-values]).astype(np.int32) - # + # Make a linked list of pixels sorted by value. -1 is the list terminator. - # prev = -np.ones(len(values), np.int32) next = -np.ones(len(values), np.int32) prev[value_sort[1:]] = value_sort[:-1] next[value_sort[:-1]] = value_sort[1:] - # + # Create a rank-order value array so that the Cython inner-loop # can operate on a uniform data type - # values, value_map = rank_order(values) current = value_sort[0] reconstruction_loop(values, prev, next, strides, current, image_stride) - # + # Reshape the values array to the shape of the padded image # and return the unpadded portion of that result - # values = value_map[values[:image_stride]] - values.shape = np.array(image.shape)+2*padding + values.shape = np.array(image.shape) + 2 * padding return values[inside_slices] diff --git a/skimage/morphology/tests/test_reconstruction.py b/skimage/morphology/tests/test_reconstruction.py index 61b72033..a3461f31 100644 --- a/skimage/morphology/tests/test_reconstruction.py +++ b/skimage/morphology/tests/test_reconstruction.py @@ -26,7 +26,7 @@ def test_image_less_than_mask(): """Test reconstruction where the image is uniform and less than mask""" image = np.ones((5, 5)) mask = np.ones((5, 5)) * 2 - assert np.all(reconstruction(image,mask) == 1) + assert np.all(reconstruction(image, mask) == 1) def test_one_image_peak(): @@ -34,7 +34,7 @@ def test_one_image_peak(): image = np.ones((5, 5)) image[2, 2] = 2 mask = np.ones((5, 5)) * 3 - assert np.all(reconstruction(image,mask) == 2) + assert np.all(reconstruction(image, mask) == 2) def test_two_image_peaks(): @@ -59,7 +59,7 @@ def test_two_image_peaks(): [1, 1, 1, 1, 1, 3, 3, 3], [1, 1, 1, 1, 1, 3, 3, 3], [1, 1, 1, 1, 1, 3, 3, 3]]) - assert np.all(reconstruction(image,mask) == expected) + assert np.all(reconstruction(image, mask) == expected) def test_zero_image_one_mask(): From e1caa9d4cdc12a8f1e4c8d0747a3a51a27bffcd9 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 9 Aug 2012 00:25:26 -0400 Subject: [PATCH 04/21] Update function names that were changed since original PR. --- doc/examples/applications/plot_peak_detection_comparison.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/examples/applications/plot_peak_detection_comparison.py b/doc/examples/applications/plot_peak_detection_comparison.py index 2cf7f933..5c0c7669 100644 --- a/doc/examples/applications/plot_peak_detection_comparison.py +++ b/doc/examples/applications/plot_peak_detection_comparison.py @@ -158,7 +158,7 @@ White tophat """ selem = morph.disk(10) img_t = np.uint8(img_smooth) -opening = morph.greyscale_open(img_t, selem) +opening = morph.opening(img_t, selem) top_hat = img_t - opening imshow(opening, vmin=0, vmax=255) @@ -177,7 +177,7 @@ plt.title("Tophat with disk of r = 10") """ selem = morph.disk(5) -top_hat = morph.greyscale_white_top_hat(img_t, selem) +top_hat = morph.white_tophat(img_t, selem) imshow(top_hat) plt.title("Tophat with disk of r = 5") @@ -187,7 +187,7 @@ plt.title("Tophat with disk of r = 5") """ selem = morph.square(20) -opening = morph.greyscale_open(img_t, selem) +opening = morph.opening(img_t, selem) # scikit's top hat filter uses uint8 and doesn't check for over(under)flow. mask = opening > img_t opening[mask] = img_t[mask] From 9cea60d817670b7439ecba97a455c319c77f565a Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 9 Aug 2012 09:49:22 -0400 Subject: [PATCH 05/21] STY: minor PEP8 change --- skimage/morphology/greyreconstruct.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index 8556e444..145d007a 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -103,9 +103,9 @@ def reconstruction(image, mask, selem=None, offset=None): # image and mask pixels when sorting which makes list manipulations easier padding = (np.array(selem.shape) / 2).astype(int) dims = np.zeros(image.ndim + 1, dtype=int) - dims[1:] = np.array(image.shape)+2*padding + dims[1:] = np.array(image.shape) + 2 * padding dims[0] = 2 - inside_slices = [slice(p,-p) for p in padding] + inside_slices = [slice(p, -p) for p in padding] values = np.ones(dims) * np.min(image) values[[0] + inside_slices] = image values[[1] + inside_slices] = mask From e6d03eaebc54baf0f61ba65fdeb0c43f7fe1b379 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 17 Aug 2012 00:42:56 -0400 Subject: [PATCH 06/21] STY: Use standard skimage data type conversion. --- .../plot_peak_detection_comparison.py | 39 ++++++++++++------- 1 file changed, 26 insertions(+), 13 deletions(-) diff --git a/doc/examples/applications/plot_peak_detection_comparison.py b/doc/examples/applications/plot_peak_detection_comparison.py index 5c0c7669..8fe26c02 100644 --- a/doc/examples/applications/plot_peak_detection_comparison.py +++ b/doc/examples/applications/plot_peak_detection_comparison.py @@ -84,23 +84,36 @@ peaks on the left to go undetected. Morphological reconstruction ============================ + +Morphological reconstruction uses two images: a seed image and a mask image. +Initially, all values of the reconstructed image start off pixel at the values +of the seed image. A seed pixels of high intensity spread outwards until it +hits the mask image (i.e. the mask value for a pixel is lower than the +high-intensity value). Note that the mask is a gray-scale image that limits the +maximum intensity at a pixel. This algorithm is clearer with pictures, which +we generate below. + +One common case uses mask and seed images derived from the same image but +shifted in intensity. Note: be careful when shifting images integer values, +since this can lead to under/overflow of values. To prevent the uint8 image we +started with from underflowing during subtraction, we first convert to float: """ -import numpy as np -img_r = np.int32(img_smooth) +from skimage import img_as_float +img_r = img_as_float(img_smooth) import skimage.morphology as morph -h = 20 +h = 0.1 rec = morph.reconstruction(img_r-h, img_r) -imshow(img_r, vmin=0, vmax=255) +imshow(img_r, vmin=0, vmax=1) plt.title("original (smoothed) image") """ .. image:: PLOT2RST.current_figure """ -imshow(rec, vmin=0, vmax=255) +imshow(rec, vmin=0, vmax=1) plt.title("background image (reconstruction)") """ @@ -156,10 +169,10 @@ features. White tophat ============ """ + selem = morph.disk(10) -img_t = np.uint8(img_smooth) -opening = morph.opening(img_t, selem) -top_hat = img_t - opening +opening = morph.opening(img_smooth, selem) +top_hat = img_smooth - opening imshow(opening, vmin=0, vmax=255) plt.title("Greyscale opening of image") @@ -177,7 +190,7 @@ plt.title("Tophat with disk of r = 10") """ selem = morph.disk(5) -top_hat = morph.white_tophat(img_t, selem) +top_hat = morph.white_tophat(img_smooth, selem) imshow(top_hat) plt.title("Tophat with disk of r = 5") @@ -187,11 +200,11 @@ plt.title("Tophat with disk of r = 5") """ selem = morph.square(20) -opening = morph.opening(img_t, selem) +opening = morph.opening(img_smooth, selem) # scikit's top hat filter uses uint8 and doesn't check for over(under)flow. -mask = opening > img_t -opening[mask] = img_t[mask] -top_hat = img_t - opening +mask = opening > img_smooth +opening[mask] = img_smooth[mask] +top_hat = img_smooth - opening imshow(opening, vmin=0, vmax=255) plt.title("Greyscale opening of image") From 33a19a8c7bb8560af53cc7b0d91e675862c6a8ca Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 17 Aug 2012 22:48:02 -0400 Subject: [PATCH 07/21] DOC: add comments to clarify algorithm --- skimage/morphology/greyreconstruct.py | 39 ++++++++++++++++----------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index 145d007a..8c12a9c5 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -1,5 +1,6 @@ """ -`reconstruction` originally part of CellProfiler, code licensed under both GPL and BSD licenses. +`reconstruction` originally part of CellProfiler, code licensed under both GPL +and BSD licenses. Website: http://www.cellprofiler.org Copyright (c) 2003-2009 Massachusetts Institute of Technology @@ -14,17 +15,18 @@ from skimage.filter.rank_order import rank_order def reconstruction(image, mask, selem=None, offset=None): - """Perform a morphological reconstruction of the image. + """Perform a morphological reconstruction of an image. - Reconstruction requires a "seed" image and a "mask" image. The seed image - gets dilated until it is constrained by the mask. The "seed" and "mask" + Reconstruction requires a "seed" image and a "mask" image. Currently, this + only implements reconstruction by dilation, such that the seed image is + dilated until it is constrained by the mask. Thus, he "seed" and "mask" images will be the minimum and maximum possible values of the reconstructed image, respectively. Parameters ---------- image : ndarray - The seed image. + The seed image; a.k.a. marker image. mask : ndarray The maximum allowed value at each point. selem : ndarray @@ -42,9 +44,13 @@ def reconstruction(image, mask, selem=None, offset=None): Pattern Recognition Letters 25 (2004) 1759-1767. Applications for greyscale reconstruction are discussed in: - Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: - Applications and Efficient Algorithms", IEEE Transactions on Image - Processing (1993) + + [1] Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: + Applications and Efficient Algorithms", IEEE Transactions on Image + Processing (1993) + + [2] Soille, P., "Morphological Image Analysis: Principles and Applications", + Chapter 6, 2nd edition (2003), ISBN 3540429883. Examples -------- @@ -98,27 +104,27 @@ def reconstruction(image, mask, selem=None, offset=None): # Cross out the center of the selem selem[[slice(d, d + 1) for d in offset]] = False - # Construct an array that's padded on the edges so we can ignore boundaries - # The array is a dstack of the image and the mask; this lets us interleave - # image and mask pixels when sorting which makes list manipulations easier + # Make padding for edges of reconstructed image so we can ignore boundaries padding = (np.array(selem.shape) / 2).astype(int) dims = np.zeros(image.ndim + 1, dtype=int) dims[1:] = np.array(image.shape) + 2 * padding dims[0] = 2 inside_slices = [slice(p, -p) for p in padding] + # Set padded region to minimum image intensity and mask along first axis so + # we can interleave image and mask pixels when sorting. values = np.ones(dims) * np.min(image) values[[0] + inside_slices] = image values[[1] + inside_slices] = mask - # Create a list of strides across the array to get the neighbors - # within a flattened array + # Create a list of strides across the array to get the neighbors within + # a flattened array value_stride = np.array(values.strides[1:]) / values.dtype.itemsize image_stride = values.strides[0] / values.dtype.itemsize selem_mgrid = np.mgrid[[slice(-o, d - o) for d, o in zip(selem.shape, offset)]] selem_offsets = selem_mgrid[:, selem].transpose() - strides = np.array([np.sum(value_stride * selem_offset) - for selem_offset in selem_offsets], np.int32) + nb_strides = np.array([np.sum(value_stride * selem_offset) + for selem_offset in selem_offsets], np.int32) values = values.flatten() value_sort = np.lexsort([-values]).astype(np.int32) @@ -130,10 +136,11 @@ def reconstruction(image, mask, selem=None, offset=None): # Create a rank-order value array so that the Cython inner-loop # can operate on a uniform data type + # fragile: `reconstruction_loop` needs 'uint32' conversion by `rank_order` values, value_map = rank_order(values) current = value_sort[0] - reconstruction_loop(values, prev, next, strides, current, image_stride) + reconstruction_loop(values, prev, next, nb_strides, current, image_stride) # Reshape the values array to the shape of the padded image # and return the unpadded portion of that result From 29c84a8a7bea09dae0dc4fa419ed0bd89436e127 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 17 Aug 2012 23:02:08 -0400 Subject: [PATCH 08/21] STY: Rename returned image to distinguish input from output. --- skimage/morphology/greyreconstruct.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index 8c12a9c5..c68e6e24 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -136,15 +136,13 @@ def reconstruction(image, mask, selem=None, offset=None): # Create a rank-order value array so that the Cython inner-loop # can operate on a uniform data type - # fragile: `reconstruction_loop` needs 'uint32' conversion by `rank_order` - values, value_map = rank_order(values) + rec_img, value_map = rank_order(values) current = value_sort[0] - reconstruction_loop(values, prev, next, nb_strides, current, image_stride) + reconstruction_loop(rec_img, prev, next, nb_strides, current, image_stride) - # Reshape the values array to the shape of the padded image - # and return the unpadded portion of that result - values = value_map[values[:image_stride]] - values.shape = np.array(image.shape) + 2 * padding - return values[inside_slices] + # Reshape reconstructed image to original image shape and remove padding. + rec_img = value_map[rec_img[:image_stride]] + rec_img.shape = np.array(image.shape) + 2 * padding + return rec_img[inside_slices] From 969772c036e80930d5f36ee46e52c6e1085bc908 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 18 Aug 2012 17:59:44 -0400 Subject: [PATCH 09/21] ENH: Add reconstruction by erosion. --- skimage/morphology/greyreconstruct.py | 54 ++++++++++++++++++--------- 1 file changed, 36 insertions(+), 18 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index c68e6e24..d1434246 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -14,14 +14,12 @@ import numpy as np from skimage.filter.rank_order import rank_order -def reconstruction(image, mask, selem=None, offset=None): +def reconstruction(image, mask, selem=None, offset=None, method='dilation'): """Perform a morphological reconstruction of an image. - Reconstruction requires a "seed" image and a "mask" image. Currently, this - only implements reconstruction by dilation, such that the seed image is - dilated until it is constrained by the mask. Thus, he "seed" and "mask" - images will be the minimum and maximum possible values of the reconstructed - image, respectively. + Reconstruction requires a "seed" image and a "mask" image of equal shape. + These images set the minimum and maximum possible values of the + reconstructed image. Parameters ---------- @@ -31,6 +29,11 @@ def reconstruction(image, mask, selem=None, offset=None): The maximum allowed value at each point. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. + method : {'dilation'|'erosion'} + Perform reconstruction by dilation or erosion. In dilation (erosion), + the seed image is dilated (eroded) until limited by the mask image. + For dilation, each seed value must be less than or equal to the + corresponding mask value; for erosion, the reverse is true. Returns ------- @@ -48,7 +51,6 @@ def reconstruction(image, mask, selem=None, offset=None): [1] Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: Applications and Efficient Algorithms", IEEE Transactions on Image Processing (1993) - [2] Soille, P., "Morphological Image Analysis: Principles and Applications", Chapter 6, 2nd edition (2003), ISBN 3540429883. @@ -62,7 +64,7 @@ def reconstruction(image, mask, selem=None, offset=None): we want to extract: >>> import numpy as np - >>> from scikits.image.morphology.grey import grey_reconstruction + >>> from skimage.morphology import reconstruction >>> y, x = np.mgrid[:20:0.5, :20:0.5] >>> bumps = np.sin(x) + np.sin(y) @@ -71,7 +73,7 @@ def reconstruction(image, mask, selem=None, offset=None): >>> h = 0.3 >>> seed = bumps - h - >>> rec = grey_reconstruction(seed, bumps) + >>> rec = reconstruction(seed, bumps) The resulting reconstructed image looks exactly like the original image, but with the peaks of the bumps cut off. Subtracting this reconstructed @@ -86,7 +88,12 @@ def reconstruction(image, mask, selem=None, offset=None): """ assert tuple(image.shape) == tuple(mask.shape) - assert np.all(image <= mask) + if method == 'dilation' and np.any(image > mask): + raise ValueError("Intensity of seed image must be less than that " + "of the mask image for reconstruction by dilation.") + elif method == 'erosion' and np.any(image < mask): + raise ValueError("Intensity of seed image must be greater than that " + "of the mask image for reconstruction by erosion.") try: from ._greyreconstruct import reconstruction_loop except ImportError: @@ -112,7 +119,11 @@ def reconstruction(image, mask, selem=None, offset=None): inside_slices = [slice(p, -p) for p in padding] # Set padded region to minimum image intensity and mask along first axis so # we can interleave image and mask pixels when sorting. - values = np.ones(dims) * np.min(image) + if method == 'dilation': + pad_value = np.min(image) + elif method == 'erosion': + pad_value = np.max(image) + values = np.ones(dims) * pad_value values[[0] + inside_slices] = image values[[1] + inside_slices] = mask @@ -126,23 +137,30 @@ def reconstruction(image, mask, selem=None, offset=None): nb_strides = np.array([np.sum(value_stride * selem_offset) for selem_offset in selem_offsets], np.int32) values = values.flatten() - value_sort = np.lexsort([-values]).astype(np.int32) + index_sorted = np.argsort(-values).astype(np.int32) + if method == 'erosion': + index_sorted = index_sorted[::-1] # Make a linked list of pixels sorted by value. -1 is the list terminator. prev = -np.ones(len(values), np.int32) next = -np.ones(len(values), np.int32) - prev[value_sort[1:]] = value_sort[:-1] - next[value_sort[:-1]] = value_sort[1:] + prev[index_sorted[1:]] = index_sorted[:-1] + next[index_sorted[:-1]] = index_sorted[1:] # Create a rank-order value array so that the Cython inner-loop # can operate on a uniform data type - rec_img, value_map = rank_order(values) - current = value_sort[0] + if method == 'dilation': + value_rank, value_map = rank_order(values) + elif method == 'erosion': + value_rank, value_map = rank_order(-values) + value_map = -value_map + current = index_sorted[0] - reconstruction_loop(rec_img, prev, next, nb_strides, current, image_stride) + reconstruction_loop(value_rank, prev, next, nb_strides, current, + image_stride) # Reshape reconstructed image to original image shape and remove padding. - rec_img = value_map[rec_img[:image_stride]] + rec_img = value_map[value_rank[:image_stride]] rec_img.shape = np.array(image.shape) + 2 * padding return rec_img[inside_slices] From 7a56a7f35ee41813d42d6f1bfb3447fa0c60d569 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 18 Aug 2012 18:04:01 -0400 Subject: [PATCH 10/21] DOC: Clarify code comments and docstring --- skimage/morphology/_greyreconstruct.pyx | 61 +++++++++++++++++-------- 1 file changed, 41 insertions(+), 20 deletions(-) diff --git a/skimage/morphology/_greyreconstruct.pyx b/skimage/morphology/_greyreconstruct.pyx index 18197812..a59143f5 100644 --- a/skimage/morphology/_greyreconstruct.pyx +++ b/skimage/morphology/_greyreconstruct.pyx @@ -18,26 +18,47 @@ cimport cython @cython.boundscheck(False) def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, - negative_indices = False, - mode = 'c'] avalues, + negative_indices=False, mode='c'] avalues, np.ndarray[dtype=np.int32_t, ndim=1, - negative_indices = False, - mode = 'c'] aprev, + negative_indices=False, mode='c'] aprev, np.ndarray[dtype=np.int32_t, ndim=1, - negative_indices = False, - mode = 'c'] anext, + negative_indices=False, mode='c'] anext, np.ndarray[dtype=np.int32_t, ndim=1, - negative_indices = False, - mode = 'c'] astrides, + negative_indices=False, mode='c'] astrides, np.int32_t current, int image_stride): - """The inner loop for reconstruction""" + """The inner loop for reconstruction. + + This algorithm uses the rank-order of pixels. If low intensity pixels have + a low rank and high intensity pixels have a high rank, then this loop + performs reconstruction by dilation. If this ranking is reversed, the + result is reconstruction by erosion. + + For each pixel in the seed image, check its neighbors. If its neighbor's + rank is below that of the current pixel, replace the neighbor's rank with + the rank of the current pixel. This dilation is limited by the mask, i.e. + the rank at each pixel cannot exceed the mask as that pixel. + + Parameters + ---------- + avalues : array + The rank order of the flattened seed and mask images. + aprev, anext: arrays + Indices of previous and next pixels in rank sorted order. + astrides : array + Strides to neighbors of the current pixel. + current : int + Index of lowest-ranked pixel used as starting point in reconstruction + loop. + image_stride : int + Stride between seed image and mask image in `avalues`. + """ cdef: np.int32_t neighbor np.uint32_t neighbor_value np.uint32_t current_value np.uint32_t mask_value - np.int32_t link + np.int32_t current_link int i np.int32_t nprev np.int32_t nnext @@ -55,18 +76,18 @@ def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, for i in range(nstrides): neighbor = current + strides[i] neighbor_value = values[neighbor] - # Only do neighbors less than the current value + # Only propagate neighbors ranked below the current rank if neighbor_value < current_value: mask_value = values[neighbor + image_stride] - # Only do neighbors less than the mask value + # Only propagate neighbors ranked below the mask rank if neighbor_value < mask_value: - # Raise the neighbor to the mask value if - # the mask is less than current + # Raise the neighbor to the mask rank if + # the mask ranked below the current rank if mask_value < current_value: - link = neighbor + image_stride + current_link = neighbor + image_stride values[neighbor] = mask_value else: - link = current + current_link = current values[neighbor] = current_value # unlink the neighbor nprev = prev[neighbor] @@ -74,12 +95,12 @@ def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, next[nprev] = nnext if nnext != -1: prev[nnext] = nprev - # link the neighbor after the link - nnext = next[link] + # link to the neighbor after the current link + nnext = next[current_link] next[neighbor] = nnext - prev[neighbor] = link + prev[neighbor] = current_link if nnext >= 0: prev[nnext] = neighbor - next[link] = neighbor + next[current_link] = neighbor current = next[current] From ab7626da3deae657e3d03af861cfd9187fea2e20 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 18 Aug 2012 19:21:55 -0400 Subject: [PATCH 11/21] STY: Rename variables for clarity. In particular, it wasn't clear whether `image` was the seed image or the mask image rnd `values` was used to refer to both image intensity values and their rank-order. --- skimage/morphology/_greyreconstruct.pyx | 62 ++++++++++++------------- skimage/morphology/greyreconstruct.py | 61 ++++++++++++------------ 2 files changed, 62 insertions(+), 61 deletions(-) diff --git a/skimage/morphology/_greyreconstruct.pyx b/skimage/morphology/_greyreconstruct.pyx index a59143f5..dff670bb 100644 --- a/skimage/morphology/_greyreconstruct.pyx +++ b/skimage/morphology/_greyreconstruct.pyx @@ -18,14 +18,14 @@ cimport cython @cython.boundscheck(False) def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, - negative_indices=False, mode='c'] avalues, + negative_indices=False, mode='c'] rank_array, np.ndarray[dtype=np.int32_t, ndim=1, negative_indices=False, mode='c'] aprev, np.ndarray[dtype=np.int32_t, ndim=1, negative_indices=False, mode='c'] anext, np.ndarray[dtype=np.int32_t, ndim=1, negative_indices=False, mode='c'] astrides, - np.int32_t current, + np.int32_t current_idx, int image_stride): """The inner loop for reconstruction. @@ -41,66 +41,66 @@ def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, Parameters ---------- - avalues : array + rank_array : array The rank order of the flattened seed and mask images. aprev, anext: arrays Indices of previous and next pixels in rank sorted order. astrides : array Strides to neighbors of the current pixel. - current : int + current_idx : int Index of lowest-ranked pixel used as starting point in reconstruction loop. image_stride : int - Stride between seed image and mask image in `avalues`. + Stride between seed image and mask image in `rank_array`. """ cdef: - np.int32_t neighbor - np.uint32_t neighbor_value - np.uint32_t current_value - np.uint32_t mask_value + np.int32_t neighbor_idx + np.uint32_t neighbor_rank + np.uint32_t current_rank + np.uint32_t mask_rank np.int32_t current_link int i np.int32_t nprev np.int32_t nnext int nstrides = astrides.shape[0] - np.uint32_t *values = (avalues.data) + np.uint32_t *ranks = (rank_array.data) np.int32_t *prev = (aprev.data) np.int32_t *next = (anext.data) np.int32_t *strides = (astrides.data) - while current != -1: - if current < image_stride: - current_value = values[current] - if current_value == 0: + while current_idx != -1: + if current_idx < image_stride: + current_rank = ranks[current_idx] + if current_rank == 0: break for i in range(nstrides): - neighbor = current + strides[i] - neighbor_value = values[neighbor] + neighbor_idx = current_idx + strides[i] + neighbor_rank = ranks[neighbor_idx] # Only propagate neighbors ranked below the current rank - if neighbor_value < current_value: - mask_value = values[neighbor + image_stride] + if neighbor_rank < current_rank: + mask_rank = ranks[neighbor_idx + image_stride] # Only propagate neighbors ranked below the mask rank - if neighbor_value < mask_value: + if neighbor_rank < mask_rank: # Raise the neighbor to the mask rank if # the mask ranked below the current rank - if mask_value < current_value: - current_link = neighbor + image_stride - values[neighbor] = mask_value + if mask_rank < current_rank: + current_link = neighbor_idx + image_stride + ranks[neighbor_idx] = mask_rank else: - current_link = current - values[neighbor] = current_value + current_link = current_idx + ranks[neighbor_idx] = current_rank # unlink the neighbor - nprev = prev[neighbor] - nnext = next[neighbor] + nprev = prev[neighbor_idx] + nnext = next[neighbor_idx] next[nprev] = nnext if nnext != -1: prev[nnext] = nprev # link to the neighbor after the current link nnext = next[current_link] - next[neighbor] = nnext - prev[neighbor] = current_link + next[neighbor_idx] = nnext + prev[neighbor_idx] = current_link if nnext >= 0: - prev[nnext] = neighbor - next[current_link] = neighbor - current = next[current] + prev[nnext] = neighbor_idx + next[current_link] = neighbor_idx + current_idx = next[current_idx] diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index d1434246..042d86f2 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -1,6 +1,6 @@ """ -`reconstruction` originally part of CellProfiler, code licensed under both GPL -and BSD licenses. +This morphological reconstruction routine was adapted from CellProfiler, code +licensed under both GPL and BSD licenses. Website: http://www.cellprofiler.org Copyright (c) 2003-2009 Massachusetts Institute of Technology @@ -14,7 +14,7 @@ import numpy as np from skimage.filter.rank_order import rank_order -def reconstruction(image, mask, selem=None, offset=None, method='dilation'): +def reconstruction(seed, mask, selem=None, offset=None, method='dilation'): """Perform a morphological reconstruction of an image. Reconstruction requires a "seed" image and a "mask" image of equal shape. @@ -23,7 +23,7 @@ def reconstruction(image, mask, selem=None, offset=None, method='dilation'): Parameters ---------- - image : ndarray + seed : ndarray The seed image; a.k.a. marker image. mask : ndarray The maximum allowed value at each point. @@ -87,11 +87,11 @@ def reconstruction(image, mask, selem=None, offset=None, method='dilation'): transforms, but don't require a structuring element. """ - assert tuple(image.shape) == tuple(mask.shape) - if method == 'dilation' and np.any(image > mask): + assert tuple(seed.shape) == tuple(mask.shape) + if method == 'dilation' and np.any(seed > mask): raise ValueError("Intensity of seed image must be less than that " "of the mask image for reconstruction by dilation.") - elif method == 'erosion' and np.any(image < mask): + elif method == 'erosion' and np.any(seed < mask): raise ValueError("Intensity of seed image must be greater than that " "of the mask image for reconstruction by erosion.") try: @@ -100,7 +100,7 @@ def reconstruction(image, mask, selem=None, offset=None, method='dilation'): raise ImportError("_greyreconstruct extension not available.") if selem is None: - selem = np.ones([3]*image.ndim, dtype=bool) + selem = np.ones([3] * seed.ndim, dtype=bool) else: selem = selem.copy() @@ -113,54 +113,55 @@ def reconstruction(image, mask, selem=None, offset=None, method='dilation'): # Make padding for edges of reconstructed image so we can ignore boundaries padding = (np.array(selem.shape) / 2).astype(int) - dims = np.zeros(image.ndim + 1, dtype=int) - dims[1:] = np.array(image.shape) + 2 * padding + dims = np.zeros(seed.ndim + 1, dtype=int) + dims[1:] = np.array(seed.shape) + 2 * padding dims[0] = 2 inside_slices = [slice(p, -p) for p in padding] # Set padded region to minimum image intensity and mask along first axis so # we can interleave image and mask pixels when sorting. if method == 'dilation': - pad_value = np.min(image) + pad_value = np.min(seed) elif method == 'erosion': - pad_value = np.max(image) - values = np.ones(dims) * pad_value - values[[0] + inside_slices] = image - values[[1] + inside_slices] = mask + pad_value = np.max(seed) + images = np.ones(dims) * pad_value + images[[0] + inside_slices] = seed + images[[1] + inside_slices] = mask # Create a list of strides across the array to get the neighbors within # a flattened array - value_stride = np.array(values.strides[1:]) / values.dtype.itemsize - image_stride = values.strides[0] / values.dtype.itemsize + value_stride = np.array(images.strides[1:]) / images.dtype.itemsize + image_stride = images.strides[0] / images.dtype.itemsize selem_mgrid = np.mgrid[[slice(-o, d - o) for d, o in zip(selem.shape, offset)]] selem_offsets = selem_mgrid[:, selem].transpose() nb_strides = np.array([np.sum(value_stride * selem_offset) for selem_offset in selem_offsets], np.int32) - values = values.flatten() - index_sorted = np.argsort(-values).astype(np.int32) - if method == 'erosion': + + images = images.flatten() + + # Erosion goes smallest to largest; dilation goes largest to smallest. + index_sorted = np.argsort(images).astype(np.int32) + if method == 'dilation': index_sorted = index_sorted[::-1] # Make a linked list of pixels sorted by value. -1 is the list terminator. - prev = -np.ones(len(values), np.int32) - next = -np.ones(len(values), np.int32) + prev = -np.ones(len(images), np.int32) + next = -np.ones(len(images), np.int32) prev[index_sorted[1:]] = index_sorted[:-1] next[index_sorted[:-1]] = index_sorted[1:] - # Create a rank-order value array so that the Cython inner-loop - # can operate on a uniform data type + # Cython inner-loop compares the rank of pixel values. if method == 'dilation': - value_rank, value_map = rank_order(values) + value_rank, value_map = rank_order(images) elif method == 'erosion': - value_rank, value_map = rank_order(-values) + value_rank, value_map = rank_order(-images) value_map = -value_map - current = index_sorted[0] - reconstruction_loop(value_rank, prev, next, nb_strides, current, - image_stride) + start = index_sorted[0] + reconstruction_loop(value_rank, prev, next, nb_strides, start, image_stride) # Reshape reconstructed image to original image shape and remove padding. rec_img = value_map[value_rank[:image_stride]] - rec_img.shape = np.array(image.shape) + 2 * padding + rec_img.shape = np.array(seed.shape) + 2 * padding return rec_img[inside_slices] From 79fca0e20d0dc10b10f2de2116f96cd064e83f36 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 18 Aug 2012 21:54:06 -0400 Subject: [PATCH 12/21] DOC: Reorder docstring sections. --- skimage/morphology/greyreconstruct.py | 43 +++++++++++---------------- 1 file changed, 18 insertions(+), 25 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index 042d86f2..04ea54a3 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -40,29 +40,11 @@ def reconstruction(seed, mask, selem=None, offset=None, method='dilation'): reconstructed : ndarray The result of morphological reconstruction. - Notes - ----- - The algorithm is taken from: - Robinson, "Efficient morphological reconstruction: a downhill filter", - Pattern Recognition Letters 25 (2004) 1759-1767. - - Applications for greyscale reconstruction are discussed in: - - [1] Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: - Applications and Efficient Algorithms", IEEE Transactions on Image - Processing (1993) - [2] Soille, P., "Morphological Image Analysis: Principles and Applications", - Chapter 6, 2nd edition (2003), ISBN 3540429883. - Examples -------- - Uses for greyscale reconstruction are described in Vincent (1993). For - example, let's try to extract the features of an image by subtracting a + Here, we try to extract the bright features of an image by subtracting a background image created by reconstruction. - First, create an image where the "bumps" are the features that - we want to extract: - >>> import numpy as np >>> from skimage.morphology import reconstruction >>> y, x = np.mgrid[:20:0.5, :20:0.5] @@ -73,19 +55,30 @@ def reconstruction(seed, mask, selem=None, offset=None, method='dilation'): >>> h = 0.3 >>> seed = bumps - h - >>> rec = reconstruction(seed, bumps) + >>> background = reconstruction(seed, bumps) The resulting reconstructed image looks exactly like the original image, but with the peaks of the bumps cut off. Subtracting this reconstructed image from the original image leaves just the peaks of the bumps - >>> hdome = bumps - rec + >>> hdome = bumps - background - This operation is known as the h-dome of the image, which leaves features - of height `h` in the subtracted image. The h-dome transform, and its - inverse h-basin, are analogous to the white top-hat and black top-hat - transforms, but don't require a structuring element. + This operation is known as the h-dome of the image and leaves features + of height `h` in the subtracted image. + Notes + ----- + The algorithm is taken from: + [1] Robinson, "Efficient morphological reconstruction: a downhill filter", + Pattern Recognition Letters 25 (2004) 1759-1767. + + Applications for greyscale reconstruction are discussed in: + + [2] Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: + Applications and Efficient Algorithms", IEEE Transactions on Image + Processing (1993) + [3] Soille, P., "Morphological Image Analysis: Principles and Applications", + Chapter 6, 2nd edition (2003), ISBN 3540429883. """ assert tuple(seed.shape) == tuple(mask.shape) if method == 'dilation' and np.any(seed > mask): From 2e87dd7a3ca00f7002e86d30d6a32d6b0f18a363 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 18 Aug 2012 23:06:08 -0400 Subject: [PATCH 13/21] ENH: Add examples of morphological reconstruction. --- doc/examples/plot_fill_holes.py | 70 +++++++++++++++++++++++++++++++++ doc/examples/plot_find_spots.py | 69 ++++++++++++++++++++++++++++++++ 2 files changed, 139 insertions(+) create mode 100644 doc/examples/plot_fill_holes.py create mode 100644 doc/examples/plot_find_spots.py diff --git a/doc/examples/plot_fill_holes.py b/doc/examples/plot_fill_holes.py new file mode 100644 index 00000000..7257d2bb --- /dev/null +++ b/doc/examples/plot_fill_holes.py @@ -0,0 +1,70 @@ +""" +========== +Fill holes +========== + +In this example, we fill holes (i.e. isolated, dark spots) in an image using +morphological reconstruction by erosion. Erosion expands the minimal values of +the seed image until it encounters a mask image. Thus, the seed image and mask +image represent the maximum and minimum possible values of the reconstructed +image. + +We start with an image containing both peaks and holes: +""" +import matplotlib.pyplot as plt + +from skimage import data +from skimage.exposure import rescale_intensity + +image = data.moon() +# Rescale image intensity so that we can see dim features. +image = rescale_intensity(image, in_range=(50, 200)) + +# convenience function for plotting images +def imshow(image, **kwargs): + plt.figure(figsize=(5, 4)) + plt.imshow(image, **kwargs) + plt.axis('off') + +imshow(image) +plt.title('original image') + +""" +.. image:: PLOT2RST.current_figure + +Now we need to create the seed image, where the minima represent the starting +points for erosion. To fill holes, we initialize the seed image to the maximum +value of the original image. Along the borders, however, we use the original +values of the image. These border pixels will be the starting points for the +erosion process. We then limit the erosion by setting the mask to the values +of the original image. +""" + +import numpy as np +from skimage.morphology import reconstruction + +seed = np.copy(image) +seed[1:-1, 1:-1] = image.max() +mask = image + +filled = reconstruction(seed, mask, method='erosion') + +imshow(filled, vmin=image.min(), vmax=image.max()) +plt.title('after filling holes') + +""" +.. image:: PLOT2RST.current_figure + +As shown above, eroding inward from the edges removes holes, since (by +definition) holes are surrounded by pixels of brighter value. Finally, we can +isolate the dark regions by subtracting the reconstructed image from the +original image. +""" + +imshow(image - filled) +plt.title('dark holes') +plt.show() + +""" +.. image:: PLOT2RST.current_figure +""" diff --git a/doc/examples/plot_find_spots.py b/doc/examples/plot_find_spots.py new file mode 100644 index 00000000..ff6e2196 --- /dev/null +++ b/doc/examples/plot_find_spots.py @@ -0,0 +1,69 @@ +""" +========== +Find spots +========== + +In this example, we find bright spots in an image using morphological +reconstruction by dilation. Dilation expands the maximal values of the seed +image until it encounters a mask image. Thus, the seed image and mask image +represent the minimum and maximum possible values of the reconstructed image. + +We start with an image containing both peaks and holes: +""" +import matplotlib.pyplot as plt + +from skimage import data +from skimage.exposure import rescale_intensity + +image = data.moon() +# Rescale image intensity so that we can see dim features. +image = rescale_intensity(image, in_range=(50, 200)) + +# convenience function for plotting images +def imshow(image, **kwargs): + plt.figure(figsize=(5, 4)) + plt.imshow(image, **kwargs) + plt.axis('off') + +imshow(image) +plt.title('original image') + +""" +.. image:: PLOT2RST.current_figure + +Now we need to create the seed image, where the maxima represent the starting +points for dilation. To find spots, we initialize the seed image to the minimum +value of the original image. Along the borders, however, we use the original +values of the image. These border pixels will be the starting points for the +dilation process. We then limit the dilation by setting the mask to the values +of the original image. +""" + +import numpy as np +from skimage.morphology import reconstruction + +seed = np.copy(image) +seed[1:-1, 1:-1] = image.min() +mask = image + +rec = reconstruction(seed, mask, method='dilation') + +imshow(rec, vmin=image.min(), vmax=image.max()) +plt.title('') + +""" +.. image:: PLOT2RST.current_figure + +As shown above, dilating inward from the edges removes peaks, since (by +definition) peaks are surrounded by pixels of darker value. Finally, we can +isolate the bright spots by subtracting the reconstructed image from the +original image. +""" + +imshow(image - rec) +plt.title('"holes"') +plt.show() + +""" +.. image:: PLOT2RST.current_figure +""" From 3ad1ed3a285cff9abfe29a50934ef0ffc2bf6d13 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 18 Aug 2012 23:09:47 -0400 Subject: [PATCH 14/21] DOC: Remove peak detection tutorial. The tutorial needs a lot of work and isn't a crucial part of this PR. Note: tutorial saved in a separate branch. --- .../plot_peak_detection_comparison.py | 231 ------------------ skimage/data/noisy_circles.jpg | Bin 26206 -> 0 bytes 2 files changed, 231 deletions(-) delete mode 100644 doc/examples/applications/plot_peak_detection_comparison.py delete mode 100644 skimage/data/noisy_circles.jpg diff --git a/doc/examples/applications/plot_peak_detection_comparison.py b/doc/examples/applications/plot_peak_detection_comparison.py deleted file mode 100644 index 8fe26c02..00000000 --- a/doc/examples/applications/plot_peak_detection_comparison.py +++ /dev/null @@ -1,231 +0,0 @@ -""" -============== -Peak detection -============== - -Peak detection (a.k.a. spot detection or particle detection) is a common -image analysis step. For example, it's used to detect tracer particles in a -flow for particle image velocimetry (i.e. PIV) and to identify features in -the Hough transform. - -To simplify plotting code, let's define a simple function that creates a new -figure on each call, and removes tick labels. -""" - -import matplotlib.pyplot as plt - -plt.rcParams['axes.titlesize'] = 10 -plt.rcParams['font.size'] = 10 -def imshow(image, **kwargs): - plt.figure(figsize=(2.5, 2.5)) - plt.imshow(image, **kwargs) - plt.axis('off') - -""" -To explore different peak detection techniques, we use an image of circles with -added noise: -""" - -from skimage import data -img = data.load('noisy_circles.jpg') -imshow(img) - -""" - -.. image:: PLOT2RST.current_figure - -This image is noisy and has uneven background illumination. The peaks in the -image, while readily identified by eye, can be tricky to find algorithmically. -The first thing we need to do is remove the high-frequency noise; this can -be done with a simple Gaussian filter. -""" - -import scipy.ndimage as ndimg -img_smooth = ndimg.gaussian_filter(img, 3) - -imshow(img_smooth) - -""" - -.. image:: PLOT2RST.current_figure - -Thresholding -============ - -One way to extract the background is to threshold the image. -""" - -thresh_value = 100 -background = img_smooth.copy() -background[img_smooth > thresh_value] = 0 -peaks = img_smooth - background - -""" -Here, all pixels values below the threshold value are subtracted from the -image. The resulting background image and the extracted peaks are shown below. -""" - -imshow(background, vmin=0, vmax=255) -plt.title("background image (thresholding)") - -""" -.. image:: PLOT2RST.current_figure -""" - -imshow(peaks, vmin=0, vmax=255) -plt.title("peaks (thresholding)") - -""" -.. image:: PLOT2RST.current_figure - -Because of uneven illumination, peaks on the right bleed into each other. -Increasing the threshold will fix this problem, but it will also cause some -peaks on the left to go undetected. - -Morphological reconstruction -============================ - -Morphological reconstruction uses two images: a seed image and a mask image. -Initially, all values of the reconstructed image start off pixel at the values -of the seed image. A seed pixels of high intensity spread outwards until it -hits the mask image (i.e. the mask value for a pixel is lower than the -high-intensity value). Note that the mask is a gray-scale image that limits the -maximum intensity at a pixel. This algorithm is clearer with pictures, which -we generate below. - -One common case uses mask and seed images derived from the same image but -shifted in intensity. Note: be careful when shifting images integer values, -since this can lead to under/overflow of values. To prevent the uint8 image we -started with from underflowing during subtraction, we first convert to float: -""" - -from skimage import img_as_float -img_r = img_as_float(img_smooth) - -import skimage.morphology as morph -h = 0.1 -rec = morph.reconstruction(img_r-h, img_r) - -imshow(img_r, vmin=0, vmax=1) -plt.title("original (smoothed) image") - -""" -.. image:: PLOT2RST.current_figure -""" - -imshow(rec, vmin=0, vmax=1) -plt.title("background image (reconstruction)") - -""" -.. image:: PLOT2RST.current_figure - -This reconstructed image looks pretty much like the original, except that the -peaks in the image are truncated. The reconstructed image can then be -subtracted from the original image to reveal the peaks of the image. -""" - -imshow(img_r-rec) -plt.title("h-dome of image") - -""" -.. image:: PLOT2RST.current_figure - -The result is known as the h-dome transformation [2]_, which extracts peaks of -height `h` from the original image. To better understand what's going on, -let's take a 1D slice along the middle of the image (cutting through peaks in -the image). -""" - -img_slice = img_r[99:100, :] -rec_slice = morph.reconstruction(img_slice-h, img_slice) - -""" -Plotting the reconstructed image (slice) next to the original image and the -seed image shed light on the reconstruction process -""" -plt.figure(figsize=(4, 3)) -plt.plot(img_slice[0], 'k', label='original image') -plt.plot(img_slice[0]-h, '0.5', label='seed image') -plt.plot(rec_slice[0], 'r', label='reconstructed') -plt.title("image slice") -plt.xlabel('x') -plt.ylabel('intensity') -plt.legend() - -""" -.. image:: PLOT2RST.current_figure - -Here, you see that morphological reconstruction dilates the seed image (i.e. -the `h`-shifted image) until it intersects the mask (original image). Note that -the peaks in the original image have very different intensity values (e.g. the -peak at x=200 and x=100 differ by about 80). Subtracting the reconstructed -image from the original image gives peaks of roughly equal intensity. Thus, the -h-dome transformation is quiet effective at removing uneven, dark backgrounds -from bright features. The inverse operation---the h-basin -transformation---should be used when removing bright backgrounds from dark -features. - - -White tophat -============ -""" - -selem = morph.disk(10) -opening = morph.opening(img_smooth, selem) -top_hat = img_smooth - opening - -imshow(opening, vmin=0, vmax=255) -plt.title("Greyscale opening of image") - -""" -.. image:: PLOT2RST.current_figure -""" - - -imshow(top_hat) -plt.title("Tophat with disk of r = 10") - -""" -.. image:: PLOT2RST.current_figure -""" - -selem = morph.disk(5) -top_hat = morph.white_tophat(img_smooth, selem) - -imshow(top_hat) -plt.title("Tophat with disk of r = 5") - -""" -.. image:: PLOT2RST.current_figure -""" - -selem = morph.square(20) -opening = morph.opening(img_smooth, selem) -# scikit's top hat filter uses uint8 and doesn't check for over(under)flow. -mask = opening > img_smooth -opening[mask] = img_smooth[mask] -top_hat = img_smooth - opening - -imshow(opening, vmin=0, vmax=255) -plt.title("Greyscale opening of image") - -""" -.. image:: PLOT2RST.current_figure -""" - -imshow(top_hat) -plt.title("Tophat with square of w = 10") - -plt.show() - -""" -.. image:: PLOT2RST.current_figure - - -References -========== - -.. [1] Crocker and Grier, Journal of Colloid and Interface Science (1996) -.. [2] Vincent, L., IEEE Transactions on Image Processing (1993) - -""" diff --git a/skimage/data/noisy_circles.jpg b/skimage/data/noisy_circles.jpg deleted file mode 100644 index a0f49d6767ca4ca3cfbef6c54c36039223d243ad..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 26206 zcmeIa2UJr_*FSthBIOF9T@^x$B*3+R^b$&lAP{LPC@7*-A#_M+B81{axJp+6m7*X; znhHo$1BwDFy(!WJq?&+$L|XDc!E*2OywAVA-}=`2*7`hI=bR~f_UxJ2)ABpBnQxgN zgclA45zYaCsVT4*002&a8v+1O&;X?uv z;rsE%Ie9xcBagUxle~$p-abe@3l*fPy(7VkP7TRy-koJGfc4_pMeIXOAGIJvpF zxcRsLww-?mFE=;u4k1Co9fE>F{M=jD_t!Vizds>7+qd(;dHCRPzTe?+`0uO-{P%B3 z`2SBBFzW$+P9Ph|hC(C(Hhu_{AHp003;7$f5rPG!+y>dYK!I5o6o9b7wy|??a{aa) z^6!xlfDQU1k`I8uAW$|K+cpk%&TUX`B`}g73fnE9unl+AUQp6ELXrLVq}vY;qR=~z zH8?0K`$cx@aY!XwEPVLvNL1N*=V6qzeu|~jp=UWH)m>f3$D9|-FFx<~A73(%xtjZD z#kIzsiRHVMFM2{bwfffD#i3T{)3Ncx>6LBp}m1tqk@2fxT?{E|q*b1EfkNc$ayoQe zmE=73Q-uGy1Tue%z#IU$p<8+K16W`=Js*+RttVo8T&1M+5zto|xANeMOQ++hzIE$| z57%;T=9uu`R%*{$y=QBVa%$+AtE#_{s+5{`dST=uAMC=Hjh?F`<)t=?{Ta3mQ+m`u z@7X(!{+B&2MCILttpqx^Vy3E7#acf{4T>)okJ}Vvp!J)3rcvwoJyeDA!4EVfeeqN2 zl@ev-%Rha;`13tOCZ6LJP53f+EyKm#s7pj
  • F%n8k9it%ULC>PP!Y|GKDK z=-BzU=3Ug)BAIIQd$JKV4uQrd9rsVD6s#b~1}=xW{d zK8fSL-cOctBV8g*$P=ktW;j>4llYRtvR){q_r%TKGddAb1CV}7rIdE5m9;E^o0P1K zI=MMN!~}BhE3OUQUhJnX!nh>*Ha^C>Sk%tv23!&&=_u8uhVHPvZaw)rHVzuo1K1gp zuO?@7_AADswW8i8De8V@%yu-+X$Yuclm+mN9y!9Fd7g;t#i!Psk}9`IlKdDx32gJbV&__sH0(SiCtd zwnsFy)Ca9khBE;}x4t_s-*ze*^}q159Gyi(q^fF~W_G!nc{VhOEu(81_#yrG^ih6B^qI?N0@D71i>0gfwqK>t@ z9X|Z2PSx19R+6hnylIFER(#RPw9yX;`T8DcVK+VL-HD|}y3&Yw!2ISPdGZkP(J4Q4 z7XLcs-AZm=;%fowXog*A$H>~AO&n4Xjn}*q_3~XOr8dE9pLL+HFAu*$E}g!&^wFbB zT;XW=xc=3*U$k!}<(g$AtQKzCwhKxk`82Xx-!ETNeX{;qFZ8~hkw9Qb<<)?wVs}hs zCt%L=p20g4B4FWFRgnPTfH$*YNwhpMxCm+d%;1LF-U*vJ99*m|AHRmzpb|kzQuT$3D2D*3``RbLvINaLl_8`#$$MLtp1WaR-q;@TBqsPj^kwf07_wgHQ za&mg*a}tl)jre*!2lCbUXI$u%k7K;Uh5fSOCWg^b<{K@ogIA;_va`G|jV0go?BKCm z_U=b{pi37cC`zCXzT4;zxt>vT{I|`hcq2r_FgaLC zT@2q4hX%%!Q$zX4YOZT?E7zhw-OCGPYwixWixYc4qw;>oaJhZRFrE*R3K>Uv_0djB zp0;las56VeSK&VG6ErrXXj zn$FWmP4D8GV7OFmeCZ~I-96nnrzOUAIKZ|`Xs^IquN!7r3QQo`TKjPHu|;K?&HJa6 z)YZx-fv1YP2ATpi&4a$m5er-%%qQz2<8;&evj|K;LrUddSPAL@Y^LBdec#)*kW2_7 zwr7qb(*LAs3n8pS_cQ5hMmHmHvZ%U7Uf-Rc=CdGSL_XQFXRbB?b(SW}=dYdXnvfCG z(amSuu7m(P&Yg0f^F@ey_V_;BdOl+#_F8XT+|&aFe#sCODdo(*`{D2|8iY?~x}js8 zlAQ5pU3iycq$c^~NNv%3-_UG>>`0@IsIP?&Uxs&`GD@kkalP^gkwd+Nz3uLWc4$;LVv`)8}c%X|^&Mx?>mmcc&&! zRY%bjid(vLQolSZy)-A%Ow8171yVJW9q+IA?+9H{SeTZR7AcFqw3c;9O1_2odSF!7 zbwLpQHTJl$=d~vfyK8kN14E+k7W*91@T|O>hf}J_HD4i5*YFw=q;lWrr@=o`6WYq& zF#)gKycG#cf15zDQ)75>R02V#ua1tTR=N#W+PlIU+5E*W`N!MZ__INc2K2)DNaxDe z8)j=#g}ghZ)x`A8>ZYb^D5g?sv{~MYHv&${V%a3Rw@5$@@OHca*1o$@_`{2wvv)Tt z_>b2n>cT|=iI9BvQ$yXW&2d2DTdi0mKJ5!xZ8SO6Nj>C?J6=|#LMf-(t^D}p9TDrk zkh$?sm~73!%0{vJIqVbctg z`Uyiw-aEQ5Zzj;<;Z|j1d_%_I<-CqgB4O=tj@!UAfsP4i&YCSZe-k!;bZ3>iTG+I^ z^0t>GtR9h~xxQylr+Z5`siJ{r9(&`%AH|hTu0GqE(&xuF&r$&WmXx9itf`nrFQAW42E+~B6}CUqpE=ALueeuPQFopXSW zgz@__hF!uny;H`zX%p$=P4vtuje}PVYwBG&9Jm|iMZ?@F*Bh@kRMQ?TKzzj7lBco8 zt{WrSU))O5Ltf5~Y94(|5WuP}1u%grclSKfeQLyjXB;nKb4aKbbKX-j8;Fz13}FHU zDXkq}ReL(oR8g;<>CpFN4aU%`YZZ@$i1v+HA8|4xq9LKfy#rqYRN?~Od51#crLqj^ zR;X+Cj<0Wwcz33zCFQ*Zn--nR<9WNYJ?e&2dkPbnbtg|1?HJdlDX^Jjj^T7@Se+f` zqXri=q<24mk6Eu^m@}>{_nUSsh}nq@VS+gsMFJL=S2X8h?0TkcbA_hpfmN%Ib}ty? z>1>V{+d z5eJ)97V(!zqaNQoD1>n0748T}EnO^uXO2!q^Py+7QY30tuY}EiMCbMCpa#)wxLCHJ z@z;YR%YU`M^E@}c`h2NgX?1`J9Q7u>!Hk^kghq)~Z?+E|YJXl@`}zum9Y1Q8HGntA zJew4X?9Po>#@>jtQHONsE5nU8EM*CBpA8xf8nya9V(pJ>G+1T1xG3+^Z7J2lf)*KT zx`2PU7I|@>W4=SA6BKG266U?RdbHmDK^KD5#le1f)<+>t;e!bBI2j8$y1_`!ll}Xk;+;Gxh(J=PWM9mC>PF*Ev3iY*ZffzIJuTokoq`xe zY7UOyj$J!5;<_>aDL4P#Oie^wI@%q%h*j2P{M8p|vF#!)24RxSwRHNB)ze3~fDRZW zANF$oW0>FW*#ok&GDq9Gl%X@|eB*ml+iji1OQ z2raXeT*lel33-=B? zC#8&vQs1(3o2RAiq9X>4Lkm8BOo*;kb@k6$qLecM&)oddSF@x2+}e-mm=sKn-|D-t zlDnln8Ab}@kCZoiEkB;23y?g#qTYP2piD%2$;=2)hr2%TbBxm%c_SZl@6sQ)7h+sz zVtdw^!2G=tb4vFk^DZ<$G@;46=E$QVY9b_E>b_I^lo1o4Fvc>Z>y7nsffy)Fe2UzK zH}~(xl%KSQh3n%4H@(kik}W+;Y+i`~g3`JJX+_Fa?)7LPjIdb6+TLrX#7E7EPEq67TmDvV(STkHl53Sf zOQGV;uFDSJV&BL=bAbswtkt=h+2qk#OUS$HcP+y?1*1>CCBhf5-j{D$;T@}q>mjr( ziI3!hs?^)F*4*#qixlujy4w5HR0?A7lF*tIT3j0{plSZh?Y{eZ4J*xp+RF>4&5t+y zz6+0h6C>6ZyJynIT4qo`L(d3ETeEqZaJT$D)sW{BMnx`jskm)q;x55!)fbwK7tGd@ zI-56MZqSlML!I52bE#tjC(wEFi>DEIY2I@4HTdpBIm>@$2MMBi3||jd6g2E?4}^#V zL1bqHE$ROH&6=DfOehnmPQF#`YDeg7G>Akkxzu%>p4z2S8j&5TNvv|y{z`L?(vyUb zsViY;8>K5wOnLdBsy0^wXe;Xhxp|&>r$p*p(bZf6`jy?_u415f*V9Mp4?1}RatJ|9V`cR5oO7O(5N$W1q!(j3~TxccWJNs>#%tnNvhfHMDLfv$U>awvA zfV`y2seov&0;kFLd#CK5+C$2AzYSM`pri z$IqAAWR7d`O1LUQKCj{TeDo`n4f2&1X^?3;LbD!$ zksYra^E)*qje0Z_Ha-&5GAykCv24ybv(DAcC*Wq3gLdNK`48`)gE)ks!1Tz-l7Pjv zzP+x&u03skc4y1Hmw2+AcE2SaBACm|`C(+>RC^BS^B!pp9}Y@qV_~RMEB-|nN7tH7 zDz3Y>gtU|9!?ASfqQrq&vTAZ>4pNav$HM zCk&&kW0;}7(HqgD7C0^08%m5AOW zc)Gaj@|pN5YVOz4$)Bn6nt4ZWUm-94k#2jm-PKi*H*DZ5Z|nKP1Q(S*r;d%P$5N@d zyP2a-?sbc-KSkn0J|SpRXH@uf5;uDHo-6>;N3*y-Z;Z6z&Du@x=0?>cqCO(j0xsWv zS#rheqM`JKuiz%XrPOg3stZ-qeKJ$KG>aN?Vp!bktAycEol~L0;3~a;21D(bnuC!` zldA{cyZ8h==Ipz00UPX}5nATo=596;sDmiAlHq=9FrPE@)J`)p`Bq3&dz}xPq!A$Y zn6OOOU_*#@&G^V7n*Q4P3U#9x`aoPF~Y@XxfYwvBnW-8=4AXx}aW1RpSn!@LQj z&;L)66QR!kMwSOMkR<9KqLQU_E3{hbMc-wQr+`O>_KzxwX`;pG&V-}eDjG%CC z@yYE;b6#DdDzWj;=}!};c0asn@PW<6I<5a<)uGUwIPf36K(0Mfrt@!NhNNYx7=9H# zcT4${M^_mElorb(*EM-wGPLw*Ux_Ft)+gc4uzB$2@G%Q!2NjbE^&l=!BzCuN-+lll z-!=qmnArAiV zLkPO0;g?;uSzp4|#f=D)x}PAc$fI9~yYq}yY3V%tFvHclV-_h1U3%{aG=Mp}JCz=G zzLs{e!6uv2qrEZg6S?H_-G=G*0-g~g_ri%d+h;X_rq}i1FSW+7MXN1)mgMesuSDvC z8|m~7+Q`j~yYaIk^MlVI;?&Sa0(EGaqIrJ#7859{I%T4afoYm+ob+DH30s&qfhr;m zQPwLAzr#VZGZi(LPRgJ8z;{|bvho4Oc7@L@(!uedSv)akMs`5w-f%xSwV*=}aEg!Q zY+S30?Z~nTsJU;5Uc{~D=3lPM9jHI;FV9XhQsCH_&n=`#bnULQi-4&^%2up$oJcM; z`3L8tZ|x6xS~6$^xpWNQk>O}Naw9qO)D3Z|>}F9;EpVx!HeJ+XzkZ`bjGq?iXr8c% z=zBpv@lc%`Dd3=K4(}hUl&|Hv;xnjSt`$VyBc*cQ8=pL+>6=ICbL$T@b{8`8$m~$f z2tH+<=6N?<3_oQ`p3WIEo4me&g6pctzwioZHk7GWmF#H~JO5 zcHzNRb6%Ys4-b=>@0e8NvA+=$B{_LO{UF{H#zI>FKvooT3V@188UB2O=nNacCdy_2 z@E;)rIeU28lf3C2>94z~bet=m!37BKPqTU;?)aVg-F54>4u z2LNpM@dE)GbAOVLKMC}*AR<7|!_UIm%ir^doAo|scuqFkXN~spmuT`oR{I!NKpR=a>u8^ruJ9&r9!bj~D50kBNf^@n?*g3n}2I z=a8p|{!b5B%)eP5baZzycW`qA^Z!o>hV^ zA&@pIHY`7K*k>p`hoZocnG z8sx`&t7WCP;^5z?Sx;>E_crI*3S;TJ_!gz`zxts|TYj#sC-{JDX$purvaa8LGywqF zl3;l&Q~*!`-w;3+e1g~@%e@s#-*Srr0P{ayP!QSNl7MrI?z;xE>4H5DyuZh>kWmn; zWIeWM5iA$f9DIDsiRA@YNbT0m19$=600AHZ!N6I-5p?;1PY=Ku4Ed=(tOP%=?<%$B z|H0*FB+Kznf*%oocNW0Q--Ff5SRgLDgSWqz6Y<}4)sdtC#%&dV)xaL!UM@e_uoU5+ z{(~fM&+mSgD8Mln*B|jLKReNb;OI;|;bF2>n19J<+g2P%#OkLUTOPdrKX`Mv_<8&L z{7VqKw;#cU;PqWgxer;exSMZ9Ks{};Od3*G;P?*Bsff1&%o(EVTN{x5X@7rOro z-T#H||3deFq5Hqk{a@(*FLeJGy8jE^|Ap@VLihjM=>8T$u^R*)0N^-ygO~*fjDx3Q zoIqI60fYgOz!AU|bd$iFxaINzDFpvF6i7f1u=oygZsqU|d*s+cJit7dA9ylcNhBZ4 zfdgJdIeXTrOgTqy&jUgBJ_i)!nM^)f|ItSm5Qmn zsgItsE5YczpR>hzGfT(w?v5Hxl2{$2RuCr0)5p`9WRDE;^zb5Lg0v;Km}5XUEBJsU za*KrIt}VHBS`uk(dIYKG?dObCk&~BolvkBUsw&7SC@HI{sqIIC2SVi!$SWRDP?S|r z!zd|Yv>B zdy(veWW9)}9~|_ZiH?4Rt?veqERG;1?oZN|1SQ=nf|KJAxj30!Ke#HMw=m>r&;b$-FXCEjF z;%p>>mkUcQ`JZBeIhdOM=cGUK@bvsiO(Y!(0L%ZsQV${v{X1~Xndt5B=jePW0IUd< zB+|(dbI#ij{N6-clHh6Y;(Wl($Js>^`Lom*6K5BD$6)*sOFeypi!%|VQ1@(LOXxPPYJLKjWFot(8L6KYSD5O&`&d0}t;0$&jOFw_-t?W!qF?cT`$==J+ z8LzJ`36@8WKybpyE2x7~swvB>Daxy=D=6Ufa0UmJ4OCSRDJv@}C@cO*tMBdT&r-l2 zX`TKr(tg9t?LGcmJivYjHZ;-S!43R~gH`8|;E->R8!3J5vJF9;=nf^m8M^}3<7gj4CklZRSYqR|Uf|Lq5Zp-pn=u&xG$y0>|JaxeXbWaL zC-~86GX%m9Y*~R3W;OzMz&+>x+zQP92MoP{hE_Ns;I0@nvQkvWO-PbSo7%LYo(v8( zErAQdQ^uRRP6~?=CK4_PqD@@GQjYuf7t8s`MJaZ0<9_Rtu_k868CNE3dl1*NJpw)| zj7n5lJG*>s(rKL};21j^HY>Jz{~bGLzaEMSD6Z|3Q7YbfiH{bBFsvJRpopP9R_O>% zp|Wv&NwMgo5noQ}W)byQ8E$W%&GzkLuj%57P9&wBk}3RDkZ8w)FcEWP6gp3RDVOIC z*M3k;KIw3>kLo_EMdi>{m$*qg;(g}Hix;ISZMjSL>Mr+tN!`3s>%_$b&bkdY%$5cL zwc>_iOkld?)Z)No@)RyMZnxFzSsqU>WHU59VJ)kT?rr_f_GGK1C>OgehtctOy(K$^ zeLGr2!4sVtCtuD6q_UZ(qVW>+`3czz7oLqlI!qGLP4jBgCF~ldmuM+s3}f)X3Gv~G zHF~BP`gV>9q=z|7DU}DB2#gnW<3Bn&O|+!bAuSzQvy53mFUX#H z+E^ucIxJBGG|Mr$Pw`W_j5RCv7tLiUAHmb_cDn)HNnB z8j~;=HpB&wlCdRMl>42(nbp*piNZst!BetUs;NSXMoNgrFybSI+rB^(L9~FOFLts> zXhr9gBr=xcK6bUg`;>Lfd_g2Pu0wqqtG4kmyRF>L9jD0E@g|%^N80w>n&P9PjgS|x zi-&LDud{NU!Nf+AH8>gYZS^T|Ayr<=6nZjy4c)FNuL%&Lb zM@r)Pa9Lkga{80|*o-=2edeVtpAN!sNmw4NTHxA~--dI~aYYp4`x&eX7@2!_Ys@GO zJtJ1JIX~g4`)~|D)WU~i#81etCTjYbLzj{cB}xs{PteoW@U8KHSkuu(hI!>~n&z9& zg-W*>W)tfs6@p5I?AUuT^_?&wS#}4!mfFzjtZ8MdCOkZTAi&T$C-!>A1x0d;nl&A>g8AF2hjJ2K80t|Uc38X+yCgnRqU z9qgdXhDytis2daKJqC0~-+UFa5fX;;@;NoSzc-~EA~U-NW+#nfr2U5^uuE35$!7FS z5hT~=&5?}{$3x?jQBljb)P%V|&X;*bjdj0ecunoTl#(wOXV}TGAgdf|C*Ckr42t!^ z;vpX=xH{mF7-}WlV?Re^CyiaXdx&bm`(Oa$VCVOHRWd#xH#0b4?QDz1PPTVxNd=fd z>~e~?_N3$m&Ocg1?%?B)O`FCP?Oi<`f&=_HcK0u#s=wwTTKxFTfzhWGPC^2?*e37ojw@KGL^xClcxuCz9R$D62N36#bR6 zLmY4DUVk47i^@&l4XhEq&MEDHgzy=z9JPifbjG3+b&Rcm$2HmdY%1d2>#v#JBeNAX zO+(oqA0iE>G>97eiDz)&StxUf42ECpeJUnE9`}{u{wVwW?Z&kJL=A1aie69;d$Tms zgm-;R#Xq$c=uZSsioW1;;4$lPGJC=VGzWV6#wys1#$*~StAG1~WCKley;Vj;A^(tn zYqBC69o4k>^fGxSr$YF8tjr8}kg<-+I4titc1j~{&+{|f*ONIhpC#s&f0ufG1|qLK zUp9G1WkK3bKlV^bd{Ab06RM6}T5Mr)Yt*KC4klPA0ooNPvN{=O?<#h2mb|oM`yUB+ zOjHsztgd_QBukTn)iwf zvv(HPcV~pfk9^#>&awwSxltHDeLOtyewXoqwIPx`QVoMHf|@Nkm+Q~>hG)4$r=U%0bJ>2gYw zfKaI$v#(7ch5fGmjLl{8(wHgGSU;s#7CmuJS$(c0j{PEekDJ{X&1cYw*rceF+G#5X z-*-xlJz~^+^T66&X74;vTVB$s%GO+>1vK(h<5l7ai=MSh#`iCgXF{w{J*9*P^*-Eu zB+$0g@^y<+_gmYqq>necCM5U_Dn~TO88ND8r78J|x-@)ryDt+6U9;(X8u^ebwueVW zqHpu#pS@IfShzuo%n2qSGb3lSw_jtCYVdd(8yGM$n;n5+%qYT=!2VM7J3lf5!Y|ge z`I)3G_mU0frvjeTTf5uyVZu^~E8tS3`ZkN%lh?%9qtWW_61}G%+1uZzBbKzECT#A! z*GfEPXE0422m99v*9Uh*L!T^9=u>Tm_O5V7@!!V4UuTyUoah6b2z<&hI=y2ev?ZTkGXopq+#LYkoCQunpXrMFg_Zm(lU8wpLLz_&FW_8nG4vZ9i7Yd&9GX? zV*19=NbgE%l!RAj`&DuosC<3XYo)$+VE$Z)!11NjJw3$CC0TSpNATvZ{oHrXfJ2Lz z?W;E8MUyufj+iBvQ5Zj6>wdl5sM~xUc7Fo8iri^hOYeuMN563MPy{&Gv|rKWoQs3eFCQ8bz|)VC zl@t$+6257Ez>KomRYE?l3Y9i#;Ls)ayWRGr$?;LyFJ`b)o^0M_uW4wENqQTHW88aS zdvSfxq%moH<64M>LMweiKNb;QrdC+JZkC)Ku44-1FDpxs-07>*u7a!3E!CXIb|MZ?nUUtqg4NRysgWh$> zBUOOP3$*25a}JNmG>5HVjhB+|p}Bu$YwWaYiz;Bh^-JaBL(SBG3IM^uPwyPnb94^l5*=hKaDi&DG-p#y37t zjtX!dB|Z<2=S8nd$zn2)uXD;TGJ(#4^ma&YJ)cuDsUKX}q*rn^`N}J2=uY^h=1u|9 zCW4}KdreEZt6LK4QdRH|X^qpI&?xLH#;cUv@ZcE1@m0#9!l>|$^znIy)oM#qwceG? z(=dS``iZ_8{Sy#avsQQz!0_9*mj{^v^Ucj#v*~$n%{MG_zfWxWs8APPloyhnv_4Wt zS&j*hDj?13gR71}Nj97uX6CaJ*ghLC{9V8nzmc-^pKO)kO=7#*5L0Hp(2} zUFH=b2`@>Tz0#>B{4fbBdD7GYY46LCObSHI)R-K#o51rcFagUXIMjgJPni;?%S^;M z1cOKXL&*1XUvqWl@kWeF_DzYML3+map&9hmf=^|>9lHET!zEExj|ue&0#Wkc)3@>{X5*aHYZCpi|lp^^2Bx!ij!#GU#w(P zalA-QQO*faF+(klWJwjV5`hfASnH7Xs*?-Ujgt2PvRWo=iyb0`qdYFxp2T({5T|u@ z&KV`l?K#^m^s_V*rN_F0)lC8(#-!DX~z!E(cq-L+gDyg3ZYT_T|}A^lQa|s=B@d(nd7wz25Tf`!^UQSh(I}-py4nt8RqNOG+)ae`arbP-IfU z8roF_15SUFFQB8g^GNoID_Fi{sZ!_5=GH9x3+4W)*~eC{%Jk5#fhF=ofhGfw*QkdnFh z8}Gr9e5WMwJREu>gCie9jz2z8F>AwyjzEOb!DPLqm-GG_%BF?IE=cq;jIUEqSoX1X zY`^0auCe>rZ2Ep~U&IEyHpCXxm>TiB#50R3sFNJeO({3QQDT1;;E&Kv69u<7dn_SA z{4#i0FwfJG^|9$%bpj~{DNF}D*3{YMUB?eqtfBezBaN&1x62gS9u-hl&>oLiw%V60 zbXq@H=o9wdNNy3W1D7PAkr`U>woPgj3ZFC>;SKDvwK5{Y?{J0Vv4)R6_8j+j&XlRv z?o+dDDNaTN=`aDo&9YYOM|@{FtlVG5z8#&GsUGzqw{G+X3YQ#k#RNvY-8`_Aitl=@ zHz3>+2az!!vQYsn5AFuAXKIzwz2EOZQx!?N3wipFtx*o~LPt zsFj^sDP{Jn!ym0_Q&kSipHC=Ixs{3>kC>IqFbdtr&P7+&Hn8ajPJu(%jSAtO1}HHX zK6yCjW=~)u!k_&@^~s%YZKawHs>poUxNAFxmjtA7*yYedU`j~Q#Uhul%-*>+O^L(Q zunA<)d8%Y~p*44!)NQ>npR;+*aGz3d3Tjj{fx5qH7sb$(D^|;W*-`~a4yFZlniX*d zOfIUR;|xLLX=|VQpS^iIBXblHhgzy4A9s&09eYP3fz9H4f}+yLk4RxRizJ7v-K!3s zJVz2NjM?ZV+%)2EK{VnbKqD);r3xQQv@e=H76YdV@_RRis$LpSKqZ_bt8SCrS(=c4 z?hPDT2M-zob*@*fZ{O*DKKPx`oP$<7uF`9!e3MO-gXXC8P;BdQH4eejY@c93D&v*k zJ0zT&Ph)&_&+~J@C-tv}p!J;2n37;SbQ;oKN8&~n)r^2lj{_DrMKQiK_EeC+Yu0+irl<2qn4t(<3fLW!S7glE1Rl za=UGHxHNYo7o(h#9P?_-PzEPEU&e6Ag((|66&%@^Z#}9AjZi3i&B$j&#zt3*yVT?U zxZ;X11l%rAQkDhPsmQ7$$_o;pnV8&Y_>xvomN0fTIj!VHqB5VqRz{%LF_e33hd4xg zQe>GDIFKV7$DWaMw3rF(n+dVBQ;$RYq9V|L)a15IMS|AXRGy9uGw%3wz;8iY=1Y4F zZ_vhE|wVEOUh-|!}%+k~32Nv1ub%Fga zZ0ynOna_CnVBpw+QaYONb!4MTdA3!Z+t)C9crwLG<4j&UaG0Kt<@LuPz5@^`dV_JK0?G?K>5?fj(~bJBs;^<`%b=K0WuQ# zk^(O3H8(vDozuRUs;s}Dbu0i?K{APj5%$6h-^41hnKKLndVn+NM#Ximqnn_@rG>XodbD)D3s>8m}%0E>}S&z*LxW2 zcvI#GZ*2FUviNYv2j1i_8!z079Tstj2&(p%wX^T6?ls=0uC#C@s&_9@W_yz^Ioc&r z1v%=QDmC4GX#FDGyz~9-tmKJ@?zh zcZ_3A?G4-+yPrA0;sZk{Z6|js%U`EPG-UC%9sGi_$oJDUT&p?!Q1pUHGjxjKM~RJ| zsuqySUCNJH0lTP8b%u~evWQ$pNQ%Ubn`*h0af7$4n*s8`sdza@h{vG|mLVPG9#7xn-nO*4%uu z4wkg~7@8YK8+t_(In@u4lnS+)7^K?xG*@>7hk|6V7C673D5Amh;Nfw|SsA*OX0Nsy zn-C;D+o3^>*%|~8e89hapgTX_fSscr_0XSzkZ*uvSdeC4ikIXT0n=P-2r&);c zB`tB#H1b4g8Rv<;=q6%a3=^=uS9auK_A>~u=+d-A(eSf9)O$%N2H`#kPKCy@We0kU zFioH#>`4>^nS_*Z6_>~Wr`9@n+z}DH0B}9hzNqNVtM_NS5|z>E*p`8N{vvY%(zc35 zT1)Fi>92eaL;x8u&_WV)yM8m>hyf3npD5ANg@$7QgG;t9bSr$gJ0u&>$1Uu$Pq6n` zsyoxgF-c{Zug)>MX4mp`?c~5Op_&9BNGhT7sBT0zeZhElkr`K zwBo_i+s+TgirrAQJN=)IBV)(3C`-9A^2!Z-T{35Qr$TPzTkr1hN8&}@H}2kzy|G8Y zyqf=IIg}VkmU8d!hgt)qvMc@<6a0Q!K`K~u)7a!^%DQiyScuT8)+W-}k zj%fOFO-5j9+KG>(6)-qITlgScC=Q>-KTdnUeviswpQ!5_RX9$k%wErwcE9;DE;fk` z7Gk^U?Kbr6)d7@G3O|SWr!r>u+qqpP4GB7qYm~`SL%FdjF(9=!y0`7+?td6`}!!{e)mYWUYI{4QMMF#AmU=a#~n+HB^ zqF=UcaAjd9vQyJin%9n*-8ND;iQ+Hwl+uhv@(M?dGwQAe?wlzQ(1IvzT2gD9IgXhn zI~cBz=chC4CV@zUl-mMil>v%>OT5^~DOR_6VB$`+ zK@RSiULp z3dL#jWoxw4Pwmjju#{R<(njiDxxAOc=$TQ}=IIDeeH?tSWIxaK6xa-4LU#8%SYDATz#W2*Q`wjz@0PKIZ8n0|ZoyrE`HZSl zV4u0*`n`kqIVFuGm8a+y(P@fkLv~aI&1@h7S%u@DAggVRxXI)~FZLS6Kim zLQPtdcwEwr2e{~lRi95Fnm8c!*c&8pg8OUVzEA#I{O#m6Z(0$li^?^T>Z+cMy_#^` zIa?rTQnS70WWn=uPJ)U~*IB;Na00K3$wk|wYC4mTM+a2K1yHs$v=su@OS7(Z#VQurez1Br(| zy(oLZ>(*3+V^BD+13CCK+0QO{`E<=5>QB|FWToVdx1R1I=yG5I%pI z^lrvqJ4jSEBxL;U+%mcC1R>ASOg}uBVGM4V0q4wa6(&vnV5I3}4^+iJjhuLgPrWJkrYj>o(+j4?WhZe0-yy zo7-n6Y{XAoAQuHd?^o&S4w6JePE1M|?wB^{3Ta$mpZ9yRu63PFUp)rt4w_sm>pmHr zv8xAGm6E*#!ge2vFEfe;Z2^;5U}J~H7qfS9wJ&c7|7IC_lHwEX25XrF*8}dn*-!67 zZ$buI_^>6*p`}&{1~5Q>ZiG<;Zg}jEyq-1>7O&S#MoWAup4vQJQJwr8V8bIMh7MDT z)O=2&jkIH|K0ht3iFU+6ao84{YFoM58=ATi!#b&zt()_HJOV3?CvcWkoJeX48`v^-vd@qMIrl+TIx Q%i7-Ox?=fb6z04C0pd5rzW@LL From b614b34eab22d8d06e87c34426229ba1fb352324 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 18 Aug 2012 23:20:49 -0400 Subject: [PATCH 15/21] ENH: Add test of reconstruction by erosion. --- skimage/morphology/tests/test_reconstruction.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/skimage/morphology/tests/test_reconstruction.py b/skimage/morphology/tests/test_reconstruction.py index a3461f31..74f61b9f 100644 --- a/skimage/morphology/tests/test_reconstruction.py +++ b/skimage/morphology/tests/test_reconstruction.py @@ -8,6 +8,7 @@ All rights reserved. Original author: Lee Kamentsky """ import numpy as np +from numpy.testing import assert_array_almost_equal as assert_close from skimage.morphology.greyreconstruct import reconstruction @@ -68,6 +69,14 @@ def test_zero_image_one_mask(): assert np.all(result == 0) +def test_fill_hole(): + """Test reconstruction by erosion, which should fill holes in mask.""" + seed = np.array([0, 8, 8, 8, 8, 8, 8, 8, 8, 0]) + mask = np.array([0, 3, 6, 2, 1, 1, 1, 4, 2, 0]) + result = reconstruction(seed, mask, method='erosion') + assert_close(result, np.array([0, 3, 6, 4, 4, 4, 4, 4, 2, 0])) + + if __name__ == '__main__': from numpy import testing testing.run_module_suite() From b5d91069664b4f323df9b2a201966a4e8ac8160d Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 19 Aug 2012 11:56:16 -0400 Subject: [PATCH 16/21] ENH: Use Cython data types instead of Numpy dtypes. Conversion to Memoryviews didn't improve performance, unfortunately. Minor slow-down of 5--10%. --- skimage/morphology/_greyreconstruct.pyx | 39 ++++++------------------- 1 file changed, 9 insertions(+), 30 deletions(-) diff --git a/skimage/morphology/_greyreconstruct.pyx b/skimage/morphology/_greyreconstruct.pyx index dff670bb..c7dcab38 100644 --- a/skimage/morphology/_greyreconstruct.pyx +++ b/skimage/morphology/_greyreconstruct.pyx @@ -10,23 +10,13 @@ Original author: Lee Kamentsky """ from __future__ import division -import numpy as np -cimport numpy as np cimport cython @cython.boundscheck(False) -def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, - negative_indices=False, mode='c'] rank_array, - np.ndarray[dtype=np.int32_t, ndim=1, - negative_indices=False, mode='c'] aprev, - np.ndarray[dtype=np.int32_t, ndim=1, - negative_indices=False, mode='c'] anext, - np.ndarray[dtype=np.int32_t, ndim=1, - negative_indices=False, mode='c'] astrides, - np.int32_t current_idx, - int image_stride): +def reconstruction_loop(unsigned int[:] ranks, int[:] prev, int[:] next, + int[:] strides, int current_idx, int image_stride): """The inner loop for reconstruction. This algorithm uses the rank-order of pixels. If low intensity pixels have @@ -41,32 +31,21 @@ def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, Parameters ---------- - rank_array : array + ranks : array The rank order of the flattened seed and mask images. - aprev, anext: arrays + prev, next: arrays Indices of previous and next pixels in rank sorted order. - astrides : array + strides : array Strides to neighbors of the current pixel. current_idx : int Index of lowest-ranked pixel used as starting point in reconstruction loop. image_stride : int - Stride between seed image and mask image in `rank_array`. + Stride between seed image and mask image in `ranks`. """ - cdef: - np.int32_t neighbor_idx - np.uint32_t neighbor_rank - np.uint32_t current_rank - np.uint32_t mask_rank - np.int32_t current_link - int i - np.int32_t nprev - np.int32_t nnext - int nstrides = astrides.shape[0] - np.uint32_t *ranks = (rank_array.data) - np.int32_t *prev = (aprev.data) - np.int32_t *next = (anext.data) - np.int32_t *strides = (astrides.data) + cdef unsigned int neighbor_rank, current_rank, mask_rank + cdef int i, current_link, neighbor_idx, nprev, nnext + cdef int nstrides = strides.shape[0] while current_idx != -1: if current_idx < image_stride: From 682f064f86b9e3052fe48dc6bdd81f5479dc21c1 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 19 Aug 2012 16:23:14 -0400 Subject: [PATCH 17/21] DOC: Fix docstring to match rank order. --- skimage/morphology/_greyreconstruct.pyx | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/skimage/morphology/_greyreconstruct.pyx b/skimage/morphology/_greyreconstruct.pyx index c7dcab38..5d2cb0c5 100644 --- a/skimage/morphology/_greyreconstruct.pyx +++ b/skimage/morphology/_greyreconstruct.pyx @@ -8,9 +8,6 @@ All rights reserved. Original author: Lee Kamentsky """ - -from __future__ import division - cimport cython @@ -38,8 +35,7 @@ def reconstruction_loop(unsigned int[:] ranks, int[:] prev, int[:] next, strides : array Strides to neighbors of the current pixel. current_idx : int - Index of lowest-ranked pixel used as starting point in reconstruction - loop. + Index of highest-ranked pixel used as starting point in loop. image_stride : int Stride between seed image and mask image in `ranks`. """ From 3fe03259d09ee8bcf4ffd31de066bde285017df1 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 19 Aug 2012 16:24:12 -0400 Subject: [PATCH 18/21] DOC: Combine examples for finding spots and filling holes. --- doc/examples/plot_find_spots.py | 69 ------------------- ..._fill_holes.py => plot_holes_and_peaks.py} | 24 +++++-- 2 files changed, 20 insertions(+), 73 deletions(-) delete mode 100644 doc/examples/plot_find_spots.py rename doc/examples/{plot_fill_holes.py => plot_holes_and_peaks.py} (73%) diff --git a/doc/examples/plot_find_spots.py b/doc/examples/plot_find_spots.py deleted file mode 100644 index ff6e2196..00000000 --- a/doc/examples/plot_find_spots.py +++ /dev/null @@ -1,69 +0,0 @@ -""" -========== -Find spots -========== - -In this example, we find bright spots in an image using morphological -reconstruction by dilation. Dilation expands the maximal values of the seed -image until it encounters a mask image. Thus, the seed image and mask image -represent the minimum and maximum possible values of the reconstructed image. - -We start with an image containing both peaks and holes: -""" -import matplotlib.pyplot as plt - -from skimage import data -from skimage.exposure import rescale_intensity - -image = data.moon() -# Rescale image intensity so that we can see dim features. -image = rescale_intensity(image, in_range=(50, 200)) - -# convenience function for plotting images -def imshow(image, **kwargs): - plt.figure(figsize=(5, 4)) - plt.imshow(image, **kwargs) - plt.axis('off') - -imshow(image) -plt.title('original image') - -""" -.. image:: PLOT2RST.current_figure - -Now we need to create the seed image, where the maxima represent the starting -points for dilation. To find spots, we initialize the seed image to the minimum -value of the original image. Along the borders, however, we use the original -values of the image. These border pixels will be the starting points for the -dilation process. We then limit the dilation by setting the mask to the values -of the original image. -""" - -import numpy as np -from skimage.morphology import reconstruction - -seed = np.copy(image) -seed[1:-1, 1:-1] = image.min() -mask = image - -rec = reconstruction(seed, mask, method='dilation') - -imshow(rec, vmin=image.min(), vmax=image.max()) -plt.title('') - -""" -.. image:: PLOT2RST.current_figure - -As shown above, dilating inward from the edges removes peaks, since (by -definition) peaks are surrounded by pixels of darker value. Finally, we can -isolate the bright spots by subtracting the reconstructed image from the -original image. -""" - -imshow(image - rec) -plt.title('"holes"') -plt.show() - -""" -.. image:: PLOT2RST.current_figure -""" diff --git a/doc/examples/plot_fill_holes.py b/doc/examples/plot_holes_and_peaks.py similarity index 73% rename from doc/examples/plot_fill_holes.py rename to doc/examples/plot_holes_and_peaks.py index 7257d2bb..a9c3a7fe 100644 --- a/doc/examples/plot_fill_holes.py +++ b/doc/examples/plot_holes_and_peaks.py @@ -1,7 +1,7 @@ """ -========== -Fill holes -========== +=============================== +Filling holes and finding peaks +=============================== In this example, we fill holes (i.e. isolated, dark spots) in an image using morphological reconstruction by erosion. Erosion expands the minimal values of @@ -62,7 +62,23 @@ original image. """ imshow(image - filled) -plt.title('dark holes') +plt.title('holes') + +""" +.. image:: PLOT2RST.current_figure + +Alternatively, we can find bright spots in an image using morphological +reconstruction by dilation. Dilation is the inverse of erosion and expands the +*maximal* values of the seed image until it encounters a mask image. Since this +is an inverse operation, we initialize the seed image to the minimum image +intensity instead of the maximum. The remainder of the process is the same. +""" + +seed = np.copy(image) +seed[1:-1, 1:-1] = image.min() +rec = reconstruction(seed, mask, method='dilation') +imshow(image - rec) +plt.title('peaks') plt.show() """ From b1007f019675621a372066433b2a374e62c3b1b1 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 19 Aug 2012 17:46:22 -0400 Subject: [PATCH 19/21] ENH: Add regional maxima example --- doc/examples/plot_regional_maxima.py | 113 +++++++++++++++++++++++++++ 1 file changed, 113 insertions(+) create mode 100644 doc/examples/plot_regional_maxima.py diff --git a/doc/examples/plot_regional_maxima.py b/doc/examples/plot_regional_maxima.py new file mode 100644 index 00000000..9d4de9b1 --- /dev/null +++ b/doc/examples/plot_regional_maxima.py @@ -0,0 +1,113 @@ +""" +========================= +Filtering regional maxima +========================= + +Here, we use morphological reconstruction to create a background image, which +we can subtract from the original image to isolate bright features (regional +maxima). + +First we try reconstruction by dilation starting at the edges of the image. We +initialize a seed image to the minimum intensity of the image, and set its +border to be the pixel values in the original image. These maximal pixels will +get dilated in order to reconstruct the background image. +""" +import numpy as np + +from skimage import data +from skimage import img_as_float +from skimage.morphology import reconstruction +from scipy.ndimage import gaussian_filter +import matplotlib.pyplot as plt + +# Convert to float: Important for subtraction later which won't work with uint8 +image = img_as_float(data.coins()) +image = gaussian_filter(image, 1) + +seed = np.copy(image) +seed[1:-1, 1:-1] = image.min() +mask = image + +dilated = reconstruction(seed, mask, method='dilation') + +""" +Subtracting the dilated image leaves an image with just the coins and a flat, +black background, as shown below. +""" + +fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 2.5)) + +ax1.imshow(image) +ax1.set_title('original image') +ax1.axis('off') + +ax2.imshow(dilated, vmin=image.min(), vmax=image.max()) +ax2.set_title('dilated') +ax2.axis('off') + +ax3.imshow(image - dilated) +ax3.set_title('image - dilated') +ax3.axis('off') + +plt.tight_layout() + +""" + +.. image:: PLOT2RST.current_figure + +Although the features (i.e. the coins) are clearly isolated, the coins +surrounded by a bright background in the original image are dimmer in the +subtracted image. We can attempt to correct this using a different seed image. + +Instead of creating a seed image with maxima along the image border, we can use +the features of the image itself to seed the reconstruction process. Here, the +seed image is the original image minus a fixed value, ``h``. +""" + +h = 0.4 +seed = image - h +dilated = reconstruction(seed, mask, method='dilation') +hdome = image - dilated + +""" +To get a feel for the reconstruction process, we plot the intensity of the +mask, seed, and dilated images along a slice of the image (indicated by red +line). +""" + +fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 2.5)) + +yslice = 197 + +ax1.plot(mask[yslice], '0.5', label='mask') +ax1.plot(seed[yslice], 'k', label='seed') +ax1.plot(dilated[yslice], 'r', label='dilated') +ax1.set_ylim(-0.2, 2) +ax1.set_title('image slice') +ax1.set_xticks([]) +ax1.legend() + +ax2.imshow(dilated, vmin=image.min(), vmax=image.max()) +ax2.axhline(yslice, color='r', alpha=0.4) +ax2.set_title('dilated') +ax2.axis('off') + +ax3.imshow(hdome) +ax3.axhline(yslice, color='r', alpha=0.4) +ax3.set_title('image - dilated') +ax3.axis('off') + +plt.tight_layout() +plt.show() + +""" +.. image:: PLOT2RST.current_figure + +As you can see in the image slice, each coin is given a different baseline +intensity in the reconstructed image; this is because we used the local +intensity (shifted by ``h``) as a seed value. As a result, the coins in the +subtracted image have similar pixel intensities. The final result is known as +the h-dome of an image since this tends to isolate regional maxima of height +``h``. This operation is particularly useful when your images are unevenly +illuminated. +""" From 945349f963b51e614f737c78825c3eea3d9e6c41 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 25 Aug 2012 08:48:03 -0400 Subject: [PATCH 20/21] DOC: Improve description of reconstruction --- skimage/morphology/greyreconstruct.py | 61 +++++++++++++++++++++------ 1 file changed, 47 insertions(+), 14 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index 04ea54a3..9362e3fb 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -14,26 +14,41 @@ import numpy as np from skimage.filter.rank_order import rank_order -def reconstruction(seed, mask, selem=None, offset=None, method='dilation'): +def reconstruction(seed, mask, method='dilation', selem=None, offset=None): """Perform a morphological reconstruction of an image. - Reconstruction requires a "seed" image and a "mask" image of equal shape. - These images set the minimum and maximum possible values of the - reconstructed image. + Morphological reconstruction by dilation is similar to basic morphological + dilation: high-intensity values will replace nearby low-intensity values. + The basic dilation operator, however, uses a structuring element to + determine how far a value in the input image can spread. In contrast, + reconstruction uses two images: a "seed" image, which specifies the values + that spread, and a "mask" image, which gives the maximum allowed value at + each pixel. The mask image, like the structuring element, limits the spread + of high-intensity values. Reconstruction by erosion is simply the inverse: + low-intensity values spread from the seed image and are limited by the mask + image, which represents the minimum allowed value. + + Alternatively, you can think of reconstruction as a way to isolate the + connected regions of an image. For dilation, reconstruction connects + regions marked by local maxima in the seed image: neighboring pixels + less-than-or-equal-to those seeds are connected to the seeded region. + Local maxima with values larger than the seed image will get truncated to + the seed value. Parameters ---------- seed : ndarray - The seed image; a.k.a. marker image. + The seed image (a.k.a. marker image), which specifies the values that + are dilated or eroded. mask : ndarray - The maximum allowed value at each point. + The maximum (dilation) / minimum (erosion) allowed value at each pixel. + method : {'dilation'|'erosion'} + Perform reconstruction by dilation or erosion. In dilation (or + erosion), the seed image is dilated (or eroded) until limited by the + mask image. For dilation, each seed value must be less than or equal + to the corresponding mask value; for erosion, the reverse is true. selem : ndarray The neighborhood expressed as a 2-D array of 1's and 0's. - method : {'dilation'|'erosion'} - Perform reconstruction by dilation or erosion. In dilation (erosion), - the seed image is dilated (eroded) until limited by the mask image. - For dilation, each seed value must be less than or equal to the - corresponding mask value; for erosion, the reverse is true. Returns ------- @@ -42,11 +57,29 @@ def reconstruction(seed, mask, selem=None, offset=None, method='dilation'): Examples -------- - Here, we try to extract the bright features of an image by subtracting a - background image created by reconstruction. - >>> import numpy as np >>> from skimage.morphology import reconstruction + + First, we create a sinusoidal mask image w/ peaks at middle and ends. + >>> x = np.linspace(0, 4 * np.pi) + >>> y_mask = np.cos(x) + + Then, we create a seed image initialized to the minimum mask value (for + reconstruction by dilation, min-intensity values don't spread) and add + "seeds" to the left and right peak, but at a fraction of peak value (1). + >>> y_seed = y_mask.min() * np.ones_like(x) + >>> y_seed[0] = 0.5 + >>> y_seed[-1] = 0 + >>> y_rec = reconstruction(y_seed, y_mask) + + The reconstructed image (or curve, in this case) is exactly the same as the + mask image, except that the peaks are truncated to 0.5 and 0. The middle + peak disappears completely: Since there were no seed values in this peak + region, its reconstructed value is truncated to the surrounding value (-1). + + As a more practical example, we try to extract the bright features of an + image by subtracting a background image created by reconstruction. + >>> y, x = np.mgrid[:20:0.5, :20:0.5] >>> bumps = np.sin(x) + np.sin(y) From 666d6f4f095adedb7f64225a609b2ab2bff11b66 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 25 Aug 2012 08:58:54 -0400 Subject: [PATCH 21/21] DOC: Bump up Cython version up to 0.16 This PR uses typed memoryviews, which were introduced in 0.16. --- DEPENDS.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DEPENDS.txt b/DEPENDS.txt index 2792870b..f06ceca3 100644 --- a/DEPENDS.txt +++ b/DEPENDS.txt @@ -2,7 +2,7 @@ Build Requirements ------------------ * `Python >= 2.5 `__ * `Numpy >= 1.6 `__ -* `Cython >= 0.15 `__ +* `Cython >= 0.16 `__ `Matplotlib >= 1.0 `__ is needed to generate the