;-----------------------------------------------------------------
; ingen.pro
;
; written by Andreas Reigber
;------------------------------------------------------------------
; Generates interferogram / coherence maps of coregistrated SAR
; images. By use of block processing memory consumption can be
; minimised. Includes also presumming and region selection.
;
; Usage: 
; ingen,file1,file2
;    file1 = filename of SLC 1 (in sarr-format)
;    file2 = filename of SLC 2 (in sarr-format)
; 
; Output goes to files in the local directory, in sarr format 
; as well as jpeg (see below for filenames)
;
; Keywords:
;    PRE_RG : Presumming in range (default = 1)
;    PRE_AZ : Presumming in azimith (default = 1)
;    KOHSX  : Window size in range for coherence calculation (default = 7)
;    KOHSY  : Window size in azimuth for coherence calculation (default = 7)
;    FE     : Array containing flat earth phase
;    STARTX : First range bin
;    ENDX   : Last range bin
;    STARTY : First azimuth bin
;    ENDY   : Last azimith bin
;    BLOCKSIZE : Blocksize in azimuth for block processing (default = 2048)
;
;------------------------------------------------------------------
; This software has been released under the terms of the GNU Public
; license. See http://www.gnu.org/copyleft/gpl.html for details.
;------------------------------------------------------------------

@dot
@rarr
@sarr
@rebinc
@smooth2d
pro ingen,f1,f2,PRE_RG=pre_rg,PRE_AZ=pre_az,KOHSX=kohsx,KOHSY=kohsy,FE=fe, $
	STARTX=startx,ENDX=endx,STARTY=starty,ENDY=endy,BLOCKSIZE=blocksize
   
	print,""
	print,"---> INGEN "
	print,""
	
	file1 = f1
	file2 = f2
	file3 = "slc1.dat"
	file4 = "slc2.dat"
	file5 = "pha.dat"
	file6 = "int.dat"
	file7 = "koh.dat"  
	file8 = "int.jpg"
	file9 = "koh.jpg"
	file10= "amp1.jpg"
	file11= "amp2.jpg"

	if not keyword_set(blocksize) then blocksize=2048
	if not keyword_set(pre_rg) then pre_rg = 1
	if not keyword_set(pre_az) then pre_az = 1
	if not keyword_set(kohsx)  then kohsx = 7
	if not keyword_set(kohsy)  then kohsy = 7
	
