FUNCTION Search3d, array, xpos, ypos, zpos, min_val, max_val, $
Decrease=decrease, Increase=increase, $
Lpf_band=smooth_band, Diagonal=diagonal
; *** Test inputs
ON_ERROR, 2
size_array = Size(array)
IF (size_array[0] NE 3L) THEN BEGIN
MESSAGE, 'Array must have three dimensions'
ENDIF
x_size = size_array[1]
y_size = size_array[2]
z_size = size_array[3]
xpos = Long(xpos[0])
ypos = Long(ypos[0])
zpos = Long(zpos[0])
IF (xpos LT 0L) THEN BEGIN
MESSAGE, 'Xpos must be >= 0'
ENDIF
IF (xpos GE x_size) THEN BEGIN
MESSAGE, 'Xpos must be < array size'
ENDIF
IF (ypos LT 0L) THEN BEGIN
MESSAGE, 'Ypos must be >= 0'
ENDIF
IF (ypos GE y_size) THEN BEGIN
MESSAGE, 'Ypos must be < array size'
ENDIF
IF (zpos LT 0L) THEN BEGIN
MESSAGE, 'Zpos must be >= 0'
ENDIF
IF (zpos GE z_size) THEN BEGIN
MESSAGE, 'Zpos must be < array size'
ENDIF
min_val = min_val[0]
max_val = max_val[0]
IF (max_val LT min_val) THEN BEGIN
MESSAGE, 'Max value must be >= min value'
ENDIF
start_val = array[xpos,ypos,zpos]
IF ((start_val LT min_val) OR (start_val GT max_val)) THEN BEGIN
MESSAGE, 'Value of array at (xpos,ypos,zpos) must be >= min_val and <= max_val'
ENDIF
dec = 0.0
inc = 0.0
range = 0B
IF (N_Elements(decrease) GT 0L) THEN BEGIN
dec = Float(decrease[0])
range = 1B
ENDIF
IF (N_Elements(increase) GT 0L) THEN BEGIN
inc = Float(increase[0])
range = 1B
ENDIF
sb = 0
IF (N_Elements(smooth_band) GT 0L) THEN sb = Fix(smooth_band[0])
IF ((sb GE x_size) OR (sb GE y_size)) THEN BEGIN
MESSAGE, 'Smooth band must be < size of array'
ENDIF
diag = 0B
IF (N_Elements(diagonal) GT 0L) THEN diag = Byte(diagonal[0])
IF (range) THEN BEGIN
; *** Calculate the edge enhanced array
IF (diag) THEN BEGIN
diff_array = Float(array)
diff_array = diff_array < $
(diff_array - Shift(array, 0, 1, 0)) < $
(diff_array - Shift(array, 1, 1, 0)) < $
(diff_array - Shift(array, 1, 0, 0)) < $
(diff_array - Shift(array, 1, -1, 0)) < $
(diff_array - Shift(array, 0, -1, 0)) < $
(diff_array - Shift(array, -1, -1, 0)) < $
(diff_array - Shift(array, -1, 0, 0)) < $
(diff_array - Shift(array, -1, 1, 0)) < $
(diff_array - Shift(array, 0, 1, 1)) < $
(diff_array - Shift(array, 1, 1, 1)) < $
(diff_array - Shift(array, 1, 0, 1)) < $
(diff_array - Shift(array, 1, -1, 1)) < $
(diff_array - Shift(array, 0, -1, 1)) < $
(diff_array - Shift(array, -1, -1, 1)) < $
(diff_array - Shift(array, -1, 0, 1)) < $
(diff_array - Shift(array, -1, 1, 1)) < $
(diff_array - Shift(array, 0, 1, -1)) < $
(diff_array - Shift(array, 1, 1, -1)) < $
(diff_array - Shift(array, 1, 0, -1)) < $
(diff_array - Shift(array, 1, -1, -1)) < $
(diff_array - Shift(array, 0, -1, -1)) < $
(diff_array - Shift(array, -1, -1, -1)) < $
(diff_array - Shift(array, -1, 0, -1)) < $
(diff_array - Shift(array, -1, 1, -1)) < $
(diff_array - Shift(array, 0, 0, 1)) < $
(diff_array - Shift(array, 0, 0, -1))
IF (sb GT 0) THEN diff_array = Smooth(diff_array, sb)
ENDIF ELSE BEGIN
diff_array = Float(array)
diff_array = diff_array < $
(diff_array - Shift(array, 0, 1, 0)) < $
(diff_array - Shift(array, 1, 0, 0)) < $
(diff_array - Shift(array, 0, -1, 0)) < $
(diff_array - Shift(array, -1, 0, 0)) < $
(diff_array - Shift(array, 0, 0, 1)) < $
(diff_array - Shift(array, 0, 0, -1))
IF (sb GT 0) THEN diff_array = Smooth(diff_array, sb)
ENDELSE
ENDIF
; *** Set up the required variables
similar_val = 1B
connect_val = 2B
c_array = Bytarr(x_size, y_size, z_size)
c_array[Where((array GE min_val) AND (array LE max_val))] = similar_val
x_size_m1 = x_size - 1L
y_size_m1 = y_size - 1L
z_size_m1 = z_size - 1L
x_ind = xpos
y_ind = ypos
z_ind = zpos
just_found = (z_ind * y_size * x_size) + (y_ind * x_size) + x_ind
c_array[just_found] = connect_val
num_found = 1L
; *** Start the search
IF (diag EQ 0B) THEN BEGIN ; *** No diagonal mode
nsew_ind = Lonarr(6, 1)
nsew_ind[0, *] = (z_ind * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + x_ind
nsew_ind[1, *] = (z_ind * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + x_ind
nsew_ind[2, *] = (z_ind * y_size * x_size) + $
(y_ind * x_size) + ((x_ind + 1L) < x_size_m1)
nsew_ind[3, *] = (z_ind * y_size * x_size) + $
(y_ind * x_size) + ((x_ind - 1L) > 0L)
nsew_ind[4, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(y_ind * x_size) + x_ind
nsew_ind[5, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(y_ind * x_size) + x_ind
cc_array = c_array[nsew_ind[*]]
just_found = Where(cc_array EQ similar_val)
; *** Loop while cells are still being found
WHILE (just_found[0] GE 0L) DO BEGIN
cc_array[just_found] = connect_val
c_array[nsew_ind[just_found]] = cc_array[just_found]
z_ind = nsew_ind[just_found] / (y_size * x_size)
y_ind = (nsew_ind[just_found] - (z_ind * y_size * x_size)) / (x_size)
x_ind = (nsew_ind[just_found] - (z_ind * y_size * x_size)) - $
(y_ind * x_size)
num_found = N_Elements(just_found)
nsew_ind = Lonarr(6, num_found, /Nozero)
nsew_ind[0, *] = (z_ind * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + x_ind
nsew_ind[1, *] = (z_ind * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + x_ind
nsew_ind[2, *] = (z_ind * y_size * x_size) + $
(y_ind * x_size) + ((x_ind + 1L) < x_size_m1)
nsew_ind[3, *] = (z_ind * y_size * x_size) + $
(y_ind * x_size) + ((x_ind - 1L) > 0L)
nsew_ind[4, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(y_ind * x_size) + x_ind
nsew_ind[5, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(y_ind * x_size) + x_ind
nsew_ind = nsew_ind[Sort(nsew_ind[*])]
nsew_ind = nsew_ind[Uniq(nsew_ind)]
cc_array = c_array[nsew_ind[*]]
IF (range) THEN BEGIN
t_array = diff_array[nsew_ind[*]]
just_found = Where((cc_array EQ similar_val) AND $
((t_array GE dec) AND $
(t_array LE inc)))
ENDIF ELSE BEGIN
just_found = Where(cc_array EQ similar_val)
ENDELSE
ENDWHILE
ENDIF ELSE BEGIN ; *** Diagonal mode
nsew_ind = Lonarr(26, 1)
nsew_ind[ 0, *] = (z_ind * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + x_ind
nsew_ind[ 1, *] = (z_ind * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + x_ind
nsew_ind[ 2, *] = (z_ind * y_size * x_size) + $
(y_ind * x_size) + ((x_ind + 1L) < x_size_m1)
nsew_ind[ 3, *] = (z_ind * y_size * x_size) + $
(y_ind * x_size) + ((x_ind - 1L) > 0L)
nsew_ind[ 4, *] = (z_ind * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[ 5, *] = (z_ind * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[ 6, *] = (z_ind * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[ 7, *] = (z_ind * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[ 8, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + x_ind
nsew_ind[ 9, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + x_ind
nsew_ind[10, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(y_ind * x_size) + ((x_ind + 1L) < x_size_m1)
nsew_ind[11, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(y_ind * x_size) + ((x_ind - 1L) > 0L)
nsew_ind[12, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[13, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[14, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[15, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[16, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + x_ind
nsew_ind[17, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + x_ind
nsew_ind[18, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(y_ind * x_size) + ((x_ind + 1L) < x_size_m1)
nsew_ind[19, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(y_ind * x_size) + ((x_ind - 1L) > 0L)
nsew_ind[20, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[21, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[22, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[23, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[24, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(y_ind * x_size) + x_ind
nsew_ind[25, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(y_ind * x_size) + x_ind
cc_array = c_array[nsew_ind[*]]
just_found = Where(cc_array EQ similar_val)
; *** Loop while cells are still being found
WHILE (just_found[0] GE 0L) DO BEGIN
cc_array[just_found] = connect_val
c_array[nsew_ind[just_found]] = cc_array[just_found]
z_ind = nsew_ind[just_found] / (y_size * x_size)
y_ind = (nsew_ind[just_found] - (z_ind * y_size * x_size)) / (x_size)
x_ind = (nsew_ind[just_found] - (z_ind * y_size * x_size)) - $
(y_ind * x_size)
num_found = N_Elements(just_found)
nsew_ind = Lonarr(26, num_found, /Nozero)
nsew_ind[ 0, *] = (z_ind * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + x_ind
nsew_ind[ 1, *] = (z_ind * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + x_ind
nsew_ind[ 2, *] = (z_ind * y_size * x_size) + $
(y_ind * x_size) + ((x_ind + 1L) < x_size_m1)
nsew_ind[ 3, *] = (z_ind * y_size * x_size) + $
(y_ind * x_size) + ((x_ind - 1L) > 0L)
nsew_ind[ 4, *] = (z_ind * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[ 5, *] = (z_ind * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[ 6, *] = (z_ind * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[ 7, *] = (z_ind * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[ 8, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + x_ind
nsew_ind[ 9, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + x_ind
nsew_ind[10, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(y_ind * x_size) + ((x_ind + 1L) < x_size_m1)
nsew_ind[11, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(y_ind * x_size) + ((x_ind - 1L) > 0L)
nsew_ind[12, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[13, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[14, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[15, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[16, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + x_ind
nsew_ind[17, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + x_ind
nsew_ind[18, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(y_ind * x_size) + ((x_ind + 1L) < x_size_m1)
nsew_ind[19, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(y_ind * x_size) + ((x_ind - 1L) > 0L)
nsew_ind[20, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[21, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[22, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind + 1L) < y_size_m1) * x_size) + $
((x_ind - 1L) > 0L)
nsew_ind[23, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(((y_ind - 1L) > 0L) * x_size) + $
((x_ind + 1L) < x_size_m1)
nsew_ind[24, *] = (((z_ind + 1L) < z_size_m1) * y_size * x_size) + $
(y_ind * x_size) + x_ind
nsew_ind[25, *] = (((z_ind - 1L) > 0L) * y_size * x_size) + $
(y_ind * x_size) + x_ind
nsew_ind = nsew_ind[Sort(nsew_ind[*])]
nsew_ind = nsew_ind[Uniq(nsew_ind)]
cc_array = c_array[nsew_ind[*]]
IF (range) THEN BEGIN
t_array = diff_array[nsew_ind[*]]
just_found = Where((cc_array EQ similar_val) AND $
((t_array GE dec) AND $
(t_array LE inc)))
ENDIF ELSE BEGIN
just_found = Where(cc_array EQ similar_val)
ENDELSE
ENDWHILE
ENDELSE
; *** Clean up and return
x_ind = 0
y_ind = 0
nsew_ind = 0
cc_array = 0
t_array = 0
diff_array = 0
index = Where(c_array EQ connect_val)
RETURN, index
END