;OBJEKT1
;This program labels objects in a segmented image.
;(Segmentation can be performed using e.g. UDISOL.PRO.)
;Objects are defined as subsets of non-zero pixels which
;are neighboring side-by-side. Single-pixel objects are
;considered as a noise and excluded, i.e. their object
;number is set to zero. Pixels forming an object are
;labelled by the same object number, i.e. all pixels of
;the n-th object are labelled by object number 'n'.
;
;CALLING SEQUENCE: objekt1,image,table
;
;INPUT: 2-D array IMAGE - a segmented image, where the
;	potential object pixels have a non-zero value,
;	other (background) pixels have a zero value.
;OUTPUT: 2-D array TABLE - a two-column table, where
;	the 0-th column contains object numbers and
;	the 1-st column contains subscripts of non-zero
;	pixels in 1-D representation (in ascending order).
;
;The program consists of main module OBJEKT1 and
;a subroutine-function module NEIGHBOR.
;
;7 March 1996, Michal.

;--------------------------------------------------------------------

FUNCTION neighbor,a,p,nx

;SUBROUTINE FOR OBJEKT1
;Search in a 1-D list of subscripts A of "existing" (non-zero)
;image pixels for 4 side-by-side neighbors of a given pixel
;with position P.
;NX is the x-size of the image (input).
;RESULT: 2 x 4 array NA containing in the 0-th column the positions
;  of existing neighbors and in the 1-st column the existence flag
;  of the neighbor (0/1).

na=intarr(2,4)

   p1=p-1			;neighbor on left
   w=where(p1 eq a,c)	        ;is in the list? and where?
   na(0,0)=w			;store position
   na(1,0)=c			;store existence flag
   p2=p+1			;neighbor on right
   w=where(p2 eq a,c)
   na(0,1)=w
   na(1,1)=c
   p3=p-nx			;neighbor below
   w=where(p3 eq a,c)
   na(0,2)=w
   na(1,2)=c
   p4=p+nx			;neighbor above
   w=where(p4 eq a,c)
   na(0,3)=w
   na(1,3)=c

RETURN,na
end

;---------------------------------------------------------------------

PRO objekt1,ima,tab
;
;This program is a bit difficult to understand - but it works!
;First, a 1-D list A of subscripts of non-zero pixels (WHERE)
;is made. Then a first pixel is checked for neighbors. If no
;neighbors are found, the pixel is not labelled by the object
;number (in array ON), only marked as 'checked' (in array TS)
;and a new pixel is investigated.
;Pixel with 1 or more neighbors is labelled by object number
;(array ON) and marked as 'checked' (value 1 in array TS).
;All its neighbors are labelled by the same object number.
;Then, all labelled and non-checked pixels (the newly found
;neighbors) are checked for their neighbors etc. So, in one
;step a neighbor is labelled and in the next step it is checked
;for its own neighbors. If all labelled pixels are checked (and
;no more neighbors appear), the object is complete.
;The search for next object begins from the next non-labelled
;and non-checked pixel; the object number is increased by 1.
;
; WARNING! Since the 1-D positions are stored in the INTEGER
; arrays, the code is limited to max. 1-D position = 32767,
; i.e. to frame size 181x181 or similar.
;

on_error,1

s=size(ima)
if s(0) ne 2 then message,'*** Input must be 2-D array'
if s(1)*s(2) gt 32767 then message,'*** Input image too large'

a=where(ima,n)		;list of coordinates of nonzero pixels (1D)
if n le 0    then message,'*** Input image is empty'
if n eq s(4) then message,'*** Input image is not segmented'
on=intarr(n)		;list of object numbers assigned to pixels
ts=intarr(n)		;list of 0/1 flags if pixel was checked
			;   for neighbors
cnt=0			;counter of object numbers


newobject:		;start search for an object

w=where(on eq 0 and ts eq 0,tslast)  ;first pixel of object
if tslast eq 0 then goto,fin	   ;no more unlabelled and
				   ; unchecked pixels - END
pos=a(w(0))		;take the first appropriate pixel

na=NEIGHBOR(a,pos,s(1))			;search for neighbors

if total(na(1,*)) le 0 then begin	;no neighbors, i.e.
  ts(w(0))=1				;no object, checked
  goto,newobject
endif

cnt=cnt+1				;this will be an object -
					;at least 1 neighbor exists
on(w(0))=cnt		;label 1st pixel by object no.

for i=0,3 do begin	;label the neighbors by object no.
  if na(1,i) gt 0 then begin
    if ts(na(0,i)) eq 1 then message, $
    '*** Something goes wrong, neighbor is a member of another object!'
    on(na(0,i))=cnt
  endif
endfor

ts(w(0))=1		;first pixel was checked for neighbors


newpix:			;loop for next pixels of the object

w=where(on eq cnt and ts eq 0,last)
if last eq 0 then goto,newobject   ;no appropriate pixels remain

pos=a(w(0))		;next non-checked and labelled pixel

na=NEIGHBOR(a,pos,s(1)) ;check for its neighbors

;label the neighbors by object no. if they exist and are unlabelled
for i=0,3 do $
 if (na(1,i) gt 0) then if (on(na(0,i)) eq 0) then on(na(0,i))=cnt

ts(w(0))=1		;this pixel was checked
goto,newpix


fin:

tab=intarr(2,n)
tab(0,*)=on
tab(1,*)=a

;print,' Number of nonzero pixels:',n
print,' Number of objects:',cnt

end