;-----------------------------------		
	head = lonarr(1)
	rarr,file1,ins1,header=head
	rarr,file2,ins2,header=head
   
   anz_az = head[2]
   anz_rg = head[1]

	if not keyword_set(fe)     then fe=dblarr(anz_rg)+1.0
	if not keyword_set(startx) then startx = 0
	if not keyword_set(endx)   then endx = anz_rg-1
	if not keyword_set(starty) then starty = 0
	if not keyword_set(endy)   then endy = anz_az-1

	print," presumming factors : "+strcompress(pre_rg)+" x"+strcompress(pre_az)
	print," coherence window   : "+strcompress(kohsx)+" x"+strcompress(kohsy)
	print,' image size         : '+strcompress(anz_rg)+" x"+strcompress(anz_az)
	print,' selected area      : '+strcompress(startx)+strcompress(starty)+strcompress(endx)+strcompress(endy)
	print,' block size         : '+strcompress(blocksize)

	x1 = startx/pre_rg
	x2 = endx/pre_rg-1
	y1 = starty/pre_az
	y2 = endy/pre_az-1
   
	anz_rg_p = long(anz_rg / pre_rg)
	anz_az_p = long(anz_az / pre_az)
	dim = 2l

	sarr,file3,outslc1,HEADER=[dim,anz_rg_p,anz_az_p,4l]
	sarr,file4,outslc2,HEADER=[dim,anz_rg_p,anz_az_p,4l]
	sarr,file5,outpha, HEADER=[dim,anz_rg_p,anz_az_p,4l]
	sarr,file6,outint, HEADER=[dim,anz_rg_p,anz_az_p,6l]
  	sarr,file7,outkoh, HEADER=[dim,anz_rg_p,anz_az_p,4l]

	big_length_p = blocksize / pre_az
	b1 = complexarr(anz_rg,blocksize)
	b2 = complexarr(anz_rg,blocksize)
	fl = (fltarr(blocksize)+1) ## exp(complex(0,-fe))

	ba = kohsy
	be = big_length_p - kohsy - 1
	anz_blocks  = (anz_az - 2*(blocksize-ba*pre_az))/(blocksize-2*ba*pre_az) + 2

	print,''
	dot,anz_blocks,1,message=" Calculating  "	
	for i=0,anz_blocks-1 do begin
		dot,1,i
		readu,ins1,b1
		readu,ins2,b2
		b2 = b2 * fl

		b3 = b1 * conj(b2)
		sb1 = rebin(abs(b1)^2,anz_rg_p,big_length_p)
		sb2 = rebin(abs(b2)^2,anz_rg_p,big_length_p)
		sb3 = rebinc(b3,anz_rg_p,big_length_p)
		
		  
		koh = float(abs(smooth2d(sb3,kohsx,kohsy))/sqrt(smooth2d(double(sb1),kohsx,kohsy)*smooth2d(double(sb2),kohsx,kohsy)))
		aux = where(koh gt 1,nn)
		if nn gt 0 then koh[aux] = 1.0

		if i eq 0 then begin
			writeu,outslc1,sb1[*,0:be]
			writeu,outslc2,sb2[*,0:be]
			writeu,outint ,sb3[*,0:be]
			writeu,outpha ,atan(sb3[*,0:be])			
			writeu,outkoh ,koh[*,0:be]
			count = be + 1
		endif else begin
			writeu,outslc1,sb1[*,ba:be]
			writeu,outslc2,sb2[*,ba:be]
			writeu,outint ,sb3[*,ba:be]
			writeu,outpha ,atan(sb3[*,ba:be])			
			writeu,outkoh ,koh[*,ba:be]
			count = count + be - ba + 1
		endelse
		point_lun,-ins1,fpos1
		point_lun,-ins2,fpos2
		point_lun,ins1,fpos1-2l*ba*anz_rg*8l*pre_az
		point_lun,ins2,fpos1-2l*ba*anz_rg*8l*pre_az	
	endfor

	point_lun,-ins1,fpos1
	fpos1 = fpos1/(anz_rg*8l)

	rest_length = anz_az - fpos1
	difmod = rest_length - rest_length / pre_az * pre_az
	rest_length = rest_length - difmod
	
	rest_length_p = rest_length / pre_az
	fl = (fltarr(rest_length)+1) ## exp(complex(0,-fe))
	b1 = complexarr(anz_rg,rest_length)
	b2 = complexarr(anz_rg,rest_length)
	dot,1,i
	readu,ins1,b1
	readu,ins2,b2
	point_lun,-ins1,fpos1

	b2 = b2 * fl
	b3 = b1 * conj(b2)
	sb1 = rebin(abs(b1)^2,anz_rg_p,rest_length_p)
	sb2 = rebin(abs(b2)^2,anz_rg_p,rest_length_p)
	sb3 = rebinc(b3,anz_rg_p,rest_length_p)

	koh = float(abs(smooth2d(sb3,kohsx,kohsy))/sqrt(smooth2d(double(sb1),kohsx,kohsy)*smooth2d(double(sb2),kohsx,kohsy)))
	aux = where(koh gt 1,nn)
	if nn gt 0 then koh[aux] = 1.0

	writeu,outslc1,sqrt(sb1[*,ba:*])
	writeu,outslc2,sqrt(sb2[*,ba:*])
	writeu,outint ,sb3[*,ba:*]
	writeu,outpha ,atan(sb3[*,ba:*])			
	writeu,outkoh ,koh[*,ba:*]
	
	free_lun,outslc1,outslc2,outint,outpha,outkoh,ins1,ins2
;-------------------------------------------------------------------------------
	print,''
	rarr,file5,bild
	write_jpeg,file8,bytscl(bild[x1:x2,y1:y2]),quality=50
	rarr,file7,bild,INFO=info
 	write_jpeg,file9,bytscl(bild[x1:x2,y1:y2],0,1),quality=50
	rarr,file3,bild
	rarr,file4,bild2
	bild  = (sqrt(bild))[x1:x2,y1:y2]
	bild2 = (sqrt(bild2))[x1:x2,y1:y2]
	write_jpeg,file10,bytscl(bild,0,2.5*avg(bild)),quality=50
	write_jpeg,file11,bytscl(bild2,0,2.5*avg(bild2)),quality=50
	print,''
end