From 40e484323b1cac99e131003a05aa80f74c6263e9 Mon Sep 17 00:00:00 2001 From: James Kerns Date: Tue, 7 Nov 2023 23:58:59 -0600 Subject: [PATCH] finish global field locator --- docs/source/changelog.rst | 2 + .../images/global_sized_field_locator.png | Bin 0 -> 15756 bytes docs/source/topics/image_metrics.rst | 64 +++++ pylinac/core/metrics.py | 237 ++++++++++++++---- tests_basic/core/test_profile_metrics.py | 137 +++++++++- 5 files changed, 384 insertions(+), 56 deletions(-) create mode 100644 docs/source/images/global_sized_field_locator.png diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 8233f8ef..9036b93b 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -13,6 +13,8 @@ Metrics without knowing where they might be. This is relatively efficient if there are multiple BBs in the image compared with using the :class:`~pylinac.core.metrics.DiskLocator` class multiple times, even when the BB locations are known. +* The metric :class:`~pylinac.core.metrics.GlobalSizedFieldLocator`` is also available. This metric + will find a number of open fields within an image. See :ref:`global_sized_field_locator` for more. Planar Imaging ^^^^^^^^^^^^^^ diff --git a/docs/source/images/global_sized_field_locator.png b/docs/source/images/global_sized_field_locator.png new file mode 100644 index 0000000000000000000000000000000000000000..bf04db6af87aa11803555c6f65ae7d82abc86481 GIT binary patch literal 15756 zcmeHucT|(>zGeUg#lp5AMT%PxQF;>uDT0*HJ0ewSB3-I{<8DAG6lXoN-062=IOH_j`Z!d7k&VuD0qidRBS_0&xtf zc0(V5pbkYKsJN;b@}DYKKF1M(xQ0jG zxNhk8c5N!i%fjd(ncRC7cL_NVFp$mdn03i0O8*vzbA`FyLLGNO*8IvLubS6Zc9%tK zLqo3#(8{Q@qgu+=C3_ZsQ)Z{TDgQdcLPPg8ozsybwMC`Oak@jw^m!j8PpD(G#%JYH zh=iH>ck|guWb(6#V~salGNW@# zHZ;#rT|hwK-rc+IRb!P!PoL6CNJ=u5+EMPupvCw6@}6BcN`w%=^< zuV~yUa$~Gxf$isBMK@XHqz3z!ChHkaM>IueYbEof_U7#F?uM?h?L#E*%E9MZd`OLYiS^-puXPx z*s){db-r%bQKg9bV*on>8SZVeKvo~RiQuIGyBy1b{=56AjlliU2N_g%e~ZryKpZ6fC;e>TK~ zuNn8rzW2&62PY@_k@s9(jg5@r&Ip<>4m;kx``Tp!R@+%!1-@;4?Z?Zh(9)AHKXD@= zRt`=$IXO)QmEP2pkdQb?AO<&OYo{y^v;151Xmg*;7;*+{h)CWJK%n}UhYA><-rJ)4w%WLk4S>-$Ev z*)A3>-%&(|KrqXkx`OC=Cu`>I?LGeC_M7SfF{_F^BjY1ypk}smJu5M@Csia08J5i-0qgEN>%C?GCLK=koz7zsvtR7j5Q^hY3ECa2osHJmo)Zc&#ZCF9_Ah;Ot^P(E zMA}Ln31<;|;blZgEQO=3uT{T~Ii#|C`QctL@`a1kUpXQJNIOM`3CGXD9cfkng=*2IB}BFV&b^>nF3!H$_821rgA5I!2E(GFqJ!T4 z0RcmUgPF_A_l%5B0CU*H*f-^K%Q?BMX zzyl7v^ohx8aBpWNPCKWl)J(ktKgbC`Q+~W z(o#LBE<*L|VH|<$*Jan6X>L{AGZGP`c=l@I;mnMTRyEP&D=x5k17l;3)DJ^Vt6Wpe z7UteqU5VE8^eh)?OSqUC*gsn7a;G%G%%M9;f`{#RvaRUV!mYOd6tpjvGQ{J&7``fO z_?btGdoLRz^))nJg8Q+uvQA1vKAl-7vn*E=8e1v6Dg zC#wS`?R*q3FStwAo!<{Apgg>`x+PcMzZC4Rf3U@kgeY#Y)-^IRI_M3q$@ucqq%^VB zDmWwQ2Ro< zdcgw1xAs4|Hs7uO_Vn!bAE$KgZq0W~6SoJZ)*xurv-Fo&SNRtYLdcs~SfrRatPB?& zM0G}<>=|LP{M{MYdld4ig2wqgHqM>@ISmqEF8{9 z{q_J;v(VM@h;His?B^3zJryICPRlDR-EFVVc0*d!5@(g@yE{0Q+HTtg#_9V0>1B5W zyG+-Bj=YPu%$j)l{ke|YRgd(VPt!k9Opn#RefopQ1$M-1^hHc@M0RaImBGASSj*jMn(IF2k?oPOw@WilS zH6Ud&OTrH6$P>M$f^bv!93A^$pLg+?!Y21Cx`^;Ue?x7gNP6U6OC|&kviH+rK*etd%EG^Elz{hTo*x zd#KK^`P3IUA|ox=dG#(#(TrMMNgf)naH2GF${Bi4(Rfp|3)0hDio_gx*oH<%3TkRh zqoHCK;4deULF+=lJTg3d`OHgX^qtb1*MD?kd$DIb{-Id+H9*z~6*1hTS}!c` z>4!-8XxDkh!gtF1R_UE55RhWo3+De&rLz~nXPLD~cO?S{T%$&p^n!OjIrr z(B#JA>hJVf`j@hooCn9!&Iq(t5iazMq^p4!MW} z`W^SVj<{R3UU?Aa&|_q!lJ^(0ERNPloHJzVlkaL;s_)6PZH$!pO zp(jk1yvSSiex>t0DFp?^wd%M@&q%Oao0G8#vSLszZAOCL|k*5BFf6j&hM?1NE&2oiZC)V&X;uxiqSD}@Jqp6>TkyjTHNdC z?Cgg6(&bV$W#7NFeYoP0qcz1N6fdRxE-u`5vyGW6Tr z3>c%FukfKfCKV4I(A9(~_<{f@DNrtTKe!3pfQLZLFLg1}xV%!}Fo9~^WbP&KaO*$w=QLP^Jro*(}-aHvXf6OGk z>%M%heTi@Q7g_(t#&>n~&S(>|u54&fLkMQUg^=J9+!@WU(dHq3DcQ*TNCwBG=YQ0p?s{RhZ>Yre-Dyx7zgPW{8!}|z(zo!c1@txL3 z)jzmLkSD9YxM7Ri>Du8$+;I@yyckmxRZ}T|C*I#o_c-JuEzs#<8z@)xE~Wjm&)dm> z$8xD`Hkp3rIr4pc>iTp}_3;!P< z2*H*el@&8(ow55(WD~}m)M8y;?k2RAu!Lp^jVf)^NbjR1vbGgZjGk!X;8`Z~v|aSR z3KI!S?!ZdXBWt%mpCC5Yn>V0X?t=AR{j?s9#jBacLNC`JPq;2wS`ngzKKLRCXY-$RVUiC>hiEsUg;IMjo;q6VA~P+beS>`U7rJ*sAt_2J zHn|*4E&d90!I*~S9V;StOnbGJDb+!#3IMezyY`*jmGS3~B^X5k@&aAt0ygBZG z-pKYW<=0paI4XT>bBpt-;;$(u=$8U*yp+$UGwg2PEt~o5U`jvFiP{a~Ny})+#PfBO zgq}yTr0QTh-u#=BGOca)qyY?DpO?2An+V?$N=7yDvz1ULTX|KdsSV(`}S|TTlBjvXU zLt~P72~SSyPcX6RZU$aHDml!Ht9a*;r+eP197`&;fhJ8vbgXb<;9huCJHzIp8EY@b%pyy( zhEX@(*VoBYs!)FSRX7jgZ+9$gauLF#j% z(QouIM6sq=vlm5f9dAh5yM0(SPL{Jav4TP-mziUQ;`AdkF?xPSp<}prL2E6U{6)=w zg@`qokl%l#c!ViGuW48&2*Qn$tMWq{(kaHh7v=M&-k6_XfTI3QmDA{fTM2cv zpgNEDtp_#VFM8P`Y@Sf+nlMcxM>(0>AE9cae9U+K`QH2`w!S`kRTcuj2C31;AxLhG zxShlJ6vow>w*)rOc<-w_>$O~tdg#1+1c4$77{X1 zR(^68L6^M0XIt!&+dI2JJaT|qv~BTo6jxb6!Or+I`9d*YY`k6gZ+ zwT`i!fGr9S2!-`I_ZT0;eH5p%sS(B=HD7{Uh$~g+{j1IH*2a)#7Sj~AwEt5C=$l-6TstavLo@6%ScT}V$ zy419M*o;;dGrKHeF7gdvqgE7JExKnLjyed;U0*p?E`VsGqxz-!u}A$)5kuXjrYUJZ z-Qi@ndEeQ_%1Vx+z&4hC7KgVdsd(&6FJY=zPX8=YFwv*Jm%#(GzoBB# zRu((hx8p~%^s%t%kM}K6@%J@y@a(FU_cbyOs{VQ!PbePic;jq59D4QWvs2VMKI@0PcY1 zH9s$}Z>FyG>RA@I4nF$ zjOckL*{_<$(y4z-6&-N&=%Gd%#%0k>kMiLWqf=8wbDeKia^#75Fy4K4jYQ4tjr9g+ zV}q~(dxJT0SqG1tStJp#BS(bl)0nIR^`f^Y&?Mu4-PJxN2HDK5UU@5+^@jj+cbN*n zHb4jasg@gYv4a_>x2TIx)x|cvrUCqhjO5Swk!WJ1pFcYjnp-V2^n{C|nd@q}-kG>Wf)H!tQE#2$u=M+5R0uj#l?&;guWB>w0Us?5YJY;fC1&D?+ z;rtEA9p1-p8JU@_44anOJ_r2J2Rki0do`v#^W^gz*Qs=`;e}%*N4dUkEPJrt%V*)L zu@gV=D$44<(wvn;xBlgm9(18^y_%QJsFd4hUmI!9Tb{&)8fYLHs7>p9iYWZdJ=iqV z>dQ}e=xHCpAbEMz(#))YGT|<9k~Rrjgr{W`q;ACw&3b!zd7j@jOg#J~as^Ly^b79LXnfo3iS=COxi2vx+sT{=vAGoA)7JujJdU(6HO6?j znV1xN`7Zc)`EpLhs~;OL_dCz`n;KWD4{h@zyYB%~wfBu7 zB-i2H0~NscQc_Y5HiZ7k+r)6g=;-(ptm^&=R%!B{2H1x6QvsFm$((+!uPviX+qpV4 zF{S)FmzhZX$heK_RW8*@;zh3py5qAPgk8xlMqOtvsuu;eG+zO8EQB%jjC0-I za)AU4_ER1x_N1$x*?<;OcptmYcx3c{M7U~KitB(6OHWE-MfL+;$$0#kgfC`tqv0Il zV$s(VAG{`2n~4GGw1mbcQY1fJf!H7+7`6Yldw$tzJ6ow8hnOF0tKLGqbN`SeqD0x- zsi~TE-(gdIJ3H42JS!3Lk9EGdd@WJ4OR-HMKTk?)l2KVaHf2 zMcU|_`%Mbp_P#Gyuzp<;*$*HYti#;Gd&sDE;a$M#m!Gy@xdkdgHrZaMy1<<^sc@)? z8?|6Cj1}^;y}G_(M;4YGXwc;$o{B#yqsmV|F|4&tY|}rT*aDq3dk3#=Q|#^X1LF$AKN{SZl}GZR)`<;z-b-0aOcpF1!Ne&?4X5GO04jC+`s> zXQU!vl7Npz6E%x%%b!pgabc;Bux2PU*))Q3%mfsK4)%{u=H7)yXhnL@F5|OUWUUKOH()?=scz;R8mu;yu1k|=2y!w~Kp=@Y?zG8=%e^9) zS6-6i&dukO3H=p2`yqYb^|a$0(&G0q-*Ja)5D0YyI*XPOGNmH z^5scYu2r`4;;L)O8;H{kfT18e=nkPgwku~Azm zTJL%79KX4@>DN{(ImpaKX{IXe9Lj6*{P$8~EEWFT&s(IFI`+zkHMT(sQ7&gZyV0ywU>M<=rO^vAq?+1{}Pa|nQEMbJiSy1v7YA8X(mb= z`$hidEW6OhOjL$mkc#K`g?bL8P<-6dOGw{TzqawlX!9j%=0p4TP~n^;lNEYfvh`Tt zMKLLE&SrzWSjds^qD!z3ECXoF}qz*oij*Se<^W3k_Va#${e3R4AgH3p|_%o3TO? z&!zW`?Iq+rldDZXo>caxmT&>w`&EkGUbJQishDnKEEwfA8K)J|$HX0!7{1y#WWe3- z5H53h9JADC8h~4_Z`xOm(sW&mfJ4G+wGTc#lRPy5c z3a)|W&QR=3kV=rOKR!g$*=AB9(NJ(-7Ep(w#=JXz2adjFg(*h&4n0rPF5b)LTpD9y zc(J(3i8)29$UPn_3_#d% z;E~;r5|2OrPD3}_@OR#)dti7oUT#q;-$*c&eHVXPT!k^nbA8T4I?9ew)1O^jMGk-; z)?h?e2+znNe7K%*&D>I01Xm}s_pSKpkNoV>lS4Y(*`Y0Yz-aTa9cQXzQC4Pvht4I* zKAeirlHWlt;luO#e{&8?bP)3|J16rQdZR6b9w@TjhCw zIkYH)#`jcP9ii~s!WUL4Vydg4L!q4-`F9$+ZtcIOp;vu?X81?#Nbvm@8|q6-s?E_K zeX*9^2_j0^5Hg_x#49i$LZ5o;R;+f;(3dY?hg1Gn<`4veml4Q0CeVXXW;r0~K0*tX z5k%NsI+gn4F$f}{nFPEA$Z4bX#PRW4pYDhtLBhm$Cf8a4bw%s- z%ZAdx^_P@5_g-5DtVzJuU5ZoBD=#lk(+&uz1Jw;T61ZzrIDHTKZHZWdE^ns0l({Z! zo%*N8^r2HzQ+u7zyMRW>=^q%XrvC;*#Q+Nyd?8D3$L(6Tc4ceQvo{eQJyCxD8$&km zeVWNvy1l%-m~FyqcE_Exfaf2p_Q)7(461vslEK6Jsy4OPR@8pF{(=3*_fPhy2=@7< zoYX5hDLzS{oS6eE=|9udCBVLk1?_pAI!fD3n8wL@equC97_kiKoW zHsQkOpWjf7Qp+6e2Gg|u5{Z%}lSq&K=|DM&zygZ@FCdUqmqRq}44du!&eWI&QWU6E zPNaSOIP6#*nmGe#W5>PiB?|hdr1WSPGiBx;{+CL!jmFbJyL#w>k_!Lms}}$iK;wstX}9KoFg#dlsZcph#=!s^=!^{QpoKht0qgGYa5TR zxJ?(91kUk&s^*4kmpZ7AvinZnmLqLMGSr|w*bGce92z10)J{mt$Y}AN7LB~N`5VYo zWCM1Yk@hg`yQWas_nYx{)BY<(w=1|fBF`d!?!T^j+8;9;F`68wRFOCYX>X)y`1dDA zc4~Tqv?;0XYuJ(PfmG7qRFfN84#HPLT)cV$VwR!-ou9Wy>Vq5x0ECb3-xRB~!7R=b z!vH?DVee2QDYP!ZC4WU%!%#r2!&2BxNtAYH8ikVDJFWaPty0wR+ldC2mjr~1MS zlHSKBklZI<+&KTq71B5JQ$Tp|BUO;Z72qM?NouBKKzK-!%43Jt$7;GAE6#>sDG5Ug#0(8E%pB-ZTayxY0Krf zehNY=6ttuP$TJAE)pG_Y8ts*?L^0K!ogL7i#&d&Go02#4i;BEq;4HfVb*0W1%Yk!( zr7<&5)FmIU?T*?&hzv3{{ZOWK1O~#)zy=##Z}_uGsM02{zd@&>C1Y(x^K3*f+92zO~m-_ zp+fdmea%3w22&Jd$0N&s3RilH#M}+tNeh7~ErtIO?^%459SJyf9zl_WK`*EWP;H&d zKePr#|7C0NdFzSduORaX0}-(9WpxreMYKmI;+#P2y-A3)1buBEMRg21Un7GYZ6lDo zfxavkBvJrK>qdVB6(eY1OiV1i<3Q=c0ZB4w<}6wloC~75^rueI`KNY87D8Sn+R?>D z*4_%HF+c%Ux??XHXaYM`vs%_k5y~A|{-3r0XA@|j{w+;t(%Af^cn!C|Vs`5L&uBi% z`G2mG-2F7=Z~+L8gf;i;v=D5J&^2q9CWF5`%_}G<7}esH4k(JlE2UR>QdgmlTz6_h0joQ{sxN_ao}GcN@+e(i<=Jch7^-qSf~k=q#yFZ%c;!6 zw|AC@_CBl1y@wMBps7=NW9gpm7kldtEKhIwTb9sCz`lX zAA9&4rVgk|(^~v|@|4RYxfEI;x1}RiKf<6;8d+KvgHFwfu)Bu&OErCq$Z>Y*OrTp} z2rS6YSEnG&wRLqUg1)*w2;nF)?6GnOiu|7QM>E~_np;Sn^)b_6x(o*jUf<;)mjNp6 z&t?3wkP!|NR97sa(1i?Nl-oEA6s{wLGzF7x+w0f1puOm_HeSKQd^CxSPbG^$O1WZ1 z+I*u0ns=zQ-P`W%sxhTSAk5!PX^M5+BjNp%l^bgJs}`@)iM2o(^1p#3^oK`M^#@FmHAe^sd>pwd2<| z{y1fiZVa-Ys`bt;E9(y=6VSdZ>oUx;I&%0-JhW@B!T3%qR^)u2DyP7xJS$F`BQhDX zDCJ~AImmpsK66DTY;NA)n&Bk-)eCOLDk9h(Mh0(ZWas6z6R#$BCxQJFAcZoO8l5sw zRecs^3$&eWxJY(+`M?2cs&dyU>8hSip(bbUL}>VIG%*Y>TqbzVysYttSM zqD?LhBXB;yK{>$?Mz8{v(~~1iwcfz@0hbRwqG&La!F!71IkrAM<`2!}1!TzHKG^4{E-pnsCF>v1I>DR~ZV*p< z;1{`IdG_*1f_$ZA2+_O9v|NQUW0N4tOey59NzX)2)5CX#zG@yG8EKUB-0{!-gmgH3C&dNwZVPWBD_XZoB1A);ZXh8-Y#cko< zau9^IZ@>t0u$cFsi?*IppaK6M-=6JkrIO7T-g-tr*r!ttgNWI8ef14J{7|*|q9X{l zK?#Bi_y1G3ZV`1$Vq?<1^_c(n=NGqE%X^o%2XuU^M*bZ*2qqJ(OiUyfu@{>dl$;zL z^&l?fV65-N;^5^3BICUlhTwBe=n}9yhRmEK7gP&MQ?n0ia5ZByv*gSS9M&%+1i%N% z3AhjPG3|?A(v>F3x2pj=+M$?>qGOc^AqD=V-R9!WJ_6y2+&|tnaCgdy#KU+V@h;O^ zI2VZOyWk|4D`1D)t+dJ+edjBk?+E_~Wj)7$OB0+r1JOtr0X;sR4_FnBq?)0(Hd#~t zHS+gAnZ`&wit&kHC(< z16tfd)Zsn}r-(FKh}eG~Xr2}Dx~N95MY%AMPAzTW!SA?(>1G^1V!T0tJKvoQ2QsE8 za{duhhm^S*R8w>Utv&Z}u?6KQ&4aHG>AK#@26UM;ZZg7@V zuPu*Uuil^i&1?bA`$#}=VWOBzYEK=ZZ<`10K&s5KE};;-hM3E^_?8ZNTgR0h3Rw;$ zzWVKcjrpuMmZH&JxJI1p7Mwat1YO1evzRol77{7+97!NU?pPhK=q$9T8#XjF6h3xT zyC;M=5|Ri7iUR4V%J1oIf%q8CRh|mVyxs|g3_A>?o(PQ44T++5FYMY_DU6qgtPC)QnnN-YA_Co8Y8Q2d+ws5zWL%q}Q!x+aui3*vOo z0WC?&LJf?K?Wf5kTyjYWapj=FZZD@RyO3)wO$ym>pEr6GU_x{LO}{TWV4NbvQb zxex`` parameter can be relatively tight if the ``max_number`` parameter is set. Without a + ``max_number`` parameter, you may have to increase the field tolerance to find all fields. Writing Custom Plugins ---------------------- @@ -315,3 +375,7 @@ API .. autoclass:: pylinac.core.metrics.GlobalDiskLocator :inherited-members: :members: + +.. autoclass:: pylinac.core.metrics.GlobalSizedFieldLocator + :inherited-members: + :members: diff --git a/pylinac/core/metrics.py b/pylinac/core/metrics.py index 68c6db11..95db96e6 100644 --- a/pylinac/core/metrics.py +++ b/pylinac/core/metrics.py @@ -76,9 +76,11 @@ def is_right_circumference(region: RegionProperties, *args, **kwargs) -> bool: def is_right_square_perimeter(region: RegionProperties, *args, **kwargs) -> bool: """Test the regionprop's perimeter attr to see if it matches - that of an equivalent square""" + that of an equivalent square. In reality, edges aren't perfectly straight, so + the real perimeter is always going to be higher than the theoretical perimeter. + We thus add a larger tolerance (20%) to the upper perimeter""" actual_perimeter = region.perimeter / kwargs["dpmm"] - upper_perimeter = 2 * ( + upper_perimeter = 1.20 * 2 * ( kwargs["field_width_mm"] + kwargs["field_tolerance_mm"] ) + 2 * (kwargs["field_height_mm"] + kwargs["field_tolerance_mm"]) lower_perimeter = 2 * ( @@ -496,37 +498,158 @@ def plot(self, axis: plt.Axes) -> None: class GlobalSizedFieldLocator(MetricBase): fields: list[Point] + boundaries: list[np.ndarray] + is_from_physical: bool = False def __init__( self, - field_with_mm: float, - field_height_mm: float, - field_tolerance_mm: float, + field_width_px: float, + field_height_px: float, + field_tolerance_px: float, min_number: int = 1, max_number: int | None = None, - min_separation_mm: float = 5, name: str = "Field Finder", detection_conditions: list[callable] = ( is_right_square_perimeter, is_right_area_square, ), + default_threshold_step_size: float = 2, ): - self.field_width_mm = field_with_mm - self.field_height_mm = field_height_mm - self.field_tolerance_mm = field_tolerance_mm + """Finds fields globally within an image. + + Parameters + ---------- + field_width_px : float + The width of the field in px. + field_height_px : float + The height of the field in px. + field_tolerance_px : float + The tolerance of the field size in px. + min_number : int + The minimum number of fields to find. If not found, an error is raised. + max_number : int, None + The maximum number of fields to find. If None, no maximum is set. + name : str + The name of the metric. + detection_conditions : list[callable] + A list of functions that take a regionprops object and return a boolean. + default_threshold_step_size : float + The default step size for the threshold iteration. This is based on the max number of fields and the field size. + """ + self.field_width_mm = field_width_px + self.field_height_mm = field_height_px + self.field_tolerance_mm = field_tolerance_px self.min_number = min_number self.max_number = max_number or 1e6 self.name = name self.detection_conditions = detection_conditions - self.min_separation_mm = min_separation_mm + self.default_threshold_step_size = default_threshold_step_size + + @classmethod + def from_physical( + cls, + field_width_mm: float, + field_height_mm: float, + field_tolerance_mm: float, + min_number: int = 1, + max_number: int | None = None, + name: str = "Field Finder", + detection_conditions: list[callable] = ( + is_right_square_perimeter, + is_right_area_square, + ), + default_threshold_step_size: float = 2, + ): + """Construct an instance using physical dimensions. + + Parameters + ---------- + field_width_mm : float + The width of the field in mm. + field_height_mm : float + The height of the field in mm. + field_tolerance_mm : float + The tolerance of the field size in mm. + min_number : int + The minimum number of fields to find. If not found, an error is raised. + max_number : int, None + The maximum number of fields to find. If None, no maximum is set. + name : str + The name of the metric. + detection_conditions : list[callable] + A list of functions that take a regionprops object and return a boolean. + default_threshold_step_size : float + The default step size for the threshold iteration. This is based on the max number of fields and the field size. + """ + instance = cls( + field_width_px=field_width_mm, + field_height_px=field_height_mm, + field_tolerance_px=field_tolerance_mm, + min_number=min_number, + max_number=max_number, + name=name, + detection_conditions=detection_conditions, + default_threshold_step_size=default_threshold_step_size, + ) + instance.is_from_physical = True + return instance + + @property + def threshold_step_size(self) -> float: + """Set the step size for the threshold. This is based on the max number of fields and the field size.""" + if not self.max_number: + return self.default_threshold_step_size + else: + # usually the threshold is actually very small + # since the field is very small compared to the + # image size. In this case, we want to increase + # the threshold much slower than the default. + # In combination with the threshold_start, + # this is actually quite sensitive and quick. + # In effect, we are shifting the threshold to whatever + # 10% of the expected total field area is or 2, whichever is smaller. + # For larger fields, this can be quite large, thus the 2 max. + calculated_step_size = ( + self.max_number + * (self.field_width_mm * self.field_height_mm) + * (self.image.dpmm**2) + / self.image.size + * 10 + ) + return min((calculated_step_size, self.default_threshold_step_size)) + + @property + def threshold_start(self) -> float: + """The starting percentile for the threshold. This is based on the max number of fields and the field size.""" + if not self.max_number: + return 5 + else: + # start at a higher threshold if we have a max number + # by using the expected total area of the fields / image size + # this offset from 100 and adds a 1.5 safety margin + # E.g. for a 10x10 field, this might result in a starting threshold of 99.6 + return ( + 100 + - 100 + * 1.5 + * self.max_number + * (self.field_width_mm * self.field_height_mm) + * (self.image.dpmm**2) + / self.image.size + ) def calculate(self) -> list[Point]: - """Find up to N BBs/disks in the image. This will look for BBs at every percentile range. - Multiple BBs may be found at different threshold levels.""" + """Find up to N fields in the image. This will look for fields at every percentile range. + Multiple fields may be found at different threshold levels.""" + if not self.is_from_physical: + self.field_width_mm /= self.image.dpmm + self.field_height_mm /= self.image.dpmm + self.field_tolerance_mm /= self.image.dpmm fields = [] + boundaries = [] sample = self.image.array # search for multiple BBs by iteratively raising the high-pass threshold value. - threshold_percentile = 5 + threshold_percentile = self.threshold_start while threshold_percentile < 100 and len(fields) < self.max_number: try: binary_array = sample > np.percentile(sample, threshold_percentile) @@ -558,32 +681,44 @@ def calculate(self) -> list[Point]: Point(region.centroid[1], region.centroid[0]) for region in fields_regions ] + # find the boundaries of the fields + # this is solely for plotting purposes + # these will be bool arrays + # we pad the boundaries to offset the ROI to the right + # position on the image. boundaries = [ - find_boundaries( - region.image, - connectivity=region.image.ndim, - mode="inner", - background=0, + np.pad( + find_boundaries( + # padding is needed as boundary edges aren't detected otherwise + np.pad( + region.image, + pad_width=1, + mode="constant", + constant_values=0, + ), + connectivity=region.image.ndim, + mode="inner", + background=0, + ), + ((region.bbox[0] - 1, 0), (region.bbox[1] - 1, 0)), + mode="constant", + constant_values=0, ) for region in fields_regions ] - - # print(boundaries) - # print(marked) - # the separation is the the minimum value + field size + # the separation is the minimum value + field size fields = deduplicate_points( original_points=fields, new_points=points, - min_separation_px=( - self.min_separation_mm - + min((self.field_height_mm, self.field_height_mm)) + min_separation_px=min( + (self.field_height_mm, self.field_width_mm) ) * self.image.dpmm, ) except (IndexError, ValueError): pass finally: - threshold_percentile += 2 + threshold_percentile += self.threshold_step_size if len(fields) < self.min_number: # didn't find the number we needed raise ValueError( @@ -593,33 +728,25 @@ def calculate(self) -> list[Point]: self.boundaries = boundaries return fields - def plot(self, axis: plt.Axes) -> None: - """Plot the BB centers""" + def plot( + self, + axis: plt.Axes, + show_boundaries: bool = True, + color: str = "red", + markersize: float = 3, + alpha: float = 0.25, + ) -> None: + """Plot the BB centers and boundary of detection.""" for point in self.fields: - axis.plot(point.x, point.y, "g+") - for boundary in self.boundaries: - axis.imshow(boundary) - - -# TODO -# class GlobalFieldLocator(MetricBase): -# def __init__( -# self, -# low_threshold_percentile: float = 5, -# high_threshold_percentile: float = 99.9, -# name: str = "Field Finder", -# ): -# self.low_threshold_percentile = low_threshold_percentile -# self.high_threshold_percentile = high_threshold_percentile -# self.name = name -# -# def calculate(self, image: BaseImage) -> Point: -# min, max = np.percentile( -# image.array, [self.low_threshold_percentile, self.high_threshold_percentile] -# ) -# threshold_img = image.as_binary((max - min) / 2 + min) -# filled_img = ndimage.binary_fill_holes(threshold_img) -# coords = ndimage.center_of_mass(filled_img) -# return Point(x=coords[-1], y=coords[0]) -# -# + axis.plot(point.x, point.y, color=color, marker="+", alpha=alpha) + if show_boundaries: + for boundary in self.boundaries: + boundary_y, boundary_x = np.nonzero(boundary) + axis.scatter( + boundary_x, + boundary_y, + c=color, + marker="s", + alpha=alpha, + s=markersize, + ) diff --git a/tests_basic/core/test_profile_metrics.py b/tests_basic/core/test_profile_metrics.py index c05baf96..6e8f60b0 100644 --- a/tests_basic/core/test_profile_metrics.py +++ b/tests_basic/core/test_profile_metrics.py @@ -10,7 +10,12 @@ GaussianFilterLayer, Layer, ) -from pylinac.core.metrics import DiskLocator, GlobalDiskLocator, deduplicate_points +from pylinac.core.metrics import ( + DiskLocator, + GlobalDiskLocator, + GlobalSizedFieldLocator, + deduplicate_points, +) from tests_basic.utils import get_file_from_cloud_test_repo @@ -34,6 +39,43 @@ def create_bb_image( return as1000.as_dicom() +def create_open_field_image( + field_size=(50, 50), + offset=(0, 0), + field: type[Layer] = FilteredFieldLayer, +) -> Dataset: + as1000 = AS1000Image( + sid=1000 + ) # this will set the pixel size and shape automatically + as1000.add_layer( + field(field_size_mm=field_size, cax_offset_mm=offset) + ) # create a 50x50mm square field + as1000.add_layer( + GaussianFilterLayer(sigma_mm=2) + ) # add an image-wide gaussian to simulate penumbra/scatter + as1000.add_layer(RandomNoiseLayer()) + return as1000.as_dicom() + + +def create_multi_open_field( + field_sizes=((50, 50), (50, 50)), + offsets=((0, 0), (40, 40)), + field: type[Layer] = FilteredFieldLayer, +) -> Dataset: + as1000 = AS1000Image( + sid=1000 + ) # this will set the pixel size and shape automatically + for field_size, offset in zip(field_sizes, offsets): + as1000.add_layer( + field(field_size_mm=field_size, cax_offset_mm=offset) + ) # create a 50x50mm square field + as1000.add_layer( + GaussianFilterLayer(sigma_mm=2) + ) # add an image-wide gaussian to simulate penumbra/scatter + as1000.add_layer(RandomNoiseLayer()) + return as1000.as_dicom() + + class TestGlobalDiskLocator(TestCase): @classmethod def setUpClass(cls) -> None: @@ -350,3 +392,96 @@ def test_shifted_bb(self): ) ] ) + + +class TestGlobalSizedFieldLocator(TestCase): + def test_perfect_image(self): + ds = create_open_field_image(field_size=(60, 60)) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=60, + field_height_mm=60, + field_tolerance_mm=2, + max_number=1, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1) + + def test_small_image(self): + ds = create_open_field_image(field_size=(10, 10)) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=10, + field_height_mm=10, + field_tolerance_mm=2, + max_number=1, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1) + + def test_large_image(self): + ds = create_open_field_image(field_size=(250, 250)) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=250, + field_height_mm=250, + field_tolerance_mm=5, + max_number=1, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1) + + def test_multiple_medium_fields(self): + ds = create_multi_open_field( + field_sizes=((30, 30), (30, 30)), offsets=((0, 0), (40, 40)) + ) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=30, + field_height_mm=30, + field_tolerance_mm=5, + max_number=2, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1) + self.assertAlmostEqual(fields[1].x, 613.6, delta=1) + self.assertAlmostEqual(fields[1].y, 485.6, delta=1) + + def test_lots_of_fields(self): + ds = create_multi_open_field( + field_sizes=((30, 30), (30, 30), (30, 30), (30, 30)), + offsets=((0, 0), (40, 40), (-40, 0), (-60, -60)), + ) + img = DicomImage.from_dataset(ds) + fields = img.compute( + metrics=GlobalSizedFieldLocator.from_physical( + field_width_mm=30, + field_height_mm=30, + field_tolerance_mm=5, + max_number=4, + ) + ) + self.assertEqual(len(fields), 4) + + def test_not_from_physical(self): + ds = create_open_field_image(field_size=(10, 10)) + img = DicomImage.from_dataset(ds) + dpmm = img.dpmm + fields = img.compute( + metrics=GlobalSizedFieldLocator( + field_width_px=10 * dpmm, + field_height_px=10 * dpmm, + field_tolerance_px=2 * dpmm, + max_number=1, + ) + ) + self.assertAlmostEqual(fields[0].x, 511.5, delta=1) + self.assertAlmostEqual(fields[0].y, 383.5, delta=1)