#! /usr/bin/env python3 # def components_1d ( A ): #*****************************************************************************80 # ## components_1d() assigns contiguous nonzero pixels to a common component. # # Discussion: # # This calculation is trivial compared to the 2D problem, and is included # primarily for comparison. # # On input, the A array contains values of 0 or 1. # # The 0 pixels are to be ignored. The 1 pixels are to be grouped # into connected components. # # The pixel A(I) is "connected" to the pixels A(I-1) and A(I+1). # # On output, COMPONENT_NUM reports the number of components of nonzero # data, and the array C contains the component assignment for # each nonzero pixel, and is 0 for zero pixels. # # Picture: # # Input A: # # 0 0 1 2 4 0 0 4 0 0 0 8 9 9 1 2 3 0 0 5 0 1 6 0 0 0 4 0 # # Output: # # C: # # 0 0 1 1 1 0 0 2 0 0 0 3 3 3 3 3 3 0 0 4 0 5 5 0 0 0 6 0 # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 July 2022 # # Author: # # John Burkardt # # Input: # # integer A(N), the pixel array. # # Output: # # integer C(N), the component array. # import numpy as np n = len ( A ) # # Initialization. # C = np.zeros ( n, dtype = np.int ) component_num = 0 # # "Read" the array one pixel at a time. If a (nonzero) pixel's west neighbor # already has a label, the current pixel inherits it. Otherwise, we have # begun a new component. # west = 0 for j in range ( 0, n ): if ( A[j] != 0 ): if ( west == 0 ): component_num = component_num + 1 C[j] = component_num west = C[j] return C def components_1d_test ( ): #*****************************************************************************80 # ## components_1d_test() tests components_1d(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 July 2022 # # Author: # # John Burkardt # import numpy as np n = 28 A = np.array ( [ \ 0, 0, 1, 2, 4, 0, 0, 4, 0, 0, \ 0, 8, 9, 9, 1, 2, 3, 0, 0, 5, \ 0, 1, 6, 0, 0, 0, 4, 0 ] ) print ( '' ) print ( 'components_1d_test():' ) print ( ' components_1d() finds and labels connected' ) print ( ' components in a 1D integer vector.' ) print ( '' ) print ( ' A:' ) print ( '' ) print ( A ) C = components_1d ( A ) component_num = np.max ( C ) print ( '' ) print ( ' Number of components =', component_num ) print ( '' ) print ( ' C:' ) print ( '' ) print ( C ) return def components_2d ( A ): #*****************************************************************************80 # ## components_2d() assigns contiguous nonzero pixels to a common component. # # Discussion: # # On input, the A array contains values of 0 or 1. # # The 0 pixels are to be ignored. The 1 pixels are to be grouped # into connected components. # # The pixel A(I,J) is "connected" to the pixels A(I-1,J), A(I+1,J), # A(I,J-1) and A(I,J+1), so most pixels have 4 neighbors. # # (Another choice would be to assume that a pixel was connected # to the other 8 pixels in the 3x3 block containing it.) # # On output, COMPONENT_NUM reports the number of components of nonzero # data, and the array C contains the component assignment for # each nonzero pixel, and is 0 for zero pixels. # # Picture: # # Input A: # # 0 2 0 0 17 0 3 # 0 0 3 0 1 0 4 # 1 0 4 8 8 0 7 # 3 0 6 45 0 0 0 # 3 17 0 5 9 2 5 # # Output: # # C: # # 0 1 0 0 2 0 3 # 0 0 2 0 2 0 3 # 4 0 2 2 2 0 3 # 4 0 2 2 0 0 0 # 4 4 0 2 2 2 2 # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 July 2022 # # Author: # # John Burkardt # # Input: # # integer A(M,N), the pixel array. # # Output: # # integer C(M,N), the component array. # import numpy as np m, n = A.shape # # Initialization. # C = np.zeros ( [ m, n ], dtype = np.int ) component_num = 0 # # P is simply used to store the component labels. The dimension used # here is, of course, usually an absurd overestimate. # p = np.arange ( 0, m * n + 1 ) # # "Read" the array one pixel at a time. If a (nonzero) pixel's north or # west neighbor already has a label, the current pixel inherits it. # In case the labels disagree, we need to adjust the P array so we can # later deal with the fact that the two labels need to be merged. # for i in range ( 0, m ): for j in range ( 0, n ): if ( i == 0 ): north = 0 else: north = C[i-1,j] if ( j == 0 ): west = 0 else: west = C[i,j-1] if ( A[i,j] != 0 ): if ( north == 0 ): if ( west == 0 ): component_num = component_num + 1 C[i,j] = component_num else: C[i,j] = west elif ( north != 0 ): if ( west == 0 or west == north ): C[i,j] = north else: C[i,j] = min ( north, west ) if ( north < west ): p[west] = north else: p[north] = west # # When a component has multiple labels, have the higher labels # point to the lowest one. # for component in range ( component_num -1, -1, -1 ): b = component while ( p[b] != b ): b = p[b] p[component] = b # # Locate the minimum label for each component. # Assign these mininum labels new consecutive indices. # q = np.zeros ( component_num, dtype = np.int ) i = 0 for component in range ( 0, component_num ): if ( p[component] == component ): i = i + 1 q[component] = i component_num = i # # Replace the labels by consecutive labels. # Need to avoid C(i,j) = 0. # i = np.where ( C != 0 ) C[i] = q [ p [ C[i] ] ] - 1 return C def components_2d_test ( ): #*****************************************************************************80 # ## components_2d_test () tests components_2d(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 July 2022 # # Author: # # John Burkardt # import numpy as np A = np.array ( [ [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ], \ [ 0, 0, 1, 0, 0, 1, 0, 0, 0 ], \ [ 0, 1, 1, 0, 1, 1, 1, 0, 0 ], \ [ 0, 1, 1, 1, 1, 1, 1, 0, 0 ], \ [ 0, 0, 1, 1, 1, 0, 0, 0, 0 ], \ [ 0, 0, 1, 1, 1, 0, 0, 0, 0 ], \ [ 0, 1, 1, 1, 0, 1, 0, 1, 0 ], \ [ 0, 1, 1, 0, 0, 1, 0, 1, 0 ], \ [ 0, 0, 1, 0, 0, 0, 0, 1, 0 ], \ [ 0, 0, 0, 0, 1, 0, 1, 1, 0 ], \ [ 0, 1, 0, 1, 1, 0, 1, 0, 0 ], \ [ 0, 1, 1, 1, 1, 1, 0, 0, 0 ], \ [ 0, 0, 1, 1, 0, 1, 0, 1, 0 ], \ [ 0, 0, 1, 1, 0, 1, 0, 1, 0 ], \ [ 0, 1, 1, 0, 1, 0, 1, 1, 0 ], \ [ 0, 1, 0, 0, 1, 0, 1, 1, 0 ], \ [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ] ] ) print ( '' ) print ( 'components_2d_test():' ) print ( ' components_2d() finds and labels connected' ) print ( ' components in a 2D integer array.' ) print ( '' ) print ( ' A:' ) print ( '' ) print ( A ) C = components_2d ( A ) component_num = np.max ( C ) print ( '' ) print ( ' Number of components = ', component_num ) print ( '' ) print ( ' C:' ) print ( '' ) print ( C ) return def components_3d ( A ): #*****************************************************************************80 # ## components_3d() assigns contiguous nonzero pixels to a common component. # # Discussion: # # On input, the A array contains values of 0 or 1. # # The 0 pixels are to be ignored. The 1 pixels are to be grouped # into connected components. # # The pixel A(I,J,K) is "connected" to the pixels: # # A(I-1,J, K ), A(I+1,J, K ), # A(I, J-1,K ), A(I, J+1,K ), # A(I, J, K-1), A(I, J, K+1), # # so most pixels have 6 neighbors. # # On output, COMPONENT_NUM reports the number of components of nonzero # data, and the array C contains the component assignment for # each nonzero pixel, and is 0 for zero pixels. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 July 2022 # # Author: # # John Burkardt # # Input: # # integer L, M, N, the order of the array. # # integer A(L,M,N), the pixel array. # # Output: # # integer C(L,M,N), the component array. # import numpy as np l, m, n = A.shape # # Initialization. # C = np.zeros ( [ l, m, n ], dtype = np.int ) component_num = 0 # # P is simply used to store the component labels. The dimension used # here is, of course, usually an absurd overestimate. # p = np.arange ( 0, l * m * n + 1 ) # # "Read" the array one pixel at a time. If a (nonzero) pixel's north or # west neighbor already has a label, the current pixel inherits it. # In case the labels disagree, we need to adjust the P array so we can # later deal with the fact that the two labels need to be merged. # for i in range ( 0, l ): for j in range ( 0, m ): for k in range ( 0, n ): if ( i == 1 ): north = 0 else: north = C[i-1,j,k] if ( j == 1 ): west = 0 else: west = C[i,j-1,k] if ( k == 1 ): up = 0 else: up = C[i,j,k-1] if ( A[i,j,k] != 0 ): # # New component? # if ( north == 0 and west == 0 and up == 0 ): component_num = component_num + 1 C[i,j,k] = component_num # # One predecessor is labeled. # elif ( north != 0 and west == 0 and up == 0 ): C[i,j,k] = north elif ( north == 0 and west != 0 and up == 0 ): C[i,j,k] = west elif ( north == 0 and west == 0 and up != 0 ): C[i,j,k] = up # # Two predecessors are labeled. # elif ( north == 0 and west != 0 and up != 0 ): C[i,j,k] = min ( west, up ) c1 = min ( p[west], p[up] ) p[west] = c1 p[up] = c1 elif ( north != 0 and west == 0 and up != 0 ): C[i,j,k] = min ( north, up ) c1 = min ( p[north], p[up] ) p[north] = c1 p[up] = c1 elif ( north != 0 and west != 0 and up == 0 ): C[i,j,k] = min ( north, west ) c1 = min ( p[north], p[west] ) p[north] = c1 p[west] = c1 # # Three predecessors are labeled. # elif ( north != 0 and west != 0 and up != 0 ): C[i,j,k] = min ( north, min ( west, up ) ) c1 = min ( p[north], min ( p[west], p[up] ) ) p[north] = c1 p[west] = c1 p[up] = c1 # # When a component has multiple labels, have the higher labels # point to the lowest one. # for component in range ( component_num -1, -1, -1 ): b = component while ( p[b] != b ): b = p[b] p[component] = b # # Locate the minimum label for each component. # Assign these mininum labels new consecutive indices. # q = np.zeros ( component_num, dtype = np.int ) i = 0 for component in range ( 0, component_num ): if ( p[component] == component ): i = i + 1 q[component] = i component_num = i # # Replace the labels by consecutive labels. # i = np.where ( C != 0 ) C[i] = q [ p [ C[i] ] ] - 1 return C def components_3d_test( ): #*****************************************************************************80 # ## components_3d_test() tests components_3d(). # # Discussion: # # Each cell is assumed to have at most 6 neighbors. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 July 2022 # # Author: # # John Burkardt # import numpy as np l = 64 m = 64 n = 26 print ( '' ) print ( 'components_3d_test():' ) print ( ' components_3d() finds and labels connected' ) print ( ' components in a 3D integer block.' ) print ( '' ) print ( ' A is a 3D block of order', l, 'x', m, 'x', n ) # # Retrieve the indices of nonzero data in A by reading a file. # filename = 'indices.txt' ijk = np.loadtxt ( filename ) ijk = ijk.astype ( int ) m1, n1 = ijk.shape # # Set the nonzero entries of A. # A = np.zeros ( [ l, m, n ], dtype = np.int ) for i1 in range ( 0, m1 ): i = ijk[i1,0] - 1 j = ijk[i1,1] - 1 k = ijk[i1,2] - 1 A[i,j,k] = 1 print ( '' ) print ( ' Number of nonzero A values is ', np.sum ( A ) ) # # Determine the components. # C = components_3d ( A ) component_num = np.max ( C ) s = np.zeros ( component_num ) print ( '' ) print ( ' Number of components = ', component_num ) for i in range ( 0, l ): for j in range ( 0, m ): for k in range ( 0, n ): if ( C[i,j,k] != 0 ): s[C[i,j,k]-1] = s[C[i,j,k]-1] + 1 print ( '' ) print ( ' Component Size' ) print ( '' ) for i in range ( 0, component_num ): print ( ' %4d %8d' % ( i, s[i] ) ) print ( '------ --------' ) print ( ' Total %8d' % np.sum ( s ) ) return def image_components_test ( ): #*****************************************************************************80 # ## image_components_test() tests image_components(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 16 July 2022 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'image_components_test():' ) print ( ' Python version: ' + platform.python_version ( ) ) print ( ' Test image_components().' ) components_1d_test ( ) components_2d_test ( ) components_3d_test ( ) # # Terminate. # print ( '' ) print ( 'image_components_test():' ) print ( ' Normal end of execution.' ) return def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 August 2019 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return if ( __name__ == '__main__' ): timestamp ( ) image_components_test ( ) timestamp ( )