2 * Translated from Jeff Siskind
's Scheme code by Stephen Weeks
4 * Here is the description from Jeff
Siskind (qobi@research
.nj
.nec
.com
)
6 * It is an implementation
of Ratio
7 * Regions
, an image segmentation
/contour finding technique due to Ingemar Cox
,
8 * Satish Rao
, and Yu Zhong
. The algorithm is a reduction to max flow
, an
9 * unpublished technique that Satish described to me
. Peter Blicher originally
10 * implemented this via a translation to Andrew Goldberg
's generic max
-flow code
.
11 * I
've reimplemented it
, specializing the max
-flow algorithm to the particular
12 * graphs that are produced by Satish
's reduction instead
of using Andrew
's code
.
13 * The max
-flow algorithm is preflow
-push
with periodic relabeling
and a
14 * wave
-based heuristic for scheduling pushes
and lifts due to Sebastien Roy
.
21 fun doo(max
: int, f
: int -> unit
): unit
=
22 let fun loop i
= if i
>= max
then () else (f i
; loop(i
+ 1))
28 val make_vector
= Array
.array
29 val vector_length
= Array
.length
30 val vector_ref
= Array
.sub
31 val vector_set
= Array
.update
32 val map_n_vector
= Array
.tabulate
33 val string_length
= String.size
34 val string_ref
= String.sub
35 fun write_char c
= () (* TextIO.output1(TextIO.stdOut
, c
) *)
37 val quotient
= Int.quot
38 fun for_each(l
, f
) = List.app f l
39 fun negative x
= x
< 0
40 fun positive x
= x
> 0
48 | x
:: l
=> loop(l
, Int.min(min
, x
))
51 | _
=> raise Fail
"min"
54 let fun loop i
= i
>= n
orelse (p i
andalso loop(i
+ 1))
58 fun some(l
, p
) = List.exists p l
61 let fun loop i
= i
< n
andalso (p i
orelse loop(i
+ 1))
65 fun some_vector(v
, p
) =
69 andalso (p(vector_ref(v
, i
))
77 datatype 'a matrix
= Matrix
of 'a array array
79 fun make_matrix(m
: int, n
: int, a
: 'a
): 'a matrix
=
80 Matrix(map_n_vector(m
, fn i
=> make_vector(n
, a
)))
82 fun matrix_rows(Matrix a
) = vector_length a
83 fun matrix_columns(Matrix a
) = vector_length(vector_ref(a
, 0))
84 fun matrix_ref(Matrix a
, i
, j
) = vector_ref(vector_ref(a
, i
), j
)
85 fun matrix_set(Matrix a
, i
, j
, x
) = vector_set(vector_ref(a
, i
), j
, x
)
87 datatype pormatValue
=
91 fun pormat(control_string
: string, values
: pormatValue list
): unit
=
93 fun loop(i
: int, values
: pormatValue list
): unit
=
94 if not(i
= string_length control_string
)
96 let val c
= string_ref(control_string
, i
)
98 then let val c2
= string_ref(control_string
, i
+ 1)
99 in case (c2
, values
) of
100 (#
"s", Int n
:: values
) =>
101 (print(Int.toString n
) ; loop(i
+ 2, values
))
102 |
(#
"a", String s
:: values
) =>
103 (print s
; loop(i
+ 2, values
))
105 (print
"\n"; loop(i
+ 2, values
))
106 | _
=> (write_char c
; loop(i
+ 1, values
))
108 else (write_char c
; loop(i
+ 1, values
))
114 (* The vertices are s
, t
, and (y
,x
).
115 * C_RIGHT
[y
,x
] is the capacity
from (y
,x
) to (y
,x
+1) which is the same
as the
116 * capacity
from (y
,x
+1) to (y
,x
).
117 * C_DOWN
[y
,x
] is the capacity
from (y
,x
) to (y
+1,x
) which is the same
as the
118 * capacity
from (y
+1,x
) to (y
,x
).
119 * The capacity from s
to (y
,0), (0,x
), (y
,Y_1
), (0,X_1
) is implicitly
121 * The capacity
from (x
,y
) to t is V
*W
[y
,x
].
122 * F_RIGHT
[y
,x
] is the preflow
from (y
,x
) to (y
,x
+1) which is the negation
of
123 * the preflow
from (y
,x
+1) to (y
,x
).
124 * F_DOWN
[y
,x
] is the preflow
from (y
,x
) to (y
+1,x
) which is the negation
of
125 * the preflow
from (y
+1,x
) to (y
,x
).
126 * We
do not record the preflow from s
to (y
,X_1
), (y
,0), (Y_1
,x
), and (0,x
)
127 * and from (y
,X_1
), (y
,0), (Y_1
,x
), and (0,x
) to s
.
128 * F_T
[y
,x
] is the preflow
from (y
,x
) to t
.
129 * We
do not record the preflow from t
to (y
,x
).
130 * {C
,F
}_RIGHT
[0:Y_1
,0:X_2
].
131 * {C
,F
}_DOWN
[0:Y_2
,0:X_1
].
133 * For now
, we will keep all
capacities (and thus all preflows
) as integers
.
134 * (CF_RIGHT y x
) is the residual capacity
from (y
,x
) to (y
,x
+1).
135 * (CF_LEFT y x
) is the residual capacity
from (y
,x
) to (y
,x_1
).
136 * (CF_DOWN y x
) is the residual capacity
from (y
,x
) to (y
+1,x
).
137 * (CF_UP y x
) is the residual capacity
from (y
,x
) to (y_1
,x
).
138 * We
do not compute the residual capacities from s
to (y
,X_1
), (y
,0),
139 * (Y_1
,x
), and (0,x
) because they are all infinite
.
140 * We
do not compute the residual capacities
from (y
,X_1
), (y
,0), (Y_1
,x
),
141 * and (0,x
) to s because they will never be used
.
142 * (CF_T y x
) is the residual capacity
from (y
,x
) to t
.
143 * We
do not compute the residual capacity from t
to (y
,x
) because it will
145 * (EF_RIGHT? y x
) is
true if there is an edge
from (y
,x
) to (y
,x
+1) in the
147 * (EF_LEFT? y x
) is
true if there is an edge
from (y
,x
) to (y
,x_1
) in the
149 * (EF_DOWN? y x
) is
true if there is an edge
from (y
,x
) to (y
+1,x
) in the
151 * (EF_UP? y x
) is
true if there is an edge
from (y
,x
) to (y_1
,x
) in the
153 * (EF_T? y x
) is
true if there is an edge
from (y
,x
) to t
in the
155 * There are always edges
in the residual network from s
to (y
,X_1
), (y
,0),
156 * (Y_1
,x
), and (0,x
).
157 * We don
't care whether there are edges
in the residual network from
158 * (y
,X_1
), (y
,0), (Y_1
,x
), and (0,x
) to s because they will never be used
.
159 * We don
't care whether there are edges
in the residual network from t to
160 * (y
,x
) because they will never be used
.
163 fun positive_min(x
, y
) = if negative x
then y
else Int.min(x
, y
)
165 fun positive_minus(x
, y
) = if negative x
then x
else x
- y
167 fun positive_plus(x
, y
) = if negative x
then x
else x
+ y
169 fun rao_ratio_region(c_right
, c_down
, w
, lg_max_v
) =
170 let val height
= matrix_rows w
171 val width
= matrix_columns w
172 val f_right
= make_matrix(height
, width
- 1, 0)
173 val f_down
= make_matrix(height
- 1, width
, 0)
174 val f_t
= make_matrix(height
, width
, 0)
175 val h
= make_matrix(height
, width
, 0)
176 val e
= make_matrix(height
, width
, 0)
177 val marked
= make_matrix(height
, width
, false)
178 val m1
= height
* width
+ 2
179 val m2
= 2 * height
* width
+ 2
180 val q
= make_vector(2 * height
* width
+ 3, [])
182 matrix_ref(c_right
, y
, x
) - matrix_ref(f_right
, y
, x
)
184 matrix_ref(c_right
, y
, x
- 1) + matrix_ref(f_right
, y
, x
- 1)
186 matrix_ref(c_down
, y
, x
) - matrix_ref(f_down
, y
, x
)
188 matrix_ref(c_down
, y
- 1, x
) + matrix_ref(f_down
, y
- 1, x
)
189 fun ef_right(y
, x
) = positive(cf_right(y
, x
))
190 fun ef_left(y
, x
) = positive(cf_left(y
, x
))
191 fun ef_down(y
, x
) = positive(cf_down(y
, x
))
192 fun ef_up(y
, x
) = positive(cf_up(y
, x
))
196 if not(matrix_ref(marked
, y
, x
))
201 vector_ref(q
, matrix_ref(h
, y
, x
)))))
202 ; matrix_set(marked
, y
, x
, true))
204 fun cf_t(y
, x
) = v
* matrix_ref(w
, y
, x
) - matrix_ref(f_t
, y
, x
)
205 fun ef_t(y
, x
) = positive(cf_t(y
, x
))
206 fun can_push_right(y
, x
) =
208 andalso not(zero(matrix_ref(e
, y
, x
)))
209 andalso ef_right(y
, x
)
210 andalso matrix_ref(h
, y
, x
) = matrix_ref(h
, y
, x
+ 1) + 1
211 fun can_push_left(y
, x
) =
213 andalso not(zero(matrix_ref(e
, y
, x
)))
214 andalso ef_left(y
, x
)
215 andalso matrix_ref(h
, y
, x
) = matrix_ref(h
, y
, x
- 1) + 1
216 fun can_push_down(y
, x
) =
218 andalso not(zero(matrix_ref(e
, y
, x
)))
219 andalso ef_down(y
, x
)
220 andalso matrix_ref(h
, y
, x
) = matrix_ref(h
, y
+ 1, x
) + 1
221 fun can_push_up(y
, x
) =
223 andalso not(zero(matrix_ref(e
, y
, x
)))
225 andalso matrix_ref(h
, y
, x
) = matrix_ref(h
, y
- 1, x
) + 1
226 fun can_push_t(y
, x
) =
227 not(zero(matrix_ref(e
, y
, x
)))
229 andalso matrix_ref(h
, y
, x
) = 1
231 not(zero(matrix_ref(e
, y
, x
)))
232 andalso (if x
= width
- 1
233 then matrix_ref(h
, y
, x
) <= m1
234 else (not(ef_right(y
, x
))
236 matrix_ref(h
, y
, x
) <= matrix_ref(h
, y
, x
+ 1)))
238 then matrix_ref(h
, y
, x
) <= m1
239 else (not(ef_left(y
, x
))
241 matrix_ref(h
, y
, x
) <= matrix_ref(h
, y
, x
- 1)))
242 andalso (if y
= height
- 1
243 then matrix_ref(h
, y
, x
) <= m1
244 else (not(ef_down(y
, x
))
246 matrix_ref(h
, y
, x
) <= matrix_ref(h
, y
+ 1, x
)))
248 then matrix_ref(h
, y
, x
) <= m1
249 else (not(ef_up(y
, x
))
251 matrix_ref(h
, y
, x
) <= matrix_ref(h
, y
- 1, x
)))
252 andalso (not(ef_t(y
, x
)) orelse matrix_ref(h
, y
, x
) = 0)
253 fun push_right(y
, x
) =
254 (* (pormat
"Push right ~s ~s~%" y x
) *)
255 let val df_u_v
= positive_min(matrix_ref(e
, y
, x
), cf_right(y
, x
))
256 in matrix_set(f_right
, y
, x
, matrix_ref(f_right
, y
, x
) + df_u_v
)
257 ; matrix_set(e
, y
, x
,
258 positive_minus(matrix_ref(e
, y
, x
), df_u_v
))
259 ; matrix_set(e
, y
, x
+ 1,
260 positive_plus(matrix_ref(e
, y
, x
+ 1), df_u_v
))
263 fun push_left(y
, x
) =
264 (* (pormat
"Push left ~s ~s~%" y x
) *)
265 let val df_u_v
= positive_min(matrix_ref(e
, y
, x
), cf_left(y
, x
))
266 in matrix_set(f_right
, y
, x
- 1,
267 matrix_ref(f_right
, y
, x
- 1) - df_u_v
)
268 ; matrix_set(e
, y
, x
,
269 positive_minus(matrix_ref(e
, y
, x
), df_u_v
))
270 ; matrix_set(e
, y
, x
- 1,
271 positive_plus(matrix_ref(e
, y
, x
- 1), df_u_v
))
275 fun push_down(y
, x
) =
276 (* (pormat
"Push down ~s ~s~%" y x
) *)
277 let val df_u_v
= positive_min(matrix_ref(e
, y
, x
), cf_down(y
, x
))
278 in matrix_set(f_down
, y
, x
, matrix_ref(f_down
, y
, x
) + df_u_v
)
279 ; matrix_set(e
, y
, x
,
280 positive_minus(matrix_ref(e
, y
, x
), df_u_v
))
281 ; matrix_set(e
, y
+ 1, x
,
282 positive_plus(matrix_ref(e
, y
+ 1, x
), df_u_v
))
286 (* ;;(pormat
"Push up ~s ~s~%" y x
) *)
287 let val df_u_v
= positive_min(matrix_ref(e
, y
, x
), cf_up(y
, x
))
288 in matrix_set(f_down
, y
- 1, x
,
289 matrix_ref(f_down
, y
- 1, x
) - df_u_v
)
290 ; matrix_set(e
, y
, x
,
291 positive_minus(matrix_ref(e
, y
, x
), df_u_v
))
292 ; matrix_set(e
, y
- 1, x
,
293 positive_plus(matrix_ref(e
, y
- 1, x
), df_u_v
))
297 (* ;;(pormat
"Push t ~s ~s~%" y x
) *)
298 let val df_u_v
= positive_min(matrix_ref(e
, y
, x
), cf_t(y
, x
))
299 in matrix_set(f_t
, y
, x
, matrix_ref(f_t
, y
, x
) + df_u_v
)
300 ; matrix_set(e
, y
, x
,
301 positive_minus(matrix_ref(e
, y
, x
), df_u_v
))
304 (* ;;(pormat
"Lift ~s ~s~%" y x
) *)
307 1 + min
[if x
= width
- 1
309 else if ef_right(y
, x
)
310 then matrix_ref(h
, y
, x
+ 1)
314 else if ef_left(y
, x
)
315 then matrix_ref(h
, y
, x
- 1)
319 else if ef_down(y
, x
)
320 then matrix_ref(h
, y
+ 1, x
)
325 then matrix_ref(h
, y
- 1, x
)
327 if ef_t(y
, x
) then 0 else m2
])
329 (* ;;(pormat
"Relabel~%") *)
333 | Cons
of 'a
* 'a queue ref
334 fun null(q
: 'q queue ref
) =
338 val q
: (int * int) queue ref
= ref Nil
339 val tail
: (int * int) queue ref
= ref Nil
340 fun enqueue(y
, x
, value
) =
341 if value
< matrix_ref(h
, y
, x
)
342 then (matrix_set(h
, y
, x
, value
)
343 ; if not(matrix_ref(marked
, y
, x
))
344 then (matrix_set(marked
, y
, x
, true)
347 (tail
:= Cons((x
, y
), ref Nil
)
350 (cdr
:= Cons((x
, y
), ref Nil
)
356 Nil
=> raise Fail
"dequeue"
358 (matrix_set(marked
, y p
, x p
, false)
360 ; if null q
then tail
:= Nil
else ()
362 in doo(height
, fn y
=>
364 (matrix_set(h
, y
, x
, m1
)
365 ; matrix_set(marked
, y
, x
, false))))
366 ; doo(height
, fn y
=>
369 andalso matrix_ref(h
, y
, x
) > 1
370 then enqueue(y
, x
, 1)
376 (let val p
= dequeue()
379 val value
= matrix_ref(h
, y
, x
) + 1
380 in if x
> 0 andalso ef_right(y
, x
- 1)
381 then enqueue(y
, x
- 1, value
)
383 ; if x
< width
- 1 andalso ef_left(y
, x
+ 1)
384 then enqueue(y
, x
+ 1, value
)
386 ; if y
> 0 andalso ef_down(y
- 1, x
)
387 then enqueue(y
- 1, x
, value
)
389 ; if y
< height
- 1 andalso ef_up(y
+ 1, x
)
390 then enqueue(y
+ 1, x
, value
)
398 in doo(height
, fn y
=>
400 (matrix_set(e
, y
, x
, 0)
401 ; matrix_set(f_t
, y
, x
, 0))))
402 ; doo(height
, fn y
=>
403 doo(width
- 1, fn x
=>
404 matrix_set(f_right
, y
, x
, 0)))
405 ; doo(height
- 1, fn y
=>
407 matrix_set(f_down
, y
, x
, 0)))
408 ; doo(height
, fn y
=>
409 (matrix_set(e
, y
, width
- 1, ~
1)
410 ; matrix_set(e
, y
, 0, ~
1)))
411 ; doo(width
- 1, fn x
=>
412 (matrix_set(e
, height
- 1, x
, ~
1)
413 ; matrix_set(e
, 0, x
, ~
1)))
414 ; let val pushes
= ref
0
418 if zero(modulo(i
, 6)) andalso not p
420 ; relabels
:= !relabels
+ 1
421 ; if every_n(height
, fn y
=>
422 every_n(width
, fn x
=>
423 zero(matrix_ref(e
, y
, x
))
425 matrix_ref(h
, y
, x
) = m1
))
427 (* Every vertex
with excess capacity is not reachable from the sink
in
428 * the inverse residual network
. So terminate early because we have
429 * already found a min cut
. In this
case, the preflows
and excess
430 * capacities will not be correct
. But the cut is indicated by the
431 * heights
. Vertices reachable from the source have height
432 * HEIGHT
* WIDTH
+ 2 while vertices reachable from the sink have
433 * smaller height
. Early termination is necessary
with relabeling to
434 * prevent an infinite loop
. The loop arises because vertices that are
435 * not reachable from the sink
in the inverse residual network have
436 * their height reset to HEIGHT
* WIDTH
+ 2 by the relabeling
437 * process
. If there are such vertices
with excess capacity
, this is
438 * not high enough for the excess capacity to be pushed back to the
439 * perimeter
. So after relabeling
, vertices get lifted to try to push
440 * excess capacity back to the perimeter but
then a relabeling happens
441 * to soon
and foils this lifting
. Terminating when all vertices
with
442 * excess capacity are not reachable from the sink
in the inverse
443 * residual network eliminates this problem
.
446 ("~s push~a, ~s lift~a, ~s relabel~a, ~s wave~a, terminated early~%",
448 String(if !pushes
= 1 then "" else "es"),
450 String(if !lifts
= 1 then "" else "s"),
452 String(if !relabels
= 1 then "" else "s"),
454 String(if i
= 1 then "" else "s")]))
456 (* We need to rebuild the priority queue after relabeling since the
457 * heights might have changed
and the priority queue is indexed by
458 * height
. This also assumes that a relabel is done
before any pushes
461 (doo(vector_length q
, fn k
=>
462 vector_set(q
, k
, []))
463 ; doo(height
, fn y
=>
465 matrix_set(marked
, y
, x
, false)))
466 ; doo(height
, fn y
=>
468 if not(zero(matrix_ref(e
, y
, x
)))
472 else if some_vector(q
, fn ps
=>
476 in can_push_right(y
, x
)
477 orelse can_push_left(y
, x
)
478 orelse can_push_down(y
, x
)
479 orelse can_push_up(y
, x
)
480 orelse can_push_t(y
, x
)
481 orelse can_lift(y
, x
)
489 let val ps
= vector_ref(q
, k
)
490 in vector_set(q
, k
, [])
493 matrix_set(marked
, y p
, x p
,
499 in if can_push_right(y
, x
)
500 then (pushes
:= !pushes
+ 1
503 ; if can_push_left(y
, x
)
504 then (pushes
:= !pushes
+ 1
507 ; if can_push_down(y
, x
)
508 then (pushes
:= !pushes
+ 1
511 ; if can_push_up(y
, x
)
512 then (pushes
:= !pushes
+ 1
515 ; if can_push_t(y
, x
)
516 then (pushes
:= !pushes
+ 1
520 then (lifts
:= !lifts
+ 1
523 ; if not(zero(matrix_ref(e
, y
, x
)))
530 in loop(vector_length q
- 1)
532 ; loop(i
+ 1, false))
534 (* This is so MIN_CUT
and MIN_CUT_INCLUDES_EVERY_EDGE_TO_T work
. *)
536 ; relabels
:= !relabels
+ 1
537 ; (pormat("~s push~a, ~s lift~a, ~s relabel~a, ~s wave~a~%",
539 String(if !pushes
= 1 then "" else "es"),
541 String(if !lifts
= 1 then "" else "s"),
543 String(if !relabels
= 1 then "" else "s"),
545 String(if i
= 1 then "" else "s")])))
549 fun min_cut_includes_every_edge_to_t() =
550 (* This requires that a relabel was done immediately
before returning from
553 every_n(height
, fn y
=>
554 every_n(width
, fn x
=>
555 matrix_ref(h
, y
, x
) = m1
))
557 (* This requires that a relabel was done immediately
before returning from
562 map_n_vector(width
, fn x
=>
563 not(matrix_ref(h
, y
, x
) = m1
)))
564 fun loop(lg_v
, v_max
) =
566 then (pormat("V-MAX=~s~%",[Int v_max
])
567 ; preflow_push(v_max
+ 1)
569 else let val v
= v_max
+ let
573 else loop(i
- 1, c
+ c
)
576 in pormat("LG-V=~s, V-MAX=~s, V=~s~%",
577 [Int lg_v
, Int v_max
, Int v
])
580 if min_cut_includes_every_edge_to_t()
593 val c_right
= make_matrix(height
, width
- 1, ~
1)
594 val c_down
= make_matrix(height
- 1, width
, ~
1)
595 in doo(height
, fn y
=>
596 doo(width
- 1, fn x
=>
599 if (y
>= quotient(height
, 4)
600 andalso y
< quotient(3 * height
, 4)
601 andalso (x
= quotient(width
, 4) - 1
602 orelse x
= quotient(3 * width
, 4) - 1))
605 ; doo(height
- 1, fn y
=>
609 if (x
>= quotient(width
, 4)
610 andalso x
< quotient(3 * width
, 4)
611 andalso (y
= quotient(height
, 4) - 1
612 orelse y
= quotient(3 * height
, 4) - 1))
615 ; rao_ratio_region(c_right
, c_down
,
616 make_matrix(height
, width
, 1),