pro conlsq2,blur,psf,result,trial=trial,gammat=gammat,gammal=gammal
;+
;			conlsq2
;
; 2-Dimensional constrainted least squares solution
;
; CALLING SEQUENCE:
;	conlsq2,blur,psf,result
;	KEYWORD PARAMETERS:  trial,gammat,gammal
;
; INPUTS:
;	blur - the blurred image
;	psf - the point spread function
; OUTPUTS:
;	result - restored image
; OPTIONAL KEYWORD INPUTS:
;	trail - trail solution (if not supplied it is set to all zeros)
;	gammat - reciprocal Lagrangian multiplier for the trial solution
;		used to control the difference from the trial solution
;		constraint (default = 0.01)
;	gammal - reciprocal Lagrangian multiplier for the laplacian constraint
;		used to control the mininimization of the Laplacian.
;		default = 0.0
; METHOD:
;	solves the system
;
;		(HtH + gammat*I + gammal*QtQ)x = Htb + gammat*trial
;
;	where:
;		H, Ht is the point spread function matrix and its
;			transpose
;		Q, Qt is the laplacian matrix and its transpose
;		I - is the identity matrix
;		gammat, gammal are the reciprocal Lagrangian Multipliers
;		trial is the trial solution
;		b - is the blurred image in vector form
;		x - is the restore image
; EXAMPLES:
;	Minimum norm solution:
;		conlsq2,blur,psf,result,gammat=0.03
;
;	Minimum difference from trial solution
;		trial = blur		;used blurred image as trial
;		conlsq2,blur,psf,result,trial=trial,gammat=0.005
;
;	Minimize Laplacian
;		conlsq2,blur,psf,result,gammat=0.0,gammal=0.01
;
;	Minimize Laplacian and difference from a trial solution
;		conlsq2,blur,psf,result,gammat=0.02,gammal=0.001,trial=trial
; HISTORY:
;	version 1  D. Lindler  Mar. 1992
;-
;----------------------------------------------------------------------------
; 
; print calling sequence if not supplied
;
	if n_params(0) lt 1 then begin
		print,'CALLING SEQUENCE: conlsq2,blur,psf,result'
		print,'OPTIONAL KEYWORD PARAMETERS:  trial,gammat,gammal'
		retall
	endif
;
; set defaults
;
	if n_elements(gammat) eq 0 then gammat = 0.01
	if n_elements(gammal) eq 0 then gammal = 0.0
;
; compute FFT of the PSF
;
	psf_fft,blur,psf,fftp
;
; construct (HtH + gammat*I + gammal*QtQ)
;
	npix = n_elements(blur)
	hth = float(fft(fftp*conj(fftp),1))*npix
	hth(0)=hth(0)+gammat
	if gammal gt 0.0 then begin
		q = [[0,-1,0],[-1,4,-1],[0,-1,0]]
		psf_fft,blur,q,fftq
		qtq = float(fft(fftq*conj(fftq),1))*npix
		hth = hth+gammal*qtq
	endif
;
; construct   Htb + gammat*trial
;
	bmod = float(fft(fft(blur,-1)*conj(fftp),1))*npix
	if (gammat gt 0.0) and (n_elements(trial) gt 0) then $
			bmod = bmod + gammat*trial
;
; solve for x = result
;
	result = FLOAT(fft(fft(bmod,-1)/fft(hth,-1),1)/npix)
return
end