From fda1d1640cbac0b63a55bd0de96a1a9983254bca Mon Sep 17 00:00:00 2001 From: Vedrana Andersen Dahl <vedranaa@gmail.com> Date: Wed, 15 Mar 2023 01:55:09 +0100 Subject: [PATCH] added chapter 7 --- .../Chapter07_mcode/GraphCut/GraphCutMex.cpp | 87 +++ .../GraphCut/GraphCutMex.mexa64 | Bin 0 -> 22225 bytes .../GraphCut/GraphCutMex.mexmaci64 | Bin 0 -> 18052 bytes .../GraphCut/GraphCutMex.mexw64 | Bin 0 -> 16384 bytes .../Chapter07_mcode/GraphCut/MaxFlow.cpp | 690 ++++++++++++++++++ Chapter07/Chapter07_mcode/GraphCut/block.h | 268 +++++++ .../GraphCut/compile_mex_functions.m | 3 + Chapter07/Chapter07_mcode/GraphCut/graph.cpp | 114 +++ Chapter07/Chapter07_mcode/GraphCut/graph.h | 506 +++++++++++++ .../Chapter07_mcode/GraphCut/instances.inc | 16 + .../Chapter07_mcode/functions/grid_cut.m | 252 +++++++ .../Chapter07_mcode/in_region_cost_example.m | 25 + .../Chapter07_mcode/on_surface_cost_example.m | 44 ++ .../w1_unwrapping_image_exercise.m | 18 + Chapter07/in_region_cost_example.py | 38 + Chapter07/on_surface_cost_example.py | 66 ++ Chapter07/w1_unwrapping_image_exercise.py | 35 + 17 files changed, 2162 insertions(+) create mode 100755 Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.cpp create mode 100755 Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexa64 create mode 100755 Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexmaci64 create mode 100755 Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexw64 create mode 100755 Chapter07/Chapter07_mcode/GraphCut/MaxFlow.cpp create mode 100755 Chapter07/Chapter07_mcode/GraphCut/block.h create mode 100755 Chapter07/Chapter07_mcode/GraphCut/compile_mex_functions.m create mode 100755 Chapter07/Chapter07_mcode/GraphCut/graph.cpp create mode 100755 Chapter07/Chapter07_mcode/GraphCut/graph.h create mode 100755 Chapter07/Chapter07_mcode/GraphCut/instances.inc create mode 100755 Chapter07/Chapter07_mcode/functions/grid_cut.m create mode 100644 Chapter07/Chapter07_mcode/in_region_cost_example.m create mode 100755 Chapter07/Chapter07_mcode/on_surface_cost_example.m create mode 100644 Chapter07/Chapter07_mcode/w1_unwrapping_image_exercise.m create mode 100755 Chapter07/in_region_cost_example.py create mode 100755 Chapter07/on_surface_cost_example.py create mode 100644 Chapter07/w1_unwrapping_image_exercise.py diff --git a/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.cpp b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.cpp new file mode 100755 index 0000000..0806d9e --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.cpp @@ -0,0 +1,87 @@ +#include <math.h> +#include "mex.h" + +#include "graph.h" +#include <vector> + +void +mexFunction(int nout, mxArray *out[], + int nin, const mxArray *in[]) +{ + + if (nin != 3) + mexErrMsgTxt("Three arguments are required (nNodes,TerminalWeights,EdgeWeights)") ; + if (nout > 2) + mexErrMsgTxt("Too many output arguments."); + + + int nNodes = (int) *mxGetPr(in[0]); + + + const int* twDim = mxGetDimensions(in[1]) ; + int twRows = twDim[0] ; + int twCols = twDim[1] ; + double* twPt = mxGetPr(in[1]) ; + if(twCols!=3) + mexErrMsgTxt("The Terminal Weight matix should have 3 columns, (Node,sink,source)."); + + + const int* ewDim = mxGetDimensions(in[2]) ; + int ewRows = ewDim[0] ; + int ewCols = ewDim[1] ; + double* ewPt = mxGetPr(in[2]) ; + if(ewCols!=4) + mexErrMsgTxt("The Terminal Weight matix should have 4 columns, (From,To,Capacity,Rev_Capacity)."); + + + typedef Graph<double,double,double> GraphType; + GraphType G(static_cast<int>(nNodes), static_cast<int>(ewRows+twRows)); + G.add_node(nNodes); + + for(int cTw=0;cTw<twRows;cTw++) + { + //Test for nodes in range + int node=(int)twPt[cTw]-1; + if(node<0 || node>=nNodes) + mexErrMsgTxt("index out of bounds in TerminalWeight Matrix."); + G.add_tweights(node,twPt[cTw+twRows],twPt[cTw+2*twRows]); + } + + for(int cEw=0;cEw<ewRows;cEw++) + { + //Test for nodes in range + int From=(int)ewPt[cEw]-1; + int To=(int)ewPt[cEw+ewRows]-1; + if(From<0 || From>=nNodes) + mexErrMsgTxt("From index out of bounds in Edge Weight Matrix."); + if(To<0 || To>=nNodes) + mexErrMsgTxt("To index out of bounds in Edge Weight Matrix."); + + G.add_edge(From,To,ewPt[cEw+2*ewRows],ewPt[cEw+3*ewRows]); + } + + double flow=G.maxflow(); + + std::vector<int> SourceCut; + for(int cNode=0;cNode<nNodes;cNode++) + { + if(G.what_segment(cNode)== GraphType::SOURCE) + SourceCut.push_back(cNode+1); + } + + out[0] = mxCreateDoubleMatrix(SourceCut.size(), 1, mxREAL) ; + double* pOut=mxGetPr(out[0]); + std::vector<int>::const_iterator Itt(SourceCut.begin()); + for(;Itt!=SourceCut.end();++Itt) + { + *pOut=*Itt; + ++pOut; + } + + if(nout==2) + { + out[1] = mxCreateDoubleMatrix(1, 1, mxREAL) ; + *mxGetPr(out[1])=flow; + } + +} \ No newline at end of file diff --git a/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexa64 b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexa64 new file mode 100755 index 0000000000000000000000000000000000000000..a9df40a4455e24c1d68dffa4196f96e4cd9f291e GIT binary patch literal 22225 zcmb<-^>JfjWMqH=W(GS35bpvPM8p9?F`RIOG8h;b92hJZ*cluc<QZhZYFQCd7_1Y< z5Fr@Nz@PvY0NKIDz`)R93+6I(*h0i%v;Z?iTmeeM^f56oFo5)Jfa==-)d!;$KnfTb z7+^F9#25w^1_m%b0kW8ZfdQRXfvSViAoU=jfTtxXAhv)Zh{wQyPM1K08DKO>El4Qv zX-Nvmogg+b80O9(sC~G?r2^_d7!CCf*l!{X3=E*)a`y{msPx{}ZGEp&Bm9N)mn=5- z?S19S;20HTU|^7dI*1!egY*eAFfa%(FffQRFfgz)Ffi~y<v1A_7=#!Y7(kl&pkg3D z@<Q48D3Du_<ze9g@+&CKh(jgC7#J8p-r{FqU;w!fDh#GT;vjjDKSBBhKnbz|l2$<C z>zF&ZO&e!!11rWsFf(v7OaOTVNq#-lUr0h==4l+_Ap0S<f|bDuc~Dw|^AL<E9On4p z5Vyi14$3~D@JE;oVTIyQ4~j!2h!PYs0f+kSIK<s?h-cssUxY*49*6iP9OAEWh&SR8 zZ^t2C$OuYHAOE8S0s{lXG_ZPJh6T{X9}Y^0AVu-<>ACrN@x>*HMJ4g^4Ds>F6^Zd_ znR$shnN_I_@m^)|L8<AP#U-gl&N+$2#i_*%xrsSB`N<5a6`3UrMX3-@d{jzeKp8`P zR7y@j8ADo8YAQo~R9<2(LvCtrZhjev7hGa!9v@$lQIuaEpOlyq57r3Ml2@L~5FeUX zo|%^tACy{Lnw!dyn_A&knwMOXnV$z@yA~Ds7N>_)lrZF0xTlr`6hUa0%-qzx;>`TK z;&^id1CUk4sURiJMX8A;sV@1YNja&$i6upu6=3=J`1Fd3_=41;;(U<rN-E>a3>b1U zlX5Hciu1vYR5&ZK1S(csl9H^gtyi3{X9iXRXQd}6#}|Wmh9L2Rl8mC%#1s%8Y-Tcu zWyauMlvt4ATw3CrTA`PlT9IgG!r<!~!Qk%f9B-^=z~CGa;ppTUZ>VR?;O^t;<Qx$Z zZ>(nmVaFTk8R?mVVhCEgfy!|f24*;f@F681D+3Dy8v_#q69Xd@K8cg%29<9z3=9l$ zpiIob0Ist@DpmwCFo5JhhJ(@zhz991V}z82p!^=7R>=*jaZ2HOnHXfD@&*hH442nu zvobKSGSouFN;D^NfLzo8SI@{00+j=qUn-T!33BTcsC<V5$QA|$P?-plhw(4iL&QO4 z9;mzp$${_<Bym{M1&KdE5{D*jkTAmpsQDo4KxI5w5JW)5L1jItj0FjSFb9%2a``QQ zBo1;DNG%9UAc=#_0*QgJ0+Kk)Zy<3EByr>-*#Jo#RM&vig0KaWII_PSki=ni6-d1Y zk~nf56@VlTtH(g<Bap;F^$$o4gcFd&L1ijP41_a~#F6W&0wi(d^11>^9G2EV<}@IQ zgX&z67zlebzu|ac&d9*<;u45BtHQwG(R_sC@PAV$GX;kKsxD><4E*vA4F6R@d{Al3 z@bba`|NsB1TA3*@WPqai<pnU`2*d|P?aKpTz7~iNiqe-Gz<ebTAC$CSE&%hTKzvY= zdN~2i7XtA?QT?(3%;y5}K~em&0L*6s@j+4hG6Bs0WvalC0gBR>0bu?o5FZqkFCD=A zS0Fwp3SSz4`Hw(+P}IFt0P}Bw_@F3zDFEhQ0`Wmn^^yV1KLz50qUhy^e<1%J0`Wmn z^YX#J|L}4i6n-AvtY1wO7(7}Jlra6j01mzhFc8{cyTDX|fx)BmsbffRX_7~??F|zJ z28I$Ik8ay?Qw4?>LjV8&_vn>vHd0{l=~XQ;Qeg1tW-T>UVDQm=>~Zinvq$r5Mi0vy z#R~lLE}+Oc!06G*8fvP*;P_wUh(|Z8m#G4Shvki;r^i^$OcfXyqa9;T9!}$zZvjOw zL#Riyt-cY&0VPr%-K@b-m-2b^+KL+~Fr57V02IF`4?|oW3UX+;W0+&8V@R+^;~Ne} z28PbPPyYS?@6oyS2ADeC{ENAS$D`MFgE9ldi^c!`|8E8>V5$QNvVa8pKnlC9UmW}U z|Nkd`9fr=<KmY#!_vuaM@a(oxd13$W|Nm}*&Z!*#|NnRC&R_ux&iecRzh~z$$L7Bf z(`SJUas?UGd|1Gv+xi73RM(Gx5M3a#RIsKaj?I5rz)o5V(j)*<I2ELf8El0|?_N+8 z?*s)%_g0WgUPQ7oFnA<`UFFd$<KS`dF^fm%4G-p42~ccUwt|vKu@F=vRL^Tx28QOn zAnlB$E*{+wrN=>rGQfEb8X!KD+mN8}7s34P(G1pN(A)}2pp388JvwjggDBq%vi8M_ zKmY%CgWbdkR^q`p7bM=i7Zgqm{Cx))7#O<2K4QGoITe&jeLKJI2AS~5m0zcp!Lj*0 zV+E&+HQ0?cpvd^9#K7=^@$dit-C*UtEW#c$J~{FWupH}d1sTeC1H$Jy=9~QAyO*Wj z<KP1pkM60U;QGWb5W=L9#;+&ppv=Hf#`cLn;?!%FPy8{5()cwle&Ub4@QFWSYQ%p~ z3=2Z(m=kc#f<8<d&4(F1n~!o>&IQ@U-yilLlvF_R!oc4LsvNspK^(>d-QdV#Jk&WA zl+`+0zx@0E-?#G_*lQlWtsvKUbhf_u_y0f0qjrv+cOg0u{yFlAUl8mg%rH;_g~8w7 zAZO};l{CL#>;y9^xVoo;qSNvMIIuK93a|VI#R13}E{rcAfx`ke92DFx%xl3;OMdI$ zt@FmCyVU}cT0jL%DUV0@UXT-C)Uq&u5(G%lqc=doqu0a1qq7&x?d}En1sv8U^&SVG zv3T^F7<wdwBiN(YrrzV=a~2QApB|k70v?@C0v_FV9-S2&9-V9+-IgAmRU96jf*ze5 z9-T%$of#6J_yt-)uKL6;=p+FS>SG+1b3q}>-+u*^Zy?Ej7bt;%)7cAsNd^X3KzST@ z-2jsDIPSUwN`ulDG=($NB{}W^34`Onr?>QsPjBu9pUwv!%#i44uHaxOVf*Hwz*wr; z>&)2cdd8!>bccs^>46eHkJfJ`Or2~V&8|CGUUU0&e)8#j|Du$GfdTBuBL<dpL7`s1 z(WAR|hlh6U29JaPnINtNRVObBfBpaO(LEJX*f4o?9)D5)6C@J4!NW3iM~T=cel4a? z{DKS~%?CI<nh&vfScmQ?eFZIpO1(Y0=Yq91ALc+9x*u$^ILPEqW{{COKS2pX0Gtj$ zhP@W{=sfJvUAw`f^?!-{3)`Rn|G$=i@s-miywLvn|3AMx!)tMnI!MUAn9mF<*zdo% z4{{$=tQ#cOcmz~#Mmxqh#yZA1#-BVq0aUJf9Ct4Og{jAJcTo3;!Q;4l1(Xl2SUrxn zH-PvI3=GHHCqP*q-RvIL?i~EhOBfg!()hzqf^wh7aR*RM#L$qS@Uwe5NNHyPC_OlU z^6SAegN6i!e;&sjK#pMW=nk*&XuVXz21-${Sv)%3Gdw!oD?B<K1dcn%fWqaqT^hgM zg))N%4+TaNba8;xf?Vm*e89lN(w&39J>~!Z|0fTF{n{D(!lyI!gHLDf50B<I3Lc%m z9XqcacL%j$89a`=TYxoo`a3lL<KF@D3%}oo=3o3S{C*ESlE3-%>XdmNe98h!IG4VI z>K<m_<R_j7pRoCKK7VllEaEKT(Jc(}7^_F?fl?WdZf6Hj4b17$Ya1@Z!0`V9B+i?` zii0_Pn6-U6vjsdlofZ6gRTO+WlQoVz2Y_91+}Q(6cgMc)=}!IN)1CVR8e_0}u=BT1 zx4v)Z8;@>x0Uym59tVH0WFA<(YSpUFLkC~UGakD5tCxj|x!l0=1HaEPj~PC_CJvql zpRxFKegpYn0@w-67|z@vjd13jZ~y;$BnKOK^oj&|9DKy$(JLb0(HYF)*K1<n(^)Lw zX?e6n)}z~5!o%8Gp;*|XTi8eQkq4uINAm#%4@+l-^0&~`26H*Q4sU+L;Q?}gz68Vr z8ek8+VD;&)S8(h+=6LX(yd&c=$BX|Q8IQSierSHc*!+XhvH9MA-{fze2cNM*J#iZ3 zz)oig3}3vGLij@L+yDQ*n(sXi{%7?}zV31G5i8Um3Q&J&fcz2R(ar2(>C90k<Iybu zia{=rUzmd>JUT)CF#ySxKY@Ez)1y0p1Eh=_t}KAVqdQmv#AOTsX@La6bx_j_8Vuk} zH(>%e-$B|k&2Kn7JOBE0r(5{8J}EKt?5<Ja@M!&CBID65529IL?D+cszh`F*NceUM zt4HVk*H1vn%@8Df+!dT+J-cfeUIc&yx@#p~xPSfsA8sc+-9PZ@Onu?gnft+~vlN{6 zJv;w+c7F9Z?#2OX9(#1hTX-~A3o!7vf(99Ux{Eowt()I7dL;h`7yr%2C0?xi0!s7F z9I(I!g?b4nW!Fj~0(&J$QEvc;M`y5rM>4m^!G|m!y)GO-JS-0t%Xl;&67c8@`rwhw z?9u#!5nQh`I}3PN9x8g`aUA3h29IuM4xdf~k8WoHP&{-73;1+ucy>=w>0n@B0EIWG z#Nc?%;n~@u0xFk4G+64m1E?9!;L$Ck;?td?!r{}MqaxtbU7{iZ2?Ka|dvyNx?EK=< zZRXMXjlU(FiGjg!=l}oz>m9)ba~gm60Y`BC;<y7;40&|B38eAsIY?Bo@wY4mxuKiy zo1+3Ff9nh|ClHhf_*=R`t)ShYGS9M7qD0W8^MFV5FUAsCkLF)YrJAqvT{>?tJIcJ) z^ym)cXg<IQk}T!vtduzJ2r29YKnhDdJi56oD<w*AzP@A&a#D%8i{%M^d4^gx$IgSE zou53q#k#w}mUu~&NOyWkbUMj^!m&Gu!=qaWE>-%}qq~d)<h@pqF}|H4DjdG8mr5kN zn?af#J6%*H9Ged^x_0`g$h=HpWMJ^=&W`YHeaqi+oq>VD@tcc^0t<gDXcW-1+f@Km zfGL0?15~Y925|5<H-o+I0G46%>~s})UE<j7$^df843A!$R#30~Mg6D$|G`Pcv)fez zWCkdy1c@V(O5Lab|2?}y!I=#ds{*fCJUd-EK<0V$27vM%$ac@<PzH~KkC;4|1w1+f zI6y+MiUzE?o5`csgT<qh%ct|9hc>sz!JjN1yFjJjs#U92dGN2l*eT%H;G?3z;@hpG z;@kSRq};Q+l;NA30%K{gN4Fa&HJ5mSOajGx8o!>CL={tsl1Fo$07I#?M|Uu+nko_i zRa3<o9-YM%9-U<Z9=!qoLD5*o>tm@<%Hh#@|HZ^V|Nr}T=W_VAzAc^gn&0!diwZbR zPIx&7)beRQ!T~M^PELT<yq&Q(d^%Gf_;lvJ@aZi5;L}+P%9#$Hoqs$!FM~2?w>yW2 zwYxxxB`Ek<y954o3P4JTUUwGHP8Su9ZfTF+@IOA9uRNL`{qfMe?$i11Mf%78{~bGH z8D1JPf|`8d9-0?C4*p;QMTMMa=Xa0J4<M1wLmmfzF?n=`9`HEs3Z4ge?d;LZ!{A|g zg1`L&D1LW<I-LLi|F`$BJm|sicF?2K^#nHkj^JjTM|157h7wf|Bwem2%2;1;egw6W zYEO6^e9rX3<|C*qJLb_Hdcd=nWg5s=EZxE&PnoFrbpH0x{Or;E_>V{DEuYTsFXn&v z|KGDWMvld^+eL-LNArhg=W);G=YM=UfBAI21BG|(hyVY57=L?q+o(XwN|_I!Cc23V zq`d47Jpr->6i@=M)v-DKxa$d!Cy%?H0r~S~2&m!Q9ecy4JN1E2ckT<H?$QrF-L>#~ z7peA=@a+8S+4=nwzkr*>Cw@WK6CT}u0v_Fh9^FnH9?d@(OB6i10|Y+t3ph!9;*WJb z0aYUalP`S$lG8xQd3KA#<Vs{d@e4TM(DwwYPk`554BYfOh%n<5e}p5*HI5RX?v%U4 zC;k}M6Q3OU1>Ge)L2ZO?1<!6!MbiA6k-vqX5!CD!068tf8Keu;ba4hb8Dup`Bfp@7 z1gIQ_wD0shB=}oTFfcIicO7D2U;rm>P>Jo^U25S8$}Rk@T_AB#mg;mp;Mi@4R(<5Z z14TYK!<EW-cDuo=4?AH*nUM7kH0}_<;n^7`;F-*az54J>W`k89p3H72)rTj@`tHyJ z9-u_!2C70pCi!$Kcy><(RTtnQq_Y*o^yyYc1V`&fP(RL*U%*9010HNHDxj7t>opMu z2L9IlV6TI^>u|4ocE=K|G%MfYs5Cu<5WX$}DWY<viN(jh-RcM*w=jW=2XKso69FjN z1wn)WC^<k9fU^WN@(=R2f;wXSU7$)59so#j-^9Sc;M-ltff?`RZ*auBn;;?tir!Ej z?>?Otp53Ix`9Dyc^MkU0BYy-a;vM+~p>h65h=Mq8eT^f|Lj({$ZU8BweVl)S`Onn@ zR1&y)cr+i809QlMJ|Jj>4b(Rxvwf7tujdAC8?6M@8PK*-%N#HV+BWJ(X&Z^c+D3|? zCQqp@Xq<0EwvXapVrd^4@gY)3+{=M&AC(~5N8q-RBfp-jL=`K)JOi2Sqr%Tbw~s)n z5`n?#6SemJ<Ix?TLB!YvsBgdkN{xtNjArmK2BSyoZT^;NYzz!O-PHm<mcPLRPCooD zKRk{*gZlmq4ao|BJ-c00DnRWNHjmB-f!C~{mPH1LsPO2F233~M8i*3|9k_~k0kW;R znuCGA<t}JI_qa2-hvnOyqhjINTr0qkHi5sj6<l1FayWL^GN9Fy?$1Hxr7Ngm3aTev z;q{~#4<c<iJ%^Q-p#q-Cz6>Z8r7WzU<(cftfT$^5!8PS$Xh+Ml+ZEK&0=K8T1w1=L z1w1+>Na$$2WMBogKlE6)^D;2-x8*W{rU9xs7)rRlIV&*ow;X3-0EOT`pKfh%&mxr( zR6&$7IDT_eVB&A}WMp6f+3jiR#=+nE4b*agc11nALsTq0Izvx*G*^N8qjN#w+3m&Q z(dm1j+X^kPq@IDAm!&5>I-v~<52S|0Q!Ye+34#=Yira4A6CRy`0v_F^Cp^#^7z(g( z^8hJIW`qTu2eTtM=pI0W4y+f^&;YCMb_6vxKqmWi8jui-p3QX}3?+7;P+=@F2esQ= zCDQoySQB{|7|PVa4XaWGcq-)Z=`0XPoA8>|x3feg14LB#c9sc1J01M3vq59j-8m|d z&YTM~1H(6GP*={184}u{luSlfZp#xaMO8Z|B6K%Dp=wd((Rs^LTSUd<;4c=BeMo(| zYn=iO-U=))v(XY;YZMa$LxZ;hBP_9jrT~4r;fc+MI*D!VV?2qif&&p8OF#+-ZDMl- zC$>^6Q1%1&G(hR?BPTiOZ3aquyTSlVZ>K={nf##xmq%FAn>afnfSn)FGrf5ncLony zJ9f9IfCf-Hoh6PtgU4-rTmP4Y`E>q%xdSxV-VN%vOaTu|d2|~<%CT-Q0r2FAlYnRQ z9<UTcDUavz7Vvl$sA2+*M3iuWhV;t5{r~@(1H|%q3-TUhfD6=bUy9af!%`pq_vpOr z&{2Bcqq*LKp;Tu-cnIdqLy$S)4)8H?=m`CB22cmcqccFlv(w$evok)zvAI5gp;W1% z_B=z0jAyqyeAN0SD<ag-J^)3_?GjaxvQky3GEjQ}G_>;{N#)`P|Nny&fjsvAh3yZB z$Ns<Qe*XVIxSoB@3Z_8ikq~@1r}L->Xa=V=)3-Z60W^rB>D!&3;oF&?;Mwh7;o0e4 z;Cb8~JSpJ`oumNQEQTzQ!p+n2XsM7VbUMQmI-LO;Q*ux6Jnjmv(mk7NIT%V7fVKT% z1`Tz&f`;-|gUg3l25^(Wn*~&_g&#=c*K-E9Y^%WoWXhh+H7Xno{H;2S3=ED9DJly8 zK%HaI=mLjhx2ptXRIS@pruhIPh%EzBFTen*!rdV)2%pXbf!8d)o#`3Ao#_?6ohd2; zzTM>!zMx@9agd?VaW`(TYuv!2Yp8`aXlT}>8$RylCh)q*v6~Hb*p1^pD3T6<>~(_= zyWM6+#3K8BP@)J!9d>i$@a%L0jcGxK-P{B`lf_U*-P{D=(bOF#;DI#i_Q+FP99&ZG zL@KG_qi(Mqz{81*{H+f`UhD=9JzBd;)N49+yQs)GTD!{B$>9%2kM3v=P{g#R|NH;{ z<yt0CGa{A2xAh5s%K|1)%*Ut*fZM0A_=2=FOjN*CD}SpB$leA&1x8RUD&X1a28}K^ z2}f9Lxye8(T4_)^1EoneP(2HZtMUxr&T?=a>)9P20dhtQ8_2wGZlG%QH@Ioyj3uVP z`3n|P9-YntuM4}GQ6nn;E;OP%x}D(>WyFMts5p><?qKAI^5}Ht0GS3I9|Z-RN3t+d zM1e)1x!eP!7b&7X;fyFh1(w%F;98-?&avA~!VwfhijLiGkO&g@?2N+}XuhCvc2GpL zmV=^32ox~LjcZBJcxY#?fJf&&pU(d;c7O%~y1}!FpfE1==#GHT*1h%u&8K*FPKN}X zfoFHPg$HukxgR9KP|EFjydBo0hn1ZF!Q+aclJoxyp%<XU4YQ+km1pOCkVcPI@Z8Cb z*WEsyFJ89&hL%p8j^Eq^?9(QcI`~3oVDus5tF8Y*^IFgu7_&Q|203Wpto1f%!t(xv z*N;J^^t#*s|9e6w4IP`^0~kTAH+S?7wj(q|UUL8b{~y|=0hQmu&_y~N{0a;n%|HH^ zNP6_j{^nC);Fo9cXnxJ{{{m=2s`<zN(vRTT9MD=Akd_zrcOdgShe5Rm1pDQeD5U1) zm!@YZ<fi847gZ`Ul%*z@<QG|4#rwu*<`t(FmBc5ORxpI*=PTqU=2a@>mzETimMA0^ zrI+TW=9Ludfek83EiF!sFDXh*EmlZQ%*)F!QAkQvC@oG+QOL|wD9K1wNXslLE>TEM z%*jzG0U1@1mXlwup{WSg&XAmun3tZK5}%U^R)k_ed1gsQerbsUk{wW83=AO|MX9L> z$0#Hgr79Gq7M5lfrKTuo<oV^Nq!#Ohq!#6ZR-%TdW~OJ96zjOAq^CmJnh@(5LNZbn zV9FFAatgVLC7Bfp#TogfIVlPmiDjt@#tO;#Ii<OI#X1TaAft7PGxM@_it|g0l2bLo z0YRKOCP?PE73JsZgyicuCl(|oXO>jz1f`b6!#EH_85kHc^HNePKw+<tpQezMUz(Rv z3<_A7lc6q#tWpIz1Z1iL4ox8M!~6)-5t2`eYG|M#C{RuT^BAR>twC$~K?~6)-1`5& zfQf<O#r^;PK~r8=9zd$_Sr7mJzrn!35c25%f6!vmgvbB?YcMh}oB_3c85tNZJo*2> zfRTYA>G}Wv6Brp7gkJvt51LW?@b>@z2aF616F&U^4_ffP;QRmo8lZt>1_lPu`c#ku z7^?z;7#J%A*k*s@Q(%zhVdnsu0FejnP4Kw?|Nk41LPtIU=Qbu@&;mUW7o^{Xfq_Bf z!T<k9K>{%O;9lk)mS$Fv-V_D~hFuT-{|AlaGQi}$n%R0-ds+HGD+3uA7(nWJ7#J9Q z9{&FiY576z320{TVe4h>W9bL01lhNRfq}u`(f|KB$oBa(v-Ys`GK0(kHHvE<{r}Gb znoWe+7v9I*%hJQz40SljoIeZ<3=xn2|DOm73z#}2b3pCQQ&0c@Uyp1~bU$++OD}5= zTQg`)BG}z7j0_Ax&;S3w30l<%GY5;ia~K&IINsp4uZNL=q2t~E{{kSBahL=0-xfv& zhJvsE|L+8?n#7?F<UdgGFjfWlfFqWDwk)p#gETB=LFRz4#Et*|L2J=jq(Re=&@fLx z4fFs186^yMGcqvDVmi*ch>PI{#}vk2oD5sIK=y+y8O5U^Fd71*Aut*OqaiRF0;3@? z8UjNm1VD>%Kp5s=5F3PHd>9RC0)yF529#ppgfgHMgD{i<rNDdEK>A>O&`=&S-wY}N zr5HeMX=HiO9ySmkv|$HC!}imH`bi-E3Q(sFy2l#C{{S`U>%af`AU<fk3Pi)2svJ-T zv>U*10@U`0?lT6dTL9J14dNi}Sq67L7#J8{K=->rwK2?qs)y~d0yW)1T48&Pe?SGG z6hi})0bZrXz`y`r_|3q;0NTI=(g)ik{pUZ#{WqWv{R`#8?EMGjgZl~$450P@VE@A0 z`2i}w0;>Nfln-?-$SS&E*#3BQ%l1J%d=5(AgVOJy^gk%g2inQSz`!5}rS+h+9hCNi z(s58a4@%cT={_hu4@$3t()*w^Y~LfyeeTZARtg%SNu_xur3xl`CVIxYCZ%AWsh*)G zM1+AsuedU|Br%CWuehWLLTA8OdHE%&dg*zkdP$|3IVrlCDGXp<Mq+UWgI-EyUU6kE zgf1z9$dsiPf%ed%aN?o5pmOm=40=VWIf)?sP*y=s34>l{US<h{UP)?234>l*W?m+P zUQvEZVo4%{UTQ{sT2W$dYJ5gY5rW5{2ii58T9KSuP!eB~n3R*spa&YWDrV42DakJ? zX3)z@E&{tQF*h@rK`%YO1Wf3GZG)JcR9wuUmz<vq8pDDG8U-|{+X>1NuyhSx7m0LE z0;rFQY&@uojm$S<U;xdEgN8N0!m#`gqt}2%K|5n13|PK{nFA_5zWx8753?UuUchMZ zJ}#v41*RWj8Uq7pH8`lOfa!;oCotLpDgi2XKunl^SUv|GD*(!~AUTkku<{8;!^$I& z90;T5U(k+)|48nK<#!kj-e-a2ewaNl8nmAuRJOqQu<`~*gZu}w7v23JF%Slg?SRT6 zn0{D!1f#(lU7&FU@+ZvyAiW@540Rh!KdgL&(G?&?P&FVXOg}9AYN2v4{jl;DM#IWy zP<(^x5On>Zef*#>g%*?W^BQ3F49p&o-5?Ay2Smd#Xb~Q={uR)23q~t|Bw>o7ZUu3{ zb7Kq)4872O;V}DQ=bymn8_FQn2n^E)6Pb*rA66g2XwXJhm|75x?*5rj{jl&m0Ll;y z3=A+j044#_%mCS+4rPL=g=qQ((9SCXjcI}vLJ4&HLFZwBf(T|!2UI_dR)I-C^~2&1 z%4GnJDI)7%13I{Zfq?-=!^Xj2>OpE@Y!D4vaSsv$rA-jO1A49t>>L#I@J6>Ebe;sV z{uxmHGoa_h5Ym4W8h$YQS3vczfa*uL51oGxdM*u2Kdiomol^rlj|P-~LH>j3hq>!H zn*Ff-{}(_9XdtzJKzczKrXNOw#)Ogm4{IO2fa+fW(u9Oz`e5`|H2Y!oaR;bi#FE}% z`d~(a<^w@uF#peh?r)z#q<+viI>`O7auQ}POus5j07`?@f|xKGv=<*{H%J_Y_d@N5 zwU1%qFn!2;v?LCe>hJ_p3=tsDf!aR|3=GZ;4Dj<BVC531xWQ)sgFuM>k7yYV-F}EF z>?AV-BmBG~WEIQ|OboF46;+&>0akCKinGAm*{I^I46u3*Rh*3hR)3+2vopZzB~)<^ z23UQAD$WUSZ=;HHF~I5vRB_DndO-5b3_J|5@*Kp6VP*zi23R=_69duA415f*@)yL1 zVP*z?c)1D_1JTS30`T*lKztZxW)Nh6m4h%b5Y5aW1TWt}d>CeC5Qd*;1rr0&%nTw7 zuyP5+hhb(0Q3hCf0}}($%nV`-uyO>%hhb(0aRylV022ez%nTCnei4Wd!^{kl46ytR z69duA3{niRd<x>jFf#+FH3;KDXi$1(!V=D)e8<9odCm@~yl2D`4qD)H1~VK$?Hfic z;ou5B8;Y0V1GJw2aSa1QG}s)3`$6sl$z?&s7eLD&Xto9S??CP^05uXp3ZWR6`$4BQ zfZRU;S`R>#GfV`VFUYV0TED@>=Yz#D<8d|A{RPna9H#ytxDA3CZWnQg-^C#gt~Zz@ z;pe=-{L9P;I#Y|6VHNZ+5>WjNvRe%*UI6VVgXiuT7#J*};!i*U!oa`)o<n3{U<d|_ zW9FksusAQn257?^JO{?Wz)-`;z`)4#@jtBo0oe({lc4JFK<x$32{JG+?1PFgfX;}+ z%8!>&@dwbuQNVM93=9m?Od#_y(=%uum61t=VF6lvnKD7bAq9H)4R}tIfq?<E`4N;J zil7BJcy5t_fnh1Q68QKZJzs2un)3r%VT0;hkiRZK)vti2H;@<zUxlhyk^=D<7{GIq z3=9n4q2dgngTojY7-0TV2Ui}L_UbW%{QL1gdO2Xu4DuIdIcdwx06wn`Vk?4l#UUQa z3_Dp5p%%o-z#(3P!yM48Ff^%ww8QaY9O}0+V?Up3Cl2*zafsh%W?)c)+lFBL#-Uz- z1$+3Y;1IXKAsz=7$2><Y3oMRs2aHvYL%auv_);9=$8d;)(igVyX8<i&K|0!uA-Sle zxTG{KO)r@tKHfRRKPcYEGdKi%;+YEy&&wG}DkVQYJtsdYF()28KmlG+lboMhkds=H znxbcJY+=Lzl8DbtsQ~dz7$Bn_DW$o&l`wl;{anD-xCJ@-y2itJAZuW3ge~CF5U?%C zM~5XB6foo_f*0QCf!OJwW5vKMmvFtJ)HFS;N1Me*`I&-E@k~icaWyo9p4FCEk{X|v zpORYan#q8w#Kq7RbRHYBE@wm6%*+6{vVde)2GFT+=vvI7&T>r(Fm=kwPtNwtbY+Nl z_wjdf^ojR(a|?D2i4Spf@^Otvv&+!5IJG1`zo;N1F)zL-HL=Jwz%SS!-UQ?_S2VRo zCQw`AkyaALg9e{4+-qc5P?VpXT3j3tS=<A)0oCyMl;ZsOjKsW@oK%MRcyR1PV*)kA zOp!t?u_(D1Vq0*D8DyoEXJ&AsXQr#Gp)tyeD_1|u`1tg^()i?x3ebsndHF@Ti8=9^ zC8<S;CHY030hz&u@vg2x-kA*X@kt<;f{xxp^Rl54%*!xursk!fsWDG1O$UvFK|E)i zSOkhAR3BKvLnAXU8*5mkf=AWjT_XJ)eLbBSz@Z8`T@Q51UTSV}YKemb_;5PtfpqBS z1!6fM5Pb3<R2TF}JIKL+sA8Z41fe2Ovq5L$K~5J$6$hOsh$<3~a!etr9Q4RRG(ph8 zfoKN?LJlp2s)IWe`#FWEs&a9iPY6vhP}NXx<yMezBB6r=(y4&)QF-OLP{UvjL^&zZ J!2x!bAOPym-T(jq literal 0 HcmV?d00001 diff --git a/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexmaci64 b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexmaci64 new file mode 100755 index 0000000000000000000000000000000000000000..75cdb2c698cc1256beeb3c1826286d56f9dceda8 GIT binary patch literal 18052 zcmX^A>+L^w1_nlE1_lNW1_lN`1_lNNHU@@P21W)DkYr$BIKjlg5Fa1n8W92#LBj^v z`Roh~46INyK`P?oOHwOJKr9GGH!pw(BEDA|LWB4$5GI6VVQ7FbAY^=eaY<=XF@(c_ z={`OUh*TPcVgT_WK7t5vK+R)=a2X)x#b>1E6r>hmF|Pz_-WjNR5T6Ok22)V;z-+Mo z<ovv1gon}n>!1lysQ{rEKztPcLd9WJe0+Lxa(rq<a%w?Id`V&wx_Jpu^J<_PKztPU zLB(NIe0*|oNl|8AI#dANya!P8^q_`-_@FcbRn7nvhH)9<<4g0(GxJj7GxO5&!Tv{0 zFMK)>rBfgj1B?&lVo>q%sTuKUMTxno5Gi!`t$>;r0d*9J4`o6$ppl?-Y{3XhCoYa5 zj!2TQ^eX@rL*s+eG00pH2C0aTPpQmF0kI%BJ|4SyFyj~u7#LuBAv#dW`1tgEL@L2_ zp8(W7pu7pv3F0H02TH$!APFdrkIzYrFRsi@%Fl@}C@Eq<cb^2*ybh>(5Fgn*kbi}s zav%y6ZXV$H^zrobc6ITDl>rVQSq26Mm`9=f7LbD&7#Kk2Am=wwdXZpYU|8Y6z@Wsy zz!2ihz+k}4z`)JGzz_g6il2dj0fa&NK*l1w1F8g0op54cP=@mu7#JKF7(n3wDn3B= z=-+YXytvj{fBUiA%;>{i=T%@%W#wXE;NW5a6*r)u0I3DJ4Q7->Q9)u!hJH?FQf`G_ zN@Y%F5=4dw^{L1jh)|zcf@A?my#)gU1HwHZ9y^16X>k$AH6UDEl9H^gt!IX$4eSq) zN`$^vHU<Vp6n&tQqgc-nSsSK4=9LWS`hqKqOHy<7oRBrDLBp}bfq?;5c0odDlpGC# z(GVC7fzc2c$synv?9(f|{PX|+j=}#`=PNTX@XIs&SDmNKz%T&}LOnX4I)*!jIfgoh z1bZ~T`N6=z;L&;7qjT$vfB*k??tSv_|NrJ+%q2V?y|xcz85mwvOEWMugN2!EK!Ut- z3=A*Qq!}1Ix~*SCNHZ{e;@4s5Y?b-{|G!6fHOC8kX$FRF4v)@Wj{pDvyL1<@fJODC z85le}k2^O1gBbl3WV8a<lxl$&+%P>K{=xMKy!b5za$2jv|Ns9ToBy$Zo%9u??*&Mr z8*G6`cPmIwuLy@n?_P`l|Nk4F?4Ap9OXs~8v1|+s9?7>odRd}94nAP<U~Y8)Inuh- z<NyEvWkOI5-CIG*U+}Smya{45mO6WM&jl+z4sMyjdCuTI2$b89tne4X{O!>UHr?QV zYsUZo|6l8PbngWlaquyVN3Tq@;Wv-YTQAC`7#O<2LB$AG=D`TDr@IvtSRS2+Jq|ur zaOpK+@I3fT!PW4|%gvw`nLN1T#SfuAx%2DjJa_zm*2Q`)$hG`)54d!`Jor+c;~+SI zbU*=gT#|voGx?}TFN*`%=iR-ap!viv;KHck)A`Bc;4cM_&U+pQ|0$&L>#^>XVqo~> z$S=V29Yl&;_vm$D)F@;7#2<0+HOD9Zm<wtA8VA$(HBNluk3IE?Kf<@+6Tcvch&l9$ zUy#9L2PnUR+Ke8IhdmjOdHz4*0r8Nd<q`fKZ%~*+V#%ZPpvS?73edFCE8_t3D#+We zhW}p%F@n-7)VnPpV;P`Ec^v$q;L+U+N(-M{_;pkcesbg2(fNrx971Fm7+&m`U|?|k zf5ydnDLkkiB7!Orq-2f+1A|ZJ2W0d41y~M#;umaX0J*&L21x1?Ozg)$P%>|g07+T5 zCj9^Zzf>5SkX~$L1!W3|T(5_NM`tfsS+9qJNB2}v1bB3M2zYcl7<eXs_2@Mz^Emj7 z#iQHC$fL7@!=sbM$TRtmN3Tu2$HC_;KHV%<KAjF49-TZ^9*oyLIxRdHuX%J9D171< z0K4lGzo3Hzhy!vEh+~0>22aLg9*l=Q7{Q6k^Z#*Z9&xlh&fl{Z<bG&QxovpK@Y{=p zf}kuc2+n36$6FQt|NjpPc#q?)8VIHVl<Bw+)DQCM-3wCY)0=wbbu=_x^2<Zw>;GYo z?pBazy_#(d7)m&tzPT7kFqRs0#$NI0-U>3(!+I;oIR2Ksj0_AOt=~#GJLO<f&Bs_g zI>D+un~xb-gL87t8jtS1V2#>fYdsFWXW9o!Zibg$oDl;BIk@QZ=w$QgJpSSsSfm%^ zzSqM1@(iE&wHU#vwD|ysNAp1z59?l#S*6b*Wg1hd!G4HYU?X2tgSA5qbmIV-SRw`r zAOUayflca$IHUD{3G0heu==?Lp!CpsyTle$*m<=6FHwYrhezvwpKkLK))%f~poCV; z;oJJnr}KO1bC1sZFU&!@kGFz!gTuTVT!g)_0||HU4fy~6{|h4!yYUF9n2vUgag240 zbBvEY?9qAnwE~Ec1`%Sf89-$Rvb*kzGBEgbZv~m`+xo4<)U&&e1LQ7AuzT!FSYK=v zWnl2_tmE)){qNKHzx26J=l2)OMHv`8Fx|5bBn)-WA`rXt{);}41)cwWIv>6GD9pg{ z(&yj*|DXy6RKIu}cP{|@=eT<Ul&*l%84$W1oW~g$7>>730EvJqp>B2$OLq?bwj~S< z3~BuD52W$SyMUbI0P2b|cyxOxcyylfXnyq5qccFEdpk&9X8|Y=cxa^Y>zybw@o0Yd z!=uwf!Q;3C2PjTGy2C3xS}&Ebd2~7nyk_y}bkFeUbg%H}bPzc1AOlkIS~rbf?_im{ zNAoF0kf{ir{uLgrCri$FbUF&WJ_Qm6xzeNgfC1Qz?P>r2|BpQk@++b)ZGOYy(dm1k z^Rj2>KgZ*4pw2CW=W#bs=a9j()6b##pMXd5F~|KN@A&Zh{_yGi0WO?-Z5&<%iGZrd zBfgz#&A$a4_kxtT@cX{@O#bU}@HyKHW3aTFgm3GU5>}7a1Eta)&9xU8O1M0FZD)%z zF#Hd_08Yg({t1I3HI2i&m#5z2;6oN4#)qDregdAIZVIl9Z#|RGdK`Sn>cM!>v(rz* z@wgkf<L`Og4K%>O;L&`911X73fcY-i1MIWkKHbJ1-QEH|njbxOtzHEI9-SvVzQ6Kd zya-|(d|~f#@Fgqb!Ojy0U&u3_xcGzd1aq;0<uQK04?fBN;Qlpn0Hw7;VNmVP4DyYP zM|ZS?N9*kpPLE#O(;^HE|1U$lqYqM(>}}xD%i{_1jtApGk4|s6?;u<OAIlFVIv(BL z5+0V`3S~+j%?A};^TYFB^FalVZe}0N4<3w;4?rr)-#~l{&PZ@yf?7kJzkRxkC49Q8 z6+F6~H9%SpzF_S<c<`k><H3tRLH=rfz{q&QC;5Tr!6&Rfo!@;rpS@Tq1oDTo1cpbL zkUWwm#K7Rmc*IxpyXV38te(l|q4D6+>8;??nXCcwLZL^uu!m(ZN0~7s<wHo1ZZ1$t z<n`!wmH@d36arix&4&y?O3NRDggv^w6+A2tmKpHNJAm>D2S@{-M|ZG<kLD+E+<0(+ zf(I0Hj2;3W%?A`bES(j~uYh73DHOtjJ-XcmJX#NUbpGaVoyq_zL5exLng5^gNd5t; zl{}h{OT6e51pAi*=3fa|#0d!_BCZUisMmwTqtjc!BU#wv|0$2wOLfK`%?BjV>I;uf zMvv}b0gq&3aCO)lBKE<f`G9~&^FaX*%TwhKJdT4hDuYM2Glx&7fk(HqfCs2`=`9fL z7>AStLqMri;1j=qn}kQVn}A2>XV1=m9-Rj#fZ8puCDZia>7hFiq|8YIRAzv}Ptc>= zi^HS&Cu8Z8Py7N-;IgEg)uZ#~YY~smn;`9wh=95q)E0UziZJgJzknymtb^rGJ-QKQ zz4*j0;H8nqAI_T%u3uwMfZ9wqKJg2Hs_h6+`vSxeaF+PQ9|Ph!OMv|5({1ru5XC;< zZW|Tf?sN{%=0A)j8K4H5iwdaF^SuFTpM7%V7j#jP@a%R`0r?7KB-k_uBq2}}@0tJu zL+O^PO`gZyIY6Urp51j8FWCS6|L@aj;oH3rR2utqR&aQBwt<+hyFm81YkcC5alHYu z!d)Yc|NTi&Sy`%4rG{e2p}$~5JiGfqIz2i)V1{^fOTT7Ko8Z|k@6%n);o1D3u_W5F z+a2ai0jLMS?HEUX0e471IPwd+gMyP)1?1vGRR=tdJ0n~i`uG2Tk4_2C?s*W)B4C!i zp7Duaz**xHe}wCePyB+;DDFiuLI`FAx_dwI3xGWv;pPGIvYUrT^AQ0oImBPUqxB>> zd6j5*Bp>!T_<+fSnbE`Ak%Paj3sjG_o-E<_==KCDD`5kbzOTz1yV*dQ=!8eFO)qHB z?u82<sOfW{RK~O0PXm-)K_%Zseg=mB7a=8|DIWuaXSW~35O6KT;@RoO;o0dX;L+<L z;nB$flJ!hxgW2Ed2o~`;?g*+pAbOK8dvr5;^omUI=oAEJ7nIVjn+=@RK^cX=l?~Ju zJMIi>-7$D}*F}KZ3ijWe6<A6wAq}HHpoU}TL66RRp4w)}6~G5@1#ry+)P@FC0Pk2C zk98V0cq_2@bX&icf+i`?&Nk49H>@W0?XCwkJosCR85tOSyVV`Pxu_^Gmb?HZWKhXc zau1vW!6`wCkAZ=I_CfwBhgvR{s(3cns0c8C0}fn8v-xx;3cP0V?M%z??M$oi?MzVt zh0uSG=ATSOGdwIMN+)|BcTqv4A;I6E1gqfL-3NA1>q-6=a5KHrk>h0;xC6QC|NsB> zpcV+IJpIHk-~n>RiBdNH?m3|N?Kb%4rohPG(huS^*KsiLw}6_yj@|n}O%hkb1Fsdq zEf+{O_UZQH@NE9cRFdr3{FAXn(6ign@&5tKI*HP!F3fH+uX#J`B#yg*TO*GD4>&gd z09jsQ=4g4rk>BTn2fxn=$L?()gFU(}UQ4A-=q#6T>2#NI{ErZ?=kVyZw=9<^eeBWg z017?+RyR<4w9^|@6Y#g_{QLjkvAYdqNT;*J|3jU@GHB@<G#(HV?17RZ#0jN{SARg6 zk>nJi0WC<O6e;PW6e)dpcKd^h5=T(nuzu$zJBdsMC6NiJN#xxRj3g2Y4p2vqm%*q> zWP(R`0mvPGAa{V%$P7@jfToewZV(5OMnEMNQW}wk`hO6l4>_XJ2Y9>))DrBD&mf|K z56WB&Am1N{xAp&n+xi}@xA|LSm_hB`fBf<co}KqWwSwnycW}Ao*pRO9*T?e5Zm?c{ z*BhSQZWSJ_2TIs{IxE2Wye`AHv#!Fovt9sPM1bmR*6*AQ45h(Ufu7Cv91JBkC>50A zcTi-4YB`XTd|Us63T1bW*X-C+9>{sE>Hq)#f0@Gs>WApD&f)|W#qJ!Q-Tn&vZRbFK zYCTZG`RyPhf6HA^CTV{7&$H8C!M9uewa^6L&b$oY&b$iW&RhZC?i3Xc-_|GmE!p5= zPDjP@o2vp->2uF+R}N2W7Zrh02hZcKh_b2l8_2C1h_VUPcX7n&R**5Eu$c}rma$|4 zC_r3PKq-;+DF*`se_Qqc|NlLkYXuni!F^0<Sn{`o{|DKVA@DK?l&oRCY2^p|<|HG` zHzNG<4A}iMn*mg86ZFs0uONHK^iSyoP;CIpTn=C_34y&-f$&l;hzUs^{4F6cFZqK- zAvKtz<2QE&mQpK76=wMZTwwWh{sUEEoga{DF;E`{y$W+|@K;~~B_&V;6k2{dcDhSA z9(M;9*Do*r{r`Ugs4KJM*Z=>o9Xz^0Z68RpqT9f;+t0$IJ5a!*Qvg(bdp3g`Sq!Dz zp2xv8At-%=>XDLtAQMu4{QvLS4Q_ma+SOYiJrsC%A5m|^b7AM@&eAI$&Gi-xo}K?n zb@qds6`xr_N!8B*R?T`GhjgBggV%0&bb3g5cDh-3cE&|GHrE9(l&W->UU6voRwCos z9S84#C9xyg@#{cCTCKNBR3QpWWqrHL96$pKpq^(XJ7{$68l=I}!pgt^QB=YO^3jps z|Np<{0I^=Z|Nr0cfX8lFpYb(6q-KPWpk86?ff5OrQ$3n1I2cMeJv)!TP6gG!zTI^R zzODaDHGR9mIX?lEh&?;q3OtXyfd>;jkGp{f7kpd)^SAtF0tJA(fRE*mQZCQqZs1`E z!vmh3ZV9jw22w_>04qPk#J~VbB8<M>>LrdZnV1+DRxt9nGJ}M>-6b3?T_x&O9lPCS z94%dC>ZHNx(X-Q)!x1Fe=_&z|`hUQ&(^UqM)sHbkMpK~~F@_aXH67rea<Jt@=}q5m z9Tm^!S`G%!ZdZlU$sWxQ{&{w~D&R2~oPA5BLB0lMbvB>Q3~(`sC}p6{HU1V*)&P|R zEc~sYtN}^%E?=RE{v9;Ye}4@gO8npf?yf^xYmN=R3M?-{qjsL1M_#fpf(kz9cplhk z;BE`3AJNAG@(rXd#@~7tS_(jY<g38n+72onn?Xg0FI9Y1k^}ZpX)44wWnYN(&1++D z0j&d?Jm~h9aJ2N5s8?|8_Lp(A^p&X-BM^(NtRQ{=!EGkP1K=V&?C<~opjrpo?Qi;s zoMZW0O#c4=4@$2--3c79?DyIy4O~9V_(o~@;PM68w$d%Woy$P&3E$QyKAjIryL~!; zyxjKZ|9@x(<Z$eE3s_-a>fqTOX94P-Yr}>mr6wT8{}@5zk=3APIC%UIG#<H}k%0j= zWCbz<l(yXhK-CJkJA^zu1REWK%p*0wk?`n_=kTzO7x3r|J>k(=dcm{vUx_2ga?Wn& zKb_2;-EkZqz0NG4UW`X?@E@P<(hH#BiRMRtJUZ`tXkPc|jy=&W?bG?~g*PJugJY-M z%UcW#44&F@9^fH(CXb!}A;8|l@<a&_sF)Pq4Jw*n3w!k1$}uxAyto5Od>4OvG(Z05 zVR^Wi$744{$AA0R!X7(dY!A!B{O!FUrM)(DLB;TkJs^uZk9&5T`E;jV@YMX`0~y#n z1R4VAeBjai;E(6QKTJN1KT0?}x_wV{b9Qog_tq?D@$Jsz@aQhR;Hml5r}Nc|Ng#6$ zzGd?4wgdH_SYL!QFfe%cx~MXG_Oi75cG`hF@7c>TAKZ;=J|Y2XhC@emUw{_+O!%+L z!ok4sUzH6^a)3!LFv$ZZ`M{(Am;}{_|5Zi6Y%wq?0Vbuuqzssp1Ct71QVC3|a4<0V z^vcd<XJB}F;qm|fKE1LtLF^+C_EZpi2ZTKl#9jko_k!39AnZ;MdkTcz3SxIa*o`1| z4TN0_Vi!QzpnBqE3WQw>5|4ne3qfoj2s;<Vc7U)mL2MHUI~By%fUpxmY#9hU7Q_~K z44JC{MIsso<qnWg!~fae8yJ{NL2NEiYJo98Ynv{B7D$5X@`E5z5C*m9Kr4YxgJhmE zF)&myFfbegahO>c7~pyiu`n=z7L8p8$(XS+Fo?m7_F;o8sk#c10h#ERU!stjmtUHm zp^%%Jn_pC^2ohyrC`v6YPK_@qN=+?RNKVYl%P&z#N>wN=PEAqB%u^`INL5J7EGjNh zNKVYjQ7Fk*$W5$B%gHa-(1hww&PdElPfdx>$pkAxF`zuNBqP7HL;=YTMTWA}<dXa% zE30_l_{_ZG)S{C3#L@}|hLDV+)KrDUqV&?-)Vz{ng~Xy%g`(8L(#)dN6a|euzx<Tc zVx5rGqTI~9#GLTd%=C<sVjb6%^i(KYlOZG@<hr~{h5XWzg3=O%S$ZI!hGe8Fz_cqs z^eE&emSk2a6ldg@=A<ZOB$lNr7%L>_=alB=73(NyfUMLh&dkf!Db6n~N>0@z!5kAL zbKHvZb9F-Ub(|9m5|cAaDs_TV%i>`ih@lJ&44HW;sTClnE99psB;}Xpr4)lA0On+< zi+vMIiZUxe4gr~}fI|}~SYUpH=?KXuMKvfH$uK}xEP^JRL3|ws2JmW6L2ed7MFu9p zixI4XpjC*V6}mnQ3=A%yRl^{~pwfzgfuVqhfgwQ?Vg?TbH-iNO3n(Jk6c`v-8`v0F za@ZM|yf_$`4LBJX%(xgB-*GW8MDQ{&%JDHU1PC!OB{DLD7V?4Y2?S|`V34aoEXJw; zpCAUt3IX=nq09^n0s=hj91|EB7;G3A7(i?M7TkxFvoX+hh#d8QAOfIP1p}Wz8xsp? zQk|EL1=QbbX76F^W$k0>XEtGAU;r)aJn$5v{0+!52FRLdFcV!rSP)__XjrD?8ARV* zs5p2e3&cj($HL3Nz<{t9<Ze*CA@m%g|0zfhK6_a>LHfaN1}O#U2NfYZp!!)sQHEDP zXk}(ITMyI(kp37(28Oy95c|QC>3H>ng0GpS4@5CQ{m;n2@C~Xz6EwVwS3hXQC#wB- z7#SE8UPJ6Z3i1?Q{UD1`?U!L<U|0s#FG$dSEdFO=V5oQtv0t7H{d1TY7$V+7^tXbX zhu8gB;-86uVcsW*{z>4>f+yW!Nk5>`j1S)+`a=oEKdSu<0te#%*E2A1Kxt_3Gl15& zgV+twoZSPZXF%yCP<jKD-UFpiK<O(``T>-F1EqgJX%^`EasenU1En>fv<Z|3xnUHK zhQMeDjE2By2#kinXb6mkz-S1JhQMeDjE2By2#kinXb4ai0-$}P%nS^m(K*mQWj0Xv z9K;7L9|qB&<-;Hvw7nNZgO=)oXwcGK5FG&-Eo1{Nv<C4(bD|&`H1i3fL35)Z8Z;md zqCtxwKs0DI5QqjXR{_x%Bt&=~@VsHX$M}q~jcEttE0$x7lUNq8w6RYD4S0b^Amtbt z7(k<rARmMF`GWW<py6@`1_sdBs2pPh!v)YFIRgWO9|Hpeh`&b|q?&;t6Uv_@0^x%u zU_kN>q7eQDsC=6kgntvtw-bl(|3LXGK=ULF3=9H{5cBh-Ao7k-egkx%vINRsAP13O z0OhNJHa{~kFdTvMcYrq3GB7X*Gchpmu|h^pK{`~K7#Kk7AwZ+4Aifb3Xm>97NCgm= z!56}ZjutU6q@(fM!2AZdd><Nr8XA8t8h;6x4;i`zS-1<rhmNW;FdRVRA4B7xL*rjZ z<KKtyp`&dK3~wQP=;###13NR=zi>WijEfO01ab`ngDONG?jBGu%QHasB+D{D_JV_w z0x0@HQ4We`P}G7(eL-m+6pf&$1EnKSl!2x-Kv7K%a}*d-Dsxi62UNr-W#*-1=B2|9 zmVg`|5T90>mmHsypPyY?zz`oF<!9<%lvt4AnUa#?YG?$VYKhOx%Pes%Lzj0pa?Q*P za4QQ)c14#mG%HF?%*n}5PAo}{&jZbzxMpIgF-207Sd?50wZY8ADJMTU+Y{L-<{&c) zN{U=vjr~%~p$d_lYMxk{4w|)b4e$##h&N6wf;kCdR3=0lrV}g+it>|Fi;Lqyvn=uX zMFkm&c@Q&9KsI9d%g7L}0X!o^j207EC?ZWl;RsQ4s8ZLY08^OD!NC?>Vg{MD^2`iQ z^vra1H8e(<y>j)ljE_&xD~(UCs4z5-kI&05%1z9P&n!tTN-W7Q@(jofHjH<54f4)p zaCde_I=mw%)sP`RH?_j8G%vX%Ge3_ZK0Xw3dPY!cacOR9e11`UP-;m=QGPkdA;Bde z)4|;Mq{NhXa2$gZMsSIdL3}(!CMPv7y(A+(wWuh+D8M@zq%kEi0IVz}2SVp1=7NOt z%0X0odPPNiL26NPeqLfuW=UmynE^w5YDH!VLws6MYAQo~F4!KBLvwTU%RmCDxy7j^ zAb|?kq9Who^pFY=ufn+~HL)bsCBHN&Clxa96mM=|z!0BX;htLJl9`*DSDcxjSBxYW NP{a@qO@|B&3;-V{=2!p# literal 0 HcmV?d00001 diff --git a/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexw64 b/Chapter07/Chapter07_mcode/GraphCut/GraphCutMex.mexw64 new file mode 100755 index 0000000000000000000000000000000000000000..d5b12fb8fa67766a0392f5226a075f06885a5b75 GIT binary patch literal 16384 zcmeZ`n!v!!z`(%5z`*eTKLf)K1_*F~P<Y7?3=9l>^1JvLws4+R+`;H`Rxu<avsj^^ zC_lX@F;^iuF)uH_L?J0vp{O)ZAu~_G#XndfH$Np+kC&HA1#0r0r~ekrXqoUChGQJ# zK2K_y@HxgY{_`v_9sPMu%Y@G<pEbe!W=WpUbHVB-e%{?O;q!i{_|LP!bkygWU^)g& z2W2K_fXswi1E&I985mO9*cdXDO+#RU9~hJrxS6;a7(kAafQWYJGB7X*fY}X<APU50 zVqjncu^A9F!vsj+DZOW4U}8vMglU56g&Tt4LsV?=VPFVkU|{g@U|^^KnOVTVumxmM z0t15y6G#Ax8zAbTJe-t)3j;#|vaACGg9FGsWHBVRUP)?22?GOzwKCL+AiskW0g@U8 z28IN^qLjpvL<R;%UIqpRP@>>tU|`TdQUnop(1R#rhyY0iFn|Jsfx!&KL0|_4h68#9 z5QShrCP2cD!2+QM!~!{8uc){vnE@ImAdi8<#0Dgbzz}t*Ir+(usAFh=+G~eU17bNa zFgSteQH&V^4#6JHZxlSbZ+diddbFM_5qhykih;rSfJY~@N2fD~N2jxZNAnSd=)-3q zz%k6D`Hg|caTgVZfB*mg_vnmKaq#FYQ3>$qtWk;Z=>Fu<{nPk1NZ7-p^-_s}$BzFH zVDDjlp<K$NmuD)-5icC17#KXdZ+LXNsCam^9w=r1FS=2mfnkDAw>yVNx44h?A<u(9 znEZNKT0zP^x=;TXT?vx+=;rt7F6Z#k{^)t|9}~#MtS^2^GB9}bvTXL~bQkb~&>S9} z@*dqFw{?Tu)(vu7^AQ8bIG8^mzBzNoqx*<QXMjMb2S;ddSh!<|W2j@8W4K3e#tM(d zHxi5t3?9c_AAo$yuUW8yU$bI{XZO$Un~u#tm`iv(dU<ON7#LpUNir}r|7NNK32Xrg z#DW~u&Hutvl7Zoq3%`zmN4G!63saD2vO#wQ3z(}T$-v;#{iFFe#H7U_jVd4|-TneE zI6;!h8W6PtFBrh$KR}^n(fpeQZ1o|KhR+fV3~s$FdalMN9lNi2^!lhcyyz5VV3^>U ze9WWMwZfxU#=+y@V-}BY9~BReUKVZ-YZsM(G9kF)J>sBv@KNzFzU0{akFiwCqf^kM z8x$#>0UBUNbArNO7~{7`vx|y{!G9N(g4aSm-G^UHlwe?RVJuN8@L-GqMdUaBz8MS* z3_iUoUqSKa-D_jx)BVr$`+ffKZ~XEu3=E%~_;qCd|36*9>S%qh1{7jJAb)5}Ffe>_ z<QEXJ`NS{4?~#1iqnF3QgZU)LVuw%s0xl;s()jg{m9c%|k2vv~<r9C*g*1Nsi)F0* z8aH1v^K0Dt#4q65@rgg;@F#vj5F_T)Cw@Vn6B-`P$2mOypDS(m=&n(*@a_KV+WMBi z?+XJ1gI6z0F363Z%tt+XO^*9^KlS{6$CvquXZKN$@8E=K;hX*T#V1F8of=b!yPy3( z{r^13<^RAgKcUaS@S;$hfuZ#_ND$(k#~{HNaR!D@`~qem=UaT@7qB@7VkLkAKckc% z5=};;ppbV_N$}~u;M41&;L&`*0TjUr9?8djdRZJi4?bY=VRlh*0EdQ%gU7)~EFR1u zDh|G#UIM<IfdZc044$1%9KM|b9^FiyosJwnom`%sLcX1hpB(uGt0X{Vp@c{CF%HlF z=lNTW|NsB*YJ9->62y9EHU<WMjSt5EK^a#8l4UiFFL`u^s2I2!Kk?}Ht?+0)$=}z* z$iU#zTh8ce{NHmgBrE^7PfO#MZvh1!1IS+vp3ODv91JCVj@>To3e67~`CHXMx;je> zJbGOVJgg7%H%|lAbzsg9{`O8Vr?kLB`@07y7hB&d&GqO$<e`1o<KPb_kLDkY#jzgR z_dO2&W&*{A_F<3^f74Ttw*QAcw6A*{{Ke$aeahqFOR)Sw59>q4EXJ2!GkY*!^sv6{ z!S8e#<V%GYzQ6zf_vj8$F#r_?8lW_w0P)%FlG<*T?VZf8%RH=am1d&Z9*JRlQKbiI zHeddQ(`IIuPGgVe!z`duQp558G0*=;Je!Xj_*g%zReZVn-~ay}y}W0185sVHF4JOQ zIOEaFyARBp0p@j|_Rv1>aquUTtMLKH-T=lIFGWBJsWii*lg*?1#EUCn?sX6A-zCy% z6Vmwem_GAAJMf8LkinzZfy1Nu5Q~TP<<jTPwFyk6@~*}Qc7a+6-KSnG0_p2?;{cg4 zK?Ib31r0n}|Cexjbf14QU4(&Q0w_owUi=aSC7dlP42&Rpj|!+2=+SrtR8~61JH|Q2 zI>yAuMaM!5T96+(!r^s{r+`QDFHoti>e0>U)625M$J&vDzpV=tv#lpf_&vHE1w2|0 zl(2y+p4VmF%>Pe#B%ko;<>^&mU|@K`1#)Be2~dS0=h5w~0aD88(aWm`3ZDx*K|c8} zx=<4oXx+{*gFFOYvv_oRb9i)m3wZQ8N_ccKdUOX1cqDUsbTfJM@^pA~dU1Gw)E{>M zwelEzK<w@q6$8(129I8m9iE-wIxrsOB_#heR|+tAbl>D}1tsThMvv}F4v1?HK4kG| zJ}B`bMwo#CT#2?GC<QrJ0_NO$9R`N~Fz0^K06DiefWsr1+2h~?7LU$g0gvVbA3T^3 zd01a4R`lp~67c9`0z3Bt*qLDSK?+WISYIr<?{OT|He>236M<U=`gE&w~5_s=dK3 zhtzmr_k${I0Z_h|09WS@{F)^y2A}u^iX?nM6?l+<NB2*D%^DR2Pyls%aCkKTWGqnu zm*1cW`ou5bAd$we|DlYTU*m%Z$m<X_r4K<;ppqIQ<=M^W(e20q%Dp9uX<+>z1;~0i zU<yi~e&QE!)JWqGXRXj?V1OIXuW^h&{QoEZ2q%e8{DMxPQcb{5!n4~=;1hq$IZzqr zC*jd;;L#n&;nDn?u_PH30Q?%qKJiC5gLHrh0SAds{4vKq@e4XgfTTb13wTJR@#~$h z(kNB)IPS&)3P6u;e~TBgpxoDK;L+U=YD;=_dT@Ajwu717#xN^N=YHZBaMJ*pn8vSn z{1d;Rn+AqKSN?(;M_3H<>^6fLR1%Xm0puKrYe8-l1Q7z>AlIG)x$hIdptl5u*%i3V z_CRvvGKi^e8o1qPfXkdu{1N3Ipb|5{qq{^Uf?uOXMd1^_AV{o6C4pa~M8&|P`G^Ck zHRS+lHGy(CsM_%CzTv{Z?UzUQPsi>rK9-3VrFxFNJ}L??0t6TsK*?XiqdVNer!#@W zr!zsoqccF_IJnjRn%NW75OcTi=!}oBEDtDEc5V4zBInWV?f}vS%Hx+bkTT6)RZymB zy<GxQSE`Jx?iI4SsUUSAqqsn&m_G{x!)p!@YbrA+Yv!mV_;lu|WcYNJs1$%QfdYT{ zH_yGRSAl^Rzq|v3N9(r|ZlBIH0UygWj#8fq9^K^$9<AR>6+ODkGd#N8D||ZBGdw!W z6FfTI3w%1$D?D1i@wcpH0u{+_0v^`KO1XSG(?K!`pi;;^0cKXo)^BbKj6U6F{H>sR z%@?HY5F?}&<OWKd$Ew))TXrxoFnD&$`B)}%@V8boFfcSHEBpf`evjU41yI@2?IzLb zEYa;I^Z!7nvy6{rvOuXc$XTFJWAo^Y0GEC786KVS6&{_@0v_G|5uoBX`I{>POKB{) z?0*pfs;ac*eYy{TYSZtpJr4e22Gy<yU)p;xUv%sgX!d1bd2R5`S%HziRgRH?q1#Qu z(pjQjuG>w<(pjcX4D6szXAV@y_;#0ZfV|bJ`Tzg_m(Lg&7(l$1y9^8ri1veLx1Rtg z<HLQt2I1q<>Sk{RMzG<O`8FEt+t6L0`r$>^&;S2D2>Mn7<l7QMs86M!J{1DhBRG6% z^#A|=*UUcMAOE+5#st8AoeTENHIMG|9-U<ZFX}-7-8~;%p0wWPZ<!CufW6Kf9^KI# zKAj3a-Qc3exAj0t7_5Ee+3mmqDw*KLkfVTSw}59SgHNZUfRE*JkXfbdKAp=u7#JA7 zIVmuDbe}I_fAJKQ(pg{Vv4i>;FjGr6d33@oFKvIR`1k*R&+acTCI9~a56aUDKHViM z1wP$5Dj6Q#IVuUDY$Mw39-wN^z|c_3{{KYti39%+v>qt6_vkLS0A&m#P^KtRNlQb- zr$;xq1>w;x@gkWAl*i&YJX&v;zVPTi_acagfx)LU-QvG{z$bnIcMlIxk1QVCTmV(A z;0ns4ThP_`q)(@3flsGrf={Psg-@qvhEHcF$fGWuouHuguny$lZw3kQho5XtQuz6a zU%*MjrQ1m%jbHzuPp6kacOyt$r;|h(n`85Xf3LaH`1MbeIXgC=V07tpQt&wL!~seR z9^H-=po}8w(HS7%(HU6b(HSIg+(`x`_gd1UGcW^Pki`eoR%kwC0CGc7^XJ13!H%G^ z&*DGWMpsaCI!7hIqdOB^s2F&3-}LFMEb!^9Oz`Qftnlfq%<$>l35qM9&Yd71_;i2u zXno7yvJq5bbe5<@fU4yLpUxDO40wUG2T>p~c=Wmocy#9q_;gx$TEFw;cX{vFoCfNk zd4rnkC2T&Og#xcxd^(FWd^(FOKwT;~4WCYfG=9BfRpF(fzLs?y{4FX>3=E#|YR}u^ zh1&Q3|9wE!UN^{TpaMd~qtlVYqq7;r_vx(Qc=-X;7w!ayI<oLhMsQtyjK6I+D1%$p zag=a7f=V|2mL-g!vdawQzJK6SF$82px7{~)1t$L1Do~TJ+fBg3(w(E!!sEC*v=a8V zcroJ}sP44z=!Uz74dj~F>>e2Q`&ia<lo)+;S70pB0oBT&rp2)`MIXz0a5?A<j$Ss8 z&I)ijSe@a~SzY1LSufzxUB&@&Q)|xu|Nmd2JH7QKINzE%Hv1_sf|Rw)WME(bw_iLU z?k}zL!0z@tUrBSj8OZIW1|ZX5O(T-s@6qiIF4WCH0nbvZzYCOpUo?FIW$$wy+T7rR z{fNi+mly>*$b&4PlmODv-3-b!j>p}=>EE%_O#+%YTiX#OL<*>Is8LA(Hw=8b!RY~1 zZWw@yi~vy46%8s5;Qd2L`Q>c^i@I*q@@qFJ50tX|bnXTvG)Vb%x`h3ODgy(<Yt|Qe zOrY`$Y>98{fznML%||R8<Kcr*XU?4Q=*&@(0Oc74NASp$PdAtGTX4<+H)Vf1?)v}# ze?5QrH<#|azO5(uTTXzhAGrHKIm8iO`OE{yJL5M81xEgsiJ%l}nZd!|($2uZ;L;5) z*1ONTSUO0QsKd%<kM01D{|7yqe=+j6^nx@r|6(eM^5|v)DJy-*>><<m6Xb$gH%IGJ zwQRoKH(k0fdUgv!OtW;7DA7d;FNk=jm(2f@olX)S-A)`H-2xCnODBoa2fp2&93X$R zI{*Lw-?uZ4!?*Qzi8@39D9RjrofsWE-DF-0f_vZ&62~3D(fP6zG+GC${X0uQ12rd0 zY$2(|43b(xK0#9pzfbo8j~%cQ2X~_BHUnpjr650p1_AlyK~c%z(JcqAIqN*Sy*WSy zNhvg2>K&_MDpB;X%m!yncUZ>6XrS~u{s;9J%4~fs4MF9CGqharw|LR;36!@aJi7Zq zi3Mahyj%cRaUPxRV2REQ4wU>2YRNi)608R}k%C46BES_MxB@~Q5~MUed3HO3+B2a3 zksi1OHybHQwf2ER3!J1tqbI0I3f8|Gr0GeLwDiRCfxPsz78F95>50F!2_ykcQ<X@? zB^A<@#s|{VRSGyeK?j)eq$_BDJ4ZzT)ajN82ahf{zfsr$8fkDf{_og*!}@|p_fMZr zPX~UV3oh0NTskWa_<e48Bp>qWthDgyHF5Cl{^)t|8H;E4x&NX*k_-&JJklP?H$A#v zcyzLR9DK|Q9vwFZ4H0*`OL(-tEdlqLtU)P0z_XV}&!f{_!GrlQBpY<*bAWPzhexNo zhDWErfJgUnm(KbCpU!$vdrSc|gcE-l)VKrnxvohwF#Nv;9ar&`0QI>#OH>3vEgl7r z?n(=f?vJ1$B6pAGdJcvX1CVD8z{MH7eDwg2=br$z2ZO;)oOoD!Fa|Uzcm9PfNc6ZH zcnZj)oA1BqEAWUks8e5}VgVk#@c6HKRh)qV(mw2rQQ-h(QfPQWMwTiiz%v-08Xn!* z3XtJ3RuAp>-3L7m{$-u1edyv3<`drzGL|s?Kj@MC5!4(%_>9$~``Ule9<XgEJUX30 zqch-P8gS@=e1Ag<smJLJa)C$kQICTUSUs2zdUQG~c=WRHd2~8!cvzn=N$}|A_OSNm zD3$?@hqN9j;quTv?!oK~9&KtqXaJHcdg9UTEa748tWa!-s6RZw#UwY#)y&=!;30=@ zZwZg)17MFZdVBzBD7p@cJ#Zy!0g6`#aKIb})qmjdKtzhAPq(9hOZO=s?c*NZH@iIz z9GOod%5%@}FMPT|Q!C7#-3L4lzOeUUKH%7W031Xv%m-Yo54!OC9Q5cu@0onSqu0d2 z<KQzEa5$xifqJwTK+{6d(JXi{ZIna`CKWMIB1v{O01x2zvWR(rMyo)xNX(}|ohJ|L z>m~B=U=W=E3O2Y`AAr632I^H<^dJ7O+9t-p@L#n8OlE<}5HNX96eNBCOilumQ^4dj zQBZP$%x&!i4OIRY-3=aYKjG1R<G-jrh~3HRp$h7tuv&ts*4w45z0QoUl|7o38D8vS z1I;M$zGh=!c(IEW%vNS+V0f{P6+AZvoe%@h9Yt_>bY=*6bXG|CbS5Z3r;nQ7C>a0u z?0(_X{jv2Df2%2|&&;}3g@K{<0DqqXBLl<#6E5BVJ6UI`FfgQgbhD@$fJ-+Ij@PUn z)(1=YjSqbA;NQsGufo7k`tCKG<9=}ajI~{bfuZ>aQwjf0dC+L63o`@5%X#2NHE)dy zXuLQMJU;ruqxArP3#cvQ(ajqTmT(mTHKgA9bU*Ov{@~gCj>E_LYf+S^rCrG<<CFhI z<wO`5JT2Kuo_d0Epo9nmgYijlNU?$>I@$hT_UPu71e^3*7?corSyVvOjsK#rg&7!p ztY4S1`dB|GS?bgM+N1mR>**fNf59WiXa9>zfld1_x(lSSll8PRIMB~1Gca`8Zc}Dp zaBcn1-+BY&|90DX$_xxI&oF=j`Q<6lXvZ%`{?<JVNP*OS-=iBnm@<^X!Nj^75=x26 zpl~kXfdrB~69dC*JC9z`Dp1(;@<xj=FdS#KRAvAr3x*eNtRTgu{T{udmJkJ6A`A>K zG+9Bagq0!l*sZrqxjed$d$hCuRbpW9u(TCZW?(2$^yuZ~Rc2s#A;APPqeRfB`{*vv zNFunP{;b5n@Sk<35(C2ykh4H(+5uAVcGsu`_;jbJM1b-zsJwNIJ$%LkQm(=Bv1j*> z?pq$+ADe$Mma@L^X9mUI@$SnW-RD3C8GwoaNZI%jw0Hy54h9XMfzo0DI4wdN!=BwA zx_|j}e*_0iXT*z6Q9&gJhSp1^LYAWPN(>C8T*kLMBVN2b&A`CW>-4|*2V=RpM>i`6 zD6h93DB<{jq5DRsi;4nhtYSH6TCe+9_jS+ib1#_w{{Ih^I>6u3gDe%z%)kJVYW~ky zYSMb3)WE0vnDHf$HhZu?pqzbwKocaavq5<gR9?bV{lC!qzl62<KXa+M56r$vfBygf zena}C^6BoIojxiGB{Ck}tcg&a0x)GT=e}k$zTFw4qHvs5PLY8DG-uTft^B${ZC_CM zLtB3TRi_AoN|1Mgpz;D#!T2OQ+JFWix-a=KA9Up3?jYjIzdeG@qxpb@C-VW<*4rg= z-2ps4)|ZO7x;;3UPx)A1EaGzn^#BfNxLThsdSrayGkh=u6fdB7h1QohJ$h|>E`WS} z{)GtxIL#=5CYnI;dj3Tds729wpoG8MLFM2J`R^B+A27aVGClx_@a7{5;8A>!?t?ph zqZk;{EICS;jSn35NcOPtU_N1dpwSYvw56117dI0~Hc*0p3dkjpGSH*hHVQO>Qfjt~ z0~FU@68}YS34n4iq(GDdu{&AK6~K8@14KazM2XjIpaRkO04S}3`dv`-K<fo~ITS!; z&MSV9CKpg${o{WrH@Gp<C&0k)x^5@Pt`~n8K~qQFytm{*s`tQEck><qi>!i*IDmCd z0r9#|`1FD%K6^t{EMA-d7jQQ~GuRgYMN2^Povgj`AjkWtIQ$n)f`}B$GcZ6boCwx% z7~*1&<~JHI+(AM6{f6{`?>GLR@Mu1w0m*L)p&s3DAPorCKXMEV9?hu|3?93U7TwpK zsm;LPVSTu?+N0N3;=kxKK9F^;5+0qd0{=xH@G&rSAMh{)v5$ac4nB~8@V0??y}k+% z-T{z0R|SuQ4-`CleFZ#>T@^fxT?IUh4|?nbrAp9xD-Y#^9y^T@HhS>89Q5d9%>_B9 z8#JZ<{enj)YZ{2#d_*DkFvvfkLP5Zz**0E?fq|iv?Z2owAFPyt`wcqh)+>8h2vmL9 z?h#^Oc+Kb0%?ujvV)}o<qxmq0hxPrUw@AqgTn0hXYOm~Uh`uRA==-ml#s|tp8DQV@ zo|I)^cnMlY^<Om(ES>-+lfYySm@MOi<u}lLvp^^$5?*<9`+B_4`SJh1@g?JvFEoGr z|3AT_w+poD=w%CNT&1@Ql!#w?-2MOG0ioxjN2jlbM>k_PN2jmK3(Y?uUEr`f16OpM zRYVq25`a=gEK&xAEAi-NG(OpA!OFmJfPc!t?t@T;aCu1DGX0M*onbhD3Ll6s(x9{t zls1FXGEn-DH$)u|ln>!CFeG_FST#@@!ed}K2j$;@((j-&3)CDwC@lu1_d(?~pnL}? zy#^|t0OjXF=_)AQ0i_o}=@n3V6O=vxrO!d>8&En3O3OiMn7iGe_QCiHP`((HW`NN! z1vuzCS`6TIgAeo}G)$Qi>tsn#8{vfuOoV}fm4SibhYJJ41E>(Fy#P~Z;|gi0__#vU zmAEo6fHDC{DM(#NMp0_2LSj*RX>Mv>NwGp=QK~{wYGG++QEG~UMxI}ON@}rANNQ1T zW?o`Wcxq;PMoF=bYf5@5l&uLeh=GA2BtKsvH!-hLA-}YwptJ;Gq8?NMlnTj6Re%|% z05L-$H?btMLZLV#zceRBAtSLYRl!&xIX|Z~H?LSnK?CFjo#M>AY@Oo#(xT*4O%xM} zGtUIcJh!6!T%C}79p}V?#N^D9N}Zt8vUnH=Zg^&1N@@kjGYa`>3Q75;c`3!9C8#i` zLtXBhSW=W(0df$?R0SNGK*0s`D@;d7J}Igh7|K$UOY)0sLhKYu^79pP^7GOm0a}z= zTAUhRQk0rnjI?03v^X^d<lK^sRE4z6qT&*T<iwmDg%VI4Rix$QmuqM$Li8tRB<7{3 zro`uDf)&A5+=2}#&n(HvFD+3(vIDA%fx$1oL?JaVzcf8VAvZNQzo-%}12z{tmdC&V zo|^>=zzEP%Oc<M*G)DkxJ|Wg@K|?+U!wP-|-3$Q+MMh-?u@52-gf0ji;49!w;4$D1 zIH$$HP~ZqD`xh8NOkD$_LG8Z(|Nq|*goufXK=PTK7=*6_>R<l<|6fv$gF&u@n?W*# zi$TPMgMlrG=>f>h8BjApszG%*h=yTMz69X}SIClCn0SK`19+*cfdIrk6F_4?3=9k# zpnQ-#NNukX#4VCsOblE`Y+yft>|^9%W8ewm_`r66rGPns$v{zrgFz&NiGlwK9|PAF zUIx|>9tP05%HL893~O8&7(h{q%xC0aX5dI-`M_krz@W&-!oc^0je#|Ui@`yjfgw<y zf#D5E56sLXybN3=JPfQU++evLc?O0bAbAFeT4X-5enveW2E7U<2B{2Y2GJ582GJC5 z22m3(2GJtH2aHl23{n|P3<4$G3<4os3<5zs4-{FM7+6y{83a@r7#!6Y7;@Yg7{FpC z9AN$wH3o(TkQf6*9+{7!Cxw&YfD!}4M`cuVKztV!28I(*bI|$7dKe|R7$j1-7$ibC z86-?N7$l5DJ_tP!C}8AbVc<$(2gOGL#C>cG1{w?uz8VY+8t#y=M&>i}voY|euru(R zurcr(aThRx!Zn41K_rBoLBxcOLBvR)fKir*LG}nYgKP;GgKP*VgRBV$gRGI{2k{4@ z7Z}-?7}$)M3?xN47(`3B7(`Pz8AMGu7({~vAH*0zx}~3W85pkVF)-{fg3xbu85r*9 zF)-W(rEdlX2B;YUoFKQGuz}qU(hEwPx%vzY3-m!fZ3rz0wGSqrZot5B)_{RQ!5EYx zLF>O+7<f%M8F-D@0vK5s8CaYc7{FlyHTQuhL>-8<@nT?5aAaV(fFu^?#lRro2pS6q zses^lUJMKl4h#$uNa~Ed85kl!3&lMkYCz-^P@4nP3IuV`@M!x5H9$yhzsRuWR=B5@ zxMb$0<`rk==M~4B8yGM+KtuzI7+$mFRyY@>CYGeS<d-Jpq(bUHup)+>%%t23y_B3B z24jZY)C$+4BH!Zl5YV;)kW4B>is3qweSAbxVoH2sPELNZLvcxpgF~caL_`$BEs#u3 zYF>ItMto{fQGO9zGQicz#UYA;!8bU}Il$1+0Bi|EA%nfWkwc_ofU9Ffyf;X`n8DuO z7-Skq$RUd1HEV8aZf<^ADuWkeymLiGNJde9xobsoYC%b6ejWou8e=?2*sUlrH`OCC zFC{0n$e4lQEl5{!Y6$~F2Wx6YW(fnsV%D^x)KmtBd#t%&SA&Lfic%p=keh<wZb}51 zSWpm=Tu|bcnNyNl#PEPIJ~6jAJsx6`J99k5i1_595{6?;@tJv<B_*jvxeUh<?D$lM z1duuL#RaL!nQ58H@fi@uF))}h$LHiHXEV%Zi7(9qv%WFLgBk<*$zT(E8RO$qa&q$X zz}DSmj)$`97~|uUb5axY;!7$EQsXo8((>c;5_40F<1_P0Qj78ub3pq;T=SCiQ&I!+ z!D0*&j4r9je0!F_($u2LfYhS2{G!~%yyR5p{8Er^20bSC)Dq{?qN3Ei5>N{?F~u{5 zf#EkuG`O<3BsDiAGdI<-*ex?B6~tnYV0KR}3CT>(1{+{z!oXm}=#p8Sn3R(W(e9I( zRFqg$>71C8Q_SG)8sz8dV{8P8b=33?;=%UrJo16mlc4epgm3sTFx>HBV0Zv3F9aDF zzJh25_oBpt4Cm4k-_#1d+|-J4GZRo^bt}zFhNKCQi6a-6&VaN<V1piKV-+r?T>*_u zjNC98Q9c)7V9*p`V3^^~zyNBOeiUF}P;g~nU=U<rh!SI9IN-v-kRZmupexS6VB^8S zU?9%G@EpqiAkM&GFU7zR;KqQ`?gX_z733Hg_JI1TatsU|atsW2Ky6{9b|$F332Q@w zXj>%)1_3t)1`wT~#K2(S#=rogw<$3&IJhyOwCj-DZ=iPFOBDu&J8ldNAo`sO1H%V5 z1_lt#0$N4l&VbUk1hp|?ZAcL9sLsGJ#f^ahM7yaoFsyK6U;xoE>I@7!KyFZH0Jrl5 zK;eVb-UGG!YBd-bOxzh5K=f)21_lpz28In93=I1;85lC$85j;|GB7OBVql1HXJELa z#lT>t4GCu(Z3c#RT?U34cLoN~KEtEB5H%-s85sWQF)%D}XJAOsXJA;X&%p4+oq^$j zJ_CcE0YtBX0RzJW0|o{T4+#C;fPulpgMs0P0RzJZLk0#74+wqP5F&oXkb!~62oiT7 z+6Br6^)*2Kp*gOQ_9Ljxc*h;mcKmL{zyM-17(?2BAU3ESXl)D?A5j@(WOx>`^^=hS zG;S+_)bEFhgT|~3u!wWAFff3c>JbGD3{i{>3{M~$7;F-Yi&Jxxaw-)na&z*EZIw!k z@~n!JGg5OCi*<7|lZ*0;^V3Rnlk;<}5{q;7$_$l2!%vxMsl_E>sYRfcsI8Kro`I4b zFPDOXO-WH{afxSMT0W7c8AJ3bDA*LICYKgvmQ;exQBY7&umO#Zr52Z@rUVpamSyIo zrl%IeC6T3FD^in7K}`go)UwnZg&YuRtCU#mnOBydom!-%P@3tOoSa%*Y^#)(m{Xjp zWM`ufHyc$m+<YVz7`EH!qqx;ZAMR;8UM?GbSa{k6I0iVnxOn=xM?^pwP>L)*XiqO_ U7$Ra}!orRvBO|Sk`VtWW0Pax%(f|Me literal 0 HcmV?d00001 diff --git a/Chapter07/Chapter07_mcode/GraphCut/MaxFlow.cpp b/Chapter07/Chapter07_mcode/GraphCut/MaxFlow.cpp new file mode 100755 index 0000000..b91856e --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/MaxFlow.cpp @@ -0,0 +1,690 @@ +/* maxflow.cpp */ + + +#include <stdio.h> +#include "graph.h" +#include "instances.inc" + + +/* + special constants for node->parent +*/ +#define TERMINAL ( (arc *) 1 ) /* to terminal */ +#define ORPHAN ( (arc *) 2 ) /* orphan */ + + +#define INFINITE_D ((int)(((unsigned)-1)/2)) /* infinite distance to the terminal */ + +/***********************************************************************/ + +/* + Functions for processing active list. + i->next points to the next node in the list + (or to i, if i is the last node in the list). + If i->next is NULL iff i is not in the list. + + There are two queues. Active nodes are added + to the end of the second queue and read from + the front of the first queue. If the first queue + is empty, it is replaced by the second queue + (and the second queue becomes empty). +*/ + + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::set_active(node *i) +{ + if (!i->next) + { + /* it's not in the list yet */ + if (queue_last[1]) queue_last[1] -> next = i; + else queue_first[1] = i; + queue_last[1] = i; + i -> next = i; + } +} + +/* + Returns the next active node. + If it is connected to the sink, it stays in the list, + otherwise it is removed from the list +*/ +template <typename captype, typename tcaptype, typename flowtype> + inline typename Graph<captype,tcaptype,flowtype>::node* Graph<captype,tcaptype,flowtype>::next_active() +{ + node *i; + + while ( 1 ) + { + if (!(i=queue_first[0])) + { + queue_first[0] = i = queue_first[1]; + queue_last[0] = queue_last[1]; + queue_first[1] = NULL; + queue_last[1] = NULL; + if (!i) return NULL; + } + + /* remove it from the active list */ + if (i->next == i) queue_first[0] = queue_last[0] = NULL; + else queue_first[0] = i -> next; + i -> next = NULL; + + /* a node in the list is active iff it has a parent */ + if (i->parent) return i; + } +} + +/***********************************************************************/ + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::set_orphan_front(node *i) +{ + nodeptr *np; + i -> parent = ORPHAN; + np = nodeptr_block -> New(); + np -> ptr = i; + np -> next = orphan_first; + orphan_first = np; +} + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::set_orphan_rear(node *i) +{ + nodeptr *np; + i -> parent = ORPHAN; + np = nodeptr_block -> New(); + np -> ptr = i; + if (orphan_last) orphan_last -> next = np; + else orphan_first = np; + orphan_last = np; + np -> next = NULL; +} + +/***********************************************************************/ + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::add_to_changed_list(node *i) +{ + if (changed_list && !i->is_in_changed_list) + { + node_id* ptr = changed_list->New(); + *ptr = (node_id)(i - nodes); + i->is_in_changed_list = true; + } +} + +/***********************************************************************/ + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::maxflow_init() +{ + node *i; + + queue_first[0] = queue_last[0] = NULL; + queue_first[1] = queue_last[1] = NULL; + orphan_first = NULL; + + TIME = 0; + + for (i=nodes; i<node_last; i++) + { + i -> next = NULL; + i -> is_marked = 0; + i -> is_in_changed_list = 0; + i -> TS = TIME; + if (i->tr_cap > 0) + { + /* i is connected to the source */ + i -> is_sink = 0; + i -> parent = TERMINAL; + set_active(i); + i -> DIST = 1; + } + else if (i->tr_cap < 0) + { + /* i is connected to the sink */ + i -> is_sink = 1; + i -> parent = TERMINAL; + set_active(i); + i -> DIST = 1; + } + else + { + i -> parent = NULL; + } + } +} + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::maxflow_reuse_trees_init() +{ + node* i; + node* j; + node* queue = queue_first[1]; + arc* a; + nodeptr* np; + + queue_first[0] = queue_last[0] = NULL; + queue_first[1] = queue_last[1] = NULL; + orphan_first = orphan_last = NULL; + + TIME ++; + + while ((i=queue)) + { + queue = i->next; + if (queue == i) queue = NULL; + i->next = NULL; + i->is_marked = 0; + set_active(i); + + if (i->tr_cap == 0) + { + if (i->parent) set_orphan_rear(i); + continue; + } + + if (i->tr_cap > 0) + { + if (!i->parent || i->is_sink) + { + i->is_sink = 0; + for (a=i->first; a; a=a->next) + { + j = a->head; + if (!j->is_marked) + { + if (j->parent == a->sister) set_orphan_rear(j); + if (j->parent && j->is_sink && a->r_cap > 0) set_active(j); + } + } + add_to_changed_list(i); + } + } + else + { + if (!i->parent || !i->is_sink) + { + i->is_sink = 1; + for (a=i->first; a; a=a->next) + { + j = a->head; + if (!j->is_marked) + { + if (j->parent == a->sister) set_orphan_rear(j); + if (j->parent && !j->is_sink && a->sister->r_cap > 0) set_active(j); + } + } + add_to_changed_list(i); + } + } + i->parent = TERMINAL; + i -> TS = TIME; + i -> DIST = 1; + } + + //test_consistency(); + + /* adoption */ + while ((np=orphan_first)) + { + orphan_first = np -> next; + i = np -> ptr; + nodeptr_block -> Delete(np); + if (!orphan_first) orphan_last = NULL; + if (i->is_sink) process_sink_orphan(i); + else process_source_orphan(i); + } + /* adoption end */ + + //test_consistency(); +} + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::augment(arc *middle_arc) +{ + node *i; + arc *a; + tcaptype bottleneck; + + + /* 1. Finding bottleneck capacity */ + /* 1a - the source tree */ + bottleneck = middle_arc -> r_cap; + for (i=middle_arc->sister->head; ; i=a->head) + { + a = i -> parent; + if (a == TERMINAL) break; + if (bottleneck > a->sister->r_cap) bottleneck = a -> sister -> r_cap; + } + if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap; + /* 1b - the sink tree */ + for (i=middle_arc->head; ; i=a->head) + { + a = i -> parent; + if (a == TERMINAL) break; + if (bottleneck > a->r_cap) bottleneck = a -> r_cap; + } + if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap; + + + /* 2. Augmenting */ + /* 2a - the source tree */ + middle_arc -> sister -> r_cap += bottleneck; + middle_arc -> r_cap -= bottleneck; + for (i=middle_arc->sister->head; ; i=a->head) + { + a = i -> parent; + if (a == TERMINAL) break; + a -> r_cap += bottleneck; + a -> sister -> r_cap -= bottleneck; + if (!a->sister->r_cap) + { + set_orphan_front(i); // add i to the beginning of the adoption list + } + } + i -> tr_cap -= bottleneck; + if (!i->tr_cap) + { + set_orphan_front(i); // add i to the beginning of the adoption list + } + /* 2b - the sink tree */ + for (i=middle_arc->head; ; i=a->head) + { + a = i -> parent; + if (a == TERMINAL) break; + a -> sister -> r_cap += bottleneck; + a -> r_cap -= bottleneck; + if (!a->r_cap) + { + set_orphan_front(i); // add i to the beginning of the adoption list + } + } + i -> tr_cap += bottleneck; + if (!i->tr_cap) + { + set_orphan_front(i); // add i to the beginning of the adoption list + } + + + flow += bottleneck; +} + +/***********************************************************************/ + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::process_source_orphan(node *i) +{ + node *j; + arc *a0, *a0_min = NULL, *a; + int d, d_min = INFINITE_D; + + /* trying to find a new parent */ + for (a0=i->first; a0; a0=a0->next) + if (a0->sister->r_cap) + { + j = a0 -> head; + if (!j->is_sink && (a=j->parent)) + { + /* checking the origin of j */ + d = 0; + while ( 1 ) + { + if (j->TS == TIME) + { + d += j -> DIST; + break; + } + a = j -> parent; + d ++; + if (a==TERMINAL) + { + j -> TS = TIME; + j -> DIST = 1; + break; + } + if (a==ORPHAN) { d = INFINITE_D; break; } + j = a -> head; + } + if (d<INFINITE_D) /* j originates from the source - done */ + { + if (d<d_min) + { + a0_min = a0; + d_min = d; + } + /* set marks along the path */ + for (j=a0->head; j->TS!=TIME; j=j->parent->head) + { + j -> TS = TIME; + j -> DIST = d --; + } + } + } + } + + if (i->parent = a0_min) + { + i -> TS = TIME; + i -> DIST = d_min + 1; + } + else + { + /* no parent is found */ + add_to_changed_list(i); + + /* process neighbors */ + for (a0=i->first; a0; a0=a0->next) + { + j = a0 -> head; + if (!j->is_sink && (a=j->parent)) + { + if (a0->sister->r_cap) set_active(j); + if (a!=TERMINAL && a!=ORPHAN && a->head==i) + { + set_orphan_rear(j); // add j to the end of the adoption list + } + } + } + } +} + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::process_sink_orphan(node *i) +{ + node *j; + arc *a0, *a0_min = NULL, *a; + int d, d_min = INFINITE_D; + + /* trying to find a new parent */ + for (a0=i->first; a0; a0=a0->next) + if (a0->r_cap) + { + j = a0 -> head; + if (j->is_sink && (a=j->parent)) + { + /* checking the origin of j */ + d = 0; + while ( 1 ) + { + if (j->TS == TIME) + { + d += j -> DIST; + break; + } + a = j -> parent; + d ++; + if (a==TERMINAL) + { + j -> TS = TIME; + j -> DIST = 1; + break; + } + if (a==ORPHAN) { d = INFINITE_D; break; } + j = a -> head; + } + if (d<INFINITE_D) /* j originates from the sink - done */ + { + if (d<d_min) + { + a0_min = a0; + d_min = d; + } + /* set marks along the path */ + for (j=a0->head; j->TS!=TIME; j=j->parent->head) + { + j -> TS = TIME; + j -> DIST = d --; + } + } + } + } + + if (i->parent = a0_min) + { + i -> TS = TIME; + i -> DIST = d_min + 1; + } + else + { + /* no parent is found */ + add_to_changed_list(i); + + /* process neighbors */ + for (a0=i->first; a0; a0=a0->next) + { + j = a0 -> head; + if (j->is_sink && (a=j->parent)) + { + if (a0->r_cap) set_active(j); + if (a!=TERMINAL && a!=ORPHAN && a->head==i) + { + set_orphan_rear(j); // add j to the end of the adoption list + } + } + } + } +} + +/***********************************************************************/ + +template <typename captype, typename tcaptype, typename flowtype> + flowtype Graph<captype,tcaptype,flowtype>::maxflow(bool reuse_trees, Block<node_id>* _changed_list) +{ + node *i, *j, *current_node = NULL; + arc *a; + nodeptr *np, *np_next; + + if (!nodeptr_block) + { + nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function); + } + + changed_list = _changed_list; + if (maxflow_iteration == 0 && reuse_trees) { if (error_function) (*error_function)("reuse_trees cannot be used in the first call to maxflow()!"); exit(1); } + if (changed_list && !reuse_trees) { if (error_function) (*error_function)("changed_list cannot be used without reuse_trees!"); exit(1); } + + if (reuse_trees) maxflow_reuse_trees_init(); + else maxflow_init(); + + // main loop + while ( 1 ) + { + // test_consistency(current_node); + + if ((i=current_node)) + { + i -> next = NULL; /* remove active flag */ + if (!i->parent) i = NULL; + } + if (!i) + { + if (!(i = next_active())) break; + } + + /* growth */ + if (!i->is_sink) + { + /* grow source tree */ + for (a=i->first; a; a=a->next) + if (a->r_cap) + { + j = a -> head; + if (!j->parent) + { + j -> is_sink = 0; + j -> parent = a -> sister; + j -> TS = i -> TS; + j -> DIST = i -> DIST + 1; + set_active(j); + add_to_changed_list(j); + } + else if (j->is_sink) break; + else if (j->TS <= i->TS && + j->DIST > i->DIST) + { + /* heuristic - trying to make the distance from j to the source shorter */ + j -> parent = a -> sister; + j -> TS = i -> TS; + j -> DIST = i -> DIST + 1; + } + } + } + else + { + /* grow sink tree */ + for (a=i->first; a; a=a->next) + if (a->sister->r_cap) + { + j = a -> head; + if (!j->parent) + { + j -> is_sink = 1; + j -> parent = a -> sister; + j -> TS = i -> TS; + j -> DIST = i -> DIST + 1; + set_active(j); + add_to_changed_list(j); + } + else if (!j->is_sink) { a = a -> sister; break; } + else if (j->TS <= i->TS && + j->DIST > i->DIST) + { + /* heuristic - trying to make the distance from j to the sink shorter */ + j -> parent = a -> sister; + j -> TS = i -> TS; + j -> DIST = i -> DIST + 1; + } + } + } + + TIME ++; + + if (a) + { + i -> next = i; /* set active flag */ + current_node = i; + + /* augmentation */ + augment(a); + /* augmentation end */ + + /* adoption */ + while ((np=orphan_first)) + { + np_next = np -> next; + np -> next = NULL; + + while ((np=orphan_first)) + { + orphan_first = np -> next; + i = np -> ptr; + nodeptr_block -> Delete(np); + if (!orphan_first) orphan_last = NULL; + if (i->is_sink) process_sink_orphan(i); + else process_source_orphan(i); + } + + orphan_first = np_next; + } + /* adoption end */ + } + else current_node = NULL; + } + // test_consistency(); + + if (!reuse_trees || (maxflow_iteration % 64) == 0) + { + delete nodeptr_block; + nodeptr_block = NULL; + } + + maxflow_iteration ++; + return flow; +} + +/***********************************************************************/ + + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::test_consistency(node* current_node) +{ + node *i; + arc *a; + int r; + int num1 = 0, num2 = 0; + + // test whether all nodes i with i->next!=NULL are indeed in the queue + for (i=nodes; i<node_last; i++) + { + if (i->next || i==current_node) num1 ++; + } + for (r=0; r<3; r++) + { + i = (r == 2) ? current_node : queue_first[r]; + if (i) + for ( ; ; i=i->next) + { + num2 ++; + if (i->next == i) + { + if (r<2) assert(i == queue_last[r]); + else assert(i == current_node); + break; + } + } + } + assert(num1 == num2); + + for (i=nodes; i<node_last; i++) + { + // test whether all edges in seach trees are non-saturated + if (i->parent == NULL) {} + else if (i->parent == ORPHAN) {} + else if (i->parent == TERMINAL) + { + if (!i->is_sink) assert(i->tr_cap > 0); + else assert(i->tr_cap < 0); + } + else + { + if (!i->is_sink) assert (i->parent->sister->r_cap > 0); + else assert (i->parent->r_cap > 0); + } + // test whether passive nodes in search trees have neighbors in + // a different tree through non-saturated edges + if (i->parent && !i->next) + { + if (!i->is_sink) + { + assert(i->tr_cap >= 0); + for (a=i->first; a; a=a->next) + { + if (a->r_cap > 0) assert(a->head->parent && !a->head->is_sink); + } + } + else + { + assert(i->tr_cap <= 0); + for (a=i->first; a; a=a->next) + { + if (a->sister->r_cap > 0) assert(a->head->parent && a->head->is_sink); + } + } + } + // test marking invariants + if (i->parent && i->parent!=ORPHAN && i->parent!=TERMINAL) + { + assert(i->TS <= i->parent->head->TS); + if (i->TS == i->parent->head->TS) assert(i->DIST > i->parent->head->DIST); + } + } +} + + + void Inst() + { + Graph<int,int,int> G(1,1); + G.maxflow(); + } \ No newline at end of file diff --git a/Chapter07/Chapter07_mcode/GraphCut/block.h b/Chapter07/Chapter07_mcode/GraphCut/block.h new file mode 100755 index 0000000..14346ef --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/block.h @@ -0,0 +1,268 @@ +/* block.h */ +/* + Template classes Block and DBlock + Implement adding and deleting items of the same type in blocks. + + If there there are many items then using Block or DBlock + is more efficient than using 'new' and 'delete' both in terms + of memory and time since + (1) On some systems there is some minimum amount of memory + that 'new' can allocate (e.g., 64), so if items are + small that a lot of memory is wasted. + (2) 'new' and 'delete' are designed for items of varying size. + If all items has the same size, then an algorithm for + adding and deleting can be made more efficient. + (3) All Block and DBlock functions are inline, so there are + no extra function calls. + + Differences between Block and DBlock: + (1) DBlock allows both adding and deleting items, + whereas Block allows only adding items. + (2) Block has an additional operation of scanning + items added so far (in the order in which they were added). + (3) Block allows to allocate several consecutive + items at a time, whereas DBlock can add only a single item. + + Note that no constructors or destructors are called for items. + + Example usage for items of type 'MyType': + + /////////////////////////////////////////////////// + #include "block.h" + #define BLOCK_SIZE 1024 + typedef struct { int a, b; } MyType; + MyType *ptr, *array[10000]; + + ... + + Block<MyType> *block = new Block<MyType>(BLOCK_SIZE); + + // adding items + for (int i=0; i<sizeof(array); i++) + { + ptr = block -> New(); + ptr -> a = ptr -> b = rand(); + } + + // reading items + for (ptr=block->ScanFirst(); ptr; ptr=block->ScanNext()) + { + printf("%d %d\n", ptr->a, ptr->b); + } + + delete block; + + ... + + DBlock<MyType> *dblock = new DBlock<MyType>(BLOCK_SIZE); + + // adding items + for (int i=0; i<sizeof(array); i++) + { + array[i] = dblock -> New(); + } + + // deleting items + for (int i=0; i<sizeof(array); i+=2) + { + dblock -> Delete(array[i]); + } + + // adding items + for (int i=0; i<sizeof(array); i++) + { + array[i] = dblock -> New(); + } + + delete dblock; + + /////////////////////////////////////////////////// + + Note that DBlock deletes items by marking them as + empty (i.e., by adding them to the list of free items), + so that this memory could be used for subsequently + added items. Thus, at each moment the memory allocated + is determined by the maximum number of items allocated + simultaneously at earlier moments. All memory is + deallocated only when the destructor is called. +*/ + +#ifndef __BLOCK_H__ +#define __BLOCK_H__ + +#include <stdlib.h> + +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ + +template <class Type> class Block +{ +public: + /* Constructor. Arguments are the block size and + (optionally) the pointer to the function which + will be called if allocation failed; the message + passed to this function is "Not enough memory!" */ + Block(int size, void (*err_function)(char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; } + + /* Destructor. Deallocates all items added so far */ + ~Block() { while (first) { block *next = first -> next; delete first; first = next; } } + + /* Allocates 'num' consecutive items; returns pointer + to the first item. 'num' cannot be greater than the + block size since items must fit in one block */ + Type *New(int num = 1) + { + Type *t; + + if (!last || last->current + num > last->last) + { + if (last && last->next) last = last -> next; + else + { + block *next = (block *) new char [sizeof(block) + (block_size-1)*sizeof(Type)]; + if (!next) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + if (last) last -> next = next; + else first = next; + last = next; + last -> current = & ( last -> data[0] ); + last -> last = last -> current + block_size; + last -> next = NULL; + } + } + + t = last -> current; + last -> current += num; + return t; + } + + /* Returns the first item (or NULL, if no items were added) */ + Type *ScanFirst() + { + for (scan_current_block=first; scan_current_block; scan_current_block = scan_current_block->next) + { + scan_current_data = & ( scan_current_block -> data[0] ); + if (scan_current_data < scan_current_block -> current) return scan_current_data ++; + } + return NULL; + } + + /* Returns the next item (or NULL, if all items have been read) + Can be called only if previous ScanFirst() or ScanNext() + call returned not NULL. */ + Type *ScanNext() + { + while (scan_current_data >= scan_current_block -> current) + { + scan_current_block = scan_current_block -> next; + if (!scan_current_block) return NULL; + scan_current_data = & ( scan_current_block -> data[0] ); + } + return scan_current_data ++; + } + + /* Marks all elements as empty */ + void Reset() + { + block *b; + if (!first) return; + for (b=first; ; b=b->next) + { + b -> current = & ( b -> data[0] ); + if (b == last) break; + } + last = first; + } + +/***********************************************************************/ + +private: + + typedef struct block_st + { + Type *current, *last; + struct block_st *next; + Type data[1]; + } block; + + int block_size; + block *first; + block *last; + + block *scan_current_block; + Type *scan_current_data; + + void (*error_function)(char *); +}; + +/***********************************************************************/ +/***********************************************************************/ +/***********************************************************************/ + +template <class Type> class DBlock +{ +public: + /* Constructor. Arguments are the block size and + (optionally) the pointer to the function which + will be called if allocation failed; the message + passed to this function is "Not enough memory!" */ + DBlock(int size, void (*err_function)(char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; } + + /* Destructor. Deallocates all items added so far */ + ~DBlock() { while (first) { block *next = first -> next; delete first; first = next; } } + + /* Allocates one item */ + Type *New() + { + block_item *item; + + if (!first_free) + { + block *next = first; + first = (block *) new char [sizeof(block) + (block_size-1)*sizeof(block_item)]; + if (!first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + first_free = & (first -> data[0] ); + for (item=first_free; item<first_free+block_size-1; item++) + item -> next_free = item + 1; + item -> next_free = NULL; + first -> next = next; + } + + item = first_free; + first_free = item -> next_free; + return (Type *) item; + } + + /* Deletes an item allocated previously */ + void Delete(Type *t) + { + ((block_item *) t) -> next_free = first_free; + first_free = (block_item *) t; + } + +/***********************************************************************/ + +private: + + typedef union block_item_st + { + Type t; + block_item_st *next_free; + } block_item; + + typedef struct block_st + { + struct block_st *next; + block_item data[1]; + } block; + + int block_size; + block *first; + block_item *first_free; + + void (*error_function)(char *); +}; + + +#endif + diff --git a/Chapter07/Chapter07_mcode/GraphCut/compile_mex_functions.m b/Chapter07/Chapter07_mcode/GraphCut/compile_mex_functions.m new file mode 100755 index 0000000..7bcfafe --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/compile_mex_functions.m @@ -0,0 +1,3 @@ +% compiling mex functions + +mex GraphCutMex.cpp MaxFlow.cpp graph.cpp diff --git a/Chapter07/Chapter07_mcode/GraphCut/graph.cpp b/Chapter07/Chapter07_mcode/GraphCut/graph.cpp new file mode 100755 index 0000000..8ba3d5a --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/graph.cpp @@ -0,0 +1,114 @@ +/* graph.cpp */ + + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include "graph.h" +#include "instances.inc" + + +template <typename captype, typename tcaptype, typename flowtype> + Graph<captype, tcaptype, flowtype>::Graph(int node_num_max, int edge_num_max, void (*err_function)(char *)) + : node_num(0), + nodeptr_block(NULL), + error_function(err_function) +{ + if (node_num_max < 16) node_num_max = 16; + if (edge_num_max < 16) edge_num_max = 16; + + nodes = (node*) malloc(node_num_max*sizeof(node)); + arcs = (arc*) malloc(2*edge_num_max*sizeof(arc)); + if (!nodes || !arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + + node_last = nodes; + node_max = nodes + node_num_max; + arc_last = arcs; + arc_max = arcs + 2*edge_num_max; + + maxflow_iteration = 0; + flow = 0; +} + +template <typename captype, typename tcaptype, typename flowtype> + Graph<captype,tcaptype,flowtype>::~Graph() +{ + if (nodeptr_block) + { + delete nodeptr_block; + nodeptr_block = NULL; + } + free(nodes); + free(arcs); +} + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::reset() +{ + node_last = nodes; + arc_last = arcs; + node_num = 0; + + if (nodeptr_block) + { + delete nodeptr_block; + nodeptr_block = NULL; + } + + maxflow_iteration = 0; + flow = 0; +} + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::reallocate_nodes(int num) +{ + int node_num_max = (int)(node_max - nodes); + node* nodes_old = nodes; + + node_num_max += node_num_max / 2; + if (node_num_max < node_num + num) node_num_max = node_num + num; + nodes = (node*) realloc(nodes_old, node_num_max*sizeof(node)); + if (!nodes) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + + node_last = nodes + node_num; + node_max = nodes + node_num_max; + + if (nodes != nodes_old) + { + arc* a; + for (a=arcs; a<arc_last; a++) + { + a->head = (node*) ((char*)a->head + (((char*) nodes) - ((char*) nodes_old))); + } + } +} + +template <typename captype, typename tcaptype, typename flowtype> + void Graph<captype,tcaptype,flowtype>::reallocate_arcs() +{ + int arc_num_max = (int)(arc_max - arcs); + int arc_num = (int)(arc_last - arcs); + arc* arcs_old = arcs; + + arc_num_max += arc_num_max / 2; if (arc_num_max & 1) arc_num_max ++; + arcs = (arc*) realloc(arcs_old, arc_num_max*sizeof(arc)); + if (!arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); } + + arc_last = arcs + arc_num; + arc_max = arcs + arc_num_max; + + if (arcs != arcs_old) + { + node* i; + arc* a; + for (i=nodes; i<node_last; i++) + { + if (i->first) i->first = (arc*) ((char*)i->first + (((char*) arcs) - ((char*) arcs_old))); + } + for (a=arcs; a<arc_last; a++) + { + if (a->next) a->next = (arc*) ((char*)a->next + (((char*) arcs) - ((char*) arcs_old))); + a->sister = (arc*) ((char*)a->sister + (((char*) arcs) - ((char*) arcs_old))); + } + } +} diff --git a/Chapter07/Chapter07_mcode/GraphCut/graph.h b/Chapter07/Chapter07_mcode/GraphCut/graph.h new file mode 100755 index 0000000..b45dd36 --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/graph.h @@ -0,0 +1,506 @@ +/* graph.h */ +/* + This software library implements the maxflow algorithm + described in + + "An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision." + Yuri Boykov and Vladimir Kolmogorov. + In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), + September 2004 + + This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov + at Siemens Corporate Research. To make it available for public use, + it was later reimplemented by Vladimir Kolmogorov based on open publications. + + If you use this software for research purposes, you should cite + the aforementioned paper in any resulting publication. + + ---------------------------------------------------------------------- + + REUSING TREES: + + Starting with version 3.0, there is a also an option of reusing search + trees from one maxflow computation to the next, as described in + + "Efficiently Solving Dynamic Markov Random Fields Using Graph Cuts." + Pushmeet Kohli and Philip H.S. Torr + International Conference on Computer Vision (ICCV), 2005 + + If you use this option, you should cite + the aforementioned paper in any resulting publication. +*/ + + + +/* + For description, license, example usage see README.TXT. +*/ + +#ifndef __GRAPH_H__ +#define __GRAPH_H__ + +#include <string.h> +#include "block.h" + +#include <assert.h> +// NOTE: in UNIX you need to use -DNDEBUG preprocessor option to supress assert's!!! + + + +// captype: type of edge capacities (excluding t-links) +// tcaptype: type of t-links (edges between nodes and terminals) +// flowtype: type of total flow +// +// Current instantiations are in instances.inc +template <typename captype, typename tcaptype, typename flowtype> class Graph +{ +public: + typedef enum + { + SOURCE = 0, + SINK = 1 + } termtype; // terminals + typedef int node_id; + + ///////////////////////////////////////////////////////////////////////// + // BASIC INTERFACE FUNCTIONS // + // (should be enough for most applications) // + ///////////////////////////////////////////////////////////////////////// + + // Constructor. + // The first argument gives an estimate of the maximum number of nodes that can be added + // to the graph, and the second argument is an estimate of the maximum number of edges. + // The last (optional) argument is the pointer to the function which will be called + // if an error occurs; an error message is passed to this function. + // If this argument is omitted, exit(1) will be called. + // + // IMPORTANT: It is possible to add more nodes to the graph than node_num_max + // (and node_num_max can be zero). However, if the count is exceeded, then + // the internal memory is reallocated (increased by 50%) which is expensive. + // Also, temporarily the amount of allocated memory would be more than twice than needed. + // Similarly for edges. + // If you wish to avoid this overhead, you can download version 2.2, where nodes and edges are stored in blocks. + Graph(int node_num_max, int edge_num_max, void (*err_function)(char *) = NULL); + + // Destructor + ~Graph(); + + // Adds node(s) to the graph. By default, one node is added (num=1); then first call returns 0, second call returns 1, and so on. + // If num>1, then several nodes are added, and node_id of the first one is returned. + // IMPORTANT: see note about the constructor + node_id add_node(int num = 1); + + // Adds a bidirectional edge between 'i' and 'j' with the weights 'cap' and 'rev_cap'. + // IMPORTANT: see note about the constructor + void add_edge(node_id i, node_id j, captype cap, captype rev_cap); + + // Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights. + // Can be called multiple times for each node. + // Weights can be negative. + // NOTE: the number of such edges is not counted in edge_num_max. + // No internal memory is allocated by this call. + void add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink); + + + // Computes the maxflow. Can be called several times. + // FOR DESCRIPTION OF reuse_trees, SEE mark_node(). + // FOR DESCRIPTION OF changed_list, SEE remove_from_changed_list(). + flowtype maxflow(bool reuse_trees = false, Block<node_id>* changed_list = NULL); + + // After the maxflow is computed, this function returns to which + // segment the node 'i' belongs (Graph<captype,tcaptype,flowtype>::SOURCE or Graph<captype,tcaptype,flowtype>::SINK). + // + // Occasionally there may be several minimum cuts. If a node can be assigned + // to both the source and the sink, then default_segm is returned. + termtype what_segment(node_id i, termtype default_segm = SOURCE); + + + + ////////////////////////////////////////////// + // ADVANCED INTERFACE FUNCTIONS // + // (provide access to the graph) // + ////////////////////////////////////////////// + +private: + struct node; + struct arc; + +public: + + //////////////////////////// + // 1. Reallocating graph. // + //////////////////////////// + + // Removes all nodes and edges. + // After that functions add_node() and add_edge() must be called again. + // + // Advantage compared to deleting Graph and allocating it again: + // no calls to delete/new (which could be quite slow). + // + // If the graph structure stays the same, then an alternative + // is to go through all nodes/edges and set new residual capacities + // (see functions below). + void reset(); + + //////////////////////////////////////////////////////////////////////////////// + // 2. Functions for getting pointers to arcs and for reading graph structure. // + // NOTE: adding new arcs may invalidate these pointers (if reallocation // + // happens). So it's best not to add arcs while reading graph structure. // + //////////////////////////////////////////////////////////////////////////////// + + // The following two functions return arcs in the same order that they + // were added to the graph. NOTE: for each call add_edge(i,j,cap,cap_rev) + // the first arc returned will be i->j, and the second j->i. + // If there are no more arcs, then the function can still be called, but + // the returned arc_id is undetermined. + typedef arc* arc_id; + arc_id get_first_arc(); + arc_id get_next_arc(arc_id a); + + // other functions for reading graph structure + int get_node_num() { return node_num; } + int get_arc_num() { return (int)(arc_last - arcs); } + void get_arc_ends(arc_id a, node_id& i, node_id& j); // returns i,j to that a = i->j + + /////////////////////////////////////////////////// + // 3. Functions for reading residual capacities. // + /////////////////////////////////////////////////// + + // returns residual capacity of SOURCE->i minus residual capacity of i->SINK + tcaptype get_trcap(node_id i); + // returns residual capacity of arc a + captype get_rcap(arc* a); + + ///////////////////////////////////////////////////////////////// + // 4. Functions for setting residual capacities. // + // NOTE: If these functions are used, the value of the flow // + // returned by maxflow() will not be valid! // + ///////////////////////////////////////////////////////////////// + + void set_trcap(node_id i, tcaptype trcap); + void set_rcap(arc* a, captype rcap); + + //////////////////////////////////////////////////////////////////// + // 5. Functions related to reusing trees & list of changed nodes. // + //////////////////////////////////////////////////////////////////// + + // If flag reuse_trees is true while calling maxflow(), then search trees + // are reused from previous maxflow computation. + // In this case before calling maxflow() the user must + // specify which parts of the graph have changed by calling mark_node(): + // add_tweights(i),set_trcap(i) => call mark_node(i) + // add_edge(i,j),set_rcap(a) => call mark_node(i); mark_node(j) + // + // This option makes sense only if a small part of the graph is changed. + // The initialization procedure goes only through marked nodes then. + // + // mark_node(i) can either be called before or after graph modification. + // Can be called more than once per node, but calls after the first one + // do not have any effect. + // + // NOTE: + // - This option cannot be used in the first call to maxflow(). + // - It is not necessary to call mark_node() if the change is ``not essential'', + // i.e. sign(trcap) is preserved for a node and zero/nonzero status is preserved for an arc. + // - To check that you marked all necessary nodes, you can call maxflow(false) after calling maxflow(true). + // If everything is correct, the two calls must return the same value of flow. (Useful for debugging). + void mark_node(node_id i); + + // If changed_list is not NULL while calling maxflow(), then the algorithm + // keeps a list of nodes which could potentially have changed their segmentation label. + // Nodes which are not in the list are guaranteed to keep their old segmentation label (SOURCE or SINK). + // Example usage: + // + // typedef Graph<int,int,int> G; + // G* g = new Graph(nodeNum, edgeNum); + // Block<G::node_id>* changed_list = new Block<G::node_id>(128); + // + // ... // add nodes and edges + // + // g->maxflow(); // first call should be without arguments + // for (int iter=0; iter<10; iter++) + // { + // ... // change graph, call mark_node() accordingly + // + // g->maxflow(true, changed_list); + // G::node_id* ptr; + // for (ptr=changed_list->ScanFirst(); ptr; ptr=changed_list->ScanNext()) + // { + // G::node_id i = *ptr; assert(i>=0 && i<nodeNum); + // g->remove_from_changed_list(i); + // // do something with node i... + // if (g->what_segment(i) == G::SOURCE) { ... } + // } + // changed_list->Reset(); + // } + // delete changed_list; + // + // NOTE: + // - If changed_list option is used, then reuse_trees must be used as well. + // - In the example above, the user may omit calls g->remove_from_changed_list(i) and changed_list->Reset() in a given iteration. + // Then during the next call to maxflow(true, &changed_list) new nodes will be added to changed_list. + // - If the next call to maxflow() does not use option reuse_trees, then calling remove_from_changed_list() + // is not necessary. ("changed_list->Reset()" or "delete changed_list" should still be called, though). + void remove_from_changed_list(node_id i) + { + assert(i>=0 && i<node_num && nodes[i].is_in_changed_list); + nodes[i].is_in_changed_list = 0; + } + + + + + + +///////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////// + +private: + // internal variables and functions + + struct node + { + arc *first; // first outcoming arc + + arc *parent; // node's parent + node *next; // pointer to the next active node + // (or to itself if it is the last node in the list) + int TS; // timestamp showing when DIST was computed + int DIST; // distance to the terminal + int is_sink : 1; // flag showing whether the node is in the source or in the sink tree (if parent!=NULL) + int is_marked : 1; // set by mark_node() + int is_in_changed_list : 1; // set by maxflow if + + tcaptype tr_cap; // if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node + // otherwise -tr_cap is residual capacity of the arc node->SINK + + }; + + struct arc + { + node *head; // node the arc points to + arc *next; // next arc with the same originating node + arc *sister; // reverse arc + + captype r_cap; // residual capacity + }; + + struct nodeptr + { + node *ptr; + nodeptr *next; + }; + static const int NODEPTR_BLOCK_SIZE = 128; + + node *nodes, *node_last, *node_max; // node_last = nodes+node_num, node_max = nodes+node_num_max; + arc *arcs, *arc_last, *arc_max; // arc_last = arcs+2*edge_num, arc_max = arcs+2*edge_num_max; + + int node_num; + + DBlock<nodeptr> *nodeptr_block; + + void (*error_function)(char *); // this function is called if a error occurs, + // with a corresponding error message + // (or exit(1) is called if it's NULL) + + flowtype flow; // total flow + + // reusing trees & list of changed pixels + int maxflow_iteration; // counter + Block<node_id> *changed_list; + + ///////////////////////////////////////////////////////////////////////// + + node *queue_first[2], *queue_last[2]; // list of active nodes + nodeptr *orphan_first, *orphan_last; // list of pointers to orphans + int TIME; // monotonically increasing global counter + + ///////////////////////////////////////////////////////////////////////// + + void reallocate_nodes(int num); // num is the number of new nodes + void reallocate_arcs(); + + // functions for processing active list + void set_active(node *i); + node *next_active(); + + // functions for processing orphans list + void set_orphan_front(node* i); // add to the beginning of the list + void set_orphan_rear(node* i); // add to the end of the list + + void add_to_changed_list(node* i); + + void maxflow_init(); // called if reuse_trees == false + void maxflow_reuse_trees_init(); // called if reuse_trees == true + void augment(arc *middle_arc); + void process_source_orphan(node *i); + void process_sink_orphan(node *i); + + void test_consistency(node* current_node=NULL); // debug function +}; + + + + + + + + + + + +/////////////////////////////////////// +// Implementation - inline functions // +/////////////////////////////////////// + + + +template <typename captype, typename tcaptype, typename flowtype> + inline typename Graph<captype,tcaptype,flowtype>::node_id Graph<captype,tcaptype,flowtype>::add_node(int num) +{ + assert(num > 0); + + if (node_last + num > node_max) reallocate_nodes(num); + + if (num == 1) + { + node_last -> first = NULL; + node_last -> tr_cap = 0; + node_last -> is_marked = 0; + node_last -> is_in_changed_list = 0; + + node_last ++; + return node_num ++; + } + else + { + memset(node_last, 0, num*sizeof(node)); + + node_id i = node_num; + node_num += num; + node_last += num; + return i; + } +} + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink) +{ + assert(i >= 0 && i < node_num); + + tcaptype delta = nodes[i].tr_cap; + if (delta > 0) cap_source += delta; + else cap_sink -= delta; + flow += (cap_source < cap_sink) ? cap_source : cap_sink; + nodes[i].tr_cap = cap_source - cap_sink; +} + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::add_edge(node_id _i, node_id _j, captype cap, captype rev_cap) +{ + assert(_i >= 0 && _i < node_num); + assert(_j >= 0 && _j < node_num); + assert(_i != _j); + assert(cap >= 0); + assert(rev_cap >= 0); + + if (arc_last == arc_max) reallocate_arcs(); + + arc *a = arc_last ++; + arc *a_rev = arc_last ++; + + node* i = nodes + _i; + node* j = nodes + _j; + + a -> sister = a_rev; + a_rev -> sister = a; + a -> next = i -> first; + i -> first = a; + a_rev -> next = j -> first; + j -> first = a_rev; + a -> head = j; + a_rev -> head = i; + a -> r_cap = cap; + a_rev -> r_cap = rev_cap; +} + +template <typename captype, typename tcaptype, typename flowtype> + inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_first_arc() +{ + return arcs; +} + +template <typename captype, typename tcaptype, typename flowtype> + inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_next_arc(arc* a) +{ + return a + 1; +} + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::get_arc_ends(arc* a, node_id& i, node_id& j) +{ + assert(a >= arcs && a < arc_last); + i = (node_id) (a->sister->head - nodes); + j = (node_id) (a->head - nodes); +} + +template <typename captype, typename tcaptype, typename flowtype> + inline tcaptype Graph<captype,tcaptype,flowtype>::get_trcap(node_id i) +{ + assert(i>=0 && i<node_num); + return nodes[i].tr_cap; +} + +template <typename captype, typename tcaptype, typename flowtype> + inline captype Graph<captype,tcaptype,flowtype>::get_rcap(arc* a) +{ + assert(a >= arcs && a < arc_last); + return a->r_cap; +} + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::set_trcap(node_id i, tcaptype trcap) +{ + assert(i>=0 && i<node_num); + nodes[i].tr_cap = trcap; +} + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::set_rcap(arc* a, captype rcap) +{ + assert(a >= arcs && a < arc_last); + a->r_cap = rcap; +} + + +template <typename captype, typename tcaptype, typename flowtype> + inline typename Graph<captype,tcaptype,flowtype>::termtype Graph<captype,tcaptype,flowtype>::what_segment(node_id i, termtype default_segm) +{ + if (nodes[i].parent) + { + return (nodes[i].is_sink) ? SINK : SOURCE; + } + else + { + return default_segm; + } +} + +template <typename captype, typename tcaptype, typename flowtype> + inline void Graph<captype,tcaptype,flowtype>::mark_node(node_id _i) +{ + node* i = nodes + _i; + if (!i->next) + { + /* it's not in the list yet */ + if (queue_last[1]) queue_last[1] -> next = i; + else queue_first[1] = i; + queue_last[1] = i; + i -> next = i; + } + i->is_marked = 1; +} + + +#endif diff --git a/Chapter07/Chapter07_mcode/GraphCut/instances.inc b/Chapter07/Chapter07_mcode/GraphCut/instances.inc new file mode 100755 index 0000000..0bcbd57 --- /dev/null +++ b/Chapter07/Chapter07_mcode/GraphCut/instances.inc @@ -0,0 +1,16 @@ +#include "graph.h" + +#ifdef _MSC_VER +#pragma warning(disable: 4661) +#endif + +// Instantiations: <captype, tcaptype, flowtype> +// IMPORTANT: +// flowtype should be 'larger' than tcaptype +// tcaptype should be 'larger' than captype + +template class Graph<int,int,int>; +template class Graph<short,int,int>; +template class Graph<float,float,float>; +template class Graph<double,double,double>; + diff --git a/Chapter07/Chapter07_mcode/functions/grid_cut.m b/Chapter07/Chapter07_mcode/functions/grid_cut.m new file mode 100755 index 0000000..cabd67a --- /dev/null +++ b/Chapter07/Chapter07_mcode/functions/grid_cut.m @@ -0,0 +1,252 @@ +function s = grid_cut(s_cost, r_cost, delta_xy, wrap_xy, delta_lu) +%GRID_CUT detects terrain-like surfaces in a volumetric data. +% +% Detected surfaces are terrain-like, i.e. height z is a function of the +% position (x,y). Only on-surface costs, only in-region costs or both type +% of costs may be used. When using in-region costs for multiple surfaces, +% surfaces are non-intersecting and ordered, so that the first surface has +% the smallest z values, and the last surface has the largest z values. +% +% S = GRID_CUT(S_COST,R_COST,DELTA_XY,WRAP_XY,DELTA_LU) +% +% S_COST is a size (X,Y,Z,K) surface cost function. Cost (x,y,z,k) is an +% inverse to likelihood of a voxel (x,y,z) containing surface k. +% K is the number of surfaces and Z is up. When using only in-region +% costs, S_COST should either contain only zeros or be empty. +% R_COST is a size (X,Y,Z,K+1) region cost function. Cost (x,y,z,k) is an +% inverse to likelihood of a voxel (x,y,z) being in a region between +% the surface k and k+1, where the first and the last region are +% boundedend with the volume boundary. When using only on-surface +% costs, R_COST should either contain only zeros or be empty. +% DELTA_XY is a size (K,2) array of stiffness parameters for x and y. If a +% size (1,2) array is given, the same stiffness parameteres are used +% for all surfaces. If size (K,1) array is given, the same stiffness +% is used for x and y. +% WRAP_XY is a length 2 array of boolean wrap options for x and y. If not +% given defaults to false. +% DELTA_LU is a size (K-1,2) array of lower and upper surface overlaps, +% [dl du]. The constraint is: dl <= surf_k - surf_k+1 <= du. When +% using in-region costs, the surfaces are non-intersecting and +% ordered which requires dl > 0 and du >= dl. If a size (1,2) array +% is given, the same overlap parameteres are used for all surfaces. +% If not given, defaults to "no constraint" (dl=-Z, du=Z) when using +% only on-surface costs and to "no overlap" (dl=1, du=Z) otherwise. +% S is a size (X,Y,K) matrix of z coordinates for K segmented surfaces. +% +% Based on: +% "Optimal surface segmentation in volumetric images-a graph-theoretic +% approach." Li, Wu, Chen and Sonka. PAMI 2006. +% "Incorporation of Regional Information in Optimal 3-D Graph Search with +% Application for Intraretinal Layer Segmentation of Optical Coherence +% Tomography Images". Haeker, Wu, Abramoff, Kardon and Sonka. IPMI 2007. +% +% Author: Vedrana Andersen Dahl, vand@dtu.dk, 2013 +% +% NOTE: needs Kolmogorovs implementation of the algorithm from +% "An Experimental Comparison of Min-Cut/Max-Flow Algorithms for +% Energy Minimization in Vision." +% Yuri Boykov and Vladimir Kolmogorov. PAMI 2004 + +% assigning an unique index to each vertex in the first 3D graph (first surface) + +if nargin<5 + delta_lu = []; +end +if nargin<4 + wrap_xy = []; +end +[s_cost, r_cost, delta_xy, wrap_xy, delta_lu] = assign_defaults... + (s_cost, r_cost, delta_xy, wrap_xy, delta_lu); +[X,Y,Z,K] = size(s_cost); +dimension = [X Y Z]; +indices = reshape(1:prod(dimension),dimension); +base = indices(:,:,1); +next_surface = prod(dimension); % adding to indices is a shift to next surface +layers = (0:(K-1))*next_surface; + +wrap = [wrap_xy 0]; % we do not wrap in z direction + +% EDGES WHICH ARE THE SAME FOR ALL SURFACES: +% intracolumn arcs pointing down, Equation (2) +Ea1 = displacement_to_edges(indices, [0 0 -1], wrap); % first surface +Ea = repmat(Ea1,[K 1]) + kron(layers(:),ones(size(Ea1))); % all surfaces +Ea = [Ea,ones(size(Ea,1),1)*[inf,0]]; % assigning [inf 0] weight +% base edges, part of intercolumn arcs, Equation (3) +erxb1 = displacement_to_edges(base, [1 0 0], wrap); +eryb1 = displacement_to_edges(base, [0 1 0], wrap); +Erb1 = [erxb1; eryb1]; % first surface +Erb = repmat(Erb1,[K 1]) + kron(layers(:),ones(size(Erb1))); % all surfaces +% assigning [inf inf] weight, base edges are in both directions +Erb = [Erb,ones(size(Erb,1),1)*[inf,inf]]; + +% EDGES WHICH DEPEND ON SURFACE STIFFNESS: +% slanted edges, part of intercolumn arcs, Equation (3) +Erpm = []; % preallocation length could be computed from X,Y,Z,delta and wrap +for k = 1:K + erxp = displacement_to_edges(indices, [1 0 -delta_xy(k,1)], wrap); + erxm = displacement_to_edges(indices, [-1 0 -delta_xy(k,1)], wrap); + eryp = displacement_to_edges(indices, [0 1 -delta_xy(k,2)], wrap); + erym = displacement_to_edges(indices, [0 -1 -delta_xy(k,2)], wrap); + erpm = [erxp; erxm; eryp; erym] + layers(k); + Erpm = [Erpm;erpm]; +end +% assigning [inf 0] weight, slanted edges are in one direction +Erpm = [Erpm,ones(size(Erpm,1),1)*[inf,0]]; +Er = [Erpm; Erb]; % all inter edges, for all surfaces +E = [Ea;Er]; % all intracolumn and intercolumn edges, for all surfaces + +% EDGES WHICH DEPEND ON DISTANCE BETWEEN SURFACES: +% intersurface arc, Equation (4) +if K>1 % only if we have more than 1 surface + Es = []; % preallocation length could be computed + for k = 1:K-1 + esl = displacement_to_edges(indices, [0 0 delta_lu(k,1)], wrap); + esl = esl + ones(size(esl,1),1)*layers([k,k+1]); + esu = displacement_to_edges(indices, [0 0 -delta_lu(k,2)], wrap); + esu = esu + ones(size(esu,1),1)*layers([k+1,k]); + es = [esl;esu]; + Es = [Es;es]; + end + Es = [Es,ones(size(Es,1),1)*[inf,0]]; % assigning [inf 0] weight + + % intersurface base edges + esb = [layers(1:end-1)',layers(2:end)']+1; %first vertex in all surfaces + Esb = [esb,ones(K-1,1)*[inf,inf]]; % assigning [inf inf] + E = [E;Es;Esb]; % all intracolumn, intercolumn and intersurface edges, for all surfaces +end + +% EDGES WHICH DEPEND ON SURFACE LIKELIHOOD: +% up to here we do not use s_cost -- consider efficient algorithm if solving +% multiple problems of the same size +% vertex s_cost, Equation (1), done simulatniously for all surfaces +w_on = -1*ones([dimension,K]); % to prevent empty solution, see second half of section 4.1 +w_on(:,:,2:end,:) = double(s_cost(:,:,2:end,:))-double(s_cost(:,:,1:end-1,:)); +% In case of layered surfaces vertices which can't be realized should be removed, +% i.e. topmost vertices of lower surface and lowest vertices of higher surface. +% Instead of removing vertices, I assign inf weight to topmost vertices. +% And this only when delta_lu are both positive (both negative can be +% avoided by proper ordering of surfaces). +if K>1 % only if we have more than 1 surface + for k=1:K-1 % assigning inf on-surface weights instead of removing vertices + if delta_lu(k,1)>0 && delta_lu(k,1)<= delta_lu(k,2) % bit clumsy + w_on(:,:,end-delta_lu(k,1)+1:end,k) = inf; + end + end +end + +% In-region cost, converting to double to avoid problems when using images +w_in = double(r_cost(:,:,:,1:end-1))-double(r_cost(:,:,:,2:end)); +% There are issues concerning in-region cost e.g. a better way of preventing +% empty solutions or topmost solutions. +w_in(:,:,1,:) = -inf; % preventing empty solution + +w = w_on + w_in; + +Vp_ix = find(w(:)>=0); % positive vertices, to be connected to sink +Vm_ix = find(w(:)<0); % negative vertices, to be connected to source +Es = [Vm_ix, -w(Vm_ix), zeros(length(Vm_ix),1)]; % source edges +Et = [Vp_ix, zeros(length(Vp_ix),1), w(Vp_ix)]; % sink edges +Est = [Es;Et]; % all terminal edges + +% FINDING GRAPH CUT USING MAGIC +[Scut, ~] = GraphCutMex(prod(dimension)*K,Est,E); + +% retreiving surfaces as the upper envelope of Scut +S = zeros(X,Y,Z,K); +S(Scut) = 1; +s = zeros(X,Y,K); +for ki = 1:K + for x = 1:X + for y = 1:Y + s(x,y,ki) = find(S(x,y,:,ki),1,'last'); + end + end +end +end + +% NOTE: when some of the intersurface edges point upwards, i.e. when +% delta_l is positive, the lowest delta_l rows of a graph for higher +% surface can not appear in any feasible solution and can be removed. +% This has NOT been implemented here. + +function edges = displacement_to_edges(indices, disp, wrap) +% indices -- 3D volume of indices +% disp -- length 3 vector of displacements, given as [x y,z] +% wrapping -- length 3 vector of boolean options, given as [wrapx wrapy wrapz] +% edges -- two columns of indices, each line is an edge. + +[x_from,x_into] = displace_1D_indices(size(indices,1),disp(1),wrap(1)); +[y_from,y_into] = displace_1D_indices(size(indices,2),disp(2),wrap(2)); +[z_from,z_into] = displace_1D_indices(size(indices,3),disp(3),wrap(3)); + +indices_from = indices(x_from,y_from,z_from); +indices_into = indices(x_into,y_into,z_into); + +edges = [indices_from(:),indices_into(:)]; +end + +function [from,into] = displace_1D_indices(dim,disp,wrap) +% dim -- length of the indices vector +% disp -- length of the displacement +% wrap -- boolean, indicating trimming or wrapping + +indices = 1:dim; +if wrap + from = indices([1-min(disp,0):end,1:-min(disp,0)]); + into = indices([1+max(disp,0):end,1:max(disp,0)]); +else + from = indices(1-min(disp,0):end-max(disp,0)); + into = indices(1+max(disp,0):end+min(disp,0)); +end +end + +function [s_cost, r_cost, delta_xy, wrap_xy, delta_lu] = assign_defaults... + (s_cost, r_cost, delta_xy, wrap_xy, delta_lu) +% either cost_s or cost_r need to be given (not empty and not all zeros) +if isempty(s_cost) + s_cost = zeros(size(r_cost)-[0,0,0,1]); +end +if isempty(r_cost) + r_cost = zeros([size(s_cost,1),size(s_cost,2),size(s_cost,3),... + size(s_cost,4)]+[0,0,0,1]); % to allow s_cost to be 3D + regional = false; +elseif all(r_cost(:)==0) + regional = false; +else + regional = true; +end + +[X,Y,Z,K] = size(s_cost); +[Xr,Yr,Zr,Krplus1] = size(r_cost); +if ~X==Xr || ~Y==Yr || ~Z==Zr || ~(K+1)==Krplus1 + error('Error using grid_cut. Dimensions of s_cost and r_cost must agree.') +end + +% smoothness constraint may be given once for all surfaces +if size(delta_xy,1)==1 + delta_xy = ones(K,1)*delta_xy; +end +% smoothness constraint may be given once for both directions +if size(delta_xy,2)==1 + delta_xy = delta_xy*ones(1,2); +end + +if isempty(wrap_xy) || (numel(wrap_xy)==1 && wrap_xy==0) + wrap_xy = [0 0]; +end + +% overlap constraint defaults to no overlap in region case +% and no constraint in surface case +if isempty(delta_lu) + if regional % no overlap + delta_lu = [1 Z]; + else % no constraint + delta_lu = [-Z Z]; + end +end + +% overlap constraint may be given once for all surface pairs +if all(size(delta_lu)==[1,2]) + delta_lu = ones(K-1,1)*delta_lu; +end +end \ No newline at end of file diff --git a/Chapter07/Chapter07_mcode/in_region_cost_example.m b/Chapter07/Chapter07_mcode/in_region_cost_example.m new file mode 100644 index 0000000..d5f725a --- /dev/null +++ b/Chapter07/Chapter07_mcode/in_region_cost_example.m @@ -0,0 +1,25 @@ +clear, close all +addpath GraphCut functions +I = imread('../../../../Data/week7/peaks_image.png'); +figure, +subplot(1,2,1) +imagesc(I), colormap gray, axis image + +% COST FUNCTION +% The region in the middle is bright compared to two darker regions. +P = permute(I, [2,3,1]); % making sure that up is the third dimension +region_cost = cat(4, P, (255-P), P); + +% GEOMETRIC CONSTRAINS +delta_xy = 1; % smoothness very constrained, try also 3 to see less smoothness +wrap_xy = 0; % 0 for terren-like surfaces +delta_ul = [1 size(I,1)]; % can be very close, but may not overlap + +% CUT +s = grid_cut([],region_cost,delta_xy,wrap_xy,delta_ul); + +% VISUALIZATION +subplot(1,2,2) +imagesc(I), axis image ij, colormap gray, hold on +plot(1:size(I,2),permute(s,[1,3,2]),'r', 'LineWidth',2) + diff --git a/Chapter07/Chapter07_mcode/on_surface_cost_example.m b/Chapter07/Chapter07_mcode/on_surface_cost_example.m new file mode 100755 index 0000000..fe8e858 --- /dev/null +++ b/Chapter07/Chapter07_mcode/on_surface_cost_example.m @@ -0,0 +1,44 @@ +clear, close all +addpath GraphCut functions + +lw = 2; + +I = imread('../../../../Data/week7/layers_A.png'); +surface_cost = permute(I,[2,3,1]); +dim = size(I); + +%% input +figure +subplot(1,4,1) +imagesc(I), axis image, colormap gray +title('input image'), drawnow + +%% one line +delta = 3; +wrap = false; +s = grid_cut(surface_cost, [], delta, wrap); +subplot(1,4,2) +imagesc(I), axis image, colormap gray, hold on +plot(1:dim(2), s, 'r', 'LineWidth', lw) +title(['delta = ',num2str(delta)]), drawnow + +%% a smoother line +delta = 1; +wrap = false; +s = grid_cut(surface_cost, [], delta, wrap); +subplot(1,4,3) +imagesc(I), axis image, colormap gray, hold on +plot(1:dim(2), s, 'r', 'LineWidth', lw) +title(['delta = ',num2str(delta)]), drawnow + +%% two lines +costs = cat(4,surface_cost,surface_cost); +delta = 3; +wrap = false; +overlap = [15, size(I,1)]; + +s = grid_cut(costs, [], delta, wrap, overlap); +subplot(1,4,4) +imagesc(I), axis image, colormap gray, hold on +plot(1:dim(2),permute(s,[1,3,2]),'r', 'LineWidth',lw) +title('two dark lines'), drawnow \ No newline at end of file diff --git a/Chapter07/Chapter07_mcode/w1_unwrapping_image_exercise.m b/Chapter07/Chapter07_mcode/w1_unwrapping_image_exercise.m new file mode 100644 index 0000000..9692cd5 --- /dev/null +++ b/Chapter07/Chapter07_mcode/w1_unwrapping_image_exercise.m @@ -0,0 +1,18 @@ +% optional exercise 1.1.5 +I = imread('../../../../Data/week1/dental/slice100.png'); + +a = 180; % number of angles for unfolding +angles = (0 : (a-1)) *2*pi/a; % angular coordinate + +center = (1+size(I))/2; +r = 1:min(size(I)/2); % radial coordinate for unwrapping + +X = center(1) + r'*cos(angles); +Y = center(2) + r'*sin(angles); + +F = griddedInterpolant(double(I)); +U = F(X,Y); + +figure +subplot(121), imagesc(I), axis image, colormap gray +subplot(122), imagesc(U), axis image, colormap gray \ No newline at end of file diff --git a/Chapter07/in_region_cost_example.py b/Chapter07/in_region_cost_example.py new file mode 100755 index 0000000..c1ac7ce --- /dev/null +++ b/Chapter07/in_region_cost_example.py @@ -0,0 +1,38 @@ +import numpy as np +import matplotlib.pyplot as plt +import skimage.io +import slgbuilder + + +I = skimage.io.imread('../../../../Data/week7/peaks_image.png').astype(np.int32) + +fig, ax = plt.subplots(1,2) +ax[0].imshow(I, cmap='gray') + + +layers = [slgbuilder.GraphObject(0*I), slgbuilder.GraphObject(0*I)] # no on-surface cost +helper = slgbuilder.MaxflowBuilder() +helper.add_objects(layers) + +# Addin regional costs, +# the region in the middle is bright compared to two darker regions. +helper.add_layered_region_cost(layers[0], I, 255-I) +helper.add_layered_region_cost(layers[1], 255-I, I) + +# Adding geometric constrains +helper.add_layered_boundary_cost() +helper.add_layered_smoothness(delta=1, wrap=False) +helper.add_layered_containment(layers[0], layers[1], min_margin=1) + +# Cut +helper.solve() +segmentations = [helper.what_segments(l).astype(np.int32) for l in layers] +segmentation_lines = [s.shape[0] - np.argmax(s[::-1,:], axis=0) - 1 for s in segmentations] + +# Visualization +ax[1].imshow(I, cmap='gray') +for line in segmentation_lines: + ax[1].plot(line, 'r') + + + diff --git a/Chapter07/on_surface_cost_example.py b/Chapter07/on_surface_cost_example.py new file mode 100755 index 0000000..eeab4b9 --- /dev/null +++ b/Chapter07/on_surface_cost_example.py @@ -0,0 +1,66 @@ +import numpy as np +import matplotlib.pyplot as plt +import skimage.io +import slgbuilder + +#%% input +I = skimage.io.imread('../../../../Data/week7/layers_A.png').astype(np.int32) + +fig, ax = plt.subplots(1,4) +ax[0].imshow(I, cmap='gray') +ax[0].set_title('input image') + +#%% one line +delta = 3 + +layer = slgbuilder.GraphObject(I) +helper = slgbuilder.MaxflowBuilder() +helper.add_object(layer) +helper.add_layered_boundary_cost() +helper.add_layered_smoothness(delta=delta, wrap=False) + +helper.solve() +segmentation = helper.what_segments(layer) +segmentation_line = segmentation.shape[0] - np.argmax(segmentation[::-1,:], axis=0) - 1 + +ax[1].imshow(I, cmap='gray') +ax[1].plot(segmentation_line, 'r') +ax[1].set_title(f'delta = {delta}') + + +#%% a smoother line +delta = 1 + +layer = slgbuilder.GraphObject(I) +helper = slgbuilder.MaxflowBuilder() +helper.add_object(layer) +helper.add_layered_boundary_cost() +helper.add_layered_smoothness(delta=delta, wrap=False) + +helper.solve() +segmentation = helper.what_segments(layer) +segmentation_line = segmentation.shape[0] - np.argmax(segmentation[::-1,:], axis=0) - 1 + +ax[2].imshow(I, cmap='gray') +ax[2].plot(segmentation_line, 'r') +ax[2].set_title(f'delta = {delta}') + + +#%% tow lines +layers = [slgbuilder.GraphObject(I), slgbuilder.GraphObject(I)] +delta = 3 + +helper = slgbuilder.MaxflowBuilder() +helper.add_objects(layers) +helper.add_layered_boundary_cost() +helper.add_layered_smoothness(delta=delta, wrap=False) +helper.add_layered_containment(layers[0], layers[1], min_margin=15) + +helper.solve() +segmentations = [helper.what_segments(l).astype(np.int32) for l in layers] +segmentation_lines = [s.shape[0] - np.argmax(s[::-1,:], axis=0) - 1 for s in segmentations] + +ax[3].imshow(I, cmap='gray') +for line in segmentation_lines: + ax[3].plot(line, 'r') +ax[3].set_title('two dark lines') diff --git a/Chapter07/w1_unwrapping_image_exercise.py b/Chapter07/w1_unwrapping_image_exercise.py new file mode 100644 index 0000000..c67fffe --- /dev/null +++ b/Chapter07/w1_unwrapping_image_exercise.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 27 09:14:57 2020 + +@author: vand +""" + +# optional exercise 1.1.5 + +import skimage.io +import numpy as np +import scipy.interpolate +import matplotlib.pyplot as plt + +I = skimage.io.imread('../../../../Data/week1/dental/slice100.png') + + +a = 180 # number of angles for unfolding +angles = np.arange(a)*2*np.pi/a # angular coordinate + +center = (np.array(I.shape)-1)/2 +r = int(min(I.shape)/2) +radii = np.arange(r) + 1 #radial coordinate for unwrapping + +X = center[0] + np.outer(radii,np.cos(angles)) +Y = center[1] + np.outer(radii,np.sin(angles)) + +F = scipy.interpolate.interp2d(np.arange(I.shape[0]), np.arange(I.shape[1]), I) +U = np.array([F(p[0],p[1]) for p in np.c_[Y.ravel(),X.ravel()]]) +U = U.reshape((r,a)).astype(np.uint8) + +fig, ax = plt.subplots(1,2) +ax[0].imshow(I, cmap='gray') +ax[1].imshow(U, cmap='gray') -- GitLab