;
;+
;  Unwraps SXT low 8 bit visible light images or other
;  low-bit wrapped images.

PRO unwrap, image, newimage, maxnum=maxnum, noshow=noshow

;
; INPUT PARAMETERS:
;	image = low-8 byte data
;
; OUTPUT PARAMETERS:
;	newimage = integer, unwrapped version
;
; OPTIONAL INPUT KEYWORDS:
;	maxnum = data value at which data wraps.  Default = 256.
;	noshow = if set, do not display intermediate images.
;		Default = display.
;
; Restrictions:
;	Has trouble with sunspots.  Output can be patched with
;	PAINT_PIXELS to finish the job.  Once a single clean
;	image is made, a data cube can be cleaned with
;	CUBE_FILL.
;	If the result looks poor, try rotating the input image.
;
; Written October 1992  Barry LaBonte
; Fixed null array bug  April 11, 1994  BJL
; Added /MAXNUM keyword April 12, 1994  BJL
;
;-

arrsize = SIZE(image)
npxlx = arrsize(1)
npxly = arrsize(2)

IF( KEYWORD_SET(maxnum) EQ 0 ) THEN maxnum = 256
mx2 = maxnum/2

gradient = INTARR(npxlx,npxly)
medimage = INTARR(npxlx,npxly)

lintot = LONARR(npxly)


IF( KEYWORD_SET(noshow) EQ 0 ) THEN TVSCL,image

; Take x "gradient" (adjacent pixel difference)
gradient = FIX(image)
gradient(1:npxlx-1,*) = gradient(1:npxlx-1,*) - gradient(0:npxlx-2,*)

IF( KEYWORD_SET(noshow) EQ 0 ) THEN TVSCL,gradient

; Remove spikes
gradient = gradient - maxnum * (gradient/mx2) 

IF( KEYWORD_SET(noshow) EQ 0 ) THEN TVSCL,gradient

; Remove excessive gradients
FOR j=0,npxly-1 DO BEGIN
	medimage(*,j) = MEDIAN(gradient(*,j),3)
ENDFOR
gradient = gradient - maxnum * ((gradient - medimage)/mx2)

IF( KEYWORD_SET(noshow) EQ 0 ) THEN TVSCL,gradient

; Reintegrate image, avoid negative values
newimage = gradient
FOR I=1,npxlx-1 DO BEGIN
	newimage(I,*) = newimage(I-1,*) + gradient(I,*)
ENDFOR
IF( MIN(newimage) LT 0 ) THEN newimage = newimage - MIN(newimage)

IF( KEYWORD_SET(noshow) EQ 0 ) THEN TVSCL,newimage

; Clean streaks, using last column and row average to count pixels
; First the row total
FOR j=0,npxly-1 DO BEGIN
	lintot(j) = TOTAL(newimage(*,j))
& ENDFOR

; Now the columns of error (streaks)
line = TRANSPOSE(newimage(npxlx-1,*))
linerr = maxnum * (line/maxnum)
linerr = linerr + maxnum * ((line - linerr)/mx2)
medlin = maxnum * ( MEDIAN(linerr)/maxnum )

; Get the smooth shape
smoothline = lintot
bad = WHERE(linerr NE medlin, bcount)
good = WHERE(linerr EQ medlin, gcount)
IF( bcount GT 0 AND gcount GT 0 ) THEN BEGIN
	smoothline(bad) = SPLINE(good, lintot(good), bad)
ENDIF

; Now get the number of bad pixels
numbad = INTARR(npxly)
nonz = WHERE(linerr NE 0, zcount)
IF( zcount GT 0 ) THEN BEGIN
	numbad(nonz) = (lintot(nonz) - smoothline(nonz))/ linerr(nonz)
	numbad = numbad < (npxlx - 1)
	numbad = numbad > 0
ENDIF

; Remove streaks
FOR j=0,npxly-1 DO BEGIN
	 newimage(npxlx-1-numbad(j):npxlx-1,j) = newimage(npxlx-1-numbad(j):npxlx-1,j) - linerr(j)
ENDFOR

IF( KEYWORD_SET(noshow) EQ 0 ) THEN TVSCL,newimage

; Now do a final cleanup
; Clean streaks - first n x +/-maxnum, then 1 x +/- maxnum
FOR I=0,npxlx-1 DO BEGIN
	 medimage(I,*) = TRANSPOSE(MEDIAN(TRANSPOSE(newimage(I,*)),5))
ENDFOR
newimage = newimage - maxnum * ((newimage - medimage)/maxnum)
IF( KEYWORD_SET(noshow) EQ 0 ) THEN TVSCL,newimage
FOR I=0,npxlx-1 DO BEGIN
	 medimage(I,*) = TRANSPOSE(MEDIAN(TRANSPOSE(newimage(I,*)),5))
ENDFOR
newimage = newimage - maxnum * ((newimage - medimage)/mx2)

IF( KEYWORD_SET(noshow) EQ 0 ) THEN TVSCL,newimage

RETURN
END