pro centroid,image,xdim,ydim,yc

;  This program calculates the centroid PIXEL position (yc)
;  for a circular image within a rectangular byte array.

   count = 0
   fac =0 

   ind = where(image ne 0)
   sum0 = n_elements(ind)
   yc0 =  fix(float(ind((float(sum0)/2.)-1))/float(ydim))
   yc = yc0
;   yc0 = (fix(ydim)/2)

 if ydim mod 2 ne 0 then begin ; Loop for odd dimension
   flag0:
   count = count + 1
   sub1 = image(*,0:yc0-1)
   sub2 = image(*,yc0+1:ydim-1)
   sum1 = n_elements(where(sub1 ne 0))
   sum2 = n_elements(where(sub2 ne 0))
   diff = sum1 - sum2
   if diff eq 0 then begin
     yc = yc0 
     goto,flag2
   endif
;  Check for local minimum
   ycn = yc0 - 1
   subn1 = image(*,0:ycn-1)
   subn2 = image(*,ycn+1:ydim-1)
   sumn1 = n_elements(where(subn1 ne 0))
   sumn2 = n_elements(where(subn2 ne 0))
   diffn = sumn1 - sumn2
   ycp = yc0 + 1
   subp1 = image(*,0:ycp-1)
   subp2 = image(*,ycp+1:ydim-1)
   sump1 = n_elements(where(subp1 ne 0))
   sump2 = n_elements(where(subp2 ne 0))
   diffp = sump1 - sump2
   if (abs(diff) le abs(diffn)) and (abs(diff) le abs(diffp)) then begin
     yc = yc0 
     goto,flag2
   endif
   if diff gt 0 then begin
     fac = float(sum1)/(float(sum2)+.0001)
     if fac gt 10 then fac = 10.
     yc0 = yc0 - 1*fac
     goto, flag0
   endif 
   if diff lt 0 then begin
     fac = float(sum2)/(float(sum1)+.0001)
     if fac gt 10 then fac = 10.
     yc0 = yc0 + 1*fac
     goto, flag0
   endif 
 endif else begin  ; Loop for even dimension
   flag1:
   count = count + 1
   sub1 = image(*,0:yc0)
   sum1 = n_elements(where(sub1 ne 0))
   diff = 2 * sum1 - sum0
   if diff eq 0 then begin
     yc = yc0 - 1
     goto,flag2
   endif
;  Check for local minimum
   ycn = yc0 - 1
   subn = image(*,0:ycn)
   sumn = n_elements(where(subn ne 0))
   diffn = 2 * sumn - sum0
   ycp = yc0 + 1
   subp = image(*,0:ycp)
   sump = n_elements(where(subp ne 0))
   diffp = 2 * sump - sum0
   if (abs(diff) le abs(diffn)) and (abs(diff) le abs(diffp)) then begin
     yc = yc0 - 1
     goto,flag2
   endif
   if diff gt 0 then begin
     fac = float(sum1)/((float(sum0)-float(sum1))+.0001) 
     if fac gt 10 then fac = 10.
     yc0 = yc0 - 1*fac
     goto, flag1
   endif 
   if diff lt 0 then begin
     fac = (float(sum0)-float(sum1))/(float(sum1)+.0001) 
     if fac gt 10 then fac = 10.
     yc0 = yc0 + 1*fac
     goto, flag1
   endif 
 endelse
flag2:
end